Homec4science
Diffusion lammps (master) Browse

lammps/lib/neuronet/ryo_force_NN.F90
775e8f3cc8f5master

Last Change
Raw File

ryo_force_NN.F90
View Options

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, rc3)
! namax = nmax
! natm = inum
! tag = type
! ra = 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),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) ! pos atom_i
is= tag(ia) ! type atom_i
do jj=1,lspr(0,ia)
ja= lspr(jj,ia) ! neighbor index j
if( ja.eq.ia ) cycle
xj(1:3)= ra(1:3,ja) ! pos atom_j
xij(1:3)= xj(1:3)-xi(1:3) ! vec from atom i to j
rij(1:3)= h(1:3,1)*xij(1) +h(1:3,2)*xij(2) +h(1:3,3)*xij(3) ! unscale
dij= sqrt(rij(1)**2 +rij(2)**2 +rij(3)**2) ! distance ij
if( dij.ge.rc ) cycle
js= tag(ja) ! type atom j
isfc=0
driji(1:3)= -rij(1:3)/dij !unit vec i->j
drijj(1:3)= -driji(1:3) !unit vec j->i
fcij= fc(dij,rc) !eq (4)
dfcij= dfc(dij,rc) !eq (5)
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:
c4science · Help