mirror of https://github.com/lammps/lammps.git
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2134 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
799bf321c9
commit
a77878d0cb
|
|
@ -54,86 +54,88 @@ c Complete the calculation of density
|
||||||
do i = 1,nlocal
|
do i = 1,nlocal
|
||||||
|
|
||||||
elti = fmap(type(i))
|
elti = fmap(type(i))
|
||||||
rho1(i) = 0.d0
|
if (elti.gt.0) then
|
||||||
rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i)
|
rho1(i) = 0.d0
|
||||||
rho3(i) = 0.d0
|
rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i)
|
||||||
do m = 1,3
|
rho3(i) = 0.d0
|
||||||
rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i)
|
do m = 1,3
|
||||||
rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i)
|
rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i)
|
||||||
enddo
|
rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i)
|
||||||
do m = 1,6
|
enddo
|
||||||
rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i)
|
do m = 1,6
|
||||||
enddo
|
rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i)
|
||||||
do m = 1,10
|
enddo
|
||||||
rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i)
|
do m = 1,10
|
||||||
enddo
|
rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i)
|
||||||
|
enddo
|
||||||
if( rho0(i) .gt. 0.0 ) then
|
|
||||||
t_ave(1,i) = t_ave(1,i)/rho0(i)
|
if( rho0(i) .gt. 0.0 ) then
|
||||||
t_ave(2,i) = t_ave(2,i)/rho0(i)
|
t_ave(1,i) = t_ave(1,i)/rho0(i)
|
||||||
t_ave(3,i) = t_ave(3,i)/rho0(i)
|
t_ave(2,i) = t_ave(2,i)/rho0(i)
|
||||||
endif
|
t_ave(3,i) = t_ave(3,i)/rho0(i)
|
||||||
|
endif
|
||||||
Gamma(i) = t_ave(1,i)*rho1(i)
|
|
||||||
$ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i)
|
Gamma(i) = t_ave(1,i)*rho1(i)
|
||||||
|
$ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i)
|
||||||
if( rho0(i) .gt. 0.0 ) then
|
|
||||||
Gamma(i) = Gamma(i)/(rho0(i)*rho0(i))
|
if( rho0(i) .gt. 0.0 ) then
|
||||||
end if
|
Gamma(i) = Gamma(i)/(rho0(i)*rho0(i))
|
||||||
|
end if
|
||||||
call G_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,errorflag)
|
|
||||||
if (errorflag.ne.0) return
|
call G_gam(Gamma(i),ibar_meam(elti),
|
||||||
if (ibar_meam(elti).le.0) then
|
$ gsmooth_factor,G,errorflag)
|
||||||
Gbar = 1.d0
|
if (errorflag.ne.0) return
|
||||||
else
|
if (ibar_meam(elti).le.0) then
|
||||||
call get_shpfcn(shp,lattce_meam(elti,elti))
|
Gbar = 1.d0
|
||||||
gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
|
else
|
||||||
$ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2)
|
call get_shpfcn(shp,lattce_meam(elti,elti))
|
||||||
call G_gam(gam,ibar_meam(elti),gsmooth_factor,
|
gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
|
||||||
$ Gbar,errorflag)
|
$ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2)
|
||||||
endif
|
call G_gam(gam,ibar_meam(elti),gsmooth_factor,
|
||||||
rho(i) = rho0(i) * G
|
$ Gbar,errorflag)
|
||||||
rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar)
|
endif
|
||||||
|
rho(i) = rho0(i) * G
|
||||||
Z = Z_meam(elti)
|
rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar)
|
||||||
call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
|
|
||||||
if (ibar_meam(elti).le.0) then
|
Z = Z_meam(elti)
|
||||||
Gbar = 1.d0
|
call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
|
||||||
dGbar = 0.d0
|
if (ibar_meam(elti).le.0) then
|
||||||
else
|
Gbar = 1.d0
|
||||||
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
dGbar = 0.d0
|
||||||
gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2)
|
else
|
||||||
$ +t_ave(3,i)*shpi(3))/(Z*Z)
|
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
||||||
call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
|
gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2)
|
||||||
$ Gbar,dGbar)
|
$ +t_ave(3,i)*shpi(3))/(Z*Z)
|
||||||
endif
|
call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
|
||||||
denom = 1.d0/(rho0_meam(elti)*Z*Gbar)
|
$ Gbar,dGbar)
|
||||||
dGamma1(i) = (G - 2*dG*Gamma(i))*denom
|
endif
|
||||||
|
denom = 1.d0/(rho0_meam(elti)*Z*Gbar)
|
||||||
if( rho0(i) .ne. 0.d0 ) then
|
dGamma1(i) = (G - 2*dG*Gamma(i))*denom
|
||||||
dGamma2(i) = (dG/rho0(i))*denom
|
|
||||||
else
|
if( rho0(i) .ne. 0.d0 ) then
|
||||||
dGamma2(i) = 0.d0
|
dGamma2(i) = (dG/rho0(i))*denom
|
||||||
end if
|
else
|
||||||
|
dGamma2(i) = 0.d0
|
||||||
dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom
|
end if
|
||||||
|
|
||||||
B = A_meam(elti)*Ec_meam(elti,elti)
|
dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom
|
||||||
|
|
||||||
if( rhob .ne. 0.d0 ) then
|
B = A_meam(elti)*Ec_meam(elti,elti)
|
||||||
fp(i) = B*(log(rhob)+1.d0)
|
|
||||||
if (eflag_either.ne.0) then
|
if( rhob .ne. 0.d0 ) then
|
||||||
if (eflag_global.ne.0) then
|
fp(i) = B*(log(rhob)+1.d0)
|
||||||
eng_vdwl = eng_vdwl + B*rhob*log(rhob)
|
if (eflag_either.ne.0) then
|
||||||
endif
|
if (eflag_global.ne.0) then
|
||||||
if (eflag_atom.ne.0) then
|
eng_vdwl = eng_vdwl + B*rhob*log(rhob)
|
||||||
eatom(i) = eatom(i) + B*rhob*log(rhob)
|
endif
|
||||||
endif
|
if (eflag_atom.ne.0) then
|
||||||
|
eatom(i) = eatom(i) + B*rhob*log(rhob)
|
||||||
|
endif
|
||||||
|
endif
|
||||||
|
else
|
||||||
|
fp(i) = B
|
||||||
endif
|
endif
|
||||||
else
|
|
||||||
fp(i) = B
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
return
|
return
|
||||||
|
|
|
||||||
|
|
@ -94,88 +94,97 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
||||||
drinv = 1.d0/delr_meam
|
drinv = 1.d0/delr_meam
|
||||||
elti = fmap(type(i))
|
elti = fmap(type(i))
|
||||||
|
|
||||||
xitmp = x(1,i)
|
if (elti.gt.0) then
|
||||||
yitmp = x(2,i)
|
|
||||||
zitmp = x(3,i)
|
|
||||||
|
|
||||||
do jn = 1,numneigh
|
xitmp = x(1,i)
|
||||||
j = firstneigh(jn)
|
yitmp = x(2,i)
|
||||||
|
zitmp = x(3,i)
|
||||||
c First compute screening function itself, sij
|
|
||||||
xjtmp = x(1,j)
|
|
||||||
yjtmp = x(2,j)
|
|
||||||
zjtmp = x(3,j)
|
|
||||||
delxij = xjtmp - xitmp
|
|
||||||
delyij = yjtmp - yitmp
|
|
||||||
delzij = zjtmp - zitmp
|
|
||||||
rij2 = delxij*delxij + delyij*delyij + delzij*delzij
|
|
||||||
rij = sqrt(rij2)
|
|
||||||
if (rij.gt.rc_meam) then
|
|
||||||
fcij = 0.0
|
|
||||||
dfcij = 0.d0
|
|
||||||
sij = 0.d0
|
|
||||||
else
|
|
||||||
rnorm = (rc_meam-rij)*drinv
|
|
||||||
call screen(i, j, nmax, x, rij2, sij,
|
|
||||||
$ numneigh_full, firstneigh_full, ntype, type, fmap)
|
|
||||||
call dfcut(rnorm,fc,dfc)
|
|
||||||
fcij = fc
|
|
||||||
dfcij = dfc*drinv
|
|
||||||
endif
|
|
||||||
|
|
||||||
c Now compute derivatives
|
|
||||||
dscrfcn(jn) = 0.d0
|
|
||||||
sfcij = sij*fcij
|
|
||||||
if (sfcij.eq.0.d0.or.sfcij.eq.1.d0) goto 100
|
|
||||||
eltj = fmap(type(j))
|
|
||||||
rbound = ebound_meam(elti,eltj) * rij2
|
|
||||||
do kn = 1,numneigh_full
|
|
||||||
k = firstneigh_full(kn)
|
|
||||||
if (k.eq.j) goto 10
|
|
||||||
eltk = fmap(type(k))
|
|
||||||
xktmp = x(1,k)
|
|
||||||
yktmp = x(2,k)
|
|
||||||
zktmp = x(3,k)
|
|
||||||
delxjk = xktmp - xjtmp
|
|
||||||
delyjk = yktmp - yjtmp
|
|
||||||
delzjk = zktmp - zjtmp
|
|
||||||
rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
|
|
||||||
if (rjk2.gt.rbound) goto 10
|
|
||||||
delxik = xktmp - xitmp
|
|
||||||
delyik = yktmp - yitmp
|
|
||||||
delzik = zktmp - zitmp
|
|
||||||
rik2 = delxik*delxik + delyik*delyik + delzik*delzik
|
|
||||||
if (rik2.gt.rbound) goto 10
|
|
||||||
xik = rik2/rij2
|
|
||||||
xjk = rjk2/rij2
|
|
||||||
a = 1 - (xik-xjk)*(xik-xjk)
|
|
||||||
if (a.eq.0.d0) goto 10
|
|
||||||
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
|
||||||
Cmax = Cmax_meam(elti,eltj,eltk)
|
|
||||||
Cmin = Cmin_meam(elti,eltj,eltk)
|
|
||||||
if (cikj.ge.Cmax.or.cikj.lt.0.d0) then
|
|
||||||
goto 10
|
|
||||||
c Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
|
|
||||||
else
|
|
||||||
delc = Cmax - Cmin
|
|
||||||
cikj = (cikj-Cmin)/delc
|
|
||||||
call dfcut(cikj,sikj,dfikj)
|
|
||||||
coef1 = dfikj/(delc*sikj)
|
|
||||||
call dCfunc(rij2,rik2,rjk2,dCikj)
|
|
||||||
dscrfcn(jn) = dscrfcn(jn) + coef1*dCikj
|
|
||||||
endif
|
|
||||||
10 continue
|
|
||||||
enddo
|
|
||||||
coef1 = sfcij
|
|
||||||
coef2 = sij*dfcij/rij
|
|
||||||
dscrfcn(jn) = dscrfcn(jn)*coef1 - coef2
|
|
||||||
100 continue
|
|
||||||
|
|
||||||
scrfcn(jn) = sij
|
do jn = 1,numneigh
|
||||||
fcpair(jn) = fcij
|
j = firstneigh(jn)
|
||||||
|
|
||||||
enddo
|
eltj = fmap(type(j))
|
||||||
|
if (eltj.gt.0) then
|
||||||
|
|
||||||
|
c First compute screening function itself, sij
|
||||||
|
xjtmp = x(1,j)
|
||||||
|
yjtmp = x(2,j)
|
||||||
|
zjtmp = x(3,j)
|
||||||
|
delxij = xjtmp - xitmp
|
||||||
|
delyij = yjtmp - yitmp
|
||||||
|
delzij = zjtmp - zitmp
|
||||||
|
rij2 = delxij*delxij + delyij*delyij + delzij*delzij
|
||||||
|
rij = sqrt(rij2)
|
||||||
|
if (rij.gt.rc_meam) then
|
||||||
|
fcij = 0.0
|
||||||
|
dfcij = 0.d0
|
||||||
|
sij = 0.d0
|
||||||
|
else
|
||||||
|
rnorm = (rc_meam-rij)*drinv
|
||||||
|
call screen(i, j, nmax, x, rij2, sij,
|
||||||
|
$ numneigh_full, firstneigh_full, ntype, type, fmap)
|
||||||
|
call dfcut(rnorm,fc,dfc)
|
||||||
|
fcij = fc
|
||||||
|
dfcij = dfc*drinv
|
||||||
|
endif
|
||||||
|
|
||||||
|
c Now compute derivatives
|
||||||
|
dscrfcn(jn) = 0.d0
|
||||||
|
sfcij = sij*fcij
|
||||||
|
if (sfcij.eq.0.d0.or.sfcij.eq.1.d0) goto 100
|
||||||
|
rbound = ebound_meam(elti,eltj) * rij2
|
||||||
|
do kn = 1,numneigh_full
|
||||||
|
k = firstneigh_full(kn)
|
||||||
|
if (k.eq.j) goto 10
|
||||||
|
eltk = fmap(type(k))
|
||||||
|
if (eltk.eq.0) goto 10
|
||||||
|
xktmp = x(1,k)
|
||||||
|
yktmp = x(2,k)
|
||||||
|
zktmp = x(3,k)
|
||||||
|
delxjk = xktmp - xjtmp
|
||||||
|
delyjk = yktmp - yjtmp
|
||||||
|
delzjk = zktmp - zjtmp
|
||||||
|
rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
|
||||||
|
if (rjk2.gt.rbound) goto 10
|
||||||
|
delxik = xktmp - xitmp
|
||||||
|
delyik = yktmp - yitmp
|
||||||
|
delzik = zktmp - zitmp
|
||||||
|
rik2 = delxik*delxik + delyik*delyik + delzik*delzik
|
||||||
|
if (rik2.gt.rbound) goto 10
|
||||||
|
xik = rik2/rij2
|
||||||
|
xjk = rjk2/rij2
|
||||||
|
a = 1 - (xik-xjk)*(xik-xjk)
|
||||||
|
if (a.eq.0.d0) goto 10
|
||||||
|
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
|
||||||
|
Cmax = Cmax_meam(elti,eltj,eltk)
|
||||||
|
Cmin = Cmin_meam(elti,eltj,eltk)
|
||||||
|
if (cikj.ge.Cmax.or.cikj.lt.0.d0) then
|
||||||
|
goto 10
|
||||||
|
c Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
|
||||||
|
else
|
||||||
|
delc = Cmax - Cmin
|
||||||
|
cikj = (cikj-Cmin)/delc
|
||||||
|
call dfcut(cikj,sikj,dfikj)
|
||||||
|
coef1 = dfikj/(delc*sikj)
|
||||||
|
call dCfunc(rij2,rik2,rjk2,dCikj)
|
||||||
|
dscrfcn(jn) = dscrfcn(jn) + coef1*dCikj
|
||||||
|
endif
|
||||||
|
10 continue
|
||||||
|
enddo
|
||||||
|
coef1 = sfcij
|
||||||
|
coef2 = sij*dfcij/rij
|
||||||
|
dscrfcn(jn) = dscrfcn(jn)*coef1 - coef2
|
||||||
|
100 continue
|
||||||
|
|
||||||
|
scrfcn(jn) = sij
|
||||||
|
fcpair(jn) = fcij
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
enddo
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -53,7 +53,7 @@ c
|
||||||
dimension t_ave(3,nmax), f(3,nmax), vatom(6,nmax)
|
dimension t_ave(3,nmax), f(3,nmax), vatom(6,nmax)
|
||||||
|
|
||||||
integer i,j,jn,k,kn,kk,m,n,p,q
|
integer i,j,jn,k,kn,kk,m,n,p,q
|
||||||
integer nv2,nv3,elti,eltj,ind
|
integer nv2,nv3,elti,eltj,eltk,ind
|
||||||
real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3
|
real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3
|
||||||
real*8 delik(3),deljk(3),v(6),fi(3),fj(3)
|
real*8 delik(3),deljk(3),v(6),fi(3),fj(3)
|
||||||
real*8 Eu,astar,astarp,third,sixth
|
real*8 Eu,astar,astarp,third,sixth
|
||||||
|
|
@ -92,400 +92,407 @@ c
|
||||||
c Compute forces atom i
|
c Compute forces atom i
|
||||||
|
|
||||||
elti = fmap(type(i))
|
elti = fmap(type(i))
|
||||||
xitmp = x(1,i)
|
|
||||||
yitmp = x(2,i)
|
if (elti.gt.0) then
|
||||||
zitmp = x(3,i)
|
xitmp = x(1,i)
|
||||||
|
yitmp = x(2,i)
|
||||||
c Treat each pair
|
zitmp = x(3,i)
|
||||||
do jn = 1,numneigh
|
|
||||||
|
|
||||||
if (scrfcn(jn).ne.0.d0) then
|
c Treat each pair
|
||||||
|
do jn = 1,numneigh
|
||||||
|
|
||||||
j = firstneigh(jn)
|
j = firstneigh(jn)
|
||||||
|
eltj = fmap(type(j))
|
||||||
|
|
||||||
sij = scrfcn(jn)*fcpair(jn)
|
if (scrfcn(jn).ne.0.d0.and.eltj.gt.0) then
|
||||||
delij(1) = x(1,j) - xitmp
|
|
||||||
delij(2) = x(2,j) - yitmp
|
sij = scrfcn(jn)*fcpair(jn)
|
||||||
delij(3) = x(3,j) - zitmp
|
delij(1) = x(1,j) - xitmp
|
||||||
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
|
delij(2) = x(2,j) - yitmp
|
||||||
$ + delij(3)*delij(3)
|
delij(3) = x(3,j) - zitmp
|
||||||
if (rij2.lt.cutforcesq) then
|
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
|
||||||
rij = sqrt(rij2)
|
$ + delij(3)*delij(3)
|
||||||
r = rij
|
if (rij2.lt.cutforcesq) then
|
||||||
eltj = fmap(type(j))
|
rij = sqrt(rij2)
|
||||||
|
r = rij
|
||||||
|
|
||||||
c Compute phi and phip
|
c Compute phi and phip
|
||||||
ind = eltind(elti,eltj)
|
ind = eltind(elti,eltj)
|
||||||
pp = rij*rdrar + 1.0D0
|
pp = rij*rdrar + 1.0D0
|
||||||
kk = pp
|
kk = pp
|
||||||
kk = min(kk,nrar-1)
|
kk = min(kk,nrar-1)
|
||||||
pp = pp - kk
|
pp = pp - kk
|
||||||
pp = min(pp,1.0D0)
|
pp = min(pp,1.0D0)
|
||||||
phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
|
phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
|
||||||
$ + phirar1(kk,ind))*pp + phirar(kk,ind)
|
$ + phirar1(kk,ind))*pp + phirar(kk,ind)
|
||||||
phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp
|
phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp
|
||||||
$ + phirar4(kk,ind)
|
$ + phirar4(kk,ind)
|
||||||
recip = 1.0d0/r
|
recip = 1.0d0/r
|
||||||
|
|
||||||
if (eflag_either.ne.0) then
|
if (eflag_either.ne.0) then
|
||||||
if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij
|
if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij
|
||||||
if (eflag_atom.ne.0) then
|
if (eflag_atom.ne.0) then
|
||||||
eatom(i) = eatom(i) + 0.5*phi*sij
|
eatom(i) = eatom(i) + 0.5*phi*sij
|
||||||
eatom(j) = eatom(j) + 0.5*phi*sij
|
eatom(j) = eatom(j) + 0.5*phi*sij
|
||||||
|
endif
|
||||||
endif
|
endif
|
||||||
endif
|
|
||||||
|
|
||||||
c write(1,*) "force_meamf: phi: ",phi
|
c write(1,*) "force_meamf: phi: ",phi
|
||||||
c write(1,*) "force_meamf: phip: ",phip
|
c write(1,*) "force_meamf: phip: ",phip
|
||||||
|
|
||||||
c Compute pair densities and derivatives
|
c Compute pair densities and derivatives
|
||||||
invrei = 1.d0/re_meam(elti,elti)
|
invrei = 1.d0/re_meam(elti,elti)
|
||||||
ai = rij*invrei - 1.d0
|
ai = rij*invrei - 1.d0
|
||||||
ro0i = rho0_meam(elti)
|
ro0i = rho0_meam(elti)
|
||||||
rhoa0i = ro0i*exp(-beta0_meam(elti)*ai)
|
rhoa0i = ro0i*exp(-beta0_meam(elti)*ai)
|
||||||
drhoa0i = -beta0_meam(elti)*invrei*rhoa0i
|
drhoa0i = -beta0_meam(elti)*invrei*rhoa0i
|
||||||
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)
|
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)
|
||||||
drhoa1i = -beta1_meam(elti)*invrei*rhoa1i
|
drhoa1i = -beta1_meam(elti)*invrei*rhoa1i
|
||||||
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)
|
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)
|
||||||
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
|
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
|
||||||
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)
|
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)
|
||||||
drhoa3i = -beta3_meam(elti)*invrei*rhoa3i
|
drhoa3i = -beta3_meam(elti)*invrei*rhoa3i
|
||||||
|
|
||||||
if (elti.ne.eltj) then
|
if (elti.ne.eltj) then
|
||||||
invrej = 1.d0/re_meam(eltj,eltj)
|
invrej = 1.d0/re_meam(eltj,eltj)
|
||||||
aj = rij*invrej - 1.d0
|
aj = rij*invrej - 1.d0
|
||||||
ro0j = rho0_meam(eltj)
|
ro0j = rho0_meam(eltj)
|
||||||
rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj)
|
rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj)
|
||||||
drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j
|
drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j
|
||||||
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)
|
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)
|
||||||
drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j
|
drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j
|
||||||
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)
|
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)
|
||||||
drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j
|
drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j
|
||||||
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)
|
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)
|
||||||
drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j
|
drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j
|
||||||
else
|
else
|
||||||
rhoa0j = rhoa0i
|
rhoa0j = rhoa0i
|
||||||
drhoa0j = drhoa0i
|
drhoa0j = drhoa0i
|
||||||
rhoa1j = rhoa1i
|
rhoa1j = rhoa1i
|
||||||
drhoa1j = drhoa1i
|
drhoa1j = drhoa1i
|
||||||
rhoa2j = rhoa2i
|
rhoa2j = rhoa2i
|
||||||
drhoa2j = drhoa2i
|
drhoa2j = drhoa2i
|
||||||
rhoa3j = rhoa3i
|
rhoa3j = rhoa3i
|
||||||
drhoa3j = drhoa3i
|
drhoa3j = drhoa3i
|
||||||
endif
|
endif
|
||||||
|
|
||||||
nv2 = 1
|
|
||||||
nv3 = 1
|
|
||||||
arg1i1 = 0.d0
|
|
||||||
arg1j1 = 0.d0
|
|
||||||
arg1i2 = 0.d0
|
|
||||||
arg1j2 = 0.d0
|
|
||||||
arg1i3 = 0.d0
|
|
||||||
arg1j3 = 0.d0
|
|
||||||
arg3i3 = 0.d0
|
|
||||||
arg3j3 = 0.d0
|
|
||||||
do n = 1,3
|
|
||||||
do p = n,3
|
|
||||||
do q = p,3
|
|
||||||
arg = delij(n)*delij(p)*delij(q)*v3D(nv3)
|
|
||||||
arg1i3 = arg1i3 + Arho3(nv3,i)*arg
|
|
||||||
arg1j3 = arg1j3 - Arho3(nv3,j)*arg
|
|
||||||
nv3 = nv3+1
|
|
||||||
enddo
|
|
||||||
arg = delij(n)*delij(p)*v2D(nv2)
|
|
||||||
arg1i2 = arg1i2 + Arho2(nv2,i)*arg
|
|
||||||
arg1j2 = arg1j2 + Arho2(nv2,j)*arg
|
|
||||||
nv2 = nv2+1
|
|
||||||
enddo
|
|
||||||
arg1i1 = arg1i1 + Arho1(n,i)*delij(n)
|
|
||||||
arg1j1 = arg1j1 - Arho1(n,j)*delij(n)
|
|
||||||
arg3i3 = arg3i3 + Arho3b(n,i)*delij(n)
|
|
||||||
arg3j3 = arg3j3 - Arho3b(n,j)*delij(n)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
c rho0 terms
|
|
||||||
drho0dr1 = drhoa0j * sij
|
|
||||||
drho0dr2 = drhoa0i * sij
|
|
||||||
|
|
||||||
c rho1 terms
|
|
||||||
a1 = 2*sij/rij
|
|
||||||
drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1
|
|
||||||
drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1
|
|
||||||
a1 = 2.d0*sij/rij
|
|
||||||
do m = 1,3
|
|
||||||
drho1drm1(m) = a1*rhoa1j*Arho1(m,i)
|
|
||||||
drho1drm2(m) = -a1*rhoa1i*Arho1(m,j)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
c rho2 terms
|
|
||||||
a2 = 2*sij/rij2
|
|
||||||
drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2
|
|
||||||
$ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij
|
|
||||||
drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2
|
|
||||||
$ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij
|
|
||||||
a2 = 4*sij/rij2
|
|
||||||
do m = 1,3
|
|
||||||
drho2drm1(m) = 0.d0
|
|
||||||
drho2drm2(m) = 0.d0
|
|
||||||
do n = 1,3
|
|
||||||
drho2drm1(m) = drho2drm1(m)
|
|
||||||
$ + Arho2(vind2D(m,n),i)*delij(n)
|
|
||||||
drho2drm2(m) = drho2drm2(m)
|
|
||||||
$ - Arho2(vind2D(m,n),j)*delij(n)
|
|
||||||
enddo
|
|
||||||
drho2drm1(m) = a2*rhoa2j*drho2drm1(m)
|
|
||||||
drho2drm2(m) = -a2*rhoa2i*drho2drm2(m)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
c rho3 terms
|
|
||||||
rij3 = rij*rij2
|
|
||||||
a3 = 2*sij/rij3
|
|
||||||
a3a = 6.d0/5.d0*sij/rij
|
|
||||||
drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3
|
|
||||||
$ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3
|
|
||||||
drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3
|
|
||||||
$ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3
|
|
||||||
a3 = 6*sij/rij3
|
|
||||||
a3a = 6*sij/(5*rij)
|
|
||||||
do m = 1,3
|
|
||||||
drho3drm1(m) = 0.d0
|
|
||||||
drho3drm2(m) = 0.d0
|
|
||||||
nv2 = 1
|
nv2 = 1
|
||||||
|
nv3 = 1
|
||||||
|
arg1i1 = 0.d0
|
||||||
|
arg1j1 = 0.d0
|
||||||
|
arg1i2 = 0.d0
|
||||||
|
arg1j2 = 0.d0
|
||||||
|
arg1i3 = 0.d0
|
||||||
|
arg1j3 = 0.d0
|
||||||
|
arg3i3 = 0.d0
|
||||||
|
arg3j3 = 0.d0
|
||||||
do n = 1,3
|
do n = 1,3
|
||||||
do p = n,3
|
do p = n,3
|
||||||
|
do q = p,3
|
||||||
|
arg = delij(n)*delij(p)*delij(q)*v3D(nv3)
|
||||||
|
arg1i3 = arg1i3 + Arho3(nv3,i)*arg
|
||||||
|
arg1j3 = arg1j3 - Arho3(nv3,j)*arg
|
||||||
|
nv3 = nv3+1
|
||||||
|
enddo
|
||||||
arg = delij(n)*delij(p)*v2D(nv2)
|
arg = delij(n)*delij(p)*v2D(nv2)
|
||||||
drho3drm1(m) = drho3drm1(m)
|
arg1i2 = arg1i2 + Arho2(nv2,i)*arg
|
||||||
$ + Arho3(vind3D(m,n,p),i)*arg
|
arg1j2 = arg1j2 + Arho2(nv2,j)*arg
|
||||||
drho3drm2(m) = drho3drm2(m)
|
nv2 = nv2+1
|
||||||
$ + Arho3(vind3D(m,n,p),j)*arg
|
|
||||||
nv2 = nv2 + 1
|
|
||||||
enddo
|
enddo
|
||||||
|
arg1i1 = arg1i1 + Arho1(n,i)*delij(n)
|
||||||
|
arg1j1 = arg1j1 - Arho1(n,j)*delij(n)
|
||||||
|
arg3i3 = arg3i3 + Arho3b(n,i)*delij(n)
|
||||||
|
arg3j3 = arg3j3 - Arho3b(n,j)*delij(n)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
c rho0 terms
|
||||||
|
drho0dr1 = drhoa0j * sij
|
||||||
|
drho0dr2 = drhoa0i * sij
|
||||||
|
|
||||||
|
c rho1 terms
|
||||||
|
a1 = 2*sij/rij
|
||||||
|
drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1
|
||||||
|
drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1
|
||||||
|
a1 = 2.d0*sij/rij
|
||||||
|
do m = 1,3
|
||||||
|
drho1drm1(m) = a1*rhoa1j*Arho1(m,i)
|
||||||
|
drho1drm2(m) = -a1*rhoa1i*Arho1(m,j)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
c rho2 terms
|
||||||
|
a2 = 2*sij/rij2
|
||||||
|
drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2
|
||||||
|
$ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij
|
||||||
|
drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2
|
||||||
|
$ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij
|
||||||
|
a2 = 4*sij/rij2
|
||||||
|
do m = 1,3
|
||||||
|
drho2drm1(m) = 0.d0
|
||||||
|
drho2drm2(m) = 0.d0
|
||||||
|
do n = 1,3
|
||||||
|
drho2drm1(m) = drho2drm1(m)
|
||||||
|
$ + Arho2(vind2D(m,n),i)*delij(n)
|
||||||
|
drho2drm2(m) = drho2drm2(m)
|
||||||
|
$ - Arho2(vind2D(m,n),j)*delij(n)
|
||||||
|
enddo
|
||||||
|
drho2drm1(m) = a2*rhoa2j*drho2drm1(m)
|
||||||
|
drho2drm2(m) = -a2*rhoa2i*drho2drm2(m)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
c rho3 terms
|
||||||
|
rij3 = rij*rij2
|
||||||
|
a3 = 2*sij/rij3
|
||||||
|
a3a = 6.d0/5.d0*sij/rij
|
||||||
|
drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3
|
||||||
|
$ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3
|
||||||
|
drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3
|
||||||
|
$ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3
|
||||||
|
a3 = 6*sij/rij3
|
||||||
|
a3a = 6*sij/(5*rij)
|
||||||
|
do m = 1,3
|
||||||
|
drho3drm1(m) = 0.d0
|
||||||
|
drho3drm2(m) = 0.d0
|
||||||
|
nv2 = 1
|
||||||
|
do n = 1,3
|
||||||
|
do p = n,3
|
||||||
|
arg = delij(n)*delij(p)*v2D(nv2)
|
||||||
|
drho3drm1(m) = drho3drm1(m)
|
||||||
|
$ + Arho3(vind3D(m,n,p),i)*arg
|
||||||
|
drho3drm2(m) = drho3drm2(m)
|
||||||
|
$ + Arho3(vind3D(m,n,p),j)*arg
|
||||||
|
nv2 = nv2 + 1
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
drho3drm1(m) = (a3*drho3drm1(m) - a3a*Arho3b(m,i))
|
||||||
|
$ *rhoa3j
|
||||||
|
drho3drm2(m) = (-a3*drho3drm2(m) + a3a*Arho3b(m,j))
|
||||||
|
$ *rhoa3i
|
||||||
enddo
|
enddo
|
||||||
drho3drm1(m) = (a3*drho3drm1(m) - a3a*Arho3b(m,i))
|
|
||||||
$ *rhoa3j
|
|
||||||
drho3drm2(m) = (-a3*drho3drm2(m) + a3a*Arho3b(m,j))
|
|
||||||
$ *rhoa3i
|
|
||||||
enddo
|
|
||||||
|
|
||||||
c Compute derivatives of weighting functions t wrt rij
|
c Compute derivatives of weighting functions t wrt rij
|
||||||
t1i = t_ave(1,i)
|
t1i = t_ave(1,i)
|
||||||
t2i = t_ave(2,i)
|
t2i = t_ave(2,i)
|
||||||
t3i = t_ave(3,i)
|
t3i = t_ave(3,i)
|
||||||
t1j = t_ave(1,j)
|
t1j = t_ave(1,j)
|
||||||
t2j = t_ave(2,j)
|
t2j = t_ave(2,j)
|
||||||
t3j = t_ave(3,j)
|
t3j = t_ave(3,j)
|
||||||
|
|
||||||
ai = 0.d0
|
|
||||||
if( rho0(i) .ne. 0.d0 ) then
|
|
||||||
ai = drhoa0j*sij/rho0(i)
|
|
||||||
end if
|
|
||||||
aj = 0.d0
|
|
||||||
if( rho0(j) .ne. 0.d0 ) then
|
|
||||||
aj = drhoa0i*sij/rho0(j)
|
|
||||||
end if
|
|
||||||
|
|
||||||
dt1dr1 = ai*(t1_meam(eltj)-t1i)
|
|
||||||
dt1dr2 = aj*(t1_meam(elti)-t1j)
|
|
||||||
dt2dr1 = ai*(t2_meam(eltj)-t2i)
|
|
||||||
dt2dr2 = aj*(t2_meam(elti)-t2j)
|
|
||||||
dt3dr1 = ai*(t3_meam(eltj)-t3i)
|
|
||||||
dt3dr2 = aj*(t3_meam(elti)-t3j)
|
|
||||||
|
|
||||||
c Compute derivatives of total density wrt rij, sij and rij(3)
|
|
||||||
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
|
||||||
call get_shpfcn(shpj,lattce_meam(eltj,eltj))
|
|
||||||
drhodr1 = dGamma1(i)*drho0dr1
|
|
||||||
$ + dGamma2(i)*
|
|
||||||
$ (dt1dr1*rho1(i)+t1i*drho1dr1
|
|
||||||
$ + dt2dr1*rho2(i)+t2i*drho2dr1
|
|
||||||
$ + dt3dr1*rho3(i)+t3i*drho3dr1)
|
|
||||||
$ - dGamma3(i)*
|
|
||||||
$ (shpi(1)*dt1dr1+shpi(2)*dt2dr1+shpi(3)*dt3dr1)
|
|
||||||
drhodr2 = dGamma1(j)*drho0dr2
|
|
||||||
$ + dGamma2(j)*
|
|
||||||
$ (dt1dr2*rho1(j)+t1j*drho1dr2
|
|
||||||
$ + dt2dr2*rho2(j)+t2j*drho2dr2
|
|
||||||
$ + dt3dr2*rho3(j)+t3j*drho3dr2)
|
|
||||||
$ - dGamma3(j)*
|
|
||||||
$ (shpj(1)*dt1dr2+shpj(2)*dt2dr2+shpj(3)*dt3dr2)
|
|
||||||
do m = 1,3
|
|
||||||
drhodrm1(m) = 0.d0
|
|
||||||
drhodrm2(m) = 0.d0
|
|
||||||
drhodrm1(m) = dGamma2(i)*
|
|
||||||
$ (t1i*drho1drm1(m)
|
|
||||||
$ + t2i*drho2drm1(m)
|
|
||||||
$ + t3i*drho3drm1(m))
|
|
||||||
drhodrm2(m) = dGamma2(j)*
|
|
||||||
$ (t1j*drho1drm2(m)
|
|
||||||
$ + t2j*drho2drm2(m)
|
|
||||||
$ + t3j*drho3drm2(m))
|
|
||||||
enddo
|
|
||||||
|
|
||||||
c Compute derivatives wrt sij, but only if necessary
|
|
||||||
if (dscrfcn(jn).ne.0.d0) then
|
|
||||||
drho0ds1 = rhoa0j
|
|
||||||
drho0ds2 = rhoa0i
|
|
||||||
a1 = 2.d0/rij
|
|
||||||
drho1ds1 = a1*rhoa1j*arg1i1
|
|
||||||
drho1ds2 = a1*rhoa1i*arg1j1
|
|
||||||
a2 = 2.d0/rij2
|
|
||||||
drho2ds1 = a2*rhoa2j*arg1i2
|
|
||||||
$ - 2.d0/3.d0*Arho2b(i)*rhoa2j
|
|
||||||
drho2ds2 = a2*rhoa2i*arg1j2
|
|
||||||
$ - 2.d0/3.d0*Arho2b(j)*rhoa2i
|
|
||||||
a3 = 2.d0/rij3
|
|
||||||
a3a = 6.d0/(5.d0*rij)
|
|
||||||
drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3
|
|
||||||
drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3
|
|
||||||
|
|
||||||
ai = 0.d0
|
ai = 0.d0
|
||||||
if( rho0(i) .ne. 0.d0 ) then
|
if( rho0(i) .ne. 0.d0 ) then
|
||||||
ai = rhoa0j/rho0(i)
|
ai = drhoa0j*sij/rho0(i)
|
||||||
end if
|
end if
|
||||||
aj = 0.d0
|
aj = 0.d0
|
||||||
if( rho0(j) .ne. 0.d0 ) then
|
if( rho0(j) .ne. 0.d0 ) then
|
||||||
aj = rhoa0i/rho0(j)
|
aj = drhoa0i*sij/rho0(j)
|
||||||
end if
|
end if
|
||||||
|
|
||||||
dt1ds1 = ai*(t1_meam(eltj)-t1i)
|
dt1dr1 = ai*(t1_meam(eltj)-t1i)
|
||||||
dt1ds2 = aj*(t1_meam(elti)-t1j)
|
dt1dr2 = aj*(t1_meam(elti)-t1j)
|
||||||
dt2ds1 = ai*(t2_meam(eltj)-t2i)
|
dt2dr1 = ai*(t2_meam(eltj)-t2i)
|
||||||
dt2ds2 = aj*(t2_meam(elti)-t2j)
|
dt2dr2 = aj*(t2_meam(elti)-t2j)
|
||||||
dt3ds1 = ai*(t3_meam(eltj)-t3i)
|
dt3dr1 = ai*(t3_meam(eltj)-t3i)
|
||||||
dt3ds2 = aj*(t3_meam(elti)-t3j)
|
dt3dr2 = aj*(t3_meam(elti)-t3j)
|
||||||
drhods1 = dGamma1(i)*drho0ds1
|
|
||||||
|
c Compute derivatives of total density wrt rij, sij and rij(3)
|
||||||
|
call get_shpfcn(shpi,lattce_meam(elti,elti))
|
||||||
|
call get_shpfcn(shpj,lattce_meam(eltj,eltj))
|
||||||
|
drhodr1 = dGamma1(i)*drho0dr1
|
||||||
$ + dGamma2(i)*
|
$ + dGamma2(i)*
|
||||||
$ (dt1ds1*rho1(i)+t1i*drho1ds1
|
$ (dt1dr1*rho1(i)+t1i*drho1dr1
|
||||||
$ + dt2ds1*rho2(i)+t2i*drho2ds1
|
$ + dt2dr1*rho2(i)+t2i*drho2dr1
|
||||||
$ + dt3ds1*rho3(i)+t3i*drho3ds1)
|
$ + dt3dr1*rho3(i)+t3i*drho3dr1)
|
||||||
$ - dGamma3(i)*
|
$ - dGamma3(i)*
|
||||||
$ (shpi(1)*dt1ds1+shpi(2)*dt2ds1+shpi(3)*dt3ds1)
|
$ (shpi(1)*dt1dr1+shpi(2)*dt2dr1+shpi(3)*dt3dr1)
|
||||||
drhods2 = dGamma1(j)*drho0ds2
|
drhodr2 = dGamma1(j)*drho0dr2
|
||||||
$ + dGamma2(j)*
|
$ + dGamma2(j)*
|
||||||
$ (dt1ds2*rho1(j)+t1j*drho1ds2
|
$ (dt1dr2*rho1(j)+t1j*drho1dr2
|
||||||
$ + dt2ds2*rho2(j)+t2j*drho2ds2
|
$ + dt2dr2*rho2(j)+t2j*drho2dr2
|
||||||
$ + dt3ds2*rho3(j)+t3j*drho3ds2)
|
$ + dt3dr2*rho3(j)+t3j*drho3dr2)
|
||||||
$ - dGamma3(j)*
|
$ - dGamma3(j)*
|
||||||
$ (shpj(1)*dt1ds2+shpj(2)*dt2ds2+shpj(3)*dt3ds2)
|
$ (shpj(1)*dt1dr2+shpj(2)*dt2dr2+shpj(3)*dt3dr2)
|
||||||
endif
|
do m = 1,3
|
||||||
|
drhodrm1(m) = 0.d0
|
||||||
|
drhodrm2(m) = 0.d0
|
||||||
|
drhodrm1(m) = dGamma2(i)*
|
||||||
|
$ (t1i*drho1drm1(m)
|
||||||
|
$ + t2i*drho2drm1(m)
|
||||||
|
$ + t3i*drho3drm1(m))
|
||||||
|
drhodrm2(m) = dGamma2(j)*
|
||||||
|
$ (t1j*drho1drm2(m)
|
||||||
|
$ + t2j*drho2drm2(m)
|
||||||
|
$ + t3j*drho3drm2(m))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
c Compute derivatives wrt sij, but only if necessary
|
||||||
|
if (dscrfcn(jn).ne.0.d0) then
|
||||||
|
drho0ds1 = rhoa0j
|
||||||
|
drho0ds2 = rhoa0i
|
||||||
|
a1 = 2.d0/rij
|
||||||
|
drho1ds1 = a1*rhoa1j*arg1i1
|
||||||
|
drho1ds2 = a1*rhoa1i*arg1j1
|
||||||
|
a2 = 2.d0/rij2
|
||||||
|
drho2ds1 = a2*rhoa2j*arg1i2
|
||||||
|
$ - 2.d0/3.d0*Arho2b(i)*rhoa2j
|
||||||
|
drho2ds2 = a2*rhoa2i*arg1j2
|
||||||
|
$ - 2.d0/3.d0*Arho2b(j)*rhoa2i
|
||||||
|
a3 = 2.d0/rij3
|
||||||
|
a3a = 6.d0/(5.d0*rij)
|
||||||
|
drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3
|
||||||
|
drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3
|
||||||
|
|
||||||
|
ai = 0.d0
|
||||||
|
if( rho0(i) .ne. 0.d0 ) then
|
||||||
|
ai = rhoa0j/rho0(i)
|
||||||
|
end if
|
||||||
|
aj = 0.d0
|
||||||
|
if( rho0(j) .ne. 0.d0 ) then
|
||||||
|
aj = rhoa0i/rho0(j)
|
||||||
|
end if
|
||||||
|
|
||||||
|
dt1ds1 = ai*(t1_meam(eltj)-t1i)
|
||||||
|
dt1ds2 = aj*(t1_meam(elti)-t1j)
|
||||||
|
dt2ds1 = ai*(t2_meam(eltj)-t2i)
|
||||||
|
dt2ds2 = aj*(t2_meam(elti)-t2j)
|
||||||
|
dt3ds1 = ai*(t3_meam(eltj)-t3i)
|
||||||
|
dt3ds2 = aj*(t3_meam(elti)-t3j)
|
||||||
|
drhods1 = dGamma1(i)*drho0ds1
|
||||||
|
$ + dGamma2(i)*
|
||||||
|
$ (dt1ds1*rho1(i)+t1i*drho1ds1
|
||||||
|
$ + dt2ds1*rho2(i)+t2i*drho2ds1
|
||||||
|
$ + dt3ds1*rho3(i)+t3i*drho3ds1)
|
||||||
|
$ - dGamma3(i)*
|
||||||
|
$ (shpi(1)*dt1ds1+shpi(2)*dt2ds1+shpi(3)*dt3ds1)
|
||||||
|
drhods2 = dGamma1(j)*drho0ds2
|
||||||
|
$ + dGamma2(j)*
|
||||||
|
$ (dt1ds2*rho1(j)+t1j*drho1ds2
|
||||||
|
$ + dt2ds2*rho2(j)+t2j*drho2ds2
|
||||||
|
$ + dt3ds2*rho3(j)+t3j*drho3ds2)
|
||||||
|
$ - dGamma3(j)*
|
||||||
|
$ (shpj(1)*dt1ds2+shpj(2)*dt2ds2+shpj(3)*dt3ds2)
|
||||||
|
endif
|
||||||
|
|
||||||
c Compute derivatives of energy wrt rij, sij and rij(3)
|
c Compute derivatives of energy wrt rij, sij and rij(3)
|
||||||
dUdrij = phip*sij
|
dUdrij = phip*sij
|
||||||
$ + fp(i)*drhodr1 + fp(j)*drhodr2
|
$ + fp(i)*drhodr1 + fp(j)*drhodr2
|
||||||
dUdsij = 0.d0
|
dUdsij = 0.d0
|
||||||
if (dscrfcn(jn).ne.0.d0) then
|
if (dscrfcn(jn).ne.0.d0) then
|
||||||
dUdsij = phi
|
dUdsij = phi
|
||||||
$ + fp(i)*drhods1 + fp(j)*drhods2
|
$ + fp(i)*drhods1 + fp(j)*drhods2
|
||||||
endif
|
endif
|
||||||
do m = 1,3
|
do m = 1,3
|
||||||
dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m)
|
dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
c Add the part of the force due to dUdrij and dUdsij
|
c Add the part of the force due to dUdrij and dUdsij
|
||||||
|
|
||||||
force = dUdrij*recip + dUdsij*dscrfcn(jn)
|
force = dUdrij*recip + dUdsij*dscrfcn(jn)
|
||||||
do m = 1,3
|
do m = 1,3
|
||||||
forcem = delij(m)*force + dUdrijm(m)
|
forcem = delij(m)*force + dUdrijm(m)
|
||||||
f(m,i) = f(m,i) + forcem
|
f(m,i) = f(m,i) + forcem
|
||||||
f(m,j) = f(m,j) - forcem
|
f(m,j) = f(m,j) - forcem
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
c Tabulate per-atom virial as symmetrized stress tensor
|
c Tabulate per-atom virial as symmetrized stress tensor
|
||||||
|
|
||||||
if (vflag_atom.ne.0) then
|
if (vflag_atom.ne.0) then
|
||||||
fi(1) = delij(1)*force + dUdrijm(1)
|
fi(1) = delij(1)*force + dUdrijm(1)
|
||||||
fi(2) = delij(2)*force + dUdrijm(2)
|
fi(2) = delij(2)*force + dUdrijm(2)
|
||||||
fi(3) = delij(3)*force + dUdrijm(3)
|
fi(3) = delij(3)*force + dUdrijm(3)
|
||||||
v(1) = -0.5 * (delij(1) * fi(1))
|
v(1) = -0.5 * (delij(1) * fi(1))
|
||||||
v(2) = -0.5 * (delij(2) * fi(2))
|
v(2) = -0.5 * (delij(2) * fi(2))
|
||||||
v(3) = -0.5 * (delij(3) * fi(3))
|
v(3) = -0.5 * (delij(3) * fi(3))
|
||||||
v(4) = -0.25 * (delij(1)*fi(2) + delij(2)*fi(1))
|
v(4) = -0.25 * (delij(1)*fi(2) + delij(2)*fi(1))
|
||||||
v(5) = -0.25 * (delij(1)*fi(3) + delij(3)*fi(1))
|
v(5) = -0.25 * (delij(1)*fi(3) + delij(3)*fi(1))
|
||||||
v(6) = -0.25 * (delij(2)*fi(3) + delij(3)*fi(2))
|
v(6) = -0.25 * (delij(2)*fi(3) + delij(3)*fi(2))
|
||||||
|
|
||||||
vatom(1,i) = vatom(1,i) + v(1)
|
vatom(1,i) = vatom(1,i) + v(1)
|
||||||
vatom(2,i) = vatom(2,i) + v(2)
|
vatom(2,i) = vatom(2,i) + v(2)
|
||||||
vatom(3,i) = vatom(3,i) + v(3)
|
vatom(3,i) = vatom(3,i) + v(3)
|
||||||
vatom(4,i) = vatom(4,i) + v(4)
|
vatom(4,i) = vatom(4,i) + v(4)
|
||||||
vatom(5,i) = vatom(5,i) + v(5)
|
vatom(5,i) = vatom(5,i) + v(5)
|
||||||
vatom(6,i) = vatom(6,i) + v(6)
|
vatom(6,i) = vatom(6,i) + v(6)
|
||||||
vatom(1,j) = vatom(1,j) + v(1)
|
vatom(1,j) = vatom(1,j) + v(1)
|
||||||
vatom(2,j) = vatom(2,j) + v(2)
|
vatom(2,j) = vatom(2,j) + v(2)
|
||||||
vatom(3,j) = vatom(3,j) + v(3)
|
vatom(3,j) = vatom(3,j) + v(3)
|
||||||
vatom(4,j) = vatom(4,j) + v(4)
|
vatom(4,j) = vatom(4,j) + v(4)
|
||||||
vatom(5,j) = vatom(5,j) + v(5)
|
vatom(5,j) = vatom(5,j) + v(5)
|
||||||
vatom(6,j) = vatom(6,j) + v(6)
|
vatom(6,j) = vatom(6,j) + v(6)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
c Now compute forces on other atoms k due to change in sij
|
c Now compute forces on other atoms k due to change in sij
|
||||||
|
|
||||||
if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100
|
if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100
|
||||||
do kn = 1,numneigh_full
|
do kn = 1,numneigh_full
|
||||||
k = firstneigh_full(kn)
|
k = firstneigh_full(kn)
|
||||||
if (k.ne.j) then
|
eltk = fmap(type(k))
|
||||||
call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
|
if (k.ne.j.and.eltk.gt.0) then
|
||||||
$ ntype,type,fmap,x,scrfcn,fcpair)
|
call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
|
||||||
if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then
|
$ ntype,type,fmap,x,scrfcn,fcpair)
|
||||||
force1 = dUdsij*dsij1
|
if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then
|
||||||
force2 = dUdsij*dsij2
|
force1 = dUdsij*dsij1
|
||||||
do m = 1,3
|
force2 = dUdsij*dsij2
|
||||||
delik(m) = x(m,k) - x(m,i)
|
do m = 1,3
|
||||||
deljk(m) = x(m,k) - x(m,j)
|
delik(m) = x(m,k) - x(m,i)
|
||||||
enddo
|
deljk(m) = x(m,k) - x(m,j)
|
||||||
do m = 1,3
|
enddo
|
||||||
f(m,i) = f(m,i) + force1*delik(m)
|
do m = 1,3
|
||||||
f(m,j) = f(m,j) + force2*deljk(m)
|
f(m,i) = f(m,i) + force1*delik(m)
|
||||||
f(m,k) = f(m,k) - force1*delik(m)
|
f(m,j) = f(m,j) + force2*deljk(m)
|
||||||
$ - force2*deljk(m)
|
f(m,k) = f(m,k) - force1*delik(m)
|
||||||
enddo
|
$ - force2*deljk(m)
|
||||||
|
enddo
|
||||||
|
|
||||||
c Tabulate per-atom virial as symmetrized stress tensor
|
c Tabulate per-atom virial as symmetrized stress tensor
|
||||||
|
|
||||||
if (vflag_atom.ne.0) then
|
if (vflag_atom.ne.0) then
|
||||||
fi(1) = force1*delik(1)
|
fi(1) = force1*delik(1)
|
||||||
fi(2) = force1*delik(2)
|
fi(2) = force1*delik(2)
|
||||||
fi(3) = force1*delik(3)
|
fi(3) = force1*delik(3)
|
||||||
fj(1) = force2*deljk(1)
|
fj(1) = force2*deljk(1)
|
||||||
fj(2) = force2*deljk(2)
|
fj(2) = force2*deljk(2)
|
||||||
fj(3) = force2*deljk(3)
|
fj(3) = force2*deljk(3)
|
||||||
v(1) = -third * (delik(1)*fi(1) + deljk(1)*fj(1))
|
v(1) = -third * (delik(1)*fi(1) + deljk(1)*fj(1))
|
||||||
v(2) = -third * (delik(2)*fi(2) + deljk(2)*fj(2))
|
v(2) = -third * (delik(2)*fi(2) + deljk(2)*fj(2))
|
||||||
v(3) = -third * (delik(3)*fi(3) + deljk(3)*fj(3))
|
v(3) = -third * (delik(3)*fi(3) + deljk(3)*fj(3))
|
||||||
v(4) = -sixth * (delik(1)*fi(2) + deljk(1)*fj(2) +
|
v(4) = -sixth * (delik(1)*fi(2) + deljk(1)*fj(2) +
|
||||||
$ delik(2)*fi(1) + deljk(2)*fj(1))
|
$ delik(2)*fi(1) + deljk(2)*fj(1))
|
||||||
v(5) = -sixth * (delik(1)*fi(3) + deljk(1)*fj(3) +
|
v(5) = -sixth * (delik(1)*fi(3) + deljk(1)*fj(3) +
|
||||||
$ delik(3)*fi(1) + deljk(3)*fj(1))
|
$ delik(3)*fi(1) + deljk(3)*fj(1))
|
||||||
v(6) = -sixth * (delik(2)*fi(3) + deljk(2)*fj(3) +
|
v(6) = -sixth * (delik(2)*fi(3) + deljk(2)*fj(3) +
|
||||||
$ delik(3)*fi(2) + deljk(3)*fj(2))
|
$ delik(3)*fi(2) + deljk(3)*fj(2))
|
||||||
|
|
||||||
|
vatom(1,i) = vatom(1,i) + v(1)
|
||||||
|
vatom(2,i) = vatom(2,i) + v(2)
|
||||||
|
vatom(3,i) = vatom(3,i) + v(3)
|
||||||
|
vatom(4,i) = vatom(4,i) + v(4)
|
||||||
|
vatom(5,i) = vatom(5,i) + v(5)
|
||||||
|
vatom(6,i) = vatom(6,i) + v(6)
|
||||||
|
vatom(1,j) = vatom(1,j) + v(1)
|
||||||
|
vatom(2,j) = vatom(2,j) + v(2)
|
||||||
|
vatom(3,j) = vatom(3,j) + v(3)
|
||||||
|
vatom(4,j) = vatom(4,j) + v(4)
|
||||||
|
vatom(5,j) = vatom(5,j) + v(5)
|
||||||
|
vatom(6,j) = vatom(6,j) + v(6)
|
||||||
|
vatom(1,k) = vatom(1,k) + v(1)
|
||||||
|
vatom(2,k) = vatom(2,k) + v(2)
|
||||||
|
vatom(3,k) = vatom(3,k) + v(3)
|
||||||
|
vatom(4,k) = vatom(4,k) + v(4)
|
||||||
|
vatom(5,k) = vatom(5,k) + v(5)
|
||||||
|
vatom(6,k) = vatom(6,k) + v(6)
|
||||||
|
endif
|
||||||
|
|
||||||
vatom(1,i) = vatom(1,i) + v(1)
|
|
||||||
vatom(2,i) = vatom(2,i) + v(2)
|
|
||||||
vatom(3,i) = vatom(3,i) + v(3)
|
|
||||||
vatom(4,i) = vatom(4,i) + v(4)
|
|
||||||
vatom(5,i) = vatom(5,i) + v(5)
|
|
||||||
vatom(6,i) = vatom(6,i) + v(6)
|
|
||||||
vatom(1,j) = vatom(1,j) + v(1)
|
|
||||||
vatom(2,j) = vatom(2,j) + v(2)
|
|
||||||
vatom(3,j) = vatom(3,j) + v(3)
|
|
||||||
vatom(4,j) = vatom(4,j) + v(4)
|
|
||||||
vatom(5,j) = vatom(5,j) + v(5)
|
|
||||||
vatom(6,j) = vatom(6,j) + v(6)
|
|
||||||
vatom(1,k) = vatom(1,k) + v(1)
|
|
||||||
vatom(2,k) = vatom(2,k) + v(2)
|
|
||||||
vatom(3,k) = vatom(3,k) + v(3)
|
|
||||||
vatom(4,k) = vatom(4,k) + v(4)
|
|
||||||
vatom(5,k) = vatom(5,k) + v(5)
|
|
||||||
vatom(6,k) = vatom(6,k) + v(6)
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
endif
|
endif
|
||||||
endif
|
|
||||||
c end of k loop
|
c end of k loop
|
||||||
enddo
|
enddo
|
||||||
|
endif
|
||||||
|
100 continue
|
||||||
endif
|
endif
|
||||||
100 continue
|
|
||||||
endif
|
|
||||||
c end of j loop
|
c end of j loop
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
c else if elti=0, this is not a meam atom
|
||||||
|
endif
|
||||||
|
|
||||||
return
|
return
|
||||||
end
|
end
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue