diff --git a/examples/neuronet_test/in.run b/examples/neuronet_test/in.run index bf06e5999..533775982 100644 --- a/examples/neuronet_test/in.run +++ b/examples/neuronet_test/in.run @@ -1,25 +1,29 @@ # 3d Lennard-Jones melt units metal atom_style atomic region box block 0 12 0 12 0 12 +boundary p p p create_box 3 box create_atoms 1 single 0 0 0 units box create_atoms 2 single 3 3 0 units box create_atoms 3 single 3 0 3 units box mass 1 1.0 mass 2 1.0 mass 3 1.0 pair_style neuronet pair_coeff * * in.const.NN in.params.NN neighbor 0.3 bin - +compute 1 all pe +variable pppp equal c_1 fix 1 all nve -dump 1 all atom 1 dump.out -run 1 +fix 2 all print 1 "potential energy ${pppp}" +dump 1 all custom 1 dump.out id xs ys zs fx fy fz +run 1 +#run 2 \ No newline at end of file diff --git a/lib/neuronet/ryo_force_NN.F90 b/lib/neuronet/ryo_force_NN.F90 index 24bb958c5..124c10126 100644 --- a/lib/neuronet/ryo_force_NN.F90 +++ b/lib/neuronet/ryo_force_NN.F90 @@ -1,530 +1,530 @@ module NN !----------------------------------------------------------------------- ! Time-stamp: <2016-09-06 16:34:13 Ryo KOBAYASHI> !----------------------------------------------------------------------- ! Parallel implementation of neural-network potential with 1 hidden ! layer. It is available for plural number of species. !----------------------------------------------------------------------- !.....cutoff region width ratio to rc real(8):: rcw = 0.9d0 contains subroutine force_NN(namax,natm,tag,ra,nnmax,force,vatom,h,hi,tcom & ,nb,nbmax,lsb,nex,lsrc,myparity,nn,sv,rc,lspr & ,mpi_world,myid,epi,epot,vflag_atom, itype, cnst, iaddr2, iaddr3, & nsp, nhl, max_ncnst, nl, nlmax, hl1, hl2, gsf, dgsf, nal, nnl, wgt11, & - wgt12, wgt21, wgt22, wgt23) + wgt12, wgt21, wgt22, wgt23, rc3) ! namax = nmax ! natm = inum ! tag = atom_type ! ra = fort_x ! no need to tranpose ! nnmax = nnmax ! force = fort_f ! no need to tranpose ! vatom = fort_vatom ! h = h ! hi = h ! tcom = dummy ! nb = nghost ! nbmax = idummy ! lsb = idummy ! nex = idummy ! lsrc = idummy ! myparity = idummy ! nn = idummy ! sv = dummy ! rc = rcin ! lspr = fort_neighbor ! mpi_world = idummy ! myid = idummy ! epi = fort_eatom ! epot = eng_ndwl ! vflag_atom = vflag_atom implicit none !include "mpif.h" !include "./params_unit.h" integer,intent(in):: nhl(0:nlmax+1), itype(nhl(0)), nl, nlmax, max_ncnst, nsp, nal, nnl real(8),intent(in):: cnst(max_ncnst,nhl(0)) integer,intent(in):: iaddr2(2,nsp,nsp),iaddr3(2,nsp,nsp,nsp) real(8),intent(in):: wgt11(nhl(0), nhl(1)),wgt12(nhl(1)) real(8),intent(in):: wgt21(nhl(0),nhl(1)),wgt22(nhl(1),nhl(2)), & wgt23(nhl(2)) integer,intent(in):: namax,natm,nnmax integer,intent(in):: nb,nbmax,lsb(0:nbmax,6),lsrc(6),myparity(3) & ,nn(6),mpi_world,myid,lspr(0:nnmax,namax),nex(3),tag(namax) real(8),intent(in):: ra(3,namax),h(3,3),hi(3,3),sv(3,6) real(8),intent(inout):: tcom,rc real(8),intent(out):: force(3,namax),epi(namax),epot,vatom(6,namax) logical:: vflag_atom !.....local integer:: i,j,k,l,m,n,ixyz,jxyz,is,js,ks,ierr,nbl,ia,ja,nexp,isf & ,icoeff,ihl0,ihl1,ihl2,jj,jsf real(8):: rcin,b_na,at(3),epotl,wgt,hl1i,hl2i,tmp2,tmp1,tmp,tmp3(3) - real(8),save:: rc3 + real(8),intent(in):: rc3 real(8),intent(inout) :: hl1(nhl(1),nal),hl2(nhl(2),nal) real(8),intent(inout) :: gsf(nhl(0),nal),dgsf(3,nhl(0),0:nnl,nal) ! real(8),allocatable:: aml(:,:,:,:),bml(:,:,:,:) !.....1st call logical,save:: l1st=.true. !.....first, calculate all the symmetry functions call eval_sf(nhl(0),namax,natm,nb,nnmax,h,tag,ra & ,lspr,rc,rc3, itype, cnst, nhl, nlmax, max_ncnst & ,iaddr2, iaddr3, nsp, gsf, dgsf, nal, nnl) !.....2nd, calculate the node values by summing contributions from !..... symmetry functions if( nl.eq.1 ) then hl1(1:nhl(1),1:natm)= 0d0 do ia=1,natm !.....debug do ihl1=1,nhl(1) tmp= 0d0 do ihl0=1,nhl(0) tmp= tmp +wgt11(ihl0,ihl1) *gsf(ihl0,ia) enddo hl1(ihl1,ia)= sigmoid(tmp) enddo enddo else if( nl.eq.2 ) then hl1(1:nhl(1),1:natm)= 0d0 hl2(1:nhl(2),1:natm)= 0d0 do ia=1,natm do ihl1=1,nhl(1) tmp= 0d0 do ihl0=1,nhl(0) tmp= tmp +wgt21(ihl0,ihl1) *gsf(ihl0,ia) enddo hl1(ihl1,ia)= sigmoid(tmp) enddo do ihl2=1,nhl(2) tmp= 0d0 do ihl1=1,nhl(1) tmp= tmp +wgt22(ihl1,ihl2) *(hl1(ihl1,ia)-0.5d0) enddo hl2(ihl2,ia)= sigmoid(tmp) enddo enddo endif !.....then calculate the energy of atom by summing up the node values epotl= 0d0 if( nl.eq.1 ) then do ia=1,natm epi(ia)= 0d0 do ihl1=1,nhl(1) epi(ia)= epi(ia) +wgt12(ihl1) *(hl1(ihl1,ia)-0.5d0) enddo epotl=epotl +epi(ia) enddo else if( nl.eq.2 ) then do ia=1,natm epi(ia)= 0d0 do ihl2=1,nhl(2) epi(ia)= epi(ia) +wgt23(ihl2) *(hl2(ihl2,ia)-0.5d0) enddo epotl=epotl +epi(ia) enddo endif !.....sum up for forces force(1:3,1:natm+nb)= 0d0 if( nl.eq.1 ) then do ia=1,natm do ihl1=1,nhl(1) hl1i= hl1(ihl1,ia) tmp= wgt12(ihl1)*hl1i*(1d0-hl1i) do jj=1,lspr(0,ia) ja= lspr(jj,ia) do ihl0=1,nhl(0) force(1:3,ja)=force(1:3,ja) & -tmp*wgt11(ihl0,ihl1)*dgsf(1:3,ihl0,jj,ia) enddo enddo !.....atom ia do ihl0= 1,nhl(0) force(1:3,ia)=force(1:3,ia) & -tmp*wgt11(ihl0,ihl1)*dgsf(1:3,ihl0,0,ia) enddo enddo enddo else if( nl.eq.2 ) then do ia=1,natm do ihl2=1,nhl(2) hl2i= hl2(ihl2,ia) tmp2= wgt23(ihl2) *hl2i*(1d0-hl2i) do ihl1=1,nhl(1) hl1i= hl1(ihl1,ia) tmp1= wgt22(ihl1,ihl2) *hl1i*(1d0-hl1i) do jj=1,lspr(0,ia) ja= lspr(jj,ia) do ihl0=1,nhl(0) force(1:3,ja)=force(1:3,ja) & -tmp2 *tmp1 & *wgt21(ihl0,ihl1)*dgsf(1:3,ihl0,jj,ia) enddo enddo !.....atom ia do ihl0= 1,nhl(0) force(1:3,ia)=force(1:3,ia) & -tmp2*tmp1*wgt21(ihl0,ihl1)*dgsf(1:3,ihl0,0,ia) enddo enddo enddo enddo endif ! call copy_dba_bk(tcom,namax,natm,nbmax,nb,lsb,nex,lsrc,myparity & ! ,nn,mpi_world,force,3) !!$ if( myid.ge.0 ) then !!$ call copy_dba_bk(tcom,namax,natm,nbmax,nb,lsb,lsrc,myparity & !!$ ,nn,mpi_world,force,3) !!$ else !!$ call reduce_dba_bk(natm,namax,tag,force,3) !!$ endif if( vflag_atom ) then call compute_stress(namax,natm,tag,ra,nnmax,vatom,h & ,tcom,nb,nbmax,lsb,nex,lsrc,myparity,nn,rc,lspr & ,mpi_world,myid, nl, wgt11, wgt12, wgt21, wgt22 & , wgt23, nlmax, nhl, hl1, hl2, nal, dgsf, nnl) endif epot= epotl return end subroutine force_NN !======================================================================= subroutine eval_sf(nsf,namax,natm,nb,nnmax,h,tag,ra,lspr,rc,rc3, itype,& cnst, nhl, nlmax, max_ncnst, iaddr2, iaddr3, nsp, gsf, dgsf, nal, nnl) ! ! Evaluate symmetry functions and derivatives for multi-species system. ! implicit none integer,intent(in):: nsf,namax,natm,nb,nnmax,lspr(0:nnmax,namax),& tag(namax), nlmax, nhl(0:nlmax+1), itype(nhl(0)), max_ncnst,& nsp, nal, nnl integer,intent(in):: iaddr2(2,nsp,nsp),iaddr3(2,nsp,nsp,nsp) real(8),intent(in):: h(3,3),ra(3,namax),rc,rc3 ! real(8),intent(out):: gsf(nsf,nal),dgsf(3,nsf,0:nnl,nal) real(8),intent(in):: cnst(max_ncnst,nhl(0)) integer:: isf,isfc,ia,jj,ja,kk,ka,is,js,ks,isfc1,isfc2 real(8):: xi(3),xj(3),xij(3),rij(3),dij,fcij,eta,rs,texp,driji(3), & dfcij,drijj(3),dgdr,xk(3),xik(3),rik(3),dik,fcik,dfcik, & driki(3),drikk(3),almbd,spijk,cs,t1,t2,dgdij,dgdik,dgcs, & dcsdj(3),dcsdk(3),dcsdi(3),tcos,tpoly,a1,a2,tmorse real(8),external:: sprod real(8),intent(out) :: gsf(nhl(0),nal), dgsf(3,nhl(0),0:nnl,nal) gsf(1:nsf,1:nal)= 0d0 dgsf(1:3,1:nsf,0:nnl,1:nal)= 0d0 do ia=1,natm xi(1:3)= ra(1:3,ia) is= tag(ia) do jj=1,lspr(0,ia) ja= lspr(jj,ia) if( ja.eq.ia ) cycle xj(1:3)= ra(1:3,ja) xij(1:3)= xj(1:3)-xi(1:3) rij(1:3)= h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3) dij= sqrt(rij(1)**2 +rij(2)**2 +rij(3)**2) if( dij.ge.rc ) cycle js= tag(ja) isfc=0 driji(1:3)= -rij(1:3)/dij drijj(1:3)= -driji(1:3) fcij= fc(dij,rc) dfcij= dfc(dij,rc) do isf=iaddr2(1,is,js),iaddr2(2,is,js) !!$ isfc= isfc+1 !!$ isf= (icmb2(is,js)-1)*nsfc1 +isfc1 if( itype(isf).eq.1 ) then ! Gaussian eta= cnst(1,isf) rs= cnst(2,isf) !.....function value texp= exp(-eta*(dij-rs)**2) gsf(isf,ia)= gsf(isf,ia) +texp*fcij !.....derivative ! dgsf(ixyz,isf,jj,ia): derivative of isf-th basis of atom-ia ! by ixyz coordinate of atom-jj. ! jj=0 means derivative by atom-ia. dgdr= -2d0*eta*(dij-rs)*texp*fcij +texp*dfcij dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr else if( itype(isf).eq.2 ) then ! cosine a1= cnst(1,isf) !.....func value tcos= (1d0+cos(dij*a1)) gsf(isf,ia)= gsf(isf,ia) +tcos*fcij !.....derivative dgdr= -a1*sin(dij*a1)*fcij +tcos*dfc(dij,rc) dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr else if( itype(isf).eq.3 ) then ! polynomial a1= cnst(1,isf) !.....func value tpoly= 1d0*dij**(-a1) gsf(isf,ia)= gsf(isf,ia) +tpoly*fcij !.....derivative dgdr= -a1*dij**(-a1-1d0)*fcij +tpoly*dfc(dij,rc) dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr else if( itype(isf).eq.4 ) then ! Morse-type a1= cnst(1,isf) a2= cnst(2,isf) !.....func value texp= exp(-a1*(dij-a2)) tmorse= ((1d0-texp)**2 -1d0) gsf(isf,ia)= gsf(isf,ia) +tmorse*fcij !.....derivative dgdr= 2d0*a1*(1d0-texp)*texp*fcij +tmorse*dfcij dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +driji(1:3)*dgdr dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +drijj(1:3)*dgdr endif enddo !!$ fcij= fc(dij,rc) !!$ dfcij= dfc(dij,rc) !!$ driji(1:3)= -rij(1:3)/dij !!$ drijj(1:3)= -driji(1:3) if( dij.gt.rc3 ) cycle do kk=1,lspr(0,ia) ka= lspr(kk,ia) ks= tag(ka) if( iaddr3(1,is,js,ks).lt.0 ) cycle if( ka.eq.ia .or. ka.le.ja ) cycle xk(1:3)= ra(1:3,ka) xik(1:3)= xk(1:3)-xi(1:3) rik(1:3)= h(1:3,1)*xik(1) +h(1:3,2)*xik(2) +h(1:3,3)*xik(3) dik= sqrt(rik(1)**2 +rik(2)**2 +rik(3)**2) if( dik.ge.rc3 ) cycle do isf=iaddr3(1,is,js,ks),iaddr3(2,is,js,ks) !!$ isf= nsfc1*nc1 +(icmb3(is,js,ks)-1)*nsfc2 +isfc2 almbd= cnst(1,isf) t2= (abs(almbd)+1d0)**2 fcik= fc(dik,rc) dfcik= dfc(dik,rc) driki(1:3)= -rik(1:3)/dik drikk(1:3)= -driki(1:3) !.....function value spijk= rij(1)*rik(1) +rij(2)*rik(2) +rij(3)*rik(3) cs= spijk/dij/dik t1= (almbd +cs)**2 gsf(isf,ia)= gsf(isf,ia) +t1/t2 *fcij*fcik !.....derivative dgdij= dfcij *fcik *t1/t2 dgdik= fcij *dfcik *t1/t2 dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) & +dgdij*driji(1:3) +dgdik*driki(1:3) dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +dgdij*drijj(1:3) dgsf(1:3,isf,kk,ia)= dgsf(1:3,isf,kk,ia) +dgdik*drikk(1:3) dgcs= 2d0*(almbd+cs)/t2 *fcij*fcik dcsdj(1:3)= rik(1:3)/dij/dik -rij(1:3)*spijk/dij**3/dik dcsdk(1:3)= rij(1:3)/dij/dik -rik(1:3)*spijk/dik**3/dij dcsdi(1:3)= -dcsdj(1:3) -dcsdk(1:3) dgsf(1:3,isf,0,ia)= dgsf(1:3,isf,0,ia) +dgcs*dcsdi(1:3) dgsf(1:3,isf,jj,ia)= dgsf(1:3,isf,jj,ia) +dgcs*dcsdj(1:3) dgsf(1:3,isf,kk,ia)= dgsf(1:3,isf,kk,ia) +dgcs*dcsdk(1:3) enddo enddo enddo enddo end subroutine eval_sf !======================================================================= function fc(r,rc) implicit none real(8),intent(in):: r,rc real(8):: fc,rs real(8),parameter:: pi= 3.14159265358979d0 rs= rc*rcw if( r.le.rs ) then fc= 1d0 else if( r.gt.rs .and. r.le.rc ) then fc= 0.5d0 *(cos((r-rs)/(rc-rs)*pi)+1d0) else fc= 0d0 endif return end function fc !======================================================================= function dfc(r,rc) implicit none real(8),intent(in):: r,rc real(8):: dfc,rs real(8),parameter:: pi= 3.14159265358979d0 rs= rc*rcw if( r.le.rs ) then dfc= 0d0 else if( r.gt.rs .and. r.le.rc ) then dfc= -pi/2/(rc-rs) *sin((r-rs)/(rc-rs)*pi) else dfc= 0d0 endif return end function dfc !======================================================================= function sigmoid(x) implicit none real(8),intent(in):: x real(8):: sigmoid sigmoid= 1d0/(1d0 +exp(-x)) return end function sigmoid !======================================================================= function dsigmoid(x) implicit none real(8),intent(in):: x real(8):: dsigmoid,sx sx= sigmoid(x) ! dsigmoid= -exp(-x)/(1d0+exp(-x))**2 dsigmoid= sx*(1d0-sx) return end function dsigmoid !======================================================================= function factorial(n,m) ! compute factorial of n, m-times. implicit none integer,intent(in):: n,m real(8):: factorial integer:: i factorial= 1 do i=0,m-1 factorial= factorial*(n-i) enddo return end function factorial !======================================================================= subroutine compute_stress(namax,natm,tag,ra,nnmax,vatom,h & ,tcom,nb,nbmax,lsb,nex,lsrc,myparity,nn,rc,lspr & ,mpi_world,myid, nl, wgt11, wgt12, wgt21, wgt22, wgt23, nlmax, nhl,& hl1, hl2, nal, dgsf, nnl) implicit none integer,intent(in):: namax,natm,nnmax,nb,nbmax,lsb(0:nbmax,6)& ,lsrc(6),myparity(3),nn(6),mpi_world,myid,lspr(0:nnmax,namax)& ,nex(3),tag(namax), nl, nlmax, nhl(0:nlmax+1), nal, nnl real(8),intent(in):: ra(3,namax),h(3,3),rc, dgsf(3,nhl(0),0:nnl,nal) real(8),intent(inout):: tcom real(8),intent(out):: vatom(6,namax) real(8),intent(in):: wgt11(nhl(0), nhl(1)),wgt12(nhl(1)) real(8),intent(in):: wgt21(nhl(0),nhl(1)),wgt22(nhl(1),nhl(2)), & wgt23(nhl(2)) real(8),intent(in) :: hl1(nhl(1),nal) real(8),intent(in):: hl2(nhl(2),nal) integer:: ia,ja,ixyz,jxyz,ihl0,ihl1,ihl2,jj real(8):: xi(3),xj(3),xji(3),rij(3),rji(3),dji,sji(6),sii& ,hl2i,hl2j,tmp2i,tmp2j,hl1i,hl1j,tmp1i,tmp1j, fi(3) vatom(1:6,1:namax) = 0d0 if( nl.eq.1 ) then do ia=1,natm xi(1:3)= ra(1:3,ia) do jj=1,lspr(0,ia) ja= lspr(jj,ia) xj(1:3)= ra(1:3,ja) xji(1:3)= xj(1:3)-xi(1:3) rji(1:3)= h(1:3,1)*xji(1) +h(1:3,2)*xji(2) +h(1:3,3)*xji(3) rij(1:3)= -rji(1:3) dji= sqrt(rji(1)**2 +rji(2)**2 +rji(3)**2) if( dji.ge.rc ) cycle do ihl1=1,nhl(1) hl1i= hl1(ihl1,ia) tmp1i= wgt12(ihl1)*hl1i*(1d0-hl1i) do ihl0=1,nhl(0) fi(1) = tmp1i*wgt11(ihl0,ihl1)*dgsf(1,ihl0,jj,ia) fi(2) = tmp1i*wgt11(ihl0,ihl1)*dgsf(2,ihl0,jj,ia) fi(3) = tmp1i*wgt11(ihl0,ihl1)*dgsf(3,ihl0,jj,ia) sji(1) = -rji(1)* 0.5 * fi(1) sji(2) = -rji(1)* 0.5 * fi(2) sji(3) = -rji(1)* 0.5 * fi(3) sji(4) = -0.25*(rji(1) * fi(2) + rji(2) * fi(1)) sji(5) = -0.25*(rji(1) * fi(3) + rji(3) * fi(1)) sji(6) = -0.25*(rji(2) * fi(3) + rji(3) * fi(2)) vatom(1,ia) = vatom(1,ia) + sji(1) vatom(2,ia) = vatom(2,ia) + sji(2) vatom(3,ia) = vatom(3,ia) + sji(3) vatom(4,ia) = vatom(4,ia) + sji(4) vatom(5,ia) = vatom(5,ia) + sji(5) vatom(6,ia) = vatom(6,ia) + sji(6) vatom(1,ja) = vatom(1,ja) + sji(1) vatom(2,ja) = vatom(2,ja) + sji(2) vatom(3,ja) = vatom(3,ja) + sji(3) vatom(4,ja) = vatom(4,ja) + sji(4) vatom(5,ja) = vatom(5,ja) + sji(5) vatom(6,ja) = vatom(6,ja) + sji(6) enddo enddo enddo enddo else if( nl.eq.2 ) then do ia=1,natm xi(1:3)= ra(1:3,ia) do jj=1,lspr(0,ia) ja= lspr(jj,ia) xj(1:3)= ra(1:3,ja) xji(1:3)= xj(1:3)-xi(1:3) rji(1:3)= h(1:3,1)*xji(1) +h(1:3,2)*xji(2) +h(1:3,3)*xji(3) rij(1:3)= -rji(1:3) dji= sqrt(rji(1)**2 +rji(2)**2 +rji(3)**2) if( dji.ge.rc ) cycle do ihl2=1,nhl(2) hl2i= hl2(ihl2,ia) tmp2i= wgt23(ihl2) *hl2i*(1d0-hl2i) do ihl1=1,nhl(1) hl1i= hl1(ihl1,ia) tmp1i= wgt22(ihl1,ihl2) *hl1i*(1d0-hl1i) do ihl0=1,nhl(0) fi(1) = tmp1i*tmp2i*wgt21(ihl0,ihl1)*dgsf(1,ihl0,jj,ia) fi(2) = tmp1i*tmp2i*wgt21(ihl0,ihl1)*dgsf(2,ihl0,jj,ia) fi(3) = tmp1i*tmp2i*wgt21(ihl0,ihl1)*dgsf(3,ihl0,jj,ia) sji(1) = -rji(1)* 0.5 * fi(1) sji(2) = -rji(1)* 0.5 * fi(2) sji(3) = -rji(1)* 0.5 * fi(3) sji(4) = -0.25*(rji(1) * fi(2) + rji(2) * fi(1)) sji(5) = -0.25*(rji(1) * fi(3) + rji(3) * fi(1)) sji(6) = -0.25*(rji(2) * fi(3) + rji(3) * fi(2)) vatom(1,ia) = vatom(1,ia) + sji(1) vatom(2,ia) = vatom(2,ia) + sji(2) vatom(3,ia) = vatom(3,ia) + sji(3) vatom(4,ia) = vatom(4,ia) + sji(4) vatom(5,ia) = vatom(5,ia) + sji(5) vatom(6,ia) = vatom(6,ia) + sji(6) vatom(1,ja) = vatom(1,ja) + sji(1) vatom(2,ja) = vatom(2,ja) + sji(2) vatom(3,ja) = vatom(3,ja) + sji(3) vatom(4,ja) = vatom(4,ja) + sji(4) vatom(5,ja) = vatom(5,ja) + sji(5) vatom(6,ja) = vatom(6,ja) + sji(6) enddo enddo enddo enddo enddo endif end subroutine compute_stress end module NN !----------------------------------------------------------------------- ! Local Variables: ! compile-command: "make pmd" ! End: diff --git a/src/NEURONET/pair_neuronet.cpp b/src/NEURONET/pair_neuronet.cpp index 18089e219..4b5142eeb 100644 --- a/src/NEURONET/pair_neuronet.cpp +++ b/src/NEURONET/pair_neuronet.cpp @@ -1,531 +1,536 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- This file is written by Till Junge ------------------------------------------------------------------------- */ #include #include #include #include #include "pair_neuronet.h" #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" #include "update.h" #include "integrate.h" //#include "respa.h" #include "math_const.h" #include "memory.h" #include "error.h" #include -#include extern "C" { void __nn_MOD_force_nn(int * namax, int * natm, int* tag, double* ra, int * nnmax, double * force, double * vatom, double *h, double *hi, double * tcom, int * nb, int * nbmax, int * lsb, int * nex, int * lsrc, int * myparity, int * nn, double * sv, double * rc, int *lspr, int * mpi_world, int * myid, double * epi, double * epot, int * vflag_atom, int * itype, double * cnst, int* iaddr2, int * iaddr3, int * nsp, int * nhl, const int * max_ncnst, int * nl, const int * nlmax, double * hl1, double *hl2, double * gsf, double * dgsf, int * nal, int * nnl, double * wgt11, double * wgt12, double * wgt21, - double * wgt22, double * wgt23); + double * wgt22, double * wgt23, double * rc3); } using namespace LAMMPS_NS; using namespace MathConst; /* ---------------------------------------------------------------------- */ PairNeuroNet::PairNeuroNet(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; writedata = 1; this->ncnst_type[0]= 2; // Gaussian this->ncnst_type[1]= 1; // cosine this->ncnst_type[2]= 1; // polynomial this->ncnst_type[3]= 2; // Morse this->ncnst_type[100]= 1; // angular for (int i = 0; i < 100; i++) { this->ncomb_type[i]= 2; // pair this->ncomb_type[100+i] = 3; // triplet } } /* ---------------------------------------------------------------------- */ PairNeuroNet::~PairNeuroNet() { if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); } } /* ---------------------------------------------------------------------- */ void PairNeuroNet::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj; int *ilist,*jlist,*numneigh,**firstneigh; evdwl = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; double **x = atom->x; double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; int nghost = atom->nghost; int nmax = atom->nlocal + atom->nghost; //atom->nmax; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; inum = list->inum; //number of atoms that have a neighbor list ilist = list->ilist; //indices of those atoms numneigh = list->numneigh; // number of their neighbors firstneigh = list->firstneigh; // pointer to first neighbor index (indices contiguous // loop over neighbors of my atoms int nnmax = 0; for (int ii = 0; ii < inum; ++ii) { i = ilist[ii]; nnmax = (nnmax > numneigh[i]) ? nnmax : numneigh[i]; - std::cout << "ii = " << ii << ", i = " << i << ", nnmax = " << nnmax << ", numneigh[" << i << "] = " << numneigh[i] << std::endl; } // ryo neighlist is matrix int *fort_neighbor = new int[(nnmax+1) * nmax] ; int *fort_type = new int[nmax] ; double * fort_x = new double [3*nmax] ; double * fort_f = new double [3*nmax] ; double * fort_vatom = new double [6*nmax] ; double * fort_eatom = new double [nmax] ; double * hl1= new double [nlocal * this->nhl[1]] ; double * hl2= new double [nlocal * this->nhl[2]] ; double * gsf= new double [nlocal * this->nhl[0]] ; double * dgsf= new double [3 * nhl[0] * (nmax+1) * nlocal] ; for (int ii = 0; ii < inum; ++ii) { i = ilist[ii]; for (int dir = 0; dir < 3; ++dir) { fort_x[dir+3*ii] = x[i][dir]; } int *neigh_ptr = firstneigh[i]; int nb_neigh = numneigh[i]; fort_neighbor[0 + (nnmax+1)*ii] = nb_neigh; for (int j = 0; j < nb_neigh; j++) { - fort_neighbor[j+1 + (nnmax+1)*ii] = neigh_ptr[j]; + fort_neighbor[j+1 + (nnmax+1)*ii] = neigh_ptr[j] + 1; } fort_type[ii] = type[i]; } double h[] = {1,0,0, 0,1,0, 0,0,1}; double dummy = 0; int idummy = 0; double delta_eng_vdwl; int fort_vflag_atom = bool(this->vflag_atom); //dummy = __nn_MOD_dsigmoid(&dummy); // __nn_MOD_force_nn(&nmax, // namax &inum, // natm fort_type, // tag fort_x, // ra &nnmax, // nnmax fort_f, // force fort_vatom, // vatom h, // h h, // hi &dummy, // tcom &nghost, // nb &idummy, // nbmax &idummy, // lsb &idummy, // nex &idummy, // lsrc &idummy, // myparity &idummy, // nn &dummy, // sv &rcin, // rc fort_neighbor, // lspr &idummy, // mpi_world &idummy, // myid fort_eatom, // epi &delta_eng_vdwl, // epot &fort_vflag_atom,// vflag_atom &itype[0], // itype &cnst[0], // cnst &iaddr2[0], // iaddr2 &iaddr3[0], // iaddr3 &nsp, // nsp &nhl[0], // nhl &max_ncnst, // max_ncnst &nl, // nl &nlmax, // nlmax hl1, // hl1 hl2, // hl2 gsf, // gsf dgsf, // dgsf &nlocal, // nal &nnmax, // nnl &wgt11[0], // wgt11 &wgt12[0], // wgt12 &wgt21[0], // wgt21 &wgt22[0], // wgt22 - &wgt23[0]); // wgt23 + &wgt23[0], // wgt23 + &rc3); // rc3 this->eng_vdwl += delta_eng_vdwl; // rewrite results into lammps form for (int ii = 0; ii < inum; ++ii) { i = ilist[ii]; if (eflag_atom) { eatom[i] += fort_eatom[ii]; } for (int dir = 0; dir < 3; ++dir) { f[i][dir] += fort_f[dir+3*ii]; } if (vflag_atom) { for (int dir = 0; dir < 6; ++dir) { vatom[i][dir] += fort_vatom[dir+6*ii]; } } } if (vflag_fdotr) virial_fdotr_compute(); delete [] fort_neighbor ; delete [] fort_type ; delete [] fort_x ; delete [] fort_f ; delete [] fort_vatom ; delete [] fort_eatom ; delete [] hl1; delete [] hl2; delete [] gsf; delete [] dgsf; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ void PairNeuroNet::settings(int narg, char **arg) // read_params { if (narg != 0) error->all(FLERR,"Illegal pair_style command"); } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ void PairNeuroNet::coeff(int narg, char **arg) { if (narg != 4) error->all(FLERR,"Not ennough arguments, should be * * constants_input_file params_input_file "); // insure I,J args are * * if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0) error->all(FLERR,"Incorrect args for pair coefficients"); // Read constants input file std::ifstream const_file(arg[2]); //TODO: check that file exist etc int dummy; const_file >> this->nl >> this->nsp >> dummy; // dummy contains number of input nodes this->nhl.push_back(dummy); for (int i = 0; i < nl; ++i) { // read the number of nodes per hidden layer const_file >> dummy; this->nhl.push_back(dummy); } int n = atom->ntypes; memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(cutsq,n+1,n+1,"pair:cutsq"); for (int i = 0; i < n+1; i++) { for (int j = 0; j < n+1; j++) { this->setflag[i][j] = 1; } } // last layer always has one node nhl.push_back(1); + while (nhl.size()nlmax+2) { + nhl.push_back(0); + } // allocate arrays for consts this->itype.resize(this->nhl[0]); this->cnst.resize(this->max_ncnst*this->nhl[0]); this->iaddr2.resize(2*this->nsp*this->nsp); this->iaddr3.resize(2*this->nsp*this->nsp*this->nsp); for (int i = 0; i < 2*this->nsp*this->nsp; i++) { this-> iaddr2[i] = -1; for (int j = 0; j < this->nsp; j++) { this->iaddr3[j+i*nsp] = -1; } } //read in consts int icmb[3] = {0}; // three is hardcoded maximum n in n-body potential - int iap = 0; - int jap = 0; - int kap = 0; + int iap = -1; + int jap = -1; + int kap = -1; int nsf1 = 0, nsf2 = 0; for (int i = 0; i < this->nhl[0]; i++) { const_file >> this->itype[i]; for (int j = 0; j < this->ncomb_type[itype[i]-1]; j++) { const_file >> icmb[j]; icmb[j]--; } double tmp; for (int j = 0; j < this->ncnst_type[itype[i]-1]; j++) { const_file >> tmp; this->cnst[j+i*this->max_ncnst] = tmp; } // distinguish n-body potentials if (itype[i] < 100) { if ((icmb[0] != iap) || (icmb[1] != jap)) { this->iaddr2[0 + 2*(icmb[0] + nsp*icmb[1])] = i+1; this->iaddr2[0 + 2*(icmb[1] + nsp*icmb[0])] = i+1; } this->iaddr2[1 + 2*(icmb[0] + nsp*icmb[1])] = i+1; this->iaddr2[1 + 2*(icmb[1] + nsp*icmb[0])] = i+1; nsf1++; iap = icmb[0]; jap = icmb[1]; } else if (itype[i] < 200){ if ((icmb[0] != iap) || (icmb[1] != jap) || (icmb[2] != kap)) { this->iaddr3[0 + 2*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1; this->iaddr3[0 + 2*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1; } this->iaddr3[1 + 2*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1; this->iaddr3[1 + 2*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1; nsf2++; iap = icmb[0]; jap = icmb[1]; kap = icmb[2]; } else { error->all(FLERR,"Unknown itype"); } } if (this-> nhl[0] != nsf1 + nsf2) { error->all(FLERR, "[Error] nsf.ne.nsf1+nsf2 !!!"); } // read in weights double nwgt[nl+1]; for (int i = 0 ; i < nl+1; i++) { nwgt[i] = nhl[i] * nhl[i+1]; } std::ifstream weights_file(arg[3]); // TODO check file existence etc int ncoeff; weights_file >> ncoeff >> this->rcin >> this->rc3; if (rc3 > rcin) { rc3 = rcin; // TODO add warning as in ryo_force_NN.F90:680 } for (int i = 0; i < atom->ntypes+1; i++) { for (int j = 0; j < atom->ntypes+1; j++) { this->cutsq[i][j] = rcin*rcin; } } int nc = 0; for (int i = 0; i < nl+1; i++) { nc += nwgt[i]; } if( ncoeff != nc ) { error->all(FLERR, "[Error] num of parameters is not correct !!!"); } // allocate them anyways this->wgt11.resize(nhl[0] * nhl[1]); this->wgt12.resize(nhl[1]); this->wgt21.resize(nhl[0] * nhl[1]); this->wgt22.resize(nhl[1] * nhl[2]); this->wgt23.resize(nhl[2]); + double fdummy1, fdummy2; if (nl == 1) { for (int ihl0 = 0; ihl0 < nhl[0]; ihl0++) { for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) { - weights_file >> this->wgt11[ihl0 + ihl1*nhl[0]]; + weights_file >> this->wgt11[ihl0 + ihl1*nhl[0]] >> fdummy1 >> fdummy2; } } for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) { - weights_file >> this->wgt12[ihl1]; + weights_file >> this->wgt12[ihl1] >> fdummy1 >> fdummy2; } } else if (nl == 2) { for (int ihl0 = 0; ihl0 < nhl[0]; ihl0++) { for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) { - weights_file >> this->wgt21[ihl0 + ihl1*nhl[0]]; + weights_file >> this->wgt21[ihl0 + ihl1*nhl[0]] >> fdummy1 >> fdummy2; } } for (int ihl0 = 0; ihl0 < nhl[1]; ihl0++) { for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) { - weights_file >> this->wgt22[ihl0 + ihl1*nhl[0]]; + weights_file >> this->wgt22[ihl0 + ihl1*nhl[0]] >> fdummy1 >> fdummy2; } } for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) { - weights_file >> this->wgt23[ihl1]; + weights_file >> this->wgt23[ihl1] >> fdummy1 >> fdummy2; } } // allocate arrays for weights this->allocated = true; } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ void PairNeuroNet::init_style() { // Don't understand this yet, but seems to be what lammps wants here int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->half=0; + neighbor->requests[irequest]->full=1; } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ double PairNeuroNet::init_one(int i, int j) { return this->rcin; } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ template void vector_write(const std::vector & vec, FILE *fp) { size_t siz = vec.size(); fwrite(&siz, sizeof(size_t), 1, fp); fwrite(&vec[0], sizeof(T), siz, fp); } template void vector_read(std::vector & vec, FILE *fp) { size_t siz; fread(&siz, sizeof(size_t), 1, fp); vec.resize(siz); fwrite(&vec[0], sizeof(T), siz, fp); } template void scalar_write(const T & scalar, FILE *fp) { fwrite(&scalar, sizeof(T), 1, fp); } template void scalar_read(T & scalar, FILE *fp) { fread(&scalar, sizeof(T), 1, fp); } void PairNeuroNet::write_restart(FILE *fp) { write_restart_settings(fp); vector_write(this->wgt11, fp); vector_write(this->wgt12, fp); vector_write(this->wgt21, fp); vector_write(this->wgt22, fp); vector_write(this->wgt23, fp); vector_write(this->nhl, fp); vector_write(this->itype, fp); vector_write(this->iaddr2, fp); vector_write(this->iaddr3, fp); vector_write(this->cnst, fp); scalar_write(this->rcin, fp); scalar_write(this->rc3, fp); } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairNeuroNet::read_restart(FILE *fp) { read_restart_settings(fp); vector_read(this->wgt11, fp); vector_read(this->wgt12, fp); vector_read(this->wgt21, fp); vector_read(this->wgt22, fp); vector_read(this->wgt23, fp); vector_read(this->nhl, fp); vector_read(this->itype, fp); vector_read(this->iaddr2, fp); vector_read(this->iaddr3, fp); vector_read(this->cnst, fp); scalar_read(this->rcin, fp); scalar_read(this->rc3, fp); } /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ void PairNeuroNet::write_restart_settings(FILE *fp) { // think this is supposed to be empty for us } /* ---------------------------------------------------------------------- proc 0 reads from restart file, bcasts ------------------------------------------------------------------------- */ void PairNeuroNet::read_restart_settings(FILE *fp) { int me = comm->me; if (me == 0) { // think this is supposed to be empty for us } //MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); //MPI_Bcast(&offset_flag,1,MPI_INT,0,world); //MPI_Bcast(&mix_flag,1,MPI_INT,0,world); //MPI_Bcast(&tail_flag,1,MPI_INT,0,world); } const int PairNeuroNet::nlmax = 2; const int PairNeuroNet::max_ncnst = 2;