diff --git a/lib/meam/meam_data.F b/lib/meam/meam_data.F
index 8562e9f03..389b2c64e 100755
--- a/lib/meam/meam_data.F
+++ b/lib/meam/meam_data.F
@@ -1,77 +1,83 @@
module meam_data
integer, parameter :: maxelt = 5
c cutforce = force cutoff
c cutforcesq = force cutoff squared
real*8 cutforce,cutforcesq
c Ec_meam = cohesive energy
c re_meam = nearest-neighbor distance
c Omega_meam = atomic volume
c B_meam = bulk modulus
c Z_meam = number of first neighbors for reference structure
c ielt_meam = atomic number of element
c A_meam = adjustable parameter
c alpha_meam = sqrt(9*Omega*B/Ec)
c rho0_meam = density scaling parameter
c delta_meam = heat of formation for alloys
c beta[0-3]_meam = electron density constants
c t[0-3]_meam = coefficients on densities in Gamma computation
+c rho_ref_meam = background density for reference structure
c ibar_meam(i) = selection parameter for Gamma function for elt i,
c lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j)
c neltypes = maximum number of element type defined
c eltind = index number of pair (similar to Voigt notation; ij = ji)
c phir = pair potential function array
c phirar[1-6] = spline coeffs
c attrac_meam = attraction parameter in Rose energy
c repuls_meam = repulsion parameter in Rose energy
c nn2_meam = 1 if second nearest neighbors are to be computed, else 0
+c zbl_meam = 1 if zbl potential for small r to be use, else 0
c Cmin_meam, Cmax_meam = min and max values in screening cutoff
c rc_meam = cutoff distance for meam
c delr_meam = cutoff region for meam
c ebound_meam = factor giving maximum boundary of sceen fcn ellipse
c augt1 = flag for whether t1 coefficient should be augmented
c ialloy = flag for newer alloy formulation (as in dynamo code)
+c mix_ref_t = flag to recover "old" way of computing t in reference config
+c erose_form = selection parameter for form of E_rose function
c gsmooth_factor = factor determining length of G smoothing region
c vind[23]D = Voight notation index maps for 2 and 3D
c v2D,v3D = array of factors to apply for Voight notation
c nr,dr = pair function discretization parameters
c nrar,rdrar = spline coeff array parameters
real*8 Ec_meam(maxelt,maxelt),re_meam(maxelt,maxelt)
real*8 Omega_meam(maxelt),Z_meam(maxelt)
real*8 A_meam(maxelt),alpha_meam(maxelt,maxelt),rho0_meam(maxelt)
real*8 delta_meam(maxelt,maxelt)
real*8 beta0_meam(maxelt),beta1_meam(maxelt)
real*8 beta2_meam(maxelt),beta3_meam(maxelt)
real*8 t0_meam(maxelt),t1_meam(maxelt)
real*8 t2_meam(maxelt),t3_meam(maxelt)
+ real*8 rho_ref_meam(maxelt)
integer ibar_meam(maxelt),ielt_meam(maxelt)
character*3 lattce_meam(maxelt,maxelt)
integer nn2_meam(maxelt,maxelt)
+ integer zbl_meam(maxelt,maxelt)
integer eltind(maxelt,maxelt)
integer neltypes
real*8, allocatable :: phir(:,:)
real*8, allocatable :: phirar(:,:),phirar1(:,:),phirar2(:,:),
$ phirar3(:,:),phirar4(:,:),phirar5(:,:),phirar6(:,:)
real*8 attrac_meam(maxelt,maxelt),repuls_meam(maxelt,maxelt)
real*8 Cmin_meam(maxelt,maxelt,maxelt)
real*8 Cmax_meam(maxelt,maxelt,maxelt)
real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt)
- integer augt1, ialloy
+ integer augt1, ialloy, mix_ref_t, erose_form
real*8 gsmooth_factor
integer vind2D(3,3),vind3D(3,3,3)
integer v2D(6),v3D(10)
integer nr,nrar
real*8 dr,rdrar
end module
diff --git a/lib/meam/meam_dens_final.F b/lib/meam/meam_dens_final.F
index d39de8ba6..35378a83d 100755
--- a/lib/meam/meam_dens_final.F
+++ b/lib/meam/meam_dens_final.F
@@ -1,222 +1,246 @@
c Extern "C" declaration has the form:
c
c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
c int *, int *, int *,
c double *, double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *, int *);
c
c Call from pair_meam.cpp has the form:
c
c meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
c &eng_vdwl,eatom,ntype,type,fmap,
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
c &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
c
subroutine meam_dens_final(nlocal, nmax,
$ eflag_either, eflag_global, eflag_atom, eng_vdwl, eatom,
$ ntype, type, fmap,
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave,
$ Gamma, dGamma1, dGamma2, dGamma3,
$ rho, rho0, rho1, rho2, rho3, fp, errorflag)
use meam_data
implicit none
integer nlocal, nmax, eflag_either, eflag_global, eflag_atom
integer ntype, type, fmap
real*8 eng_vdwl, eatom, Arho1, Arho2
real*8 Arho2b, Arho3, Arho3b
real*8 t_ave, tsq_ave
real*8 Gamma, dGamma1, dGamma2, dGamma3
real*8 rho, rho0, rho1, rho2, rho3
real*8 fp
integer errorflag
dimension eatom(nmax)
dimension type(nmax), fmap(ntype)
dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax)
dimension Arho3(10,nmax), Arho3b(3,nmax), t_ave(3,nmax)
dimension tsq_ave(3,nmax)
dimension Gamma(nmax), dGamma1(nmax), dGamma2(nmax)
dimension dGamma3(nmax), rho(nmax), rho0(nmax)
dimension rho1(nmax), rho2(nmax), rho3(nmax)
dimension fp(nmax)
integer i, elti
integer m
real*8 rhob, G, dG, Gbar, dGbar, gam, shp(3), shpi(3), Z
- real*8 B, denom
+ real*8 B, denom, rho_bkgd
c Complete the calculation of density
do i = 1,nlocal
elti = fmap(type(i))
if (elti.gt.0) then
rho1(i) = 0.d0
rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i)
rho3(i) = 0.d0
do m = 1,3
rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i)
rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i)
enddo
do m = 1,6
rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i)
enddo
do m = 1,10
rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i)
enddo
if( rho0(i) .gt. 0.0 ) then
if (ialloy.eq.1) then
t_ave(1,i) = t_ave(1,i)/tsq_ave(1,i)
t_ave(2,i) = t_ave(2,i)/tsq_ave(2,i)
t_ave(3,i) = t_ave(3,i)/tsq_ave(3,i)
+ else if (ialloy.eq.2) then
+ t_ave(1,i) = t1_meam(elti)
+ t_ave(2,i) = t2_meam(elti)
+ t_ave(3,i) = t3_meam(elti)
else
t_ave(1,i) = t_ave(1,i)/rho0(i)
t_ave(2,i) = t_ave(2,i)/rho0(i)
t_ave(3,i) = t_ave(3,i)/rho0(i)
endif
endif
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))
end if
-
+
+ Z = Z_meam(elti)
+
call G_gam(Gamma(i),ibar_meam(elti),
$ gsmooth_factor,G,errorflag)
if (errorflag.ne.0) return
if (ibar_meam(elti).le.0) then
Gbar = 1.d0
else
call get_shpfcn(shp,lattce_meam(elti,elti))
- gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2)
- $ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2)
+ if (mix_ref_t.eq.1) then
+ gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2)
+ $ +t_ave(3,i)*shpi(3))/(Z*Z)
+ else
+ gam = (t1_meam(elti)*shp(1)+t2_meam(elti)*shp(2)
+ $ +t3_meam(elti)*shp(3))/(Z*Z)
+ endif
call G_gam(gam,ibar_meam(elti),gsmooth_factor,
$ Gbar,errorflag)
endif
rho(i) = rho0(i) * G
- rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar)
-
- Z = Z_meam(elti)
- call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
- if (ibar_meam(elti).le.0) then
- Gbar = 1.d0
- dGbar = 0.d0
+
+ if (mix_ref_t.eq.1) then
+ if (ibar_meam(elti).le.0) then
+ Gbar = 1.d0
+ dGbar = 0.d0
+ else
+ call get_shpfcn(shpi,lattce_meam(elti,elti))
+ gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2)
+ $ +t_ave(3,i)*shpi(3))/(Z*Z)
+ call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
+ $ Gbar,dGbar)
+ endif
+ rho_bkgd = rho0_meam(elti)*Z*Gbar
else
- call get_shpfcn(shpi,lattce_meam(elti,elti))
- gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2)
- $ +t_ave(3,i)*shpi(3))/(Z*Z)
- call dG_gam(gam,ibar_meam(elti),gsmooth_factor,
- $ Gbar,dGbar)
+ rho_bkgd = rho_ref_meam(elti)
endif
- denom = 1.d0/(rho0_meam(elti)*Z*Gbar)
+ rhob = rho(i)/rho_bkgd
+ denom = 1.d0/rho_bkgd
+
+ call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG)
+
dGamma1(i) = (G - 2*dG*Gamma(i))*denom
if( rho0(i) .ne. 0.d0 ) then
dGamma2(i) = (dG/rho0(i))*denom
else
dGamma2(i) = 0.d0
end if
-
- dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom
-
+
+c dGamma3 is nonzero only if we are using the "mixed" rule for
+c computing t in the reference system (which is not correct, but
+c included for backward compatibility
+ if (mix_ref_t.eq.1) then
+ dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom
+ else
+ dGamma3(i) = 0.0
+ endif
+
B = A_meam(elti)*Ec_meam(elti,elti)
if( rhob .ne. 0.d0 ) then
fp(i) = B*(log(rhob)+1.d0)
if (eflag_either.ne.0) then
if (eflag_global.ne.0) then
eng_vdwl = eng_vdwl + B*rhob*log(rhob)
endif
if (eflag_atom.ne.0) then
eatom(i) = eatom(i) + B*rhob*log(rhob)
endif
endif
else
fp(i) = B
endif
endif
enddo
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine G_gam(Gamma,ibar,gsmooth_factor,G,errorflag)
c Compute G(Gamma) based on selection flag ibar:
c 0 => G = sqrt(1+Gamma)
c 1 => G = exp(Gamma/2)
c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma))
c 4 => G = sqrt(1+Gamma)
implicit none
real*8 Gamma,G
real*8 gsmooth_factor, gsmooth_switchpoint
integer ibar, errorflag
if (ibar.eq.0.or.ibar.eq.4) then
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1)
if (Gamma.lt.gsmooth_switchpoint) then
c e.g. gsmooth_factor is 99, then:
c gsmooth_switchpoint = -0.99
c G = 0.01*(-0.99/Gamma)**99
G = 1/(gsmooth_factor+1)
$ *(gsmooth_switchpoint/Gamma)**gsmooth_factor
G = sqrt(G)
else
G = sqrt(1.d0+Gamma)
endif
else if (ibar.eq.1) then
G = exp(Gamma/2.d0)
else if (ibar.eq.3) then
G = 2.d0/(1.d0+exp(-Gamma))
else
errorflag = 1
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine dG_gam(Gamma,ibar,gsmooth_factor,G,dG)
c Compute G(Gamma) and dG(gamma) based on selection flag ibar:
c 0 => G = sqrt(1+Gamma)
c 1 => G = exp(Gamma/2)
c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma))
c 4 => G = sqrt(1+Gamma)
real*8 Gamma,G,dG
real*8 gsmooth_factor, gsmooth_switchpoint
integer ibar
if (ibar.eq.0.or.ibar.eq.4) then
gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1)
if (Gamma.lt.gsmooth_switchpoint) then
c e.g. gsmooth_factor is 99, then:
c gsmooth_switchpoint = -0.99
c G = 0.01*(-0.99/Gamma)**99
G = 1/(gsmooth_factor+1)
$ *(gsmooth_switchpoint/Gamma)**gsmooth_factor
G = sqrt(G)
dG = -gsmooth_factor*G/(2.0*Gamma)
else
G = sqrt(1.d0+Gamma)
dG = 1.d0/(2.d0*G)
endif
else if (ibar.eq.1) then
G = exp(Gamma/2.d0)
dG = G/2.d0
else if (ibar.eq.3) then
G = 2.d0/(1.d0+exp(-Gamma))
dG = G*(2.d0-G)/2
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
diff --git a/lib/meam/meam_dens_init.F b/lib/meam/meam_dens_init.F
index 53ad0956d..52adb75d0 100755
--- a/lib/meam/meam_dens_init.F
+++ b/lib/meam/meam_dens_init.F
@@ -1,561 +1,564 @@
c Extern "C" declaration has the form:
c
c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *,
c int *, int *, int *, int *,
c double *, double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *, int *);
c
c
c Call from pair_meam.cpp has the form:
c
c meam_dens_init_(&i,&nmax,ntype,type,fmap,&x[0][0],
c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
c rho0,&arho1[0][0],&arho2[0][0],arho2b,
c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);
c
subroutine meam_dens_init(i, nmax,
$ ntype, type, fmap, x,
$ numneigh, firstneigh,
$ numneigh_full, firstneigh_full,
$ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b,
$ arho3, arho3b, t_ave, tsq_ave, errorflag)
use meam_data
implicit none
integer i, nmax, ntype, type, fmap
real*8 x
integer numneigh, firstneigh, numneigh_full, firstneigh_full
real*8 scrfcn, dscrfcn, fcpair
real*8 rho0, arho1, arho2
real*8 arho2b, arho3, arho3b, t_ave, tsq_ave
integer errorflag
integer j,jn
dimension x(3,nmax)
dimension type(nmax), fmap(ntype)
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh)
dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax)
dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax)
dimension t_ave(3,nmax), tsq_ave(3,nmax)
errorflag = 0
c Compute screening function and derivatives
call getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
$ numneigh, firstneigh,
$ numneigh_full, firstneigh_full,
$ ntype, type, fmap)
c Calculate intermediate density terms to be communicated
call calc_rho1(i, nmax, ntype, type, fmap, x,
$ numneigh, firstneigh,
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
$ arho3, arho3b, t_ave, tsq_ave)
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x,
$ numneigh, firstneigh,
$ numneigh_full, firstneigh_full,
$ ntype, type, fmap)
use meam_data
implicit none
integer i, nmax
real*8 scrfcn, dscrfcn, fcpair, x
integer numneigh, firstneigh, numneigh_full, firstneigh_full
integer ntype, type, fmap
dimension scrfcn(numneigh), dscrfcn(numneigh)
dimension fcpair(numneigh), x(3,nmax)
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
dimension type(nmax), fmap(ntype)
integer jn,j,kn,k
integer elti,eltj,eltk
real*8 xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij
real*8 xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2,rik
real*8 xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2,rjk
real*8 xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj
real*8 Cmin,Cmax,delc,ebound,rbound,a,coef1,coef2
real*8 coef1a,coef1b,coef2a,coef2b
real*8 dcikj
real*8 dC1a,dC1b,dC2a,dC2b
real*8 rnorm,fc,dfc,drinv
drinv = 1.d0/delr_meam
elti = fmap(type(i))
if (elti.gt.0) then
xitmp = x(1,i)
yitmp = x(2,i)
zitmp = x(3,i)
do jn = 1,numneigh
j = firstneigh(jn)
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)
c if a < 0, then ellipse equation doesn't describe this case and
c atom k can't possibly screen i-j
if (a.le.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) then
goto 10
c Note that cikj may be slightly negative (within numerical
c tolerance) if atoms are colinear, so don't reject that case here
c (other negative cikj cases were handled by the test on "a" above)
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
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine calc_rho1(i, nmax, ntype, type, fmap, x,
$ numneigh, firstneigh,
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
$ arho3, arho3b, t_ave, tsq_ave)
use meam_data
implicit none
integer i, nmax, ntype, type, fmap
real*8 x
integer numneigh, firstneigh
real*8 scrfcn, fcpair, rho0, arho1, arho2
real*8 arho2b, arho3, arho3b, t_ave, tsq_ave
dimension type(nmax), fmap(ntype), x(3,nmax)
dimension firstneigh(numneigh)
dimension scrfcn(numneigh), fcpair(numneigh)
dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax)
dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax)
dimension t_ave(3,nmax), tsq_ave(3,nmax)
integer jn,j,m,n,p,elti,eltj
integer nv2,nv3
real*8 xtmp,ytmp,ztmp,delij(3),rij2,rij,sij
real*8 ai,aj,rhoa0j,rhoa1j,rhoa2j,rhoa3j,A1j,A2j,A3j
real*8 G,Gbar,gam,shp(3)
real*8 ro0i,ro0j
real*8 rhoa0i,rhoa1i,rhoa2i,rhoa3i,A1i,A2i,A3i
elti = fmap(type(i))
xtmp = x(1,i)
ytmp = x(2,i)
ztmp = x(3,i)
do jn = 1,numneigh
if (scrfcn(jn).ne.0.d0) then
j = firstneigh(jn)
sij = scrfcn(jn)*fcpair(jn)
delij(1) = x(1,j) - xtmp
delij(2) = x(2,j) - ytmp
delij(3) = x(3,j) - ztmp
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
$ + delij(3)*delij(3)
if (rij2.lt.cutforcesq) then
eltj = fmap(type(j))
rij = sqrt(rij2)
ai = rij/re_meam(elti,elti) - 1.d0
aj = rij/re_meam(eltj,eltj) - 1.d0
ro0i = rho0_meam(elti)
ro0j = rho0_meam(eltj)
rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj)*sij
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)*sij
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)*sij
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)*sij
rhoa0i = ro0i*exp(-beta0_meam(elti)*ai)*sij
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)*sij
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)*sij
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)*sij
if (ialloy.eq.1) then
rhoa1j = rhoa1j * t1_meam(eltj)
rhoa2j = rhoa2j * t2_meam(eltj)
rhoa3j = rhoa3j * t3_meam(eltj)
rhoa1i = rhoa1i * t1_meam(elti)
rhoa2i = rhoa2i * t2_meam(elti)
rhoa3i = rhoa3i * t3_meam(elti)
endif
rho0(i) = rho0(i) + rhoa0j
rho0(j) = rho0(j) + rhoa0i
- t_ave(1,i) = t_ave(1,i) + t1_meam(eltj)*rhoa0j
- t_ave(2,i) = t_ave(2,i) + t2_meam(eltj)*rhoa0j
- t_ave(3,i) = t_ave(3,i) + t3_meam(eltj)*rhoa0j
- t_ave(1,j) = t_ave(1,j) + t1_meam(elti)*rhoa0i
- t_ave(2,j) = t_ave(2,j) + t2_meam(elti)*rhoa0i
- t_ave(3,j) = t_ave(3,j) + t3_meam(elti)*rhoa0i
+c For ialloy = 2, use single-element value (not average)
+ if (ialloy.ne.2) then
+ t_ave(1,i) = t_ave(1,i) + t1_meam(eltj)*rhoa0j
+ t_ave(2,i) = t_ave(2,i) + t2_meam(eltj)*rhoa0j
+ t_ave(3,i) = t_ave(3,i) + t3_meam(eltj)*rhoa0j
+ t_ave(1,j) = t_ave(1,j) + t1_meam(elti)*rhoa0i
+ t_ave(2,j) = t_ave(2,j) + t2_meam(elti)*rhoa0i
+ t_ave(3,j) = t_ave(3,j) + t3_meam(elti)*rhoa0i
+ endif
if (ialloy.eq.1) then
tsq_ave(1,i) = tsq_ave(1,i) +
$ t1_meam(eltj)*t1_meam(eltj)*rhoa0j
tsq_ave(2,i) = tsq_ave(2,i) +
$ t2_meam(eltj)*t2_meam(eltj)*rhoa0j
tsq_ave(3,i) = tsq_ave(3,i) +
$ t3_meam(eltj)*t3_meam(eltj)*rhoa0j
tsq_ave(1,j) = tsq_ave(1,j) +
$ t1_meam(elti)*t1_meam(elti)*rhoa0i
tsq_ave(2,j) = tsq_ave(2,j) +
$ t2_meam(elti)*t2_meam(elti)*rhoa0i
tsq_ave(3,j) = tsq_ave(3,j) +
$ t3_meam(elti)*t3_meam(elti)*rhoa0i
endif
Arho2b(i) = Arho2b(i) + rhoa2j
Arho2b(j) = Arho2b(j) + rhoa2i
A1j = rhoa1j/rij
A2j = rhoa2j/rij2
A3j = rhoa3j/(rij2*rij)
A1i = rhoa1i/rij
A2i = rhoa2i/rij2
A3i = rhoa3i/(rij2*rij)
nv2 = 1
nv3 = 1
do m = 1,3
Arho1(m,i) = Arho1(m,i) + A1j*delij(m)
Arho1(m,j) = Arho1(m,j) - A1i*delij(m)
Arho3b(m,i) = Arho3b(m,i) + rhoa3j*delij(m)/rij
Arho3b(m,j) = Arho3b(m,j) - rhoa3i*delij(m)/rij
do n = m,3
Arho2(nv2,i) = Arho2(nv2,i) + A2j*delij(m)*delij(n)
Arho2(nv2,j) = Arho2(nv2,j) + A2i*delij(m)*delij(n)
nv2 = nv2+1
do p = n,3
Arho3(nv3,i) = Arho3(nv3,i)
$ + A3j*delij(m)*delij(n)*delij(p)
Arho3(nv3,j) = Arho3(nv3,j)
$ - A3i*delij(m)*delij(n)*delij(p)
nv3 = nv3+1
enddo
enddo
enddo
endif
endif
enddo
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine screen(i, j, nmax, x, rijsq, sij,
$ numneigh_full, firstneigh_full, ntype, type, fmap)
c Screening function
c Inputs: i = atom 1 id (integer)
c j = atom 2 id (integer)
c rijsq = squared distance between i and j
c Outputs: sij = screening function
use meam_data
implicit none
integer i,j,nmax,k,nk,m
real*8 x,rijsq,sij
integer numneigh_full, firstneigh_full
integer ntype, type, fmap
dimension x(3,nmax), firstneigh_full(numneigh_full)
dimension type(nmax), fmap(ntype)
integer elti,eltj,eltk
real*8 delxik,delyik,delzik
real*8 delxjk,delyjk,delzjk
real*8 riksq,rjksq,xik,xjk,cikj,a,delc,sikj,fcij,rij
real*8 Cmax,Cmin,rbound
sij = 1.d0
elti = fmap(type(i))
eltj = fmap(type(j))
c if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
rbound = ebound_meam(elti,eltj)*rijsq
do nk = 1,numneigh_full
k = firstneigh_full(nk)
eltk = fmap(type(k))
if (k.eq.j) goto 10
delxjk = x(1,k) - x(1,j)
delyjk = x(2,k) - x(2,j)
delzjk = x(3,k) - x(3,j)
rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk
if (rjksq.gt.rbound) goto 10
delxik = x(1,k) - x(1,i)
delyik = x(2,k) - x(2,i)
delzik = x(3,k) - x(3,i)
riksq = delxik*delxik + delyik*delyik + delzik*delzik
if (riksq.gt.rbound) goto 10
xik = riksq/rijsq
xjk = rjksq/rijsq
a = 1 - (xik-xjk)*(xik-xjk)
c if a < 0, then ellipse equation doesn't describe this case and
c atom k can't possibly screen i-j
if (a.le.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) then
goto 10
c note that cikj may be slightly negative (within numerical
c tolerance) if atoms are colinear, so don't reject that case here
c (other negative cikj cases were handled by the test on "a" above)
else if (cikj.le.Cmin) then
sij = 0.d0
goto 20
else
delc = Cmax - Cmin
cikj = (cikj-Cmin)/delc
call fcut(cikj,sikj)
endif
sij = sij * sikj
10 continue
enddo
20 continue
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
$ ntype,type,fmap,x,scrfcn,fcpair)
c Inputs: i,j,k = id's of 3 atom triplet
c jn = id of i-j pair
c rij2 = squared distance between i and j
c Outputs: dsij1 = deriv. of sij w.r.t. rik
c dsij2 = deriv. of sij w.r.t. rjk
use meam_data
implicit none
integer i,j,k,jn,nmax,numneigh
integer elti,eltj,eltk
real*8 rij2,rik2,rjk2,dsij1,dsij2
integer ntype, type, fmap
real*8 x, scrfcn, fcpair
dimension type(nmax), fmap(ntype)
dimension x(3,nmax), scrfcn(numneigh), fcpair(numneigh)
real*8 dxik,dyik,dzik
real*8 dxjk,dyjk,dzjk
real*8 rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a
real*8 Cmax,Cmin,dCikj1,dCikj2
sij = scrfcn(jn)*fcpair(jn)
elti = fmap(type(i))
eltj = fmap(type(j))
eltk = fmap(type(k))
Cmax = Cmax_meam(elti,eltj,eltk)
Cmin = Cmin_meam(elti,eltj,eltk)
dsij1 = 0.d0
dsij2 = 0.d0
if ((sij.ne.0.d0).and.(sij.ne.1.d0)) then
rbound = rij2*ebound_meam(elti,eltj)
delc = Cmax-Cmin
dxjk = x(1,k) - x(1,j)
dyjk = x(2,k) - x(2,j)
dzjk = x(3,k) - x(3,j)
rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk
if (rjk2.le.rbound) then
dxik = x(1,k) - x(1,i)
dyik = x(2,k) - x(2,i)
dzik = x(3,k) - x(3,i)
rik2 = dxik*dxik + dyik*dyik + dzik*dzik
if (rik2.le.rbound) then
xik = rik2/rij2
xjk = rjk2/rij2
a = 1 - (xik-xjk)*(xik-xjk)
if (a.ne.0.d0) then
cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
if (cikj.ge.Cmin.and.cikj.le.Cmax) then
cikj = (cikj-Cmin)/delc
call dfcut(cikj,sikj,dfc)
call dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2)
a = sij/delc*dfc/sikj
dsij1 = a*dCikj1
dsij2 = a*dCikj2
endif
endif
endif
endif
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine fcut(xi,fc)
c cutoff function
implicit none
real*8 xi,fc
real*8 a
if (xi.ge.1.d0) then
fc = 1.d0
else if (xi.le.0.d0) then
fc = 0.d0
else
a = 1.d0-xi
a = a*a
a = a*a
a = 1.d0-a
fc = a*a
c fc = xi
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine dfcut(xi,fc,dfc)
c cutoff function and its derivative
implicit none
real*8 xi,fc,dfc,a,a3,a4
if (xi.ge.1.d0) then
fc = 1.d0
dfc = 0.d0
else if (xi.le.0.d0) then
fc = 0.d0
dfc = 0.d0
else
a = 1.d0-xi
a3 = a*a*a
a4 = a*a3
fc = (1.d0-a4)**2
dfc = 8*(1.d0-a4)*a3
c fc = xi
c dfc = 1.d0
endif
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine dCfunc(rij2,rik2,rjk2,dCikj)
c Inputs: rij,rij2,rik2,rjk2
c Outputs: dCikj = derivative of Cikj w.r.t. rij
implicit none
real*8 rij2,rik2,rjk2,dCikj
real*8 rij4,a,b,denom
rij4 = rij2*rij2
a = rik2-rjk2
b = rik2+rjk2
denom = rij4 - a*a
denom = denom*denom
dCikj = -4*(-2*rij2*a*a + rij4*b + a*a*b)/denom
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2)
c Inputs: rij,rij2,rik2,rjk2
c Outputs: dCikj1 = derivative of Cikj w.r.t. rik
c dCikj2 = derivative of Cikj w.r.t. rjk
implicit none
real*8 rij2,rik2,rjk2,dCikj1,dCikj2
real*8 rij4,rik4,rjk4,a,b,denom
rij4 = rij2*rij2
rik4 = rik2*rik2
rjk4 = rjk2*rjk2
a = rik2-rjk2
b = rik2+rjk2
denom = rij4 - a*a
denom = denom*denom
dCikj1 = 4*rij2*(rij4 + rik4 + 2*rik2*rjk2 - 3*rjk4 - 2*rij2*a)/
$ denom
dCikj2 = 4*rij2*(rij4 - 3*rik4 + 2*rik2*rjk2 + rjk4 + 2*rij2*a)/
$ denom
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
diff --git a/lib/meam/meam_force.F b/lib/meam/meam_force.F
index 05329a4b4..784750298 100755
--- a/lib/meam/meam_force.F
+++ b/lib/meam/meam_force.F
@@ -1,590 +1,608 @@
c Extern "C" declaration has the form:
c
c void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *,
c int *, int *, int *, int *, double *, double *,
c double *, double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *, double *, int *);
c
c Call from pair_meam.cpp has the form:
c
c meam_force_(&i,&nmax,&eflag_either,&eflag_global,&eflag_atom,&vflag_atom,
c &eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i],
c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag);
c
subroutine meam_force(i, nmax,
$ eflag_either, eflag_global, eflag_atom, vflag_atom,
$ eng_vdwl, eatom, ntype, type, fmap, x,
$ numneigh, firstneigh, numneigh_full, firstneigh_full,
$ scrfcn, dscrfcn, fcpair,
$ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp,
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f,
$ vatom, errorflag)
use meam_data
implicit none
integer eflag_either, eflag_global, eflag_atom, vflag_atom
integer nmax, ntype, type, fmap
real*8 eng_vdwl, eatom, x
integer numneigh, firstneigh, numneigh_full, firstneigh_full
real*8 scrfcn, dscrfcn, fcpair
real*8 dGamma1, dGamma2, dGamma3
real*8 rho0, rho1, rho2, rho3, fp
real*8 Arho1, Arho2, Arho2b
real*8 Arho3, Arho3b
real*8 t_ave, tsq_ave, f, vatom
integer errorflag
dimension eatom(nmax)
dimension type(nmax), fmap(ntype)
dimension x(3,nmax)
dimension firstneigh(numneigh), firstneigh_full(numneigh_full)
dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh)
dimension dGamma1(nmax), dGamma2(nmax), dGamma3(nmax)
dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax)
dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax)
dimension Arho3(10,nmax), Arho3b(3,nmax)
dimension t_ave(3,nmax), tsq_ave(3,nmax), f(3,nmax), vatom(6,nmax)
integer i,j,jn,k,kn,kk,m,n,p,q
integer nv2,nv3,elti,eltj,eltk,ind
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 Eu,astar,astarp,third,sixth
real*8 pp,phiforce,dUdrij,dUdsij,dUdrijm(3),force,forcem
real*8 B,r,recip,phi,phip,rhop,a
real*8 sij,fcij,dfcij,ds(3)
real*8 a0,a1,a1i,a1j,a2,a2i,a2j
real*8 a3i,a3j,a3i1,a3i2,a3j1,a3j2
real*8 G,dG,Gbar,dGbar,gam,shpi(3),shpj(3),Z,denom
real*8 ai,aj,ro0i,ro0j,invrei,invrej
real*8 b0,rhoa0j,drhoa0j,rhoa0i,drhoa0i
real*8 b1,rhoa1j,drhoa1j,rhoa1i,drhoa1i
real*8 b2,rhoa2j,drhoa2j,rhoa2i,drhoa2i
real*8 a3,a3a,b3,rhoa3j,drhoa3j,rhoa3i,drhoa3i
real*8 drho0dr1,drho0dr2,drho0ds1,drho0ds2
real*8 drho1dr1,drho1dr2,drho1ds1,drho1ds2
real*8 drho1drm1(3),drho1drm2(3)
real*8 drho2dr1,drho2dr2,drho2ds1,drho2ds2
real*8 drho2drm1(3),drho2drm2(3)
real*8 drho3dr1,drho3dr2,drho3ds1,drho3ds2
real*8 drho3drm1(3),drho3drm2(3)
real*8 dt1dr1,dt1dr2,dt1ds1,dt1ds2
real*8 dt2dr1,dt2dr2,dt2ds1,dt2ds2
real*8 dt3dr1,dt3dr2,dt3ds1,dt3ds2
real*8 drhodr1,drhodr2,drhods1,drhods2,drhodrm1(3),drhodrm2(3)
real*8 arg,arg1,arg2
real*8 arg1i1,arg1j1,arg1i2,arg1j2,arg2i2,arg2j2
real*8 arg1i3,arg1j3,arg2i3,arg2j3,arg3i3,arg3j3
real*8 dsij1,dsij2,force1,force2
real*8 t1i,t2i,t3i,t1j,t2j,t3j
errorflag = 0
third = 1.0/3.0
sixth = 1.0/6.0
c Compute forces atom i
elti = fmap(type(i))
if (elti.gt.0) then
xitmp = x(1,i)
yitmp = x(2,i)
zitmp = x(3,i)
c Treat each pair
do jn = 1,numneigh
j = firstneigh(jn)
eltj = fmap(type(j))
if (scrfcn(jn).ne.0.d0.and.eltj.gt.0) then
sij = scrfcn(jn)*fcpair(jn)
delij(1) = x(1,j) - xitmp
delij(2) = x(2,j) - yitmp
delij(3) = x(3,j) - zitmp
rij2 = delij(1)*delij(1) + delij(2)*delij(2)
$ + delij(3)*delij(3)
if (rij2.lt.cutforcesq) then
rij = sqrt(rij2)
r = rij
c Compute phi and phip
ind = eltind(elti,eltj)
pp = rij*rdrar + 1.0D0
kk = pp
kk = min(kk,nrar-1)
pp = pp - kk
pp = min(pp,1.0D0)
phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
$ + phirar1(kk,ind))*pp + phirar(kk,ind)
phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp
$ + phirar4(kk,ind)
recip = 1.0d0/r
if (eflag_either.ne.0) then
if (eflag_global.ne.0) eng_vdwl = eng_vdwl + phi*sij
if (eflag_atom.ne.0) then
eatom(i) = eatom(i) + 0.5*phi*sij
eatom(j) = eatom(j) + 0.5*phi*sij
endif
endif
c write(1,*) "force_meamf: phi: ",phi
c write(1,*) "force_meamf: phip: ",phip
c Compute pair densities and derivatives
invrei = 1.d0/re_meam(elti,elti)
ai = rij*invrei - 1.d0
ro0i = rho0_meam(elti)
rhoa0i = ro0i*exp(-beta0_meam(elti)*ai)
drhoa0i = -beta0_meam(elti)*invrei*rhoa0i
rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)
drhoa1i = -beta1_meam(elti)*invrei*rhoa1i
rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)
drhoa2i = -beta2_meam(elti)*invrei*rhoa2i
rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)
drhoa3i = -beta3_meam(elti)*invrei*rhoa3i
if (elti.ne.eltj) then
invrej = 1.d0/re_meam(eltj,eltj)
aj = rij*invrej - 1.d0
ro0j = rho0_meam(eltj)
rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj)
drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j
rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj)
drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j
rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj)
drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j
rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj)
drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j
else
rhoa0j = rhoa0i
drhoa0j = drhoa0i
rhoa1j = rhoa1i
drhoa1j = drhoa1i
rhoa2j = rhoa2i
drhoa2j = drhoa2i
rhoa3j = rhoa3i
drhoa3j = drhoa3i
endif
if (ialloy.eq.1) then
rhoa1j = rhoa1j * t1_meam(eltj)
rhoa2j = rhoa2j * t2_meam(eltj)
rhoa3j = rhoa3j * t3_meam(eltj)
rhoa1i = rhoa1i * t1_meam(elti)
rhoa2i = rhoa2i * t2_meam(elti)
rhoa3i = rhoa3i * t3_meam(elti)
drhoa1j = drhoa1j * t1_meam(eltj)
drhoa2j = drhoa2j * t2_meam(eltj)
drhoa3j = drhoa3j * t3_meam(eltj)
drhoa1i = drhoa1i * t1_meam(elti)
drhoa2i = drhoa2i * t2_meam(elti)
drhoa3i = drhoa3i * t3_meam(elti)
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
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
c Compute derivatives of weighting functions t wrt rij
t1i = t_ave(1,i)
t2i = t_ave(2,i)
t3i = t_ave(3,i)
t1j = t_ave(1,j)
t2j = t_ave(2,j)
t3j = t_ave(3,j)
if (ialloy.eq.1) then
a1i = 0.d0
a1j = 0.d0
a2i = 0.d0
a2j = 0.d0
a3i = 0.d0
a3j = 0.d0
if ( tsq_ave(1,i) .ne. 0.d0 ) then
a1i = drhoa0j*sij/tsq_ave(1,i)
endif
if ( tsq_ave(1,j) .ne. 0.d0 ) then
a1j = drhoa0i*sij/tsq_ave(1,j)
endif
if ( tsq_ave(2,i) .ne. 0.d0 ) then
a2i = drhoa0j*sij/tsq_ave(2,i)
endif
if ( tsq_ave(2,j) .ne. 0.d0 ) then
a2j = drhoa0i*sij/tsq_ave(2,j)
endif
if ( tsq_ave(3,i) .ne. 0.d0 ) then
a3i = drhoa0j*sij/tsq_ave(3,i)
endif
if ( tsq_ave(3,j) .ne. 0.d0 ) then
a3j = drhoa0i*sij/tsq_ave(3,j)
endif
dt1dr1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2)
dt1dr2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2)
dt2dr1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2)
dt2dr2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2)
dt3dr1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2)
dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2)
+ else if (ialloy.eq.2) then
+
+ dt1dr1 = 0.d0
+ dt1dr2 = 0.d0
+ dt2dr1 = 0.d0
+ dt2dr2 = 0.d0
+ dt3dr1 = 0.d0
+ dt3dr2 = 0.d0
+
else
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)
endif
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
if (ialloy.eq.1) then
a1i = 0.d0
a1j = 0.d0
a2i = 0.d0
a2j = 0.d0
a3i = 0.d0
a3j = 0.d0
if ( tsq_ave(1,i) .ne. 0.d0 ) then
a1i = rhoa0j/tsq_ave(1,i)
endif
if ( tsq_ave(1,j) .ne. 0.d0 ) then
a1j = rhoa0i/tsq_ave(1,j)
endif
if ( tsq_ave(2,i) .ne. 0.d0 ) then
a2i = rhoa0j/tsq_ave(2,i)
endif
if ( tsq_ave(2,j) .ne. 0.d0 ) then
a2j = rhoa0i/tsq_ave(2,j)
endif
if ( tsq_ave(3,i) .ne. 0.d0 ) then
a3i = rhoa0j/tsq_ave(3,i)
endif
if ( tsq_ave(3,j) .ne. 0.d0 ) then
a3j = rhoa0i/tsq_ave(3,j)
endif
dt1ds1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2)
dt1ds2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2)
dt2ds1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2)
dt2ds2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2)
dt3ds1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2)
dt3ds2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2)
+ else if (ialloy.eq.2) then
+
+ dt1ds1 = 0.d0
+ dt1ds2 = 0.d0
+ dt2ds1 = 0.d0
+ dt2ds2 = 0.d0
+ dt3ds1 = 0.d0
+ dt3ds2 = 0.d0
+
else
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)
endif
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)
dUdrij = phip*sij
$ + fp(i)*drhodr1 + fp(j)*drhodr2
dUdsij = 0.d0
if (dscrfcn(jn).ne.0.d0) then
dUdsij = phi
$ + fp(i)*drhods1 + fp(j)*drhods2
endif
do m = 1,3
dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m)
enddo
c Add the part of the force due to dUdrij and dUdsij
force = dUdrij*recip + dUdsij*dscrfcn(jn)
do m = 1,3
forcem = delij(m)*force + dUdrijm(m)
f(m,i) = f(m,i) + forcem
f(m,j) = f(m,j) - forcem
enddo
c Tabulate per-atom virial as symmetrized stress tensor
if (vflag_atom.ne.0) then
fi(1) = delij(1)*force + dUdrijm(1)
fi(2) = delij(2)*force + dUdrijm(2)
fi(3) = delij(3)*force + dUdrijm(3)
v(1) = -0.5 * (delij(1) * fi(1))
v(2) = -0.5 * (delij(2) * fi(2))
v(3) = -0.5 * (delij(3) * fi(3))
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(6) = -0.25 * (delij(2)*fi(3) + delij(3)*fi(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)
endif
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
do kn = 1,numneigh_full
k = firstneigh_full(kn)
eltk = fmap(type(k))
if (k.ne.j.and.eltk.gt.0) then
call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2,
$ ntype,type,fmap,x,scrfcn,fcpair)
if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then
force1 = dUdsij*dsij1
force2 = dUdsij*dsij2
do m = 1,3
delik(m) = x(m,k) - x(m,i)
deljk(m) = x(m,k) - x(m,j)
enddo
do m = 1,3
f(m,i) = f(m,i) + force1*delik(m)
f(m,j) = f(m,j) + force2*deljk(m)
f(m,k) = f(m,k) - force1*delik(m)
$ - force2*deljk(m)
enddo
c Tabulate per-atom virial as symmetrized stress tensor
if (vflag_atom.ne.0) then
fi(1) = force1*delik(1)
fi(2) = force1*delik(2)
fi(3) = force1*delik(3)
fj(1) = force2*deljk(1)
fj(2) = force2*deljk(2)
fj(3) = force2*deljk(3)
v(1) = -third * (delik(1)*fi(1) + deljk(1)*fj(1))
v(2) = -third * (delik(2)*fi(2) + deljk(2)*fj(2))
v(3) = -third * (delik(3)*fi(3) + deljk(3)*fj(3))
v(4) = -sixth * (delik(1)*fi(2) + deljk(1)*fj(2) +
$ delik(2)*fi(1) + deljk(2)*fj(1))
v(5) = -sixth * (delik(1)*fi(3) + deljk(1)*fj(3) +
$ delik(3)*fi(1) + deljk(3)*fj(1))
v(6) = -sixth * (delik(2)*fi(3) + deljk(2)*fj(3) +
$ 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
endif
endif
c end of k loop
enddo
endif
100 continue
endif
c end of j loop
enddo
c else if elti=0, this is not a meam atom
endif
return
end
diff --git a/lib/meam/meam_setup_done.F b/lib/meam/meam_setup_done.F
index a4373e9e5..597d9a6f2 100755
--- a/lib/meam/meam_setup_done.F
+++ b/lib/meam/meam_setup_done.F
@@ -1,917 +1,1005 @@
c Declaration in pair_meam.h:
c
c void meam_setup_done(double *)
c
c Call from pair_meam.cpp:
c
c meam_setup_done(&cutmax)
c
subroutine meam_setup_done(cutmax)
use meam_data
implicit none
real*8 cutmax
integer nv2, nv3, m, n, p
c Force cutoff
cutforce = rc_meam
cutforcesq = cutforce*cutforce
c Pass cutoff back to calling program
cutmax = cutforce
c Augment t1 term
t1_meam(:) = t1_meam(:) + augt1 * 3.d0/5.d0 * t3_meam(:)
c Compute off-diagonal alloy parameters
call alloyparams
c indices and factors for Voight notation
nv2 = 1
nv3 = 1
do m = 1,3
do n = m,3
vind2D(m,n) = nv2
vind2D(n,m) = nv2
nv2 = nv2+1
do p = n,3
vind3D(m,n,p) = nv3
vind3D(m,p,n) = nv3
vind3D(n,m,p) = nv3
vind3D(n,p,m) = nv3
vind3D(p,m,n) = nv3
vind3D(p,n,m) = nv3
nv3 = nv3+1
enddo
enddo
enddo
v2D(1) = 1
v2D(2) = 2
v2D(3) = 2
v2D(4) = 1
v2D(5) = 2
v2D(6) = 1
v3D(1) = 1
v3D(2) = 3
v3D(3) = 3
v3D(4) = 3
v3D(5) = 6
v3D(6) = 3
v3D(7) = 1
v3D(8) = 3
v3D(9) = 3
v3D(10) = 1
nv2 = 1
do m = 1,neltypes
do n = m,neltypes
eltind(m,n) = nv2
eltind(n,m) = nv2
nv2 = nv2+1
enddo
enddo
+c Compute background densities for reference structure
+ call compute_reference_density
+
c Compute pair potentials and setup arrays for interpolation
nr = 1000
dr = 1.1*rc_meam/nr
call compute_pair_meam
return
end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Fill off-diagonal alloy parameters
subroutine alloyparams
use meam_data
implicit none
integer i,j,k
real*8 eb
c Loop over pairs
do i = 1,neltypes
do j = 1,neltypes
c Treat off-diagonal pairs
c If i>j, set all equal to i<j case (which has aready been set,
c here or in the input file)
if (i.gt.j) then
re_meam(i,j) = re_meam(j,i)
Ec_meam(i,j) = Ec_meam(j,i)
alpha_meam(i,j) = alpha_meam(j,i)
lattce_meam(i,j) = lattce_meam(j,i)
c If i<j and term is unset, use default values (e.g. mean of i-i and j-j)
else if (j.gt.i) then
if (Ec_meam(i,j).eq.0.d0) then
if (lattce_meam(i,j).eq.'l12') then
Ec_meam(i,j) = (3*Ec_meam(i,i)+Ec_meam(j,j))/4.d0
$ - delta_meam(i,j)
else if (lattce_meam(i,j).eq.'c11') then
if (lattce_meam(i,i).eq.'dia') then
Ec_meam(i,j) = (2*Ec_meam(i,i)+Ec_meam(j,j))/3.d0
$ - delta_meam(i,j)
else
Ec_meam(i,j) = (Ec_meam(i,i)+2*Ec_meam(j,j))/3.d0
$ - delta_meam(i,j)
endif
else
Ec_meam(i,j) = (Ec_meam(i,i)+Ec_meam(j,j))/2.d0
$ - delta_meam(i,j)
endif
endif
if (alpha_meam(i,j).eq.0.d0) then
alpha_meam(i,j) = (alpha_meam(i,i)+alpha_meam(j,j))/2.d0
endif
if (re_meam(i,j).eq.0.d0) then
re_meam(i,j) = (re_meam(i,i)+re_meam(j,j))/2.d0
endif
endif
enddo
enddo
c Cmin(i,k,j) is symmetric in i-j, but not k. For all triplets
c where i>j, set equal to the i<j element. Likewise for Cmax.
do i = 2,neltypes
do j = 1,i-1
do k = 1,neltypes
Cmin_meam(i,j,k) = Cmin_meam(j,i,k)
Cmax_meam(i,j,k) = Cmax_meam(j,i,k)
enddo
enddo
enddo
c ebound gives the squared distance such that, for rik2 or rjk2>ebound,
c atom k definitely lies outside the screening function ellipse (so
c there is no need to calculate its effects). Here, compute it for all
c triplets (i,j,k) so that ebound(i,j) is the maximized over k
do i = 1,neltypes
do j = 1,neltypes
do k = 1,neltypes
eb = (Cmax_meam(i,j,k)*Cmax_meam(i,j,k))
$ /(4.d0*(Cmax_meam(i,j,k)-1.d0))
ebound_meam(i,j) = max(ebound_meam(i,j),eb)
enddo
enddo
enddo
return
end
c-----------------------------------------------------------------------
c compute MEAM pair potential for each pair of element types
c
subroutine compute_pair_meam
use meam_data
implicit none
real*8 r, temp
integer j,a,b,nv2
real*8 astar,frac,phizbl
integer n,nmax,Z1,Z2
real*8 arat,rarat,scrn,scrn2
real*8 phiaa,phibb,phitmp
real*8 C,s111,s112,s221,S11,S22
real*8, external :: phi_meam
real*8, external :: zbl
real*8, external :: compute_phi
c allocate memory for array that defines the potential
allocate(phir(nr,(neltypes*(neltypes+1))/2))
c allocate coeff memory
allocate(phirar(nr,(neltypes*(neltypes+1))/2))
allocate(phirar1(nr,(neltypes*(neltypes+1))/2))
allocate(phirar2(nr,(neltypes*(neltypes+1))/2))
allocate(phirar3(nr,(neltypes*(neltypes+1))/2))
allocate(phirar4(nr,(neltypes*(neltypes+1))/2))
allocate(phirar5(nr,(neltypes*(neltypes+1))/2))
allocate(phirar6(nr,(neltypes*(neltypes+1))/2))
c loop over pairs of element types
nv2 = 0
do a = 1,neltypes
do b = a,neltypes
nv2 = nv2 + 1
c loop over r values and compute
do j = 1,nr
r = (j-1)*dr
phir(j,nv2) = phi_meam(r,a,b)
c if using second-nearest neighbor, solve recursive problem
c (see Lee and Baskes, PRB 62(13):8564 eqn.(21))
if (nn2_meam(a,b).eq.1) then
-c if (a.ne.b) then
-c write(6,*) 'Second-neighbor currently only valid '
-c write(6,*) 'if element types are the same.'
-c call error(' ')
-c endif
call get_Zij(Z1,lattce_meam(a,b))
call get_Zij2(Z2,arat,scrn,lattce_meam(a,b),
$ Cmin_meam(a,a,b),Cmax_meam(a,a,b))
-c The B1 and L12 cases with NN2 have a trick to them; we need to
+c The B1, B2, and L12 cases with NN2 have a trick to them; we need to
c compute the contributions from second nearest neighbors, like a-a
c pairs, but need to include NN2 contributions to those pairs as
c well.
if (lattce_meam(a,b).eq.'b1'.or.
+ $ lattce_meam(a,b).eq.'b2'.or.
$ lattce_meam(a,b).eq.'l12') then
rarat = r*arat
c phi_aa
phiaa = phi_meam(rarat,a,a)
call get_Zij(Z1,lattce_meam(a,a))
call get_Zij2(Z2,arat,scrn,lattce_meam(a,a),
$ Cmin_meam(a,a,a),Cmax_meam(a,a,a))
nmax = 10
if (scrn.gt.0.0) then
do n = 1,nmax
phiaa = phiaa +
$ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,a,a)
enddo
endif
c phi_bb
phibb = phi_meam(rarat,b,b)
call get_Zij(Z1,lattce_meam(b,b))
call get_Zij2(Z2,arat,scrn,lattce_meam(b,b),
$ Cmin_meam(b,b,b),Cmax_meam(b,b,b))
nmax = 10
if (scrn.gt.0.0) then
do n = 1,nmax
phibb = phibb +
$ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,b,b)
enddo
endif
- if (lattce_meam(a,b).eq.'b1') then
-c Add contributions to the B1 potential
+ if (lattce_meam(a,b).eq.'b1'.
+ $ or.lattce_meam(a,b).eq.'b2') then
+c Add contributions to the B1 or B2 potential
call get_Zij(Z1,lattce_meam(a,b))
call get_Zij2(Z2,arat,scrn,lattce_meam(a,b),
$ Cmin_meam(a,a,b),Cmax_meam(a,a,b))
phir(j,nv2) = phir(j,nv2) -
$ Z2*scrn/(2*Z1) * phiaa
call get_Zij2(Z2,arat,scrn2,lattce_meam(a,b),
$ Cmin_meam(b,b,a),Cmax_meam(b,b,a))
phir(j,nv2) = phir(j,nv2) -
$ Z2*scrn2/(2*Z1) * phibb
else if (lattce_meam(a,b).eq.'l12') then
c The L12 case has one last trick; we have to be careful to compute
c the correct screening between 2nd-neighbor pairs. 1-1
c second-neighbor pairs are screened by 2 type 1 atoms and two type
c 2 atoms. 2-2 second-neighbor pairs are screened by 4 type 1
c atoms.
C = 1.d0
call get_sijk(C,a,a,a,s111)
call get_sijk(C,a,a,b,s112)
call get_sijk(C,b,b,a,s221)
S11 = s111 * s111 * s112 * s112
S22 = s221**4
phir(j,nv2) = phir(j,nv2) -
$ 0.75*S11*phiaa - 0.25*S22*phibb
endif
else
nmax = 10
do n = 1,nmax
phir(j,nv2) = phir(j,nv2) +
$ (-Z2*scrn/Z1)**n * phi_meam(r*arat**n,a,b)
enddo
endif
endif
+c For Zbl potential:
c if astar <= -3
c potential is zbl potential
c else if -3 < astar < -1
c potential is linear combination with zbl potential
c endif
- astar = alpha_meam(a,b) * (r/re_meam(a,b) - 1.d0)
- if (astar.le.-3.d0) then
- phir(j,nv2) = zbl(r,ielt_meam(a),ielt_meam(b))
- else if (astar.gt.-3.d0.and.astar.lt.-1.d0) then
- call fcut(1-(astar+1.d0)/(-3.d0+1.d0),frac)
- phizbl = zbl(r,ielt_meam(a),ielt_meam(b))
- phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl
+ if (zbl_meam(a,b).eq.1) then
+ astar = alpha_meam(a,b) * (r/re_meam(a,b) - 1.d0)
+ if (astar.le.-3.d0) then
+ phir(j,nv2) = zbl(r,ielt_meam(a),ielt_meam(b))
+ else if (astar.gt.-3.d0.and.astar.lt.-1.d0) then
+ call fcut(1-(astar+1.d0)/(-3.d0+1.d0),frac)
+ phizbl = zbl(r,ielt_meam(a),ielt_meam(b))
+ phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl
+ endif
endif
-
+
enddo
c call interpolation
call interpolate_meam(nv2)
-
+
enddo
enddo
return
end
c----------------------------------------------------------------------c
c Compute MEAM pair potential for distance r, element types a and b
c
real*8 recursive function phi_meam(r,a,b)result(phi_m)
use meam_data
implicit none
integer a,b
real*8 r
real*8 a1,a2,a12
real*8 t11av,t21av,t31av,t12av,t22av,t32av
real*8 G1,G2,s1(3),s2(3),s12(3),rho0_1,rho0_2
real*8 Gam1,Gam2,Z1,Z2
real*8 rhobar1,rhobar2,F1,F2
real*8 rhoa01,rhoa11,rhoa21,rhoa31
real*8 rhoa02,rhoa12,rhoa22,rhoa32
real*8 rho01,rho11,rho21,rho31
real*8 rho02,rho12,rho22,rho32
real*8 scalfac,phiaa,phibb
real*8 Eu
real*8 arat,scrn,scrn2
integer Z12, errorflag
integer n,nmax,Z1nn,Z2nn
character*3 latta,lattb
+ real*8 rho_bkgd1, rho_bkgd2
real*8, external :: erose
c Equation numbers below refer to:
c I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615
c get number of neighbors in the reference structure
c Nref(i,j) = # of i's neighbors of type j
call get_Zij(Z12,lattce_meam(a,b))
call get_densref(r,a,b,rho01,rho11,rho21,rho31,
$ rho02,rho12,rho22,rho32)
c if densities are too small, numerical problems may result; just return zero
if (rho01.le.1e-14.and.rho02.le.1e-14) then
phi_m = 0.0
return
endif
c calculate average weighting factors for the reference structure
if (lattce_meam(a,b).eq.'c11') then
- scalfac = 1.0/(rho01+rho02)
- t11av = scalfac*(t1_meam(a)*rho01 + t1_meam(b)*rho02)
- t12av = t11av
- t21av = scalfac*(t2_meam(a)*rho01 + t2_meam(b)*rho02)
- t22av = t21av
- t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02)
- t32av = t31av
+ if (ialloy.eq.2) then
+ t11av = t1_meam(a)
+ t12av = t1_meam(b)
+ t21av = t2_meam(a)
+ t22av = t2_meam(b)
+ t31av = t3_meam(a)
+ t32av = t3_meam(b)
+ else
+ scalfac = 1.0/(rho01+rho02)
+ t11av = scalfac*(t1_meam(a)*rho01 + t1_meam(b)*rho02)
+ t12av = t11av
+ t21av = scalfac*(t2_meam(a)*rho01 + t2_meam(b)*rho02)
+ t22av = t21av
+ t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02)
+ t32av = t31av
+ endif
else
c average weighting factors for the reference structure, eqn. I.8
call get_tavref(t11av,t21av,t31av,t12av,t22av,t32av,
$ t1_meam(a),t2_meam(a),t3_meam(a),
$ t1_meam(b),t2_meam(b),t3_meam(b),
$ r,a,b,lattce_meam(a,b))
endif
c for c11b structure, calculate background electron densities
if (lattce_meam(a,b).eq.'c11') then
latta = lattce_meam(a,a)
if (latta.eq.'dia') then
rhobar1 = ((Z12/2)*(rho02+rho01))**2 +
$ t11av*(rho12-rho11)**2 +
$ t21av/6.0*(rho22+rho21)**2 +
$ 121.0/40.*t31av*(rho32-rho31)**2
rhobar1 = sqrt(rhobar1)
rhobar2 = (Z12*rho01)**2 + 2.0/3.0*t21av*rho21**2
rhobar2 = sqrt(rhobar2)
else
rhobar2 = ((Z12/2)*(rho01+rho02))**2 +
$ t12av*(rho11-rho12)**2 +
$ t22av/6.0*(rho21+rho22)**2 +
$ 121.0/40.*t32av*(rho31-rho32)**2
rhobar2 = sqrt(rhobar2)
rhobar1 = (Z12*rho02)**2 + 2.0/3.0*t22av*rho22**2
rhobar1 = sqrt(rhobar1)
endif
else
c for other structures, use formalism developed in Huang's paper
c
-c composition-dependent scaling, equation I.7
-c -- note: The shape factors for the individual elements are used, since
-c this function should cancel the F(rho) terms when the atoms are
-c in the reference structure.
-c -- GJW (6/12/09): I suspect there may be a problem here... since we should be
-c using the pure element reference structure, we should also
-c be using the element "t" values (not the averages compute
-c above). It doesn't matter for reference structures with
-c s(1)=s(2)=s(3)=0, but might cause diffs otherwise. For now
-c what's here is consistent with what's done in meam_dens_final.
- Z1 = Z_meam(a)
- Z2 = Z_meam(b)
- if (ibar_meam(a).le.0) then
+c composition-dependent scaling, equation I.7
+c If using mixing rule for t, apply to reference structure; else
+c use precomputed values
+ if (mix_ref_t.eq.1) then
+ Z1 = Z_meam(a)
+ Z2 = Z_meam(b)
+ if (ibar_meam(a).le.0) then
G1 = 1.d0
- else
+ else
call get_shpfcn(s1,lattce_meam(a,a))
Gam1 = (s1(1)*t11av+s1(2)*t21av+s1(3)*t31av)/(Z1*Z1)
call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag)
- endif
- if (ibar_meam(b).le.0) then
+ endif
+ if (ibar_meam(b).le.0) then
G2 = 1.d0
- else
+ else
call get_shpfcn(s2,lattce_meam(b,b))
Gam2 = (s2(1)*t12av+s2(2)*t22av+s2(3)*t32av)/(Z2*Z2)
call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag)
- endif
- rho0_1 = rho0_meam(a)*Z1*G1
- rho0_2 = rho0_meam(b)*Z2*G2
-
-c get number of neighbors in the reference structure
-c Nref(i,j) = # of i's neighbors of type j
-c call get_Zij(Z12,lattce_meam(a,b))
-c
-c call get_densref(r,a,b,rho01,rho11,rho21,rho31,
-c $ rho02,rho12,rho22,rho32)
-
-c compute total background density, eqn I.6
- Gam1 = (t11av*rho11+t21av*rho21+t31av*rho31)
- Gam1 = Gam1/(rho01*rho01)
- Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32)
- Gam2 = Gam2/(rho02*rho02)
- call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag)
- call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag)
- rhobar1 = rho01/rho0_1*G1
- rhobar2 = rho02/rho0_2*G2
+ endif
+ rho0_1 = rho0_meam(a)*Z1*G1
+ rho0_2 = rho0_meam(b)*Z2*G2
+ endif
+ Gam1 = (t11av*rho11+t21av*rho21+t31av*rho31)
+ Gam1 = Gam1/(rho01*rho01)
+ Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32)
+ Gam2 = Gam2/(rho02*rho02)
+ call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag)
+ call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag)
+ if (mix_ref_t.eq.1) then
+ rho_bkgd1 = rho0_1
+ rho_bkgd2 = rho0_2
+ else
+ rho_bkgd1 = rho_ref_meam(a)
+ rho_bkgd2 = rho_ref_meam(b)
+ endif
+ rhobar1 = rho01/rho_bkgd1*G1
+ rhobar2 = rho02/rho_bkgd2*G2
+
endif
+
c compute embedding functions, eqn I.5
if (rhobar1.eq.0.d0) then
F1 = 0.d0
else
F1 = A_meam(a)*Ec_meam(a,a)*rhobar1*log(rhobar1)
endif
if (rhobar2.eq.0.d0) then
F2 = 0.d0
else
F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2)
endif
c compute Rose function, I.16
Eu = erose(r,re_meam(a,b),alpha_meam(a,b),
- $ Ec_meam(a,b),repuls_meam(a,b),attrac_meam(a,b))
+ $ Ec_meam(a,b),repuls_meam(a,b),attrac_meam(a,b),erose_form)
c calculate the pair energy
if (lattce_meam(a,b).eq.'c11') then
latta = lattce_meam(a,a)
if (latta.eq.'dia') then
phiaa = phi_meam(r,a,a)
phi_m = (3*Eu - F2 - 2*F1 - 5*phiaa)/Z12
else
phibb = phi_meam(r,b,b)
phi_m = (3*Eu - F1 - 2*F2 - 5*phibb)/Z12
endif
else if (lattce_meam(a,b).eq.'l12') then
phiaa = phi_meam(r,a,a)
c account for second neighbor a-a potential here...
call get_Zij(Z1nn,lattce_meam(a,a))
call get_Zij2(Z2nn,arat,scrn,lattce_meam(a,a),
$ Cmin_meam(a,a,a),Cmax_meam(a,a,a))
nmax = 10
if (scrn.gt.0.0) then
do n = 1,nmax
phiaa = phiaa +
$ (-Z2nn*scrn/Z1nn)**n * phi_meam(r*arat**n,a,a)
enddo
endif
phi_m = Eu/3. - F1/4. - F2/12. - phiaa
else
c
c potential is computed from Rose function and embedding energy
phi_m = (2*Eu - F1 - F2)/Z12
c
endif
c if r = 0, just return 0
if (r.eq.0.d0) then
phi_m = 0.d0
endif
return
end
+c----------------------------------------------------------------------c
+c Compute background density for reference structure of each element
+ subroutine compute_reference_density
+ use meam_data
+ implicit none
+
+ integer a,Z,Z2,errorflag
+ real*8 gam,Gbar,shp(3)
+ real*8 rho0,rho0_2nn,arat,scrn
+
+c loop over element types
+ do a = 1,neltypes
+
+ Z = Z_meam(a)
+ if (ibar_meam(a).le.0) then
+ Gbar = 1.d0
+ else
+ call get_shpfcn(shp,lattce_meam(a,a))
+ gam = (t1_meam(a)*shp(1)+t2_meam(a)*shp(2)
+ $ +t3_meam(a)*shp(3))/(Z*Z)
+ call G_gam(gam,ibar_meam(a),gsmooth_factor,
+ $ Gbar,errorflag)
+ endif
+
+c The zeroth order density in the reference structure, with
+c equilibrium spacing, is just the number of first neighbors times
+c the rho0_meam coefficient...
+ rho0 = rho0_meam(a)*Z
+
+c ...unless we have unscreened second neighbors, in which case we
+c add on the contribution from those (accounting for partial
+c screening)
+ if (nn2_meam(a,a).eq.1) then
+ call get_Zij2(Z2,arat,scrn,lattce_meam(a,a),
+ $ Cmin_meam(a,a,a),Cmax_meam(a,a,a))
+ rho0_2nn = rho0_meam(a)*exp(-beta0_meam(a)*(arat-1))
+ rho0 = rho0 + Z2*rho0_2nn*scrn
+ endif
+
+ rho_ref_meam(a) = rho0*Gbar
+
+ enddo
+
+ return
+ end
+
c----------------------------------------------------------------------c
c Shape factors for various configurations
subroutine get_shpfcn(s,latt)
implicit none
real*8 s(3)
character*3 latt
- if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'b1') then
+ if (latt.eq.'fcc'.or.latt.eq.'bcc'.
+ $ or.latt.eq.'b1'.or.latt.eq.'b2') then
s(1) = 0.d0
s(2) = 0.d0
s(3) = 0.d0
else if (latt.eq.'hcp') then
s(1) = 0.d0
s(2) = 0.d0
s(3) = 1.d0/3.d0
else if (latt.eq.'dia') then
s(1) = 0.d0
s(2) = 0.d0
s(3) = 32.d0/9.d0
else if (latt.eq.'dim') then
s(1) = 1.d0
s(2) = 2.d0/3.d0
c s(3) = 1.d0
s(3) = 0.4d0
else
s(1) = 0.0
c call error('Lattice not defined in get_shpfcn.')
endif
return
end
c------------------------------------------------------------------------------c
c Average weighting factors for the reference structure
subroutine get_tavref(t11av,t21av,t31av,t12av,t22av,t32av,
$ t11,t21,t31,t12,t22,t32,
$ r,a,b,latt)
use meam_data
implicit none
real*8 t11av,t21av,t31av,t12av,t22av,t32av
real*8 t11,t21,t31,t12,t22,t32,r
integer a,b
character*3 latt
real*8 rhoa01,rhoa02,a1,a2,rho01,rho02
- if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'dia'
- $ .or.latt.eq.'hcp'.or.latt.eq.'b1'
- $ .or.latt.eq.'dim') then
-c all neighbors are of the opposite type
- t11av = t12
- t21av = t22
- t31av = t32
- t12av = t11
- t22av = t21
- t32av = t31
+
+c For ialloy = 2, no averaging is done
+ if (ialloy.eq.2) then
+ t11av = t11
+ t21av = t21
+ t31av = t31
+ t12av = t12
+ t22av = t22
+ t32av = t32
else
- a1 = r/re_meam(a,a) - 1.d0
- a2 = r/re_meam(b,b) - 1.d0
- rhoa01 = rho0_meam(a)*exp(-beta0_meam(a)*a1)
- rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2)
- if (latt.eq.'l12') then
- rho01 = 8*rhoa01 + 4*rhoa02
- t11av = (8*t11*rhoa01 + 4*t12*rhoa02)/rho01
+ if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'dia'
+ $ .or.latt.eq.'hcp'.or.latt.eq.'b1'
+ $ .or.latt.eq.'dim'.or.latt.eq.'b2') then
+c all neighbors are of the opposite type
+ t11av = t12
+ t21av = t22
+ t31av = t32
t12av = t11
- t21av = (8*t21*rhoa01 + 4*t22*rhoa02)/rho01
t22av = t21
- t31av = (8*t31*rhoa01 + 4*t32*rhoa02)/rho01
t32av = t31
else
-c call error('Lattice not defined in get_tavref.')
+ a1 = r/re_meam(a,a) - 1.d0
+ a2 = r/re_meam(b,b) - 1.d0
+ rhoa01 = rho0_meam(a)*exp(-beta0_meam(a)*a1)
+ rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2)
+ if (latt.eq.'l12') then
+ rho01 = 8*rhoa01 + 4*rhoa02
+ t11av = (8*t11*rhoa01 + 4*t12*rhoa02)/rho01
+ t12av = t11
+ t21av = (8*t21*rhoa01 + 4*t22*rhoa02)/rho01
+ t22av = t21
+ t31av = (8*t31*rhoa01 + 4*t32*rhoa02)/rho01
+ t32av = t31
+ else
+c call error('Lattice not defined in get_tavref.')
+ endif
endif
endif
return
end
c------------------------------------------------------------------------------c
c Number of neighbors for the reference structure
subroutine get_Zij(Zij,latt)
implicit none
integer Zij
character*3 latt
if (latt.eq.'fcc') then
Zij = 12
else if (latt.eq.'bcc') then
Zij = 8
else if (latt.eq.'hcp') then
Zij = 12
else if (latt.eq.'b1') then
Zij = 6
else if (latt.eq.'dia') then
Zij = 4
else if (latt.eq.'dim') then
Zij = 1
else if (latt.eq.'c11') then
Zij = 10
else if (latt.eq.'l12') then
Zij = 12
+ else if (latt.eq.'b2') then
+ Zij = 8
else
c call error('Lattice not defined in get_Zij.')
endif
return
end
c------------------------------------------------------------------------------c
c Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second
c neighbor screening function for lattice type "latt"
subroutine get_Zij2(Zij2,a,S,latt,cmin,cmax)
implicit none
integer Zij2
real*8 a,S,cmin,cmax
character*3 latt
real*8 rratio,C,x,sijk
integer numscr
if (latt.eq.'bcc') then
Zij2 = 6
a = 2.d0/sqrt(3.d0)
numscr = 4
else if (latt.eq.'fcc') then
Zij2 = 6
a = sqrt(2.d0)
numscr = 4
else if (latt.eq.'dia') then
Zij2 = 0
a = sqrt(8.d0/3.d0)
numscr = 4
if (cmin.lt.0.500001) then
c call error('can not do 2NN MEAM for dia')
endif
else if (latt.eq.'hcp') then
Zij2 = 6
a = sqrt(2.d0)
numscr = 4
else if (latt.eq.'b1') then
Zij2 = 12
a = sqrt(2.d0)
numscr = 2
else if (latt.eq.'l12') then
Zij2 = 6
a = sqrt(2.d0)
numscr = 4
+ else if (latt.eq.'b2') then
+ Zij2 = 6
+ a = 2.d0/sqrt(3.d0)
+ numscr = 4
+ else if (latt.eq.'dim') then
+c this really shouldn't be allowed; make sure screening is zero
+ Zij2 = 0
+ a = 1
+ S = 0
+ return
else
c call error('Lattice not defined in get_Zij2.')
endif
c Compute screening for each first neighbor
C = 4.d0/(a*a) - 1.d0
x = (C-cmin)/(cmax-cmin)
call fcut(x,sijk)
c There are numscr first neighbors screening the second neighbors
S = sijk**numscr
return
end
c------------------------------------------------------------------------------c
subroutine get_sijk(C,i,j,k,sijk)
use meam_data
implicit none
real*8 C,sijk
integer i,j,k
real*8 x
x = (C-Cmin_meam(i,j,k))/(Cmax_meam(i,j,k)-Cmin_meam(i,j,k))
call fcut(x,sijk)
return
end
c------------------------------------------------------------------------------c
c Calculate density functions, assuming reference configuration
subroutine get_densref(r,a,b,rho01,rho11,rho21,rho31,
$ rho02,rho12,rho22,rho32)
use meam_data
implicit none
real*8 r,rho01,rho11,rho21,rho31,rho02,rho12,rho22,rho32
real*8 a1,a2
real*8 rhoa01,rhoa11,rhoa21,rhoa31,rhoa02,rhoa12,rhoa22,rhoa32
real*8 s(3)
character*3 lat
integer a,b
integer Zij1nn,Zij2nn
real*8 rhoa01nn,rhoa02nn
real*8 arat,scrn,denom
real*8 C,s111,s112,s221,S11,S22
a1 = r/re_meam(a,a) - 1.d0
a2 = r/re_meam(b,b) - 1.d0
rhoa01 = rho0_meam(a)*exp(-beta0_meam(a)*a1)
rhoa11 = rho0_meam(a)*exp(-beta1_meam(a)*a1)
rhoa21 = rho0_meam(a)*exp(-beta2_meam(a)*a1)
rhoa31 = rho0_meam(a)*exp(-beta3_meam(a)*a1)
rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2)
rhoa12 = rho0_meam(b)*exp(-beta1_meam(b)*a2)
rhoa22 = rho0_meam(b)*exp(-beta2_meam(b)*a2)
rhoa32 = rho0_meam(b)*exp(-beta3_meam(b)*a2)
lat = lattce_meam(a,b)
rho11 = 0.d0
rho21 = 0.d0
rho31 = 0.d0
rho12 = 0.d0
rho22 = 0.d0
rho32 = 0.d0
call get_Zij(Zij1nn,lat)
if (lat.eq.'fcc') then
rho01 = 12.d0*rhoa02
rho02 = 12.d0*rhoa01
else if (lat.eq.'bcc') then
rho01 = 8.d0*rhoa02
rho02 = 8.d0*rhoa01
else if (lat.eq.'b1') then
rho01 = 6*rhoa02
rho02 = 6*rhoa01
else if (lat.eq.'dia') then
rho01 = 4*rhoa02
rho02 = 4*rhoa01
rho31 = 32.d0/9.d0*rhoa32*rhoa32
rho32 = 32.d0/9.d0*rhoa31*rhoa31
else if (lat.eq.'hcp') then
rho01 = 12*rhoa02
rho02 = 12*rhoa01
rho31 = 1.d0/3.d0*rhoa32*rhoa32
rho32 = 1.d0/3.d0*rhoa31*rhoa31
else if (lat.eq.'dim') then
call get_shpfcn(s,'dim')
rho01 = rhoa02
rho02 = rhoa01
rho11 = s(1)*rhoa12*rhoa12
rho12 = s(1)*rhoa11*rhoa11
rho21 = s(2)*rhoa22*rhoa22
rho22 = s(2)*rhoa21*rhoa21
rho31 = s(3)*rhoa32*rhoa32
rho32 = s(3)*rhoa31*rhoa31
else if (lat.eq.'c11') then
rho01 = rhoa01
rho02 = rhoa02
rho11 = rhoa11
rho12 = rhoa12
rho21 = rhoa21
rho22 = rhoa22
rho31 = rhoa31
rho32 = rhoa32
else if (lat.eq.'l12') then
rho01 = 8*rhoa01 + 4*rhoa02
rho02 = 12*rhoa01
- if (ialloy.eq.0) then
- rho21 = 8./3.*(rhoa21-rhoa22)*(rhoa21-rhoa22)
- else
+ if (ialloy.eq.1) then
rho21 = 8./3.*(rhoa21*t2_meam(a)-rhoa22*t2_meam(b))**2
denom = 8*rhoa01*t2_meam(a)**2 + 4*rhoa02*t2_meam(b)**2
if (denom.gt.0.) then
rho21 = rho21/denom * rho01
endif
+ else
+ rho21 = 8./3.*(rhoa21-rhoa22)*(rhoa21-rhoa22)
endif
+ else if (lat.eq.'b2') then
+ rho01 = 8.d0*rhoa02
+ rho02 = 8.d0*rhoa01
else
c call error('Lattice not defined in get_densref.')
endif
if (nn2_meam(a,b).eq.1) then
call get_Zij2(Zij2nn,arat,scrn,lat,
$ Cmin_meam(a,a,b),Cmax_meam(a,a,b))
a1 = arat*r/re_meam(a,a) - 1.d0
a2 = arat*r/re_meam(b,b) - 1.d0
rhoa01nn = rho0_meam(a)*exp(-beta0_meam(a)*a1)
rhoa02nn = rho0_meam(b)*exp(-beta0_meam(b)*a2)
if (lat.eq.'l12') then
c As usual, L12 thinks it's special; we need to be careful computing
c the screening functions
C = 1.d0
call get_sijk(C,a,a,a,s111)
call get_sijk(C,a,a,b,s112)
call get_sijk(C,b,b,a,s221)
S11 = s111 * s111 * s112 * s112
S22 = s221**4
rho01 = rho01 + 6*S11*rhoa01nn
rho02 = rho02 + 6*S22*rhoa02nn
else
c For other cases, assume that second neighbor is of same type,
c first neighbor may be of different type
rho01 = rho01 + Zij2nn*scrn*rhoa01nn
c Assume Zij2nn and arat don't depend on order, but scrn might
call get_Zij2(Zij2nn,arat,scrn,lat,
$ Cmin_meam(b,b,a),Cmax_meam(b,b,a))
rho02 = rho02 + Zij2nn*scrn*rhoa02nn
endif
endif
return
end
c---------------------------------------------------------------------
c Compute ZBL potential
c
real*8 function zbl(r,z1,z2)
implicit none
integer i,z1,z2
real*8 r,c,d,a,azero,cc,x
dimension c(4),d(4)
data c /0.028171,0.28022,0.50986,0.18175/
data d /0.20162,0.40290,0.94229,3.1998/
data azero /0.4685/
data cc /14.3997/
c azero = (9pi^2/128)^1/3 (0.529) Angstroms
a = azero/(z1**0.23+z2**0.23)
zbl = 0.0
x = r/a
do i=1,4
zbl = zbl + c(i)*exp(-d(i)*x)
enddo
if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc
return
end
c---------------------------------------------------------------------
c Compute Rose energy function, I.16
c
- real*8 function erose(r,re,alpha,Ec,repuls,attrac)
+ real*8 function erose(r,re,alpha,Ec,repuls,attrac,form)
implicit none
real*8 r,re,alpha,Ec,repuls,attrac,astar,a3
+ integer form
erose = 0.d0
if (r.gt.0.d0) then
astar = alpha * (r/re - 1.d0)
a3 = 0.d0
if (astar.ge.0) then
a3 = attrac
else if (astar.lt.0) then
a3 = repuls
endif
- erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar)
-c erose = -Ec*(1+astar+(-attrac+repuls/r)*(astar**3))*exp(-astar)
+ if (form.eq.1) then
+ erose = -Ec*(1+astar+(-attrac+repuls/r)*
+ $ (astar**3))*exp(-astar)
+ else if (form.eq.2) then
+ erose = -Ec * (1 +astar + a3*(astar**3))*exp(-astar)
+ else
+ erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar)
+ endif
endif
return
end
c -----------------------------------------------------------------------
subroutine interpolate_meam(ind)
use meam_data
implicit none
integer j,ind
real*8 drar
c map to coefficient space
nrar = nr
drar = dr
rdrar = 1.0D0/drar
c phir interp
do j = 1,nrar
phirar(j,ind) = phir(j,ind)
enddo
phirar1(1,ind) = phirar(2,ind)-phirar(1,ind)
phirar1(2,ind) = 0.5D0*(phirar(3,ind)-phirar(1,ind))
phirar1(nrar-1,ind) = 0.5D0*(phirar(nrar,ind)
$ -phirar(nrar-2,ind))
phirar1(nrar,ind) = 0.0D0
do j = 3,nrar-2
phirar1(j,ind) = ((phirar(j-2,ind)-phirar(j+2,ind)) +
$ 8.0D0*(phirar(j+1,ind)-phirar(j-1,ind)))/12.
enddo
do j = 1,nrar-1
phirar2(j,ind) = 3.0D0*(phirar(j+1,ind)-phirar(j,ind)) -
$ 2.0D0*phirar1(j,ind) - phirar1(j+1,ind)
phirar3(j,ind) = phirar1(j,ind) + phirar1(j+1,ind) -
$ 2.0D0*(phirar(j+1,ind)-phirar(j,ind))
enddo
phirar2(nrar,ind) = 0.0D0
phirar3(nrar,ind) = 0.0D0
do j = 1,nrar
phirar4(j,ind) = phirar1(j,ind)/drar
phirar5(j,ind) = 2.0D0*phirar2(j,ind)/drar
phirar6(j,ind) = 3.0D0*phirar3(j,ind)/drar
enddo
end
c---------------------------------------------------------------------
c Compute Rose energy function, I.16
c
real*8 function compute_phi(rij, elti, eltj)
use meam_data
implicit none
real*8 rij, pp
integer elti, eltj, ind, kk
ind = eltind(elti, eltj)
pp = rij*rdrar + 1.0D0
kk = pp
kk = min(kk,nrar-1)
pp = pp - kk
pp = min(pp,1.0D0)
compute_phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp
$ + phirar1(kk,ind))*pp + phirar(kk,ind)
return
end
diff --git a/lib/meam/meam_setup_global.F b/lib/meam/meam_setup_global.F
index 9603c4ecf..548c97b25 100755
--- a/lib/meam/meam_setup_global.F
+++ b/lib/meam/meam_setup_global.F
@@ -1,106 +1,109 @@
c
c declaration in pair_meam.h:
c
c void meam_setup_global(int *, int *, double *, int *, double *, double *,
c double *, double *, double *, double *, double *,
c double *, double *, double *, double *, double *,
c double *, double *, int *);
c
c call in pair_meam.cpp:
c
c meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
c alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
c
c
subroutine meam_setup_global(nelt, lat, z, ielement, atwt, alpha,
$ b0, b1, b2, b3, alat, esub, asub,
$ t0, t1, t2, t3, rozero, ibar)
use meam_data
implicit none
integer nelt, lat, ielement, ibar
real*8 z, atwt, alpha, b0, b1, b2, b3
real*8 alat, esub, asub, t0, t1, t2, t3
real*8 rozero
dimension lat(nelt), ielement(nelt), ibar(nelt)
dimension z(nelt), atwt(nelt), alpha(nelt)
dimension b0(nelt), b1(nelt), b2(nelt), b3(nelt)
dimension alat(nelt), esub(nelt), asub(nelt)
dimension t0(nelt), t1(nelt), t2(nelt), t3(nelt), rozero(nelt)
integer i
real*8 tmplat(maxelt)
neltypes = nelt
do i = 1,nelt
if (lat(i).eq.0) then
lattce_meam(i,i) = 'fcc'
else if (lat(i).eq.1) then
lattce_meam(i,i) = 'bcc'
else if (lat(i).eq.2) then
lattce_meam(i,i) = 'hcp'
else if (lat(i).eq.3) then
lattce_meam(i,i) = 'dim'
else if (lat(i).eq.4) then
lattce_meam(i,i) = 'dia'
else
c unknown
endif
Z_meam(i) = z(i)
ielt_meam(i) = ielement(i)
alpha_meam(i,i) = alpha(i)
beta0_meam(i) = b0(i)
beta1_meam(i) = b1(i)
beta2_meam(i) = b2(i)
beta3_meam(i) = b3(i)
tmplat(i) = alat(i)
Ec_meam(i,i) = esub(i)
A_meam(i) = asub(i)
t0_meam(i) = t0(i)
t1_meam(i) = t1(i)
t2_meam(i) = t2(i)
t3_meam(i) = t3(i)
rho0_meam(i) = rozero(i)
ibar_meam(i) = ibar(i)
if (lattce_meam(i,i).eq.'fcc') then
re_meam(i,i) = tmplat(i)/sqrt(2.d0)
elseif (lattce_meam(i,i).eq.'bcc') then
re_meam(i,i) = tmplat(i)*sqrt(3.d0)/2.d0
elseif (lattce_meam(i,i).eq.'hcp') then
re_meam(i,i) = tmplat(i)
elseif (lattce_meam(i,i).eq.'dim') then
re_meam(i,i) = tmplat(i)
elseif (lattce_meam(i,i).eq.'dia') then
re_meam(i,i) = tmplat(i)*sqrt(3.d0)/4.d0
else
c error
endif
enddo
c Set some defaults
rc_meam = 4.0
delr_meam = 0.1
attrac_meam(:,:) = 0.0
repuls_meam(:,:) = 0.0
Cmax_meam(:,:,:) = 2.8
Cmin_meam(:,:,:) = 2.0
ebound_meam(:,:) = (2.8d0**2)/(4.d0*(2.8d0-1.d0))
delta_meam(:,:) = 0.0
nn2_meam(:,:) = 0
+ zbl_meam(:,:) = 1
gsmooth_factor = 99.0
augt1 = 1
ialloy = 0
+ mix_ref_t = 0
+ erose_form = 0
return
end
diff --git a/lib/meam/meam_setup_param.F b/lib/meam/meam_setup_param.F
index d70afa74c..70fccf8da 100755
--- a/lib/meam/meam_setup_param.F
+++ b/lib/meam/meam_setup_param.F
@@ -1,128 +1,150 @@
c
c Declaration in pair_meam.h:
c
c void meam_setup_param(int *, double *, int *, int *, int *);
c
c Call in pair_meam.cpp
c
c meam_setup_param(&which,&value,&nindex,index,&errorflag);
c
c
c
c The "which" argument corresponds to the index of the "keyword" array
c in pair_meam.cpp:
c
c 0 = Ec_meam
c 1 = alpha_meam
c 2 = rho0_meam
c 3 = delta_meam
c 4 = lattce_meam
c 5 = attrac_meam
c 6 = repuls_meam
c 7 = nn2_meam
c 8 = Cmin_meam
c 9 = Cmax_meam
c 10 = rc_meam
c 11 = delr_meam
c 12 = augt1
c 13 = gsmooth_factor
c 14 = re_meam
c 15 = ialloy
+c 16 = mixture_ref_t
+c 17 = erose_form
+c 18 = zbl_meam
subroutine meam_setup_param(which, value, nindex,
$ index, errorflag)
use meam_data
implicit none
integer which, nindex, index(3), errorflag
real*8 value
+ integer i1, i2
errorflag = 0
c 0 = Ec_meam
if (which.eq.0) then
Ec_meam(index(1),index(2)) = value
c 1 = alpha_meam
else if (which.eq.1) then
alpha_meam(index(1),index(2)) = value
c 2 = rho0_meam
else if (which.eq.2) then
rho0_meam(index(1)) = value
c 3 = delta_meam
else if (which.eq.3) then
delta_meam(index(1),index(2)) = value
c 4 = lattce_meam
else if (which.eq.4) then
if (value.eq.0) then
lattce_meam(index(1),index(2)) = "fcc"
else if (value.eq.1) then
lattce_meam(index(1),index(2)) = "bcc"
else if (value.eq.2) then
lattce_meam(index(1),index(2)) = "hcp"
else if (value.eq.3) then
lattce_meam(index(1),index(2)) = "dim"
else if (value.eq.4) then
lattce_meam(index(1),index(2)) = "dia"
else if (value.eq.5) then
lattce_meam(index(1),index(2)) = 'b1'
else if (value.eq.6) then
lattce_meam(index(1),index(2)) = 'c11'
else if (value.eq.7) then
lattce_meam(index(1),index(2)) = 'l12'
+ else if (value.eq.8) then
+ lattce_meam(index(1),index(2)) = 'b2'
endif
c 5 = attrac_meam
else if (which.eq.5) then
attrac_meam(index(1),index(2)) = value
c 6 = repuls_meam
else if (which.eq.6) then
repuls_meam(index(1),index(2)) = value
c 7 = nn2_meam
else if (which.eq.7) then
- nn2_meam(index(1),index(2)) = value
+ i1 = min(index(1),index(2))
+ i2 = max(index(1),index(2))
+ nn2_meam(i1,i2) = value
c 8 = Cmin_meam
else if (which.eq.8) then
Cmin_meam(index(1),index(2),index(3)) = value
c 9 = Cmax_meam
else if (which.eq.9) then
Cmax_meam(index(1),index(2),index(3)) = value
c 10 = rc_meam
else if (which.eq.10) then
rc_meam = value
c 11 = delr_meam
else if (which.eq.11) then
delr_meam = value
c 12 = augt1
else if (which.eq.12) then
augt1 = value
c 13 = gsmooth
else if (which.eq.13) then
gsmooth_factor = value
c 14 = re_meam
else if (which.eq.14) then
re_meam(index(1),index(2)) = value
c 15 = ialloy
else if (which.eq.15) then
ialloy = value
+c 16 = mixture_ref_t
+ else if (which.eq.16) then
+ mix_ref_t = value
+
+c 17 = erose_form
+ else if (which.eq.17) then
+ erose_form = value
+
+c 18 = zbl_meam
+ else if (which.eq.18) then
+ i1 = min(index(1),index(2))
+ i2 = max(index(1),index(2))
+ zbl_meam(i1,i2) = value
+
else
errorflag = 1
endif
return
end