forked from lijiext/lammps
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3084 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
sjplimp
parent
826e513fd7
commit
8a6024b3ea
|
|
@ -34,6 +34,7 @@ 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 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
|
||||
|
|
@ -65,7 +66,7 @@ c nrar,rdrar = spline coeff array parameters
|
|||
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
|
||||
integer augt1, ialloy
|
||||
real*8 gsmooth_factor
|
||||
integer vind2D(3,3),vind3D(3,3,3)
|
||||
integer v2D(6),v3D(10)
|
||||
|
|
|
|||
|
|
@ -4,21 +4,21 @@ 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 *, int *);
|
||||
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],gamma,dgamma1,
|
||||
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,
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave,
|
||||
$ Gamma, dGamma1, dGamma2, dGamma3,
|
||||
$ rho, rho0, rho1, rho2, rho3, fp, errorflag)
|
||||
|
||||
|
|
@ -29,7 +29,7 @@ c
|
|||
integer ntype, type, fmap
|
||||
real*8 eng_vdwl, eatom, Arho1, Arho2
|
||||
real*8 Arho2b, Arho3, Arho3b
|
||||
real*8 t_ave
|
||||
real*8 t_ave, tsq_ave
|
||||
real*8 Gamma, dGamma1, dGamma2, dGamma3
|
||||
real*8 rho, rho0, rho1, rho2, rho3
|
||||
real*8 fp
|
||||
|
|
@ -39,6 +39,7 @@ c
|
|||
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)
|
||||
|
|
@ -70,9 +71,15 @@ c Complete the calculation of density
|
|||
enddo
|
||||
|
||||
if( rho0(i) .gt. 0.0 ) then
|
||||
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)
|
||||
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
|
||||
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)
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@ 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 *, int *);
|
||||
c double *, double *, double *, double *, double *, int *);
|
||||
c
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
|
|
@ -12,7 +12,7 @@ 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],&errorflag);
|
||||
c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_dens_init(i, nmax,
|
||||
|
|
@ -20,7 +20,7 @@ c
|
|||
$ numneigh, firstneigh,
|
||||
$ numneigh_full, firstneigh_full,
|
||||
$ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave, errorflag)
|
||||
$ arho3, arho3b, t_ave, tsq_ave, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
|
@ -30,7 +30,7 @@ c
|
|||
integer numneigh, firstneigh, numneigh_full, firstneigh_full
|
||||
real*8 scrfcn, dscrfcn, fcpair
|
||||
real*8 rho0, arho1, arho2
|
||||
real*8 arho2b, arho3, arho3b, t_ave
|
||||
real*8 arho2b, arho3, arho3b, t_ave, tsq_ave
|
||||
integer errorflag
|
||||
integer j,jn
|
||||
|
||||
|
|
@ -40,7 +40,7 @@ c
|
|||
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)
|
||||
dimension t_ave(3,nmax), tsq_ave(3,nmax)
|
||||
|
||||
errorflag = 0
|
||||
|
||||
|
|
@ -54,7 +54,7 @@ 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)
|
||||
$ arho3, arho3b, t_ave, tsq_ave)
|
||||
|
||||
return
|
||||
end
|
||||
|
|
@ -194,7 +194,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|||
subroutine calc_rho1(i, nmax, ntype, type, fmap, x,
|
||||
$ numneigh, firstneigh,
|
||||
$ scrfcn, fcpair, rho0, arho1, arho2, arho2b,
|
||||
$ arho3, arho3b, t_ave)
|
||||
$ arho3, arho3b, t_ave, tsq_ave)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
|
@ -203,14 +203,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|||
real*8 x
|
||||
integer numneigh, firstneigh
|
||||
real*8 scrfcn, fcpair, rho0, arho1, arho2
|
||||
real*8 arho2b, arho3, arho3b, t_ave
|
||||
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)
|
||||
dimension t_ave(3,nmax), tsq_ave(3,nmax)
|
||||
|
||||
integer jn,j,m,n,p,elti,eltj
|
||||
integer nv2,nv3
|
||||
|
|
@ -248,6 +248,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|||
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
|
||||
|
|
@ -256,6 +264,20 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|||
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
|
||||
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
|
||||
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@ 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 *, int *);
|
||||
c double *, double *, double *, double *, double *, double *, int *);
|
||||
c
|
||||
c Call from pair_meam.cpp has the form:
|
||||
c
|
||||
|
|
@ -14,7 +14,7 @@ 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],&f[0][0],&vatom[0][0],&errorflag);
|
||||
c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag);
|
||||
c
|
||||
|
||||
subroutine meam_force(i, nmax,
|
||||
|
|
@ -23,8 +23,8 @@ c
|
|||
$ numneigh, firstneigh, numneigh_full, firstneigh_full,
|
||||
$ scrfcn, dscrfcn, fcpair,
|
||||
$ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp,
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, f, vatom,
|
||||
$ errorflag)
|
||||
$ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f,
|
||||
$ vatom, errorflag)
|
||||
|
||||
use meam_data
|
||||
implicit none
|
||||
|
|
@ -38,7 +38,7 @@ c
|
|||
real*8 rho0, rho1, rho2, rho3, fp
|
||||
real*8 Arho1, Arho2, Arho2b
|
||||
real*8 Arho3, Arho3b
|
||||
real*8 t_ave, f, vatom
|
||||
real*8 t_ave, tsq_ave, f, vatom
|
||||
integer errorflag
|
||||
|
||||
dimension eatom(nmax)
|
||||
|
|
@ -50,7 +50,7 @@ c
|
|||
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), f(3,nmax), vatom(6,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
|
||||
|
|
@ -176,6 +176,21 @@ c Compute pair densities and derivatives
|
|||
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
|
||||
|
|
@ -277,21 +292,59 @@ c Compute derivatives of weighting functions t wrt rij
|
|||
t2j = t_ave(2,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
|
||||
if (ialloy.eq.1) then
|
||||
|
||||
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)
|
||||
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
|
||||
|
||||
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))
|
||||
|
|
@ -340,21 +393,60 @@ c Compute derivatives wrt sij, but only if necessary
|
|||
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
|
||||
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
|
||||
|
||||
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
|
||||
|
||||
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
|
||||
|
|
|
|||
|
|
@ -158,11 +158,12 @@ c
|
|||
integer j,a,b,nv2
|
||||
real*8 astar,frac,phizbl
|
||||
integer n,nmax,Z1,Z2
|
||||
real*8 arat,scrn
|
||||
real*8 arat,scrn,scrn2
|
||||
real*8 phitmp
|
||||
|
||||
real*8, external :: phi_meam
|
||||
real*8, external :: zbl
|
||||
real*8, external :: compute_phi
|
||||
|
||||
|
||||
c allocate memory for array that defines the potential
|
||||
|
|
@ -178,6 +179,10 @@ c allocate coeff memory
|
|||
allocate(phirar5(nr,(neltypes*(neltypes+1))/2))
|
||||
allocate(phirar6(nr,(neltypes*(neltypes+1))/2))
|
||||
|
||||
c HACK for debug: compute phi_meam for Ni3Si equilibrium spacing
|
||||
r = 2.47770216
|
||||
phitmp = phi_meam(r,1,2)
|
||||
|
||||
c loop over pairs of element types
|
||||
nv2 = 0
|
||||
do a = 1,neltypes
|
||||
|
|
@ -194,20 +199,29 @@ c loop over r values and compute
|
|||
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
|
||||
if (a.ne.b) then
|
||||
write(6,*) 'Second-neighbor currently only valid '
|
||||
write(6,*) 'if element types are the same.'
|
||||
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(' ')
|
||||
endif
|
||||
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))
|
||||
|
||||
nmax = 10
|
||||
do n = 1,nmax
|
||||
phir(j,nv2) = phir(j,nv2) +
|
||||
$ (-Z2*scrn/Z1)**n * phi_meam(r*arat**n,a,b)
|
||||
enddo
|
||||
if (lattce_meam(a,b).eq.'b1') then
|
||||
phir(j,nv2) = phir(j,nv2) -
|
||||
$ Z2*scrn/(2*Z1) * phi_meam(r*arat,a,a)
|
||||
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) * phi_meam(r*arat,b,b)
|
||||
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
|
||||
|
||||
|
|
@ -258,7 +272,7 @@ c
|
|||
real*8 rho02,rho12,rho22,rho32
|
||||
real*8 scalfac,phiaa,phibb
|
||||
real*8 Eu
|
||||
integer Z12
|
||||
integer Z12, errorflag
|
||||
character*3 latta,lattb
|
||||
|
||||
real*8, external :: erose
|
||||
|
|
@ -285,8 +299,9 @@ c calculate average weighting factors for the reference structure
|
|||
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),lattce_meam(a,b))
|
||||
$ 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
|
||||
|
|
@ -316,6 +331,12 @@ 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
|
||||
|
|
@ -323,14 +344,14 @@ c in the reference structure.
|
|||
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)
|
||||
call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag)
|
||||
endif
|
||||
if (ibar_meam(b).le.0) then
|
||||
G2 = 1.d0
|
||||
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)
|
||||
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
|
||||
|
|
@ -347,8 +368,8 @@ c compute total background density, eqn I.6
|
|||
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)
|
||||
call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2)
|
||||
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
|
||||
|
|
@ -370,17 +391,20 @@ c compute Rose function, I.16
|
|||
|
||||
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
|
||||
c phi_m = 0.d0
|
||||
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
|
||||
c phi_m = 0.d0
|
||||
else if (lattce_meam(a,b).eq.'l12') then
|
||||
phiaa = phi_meam(r,a,a)
|
||||
phi_m = Eu/3. - F1/4. - F2/12. - phiaa
|
||||
else
|
||||
c
|
||||
c
|
||||
c potential is computed from Rose function and embedding energy
|
||||
phi_m = (2*Eu - F1 - F2)/Z12
|
||||
c
|
||||
|
|
@ -426,11 +450,15 @@ c call error('Lattice not defined in get_shpfcn.')
|
|||
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,latt)
|
||||
$ 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
|
||||
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
|
||||
|
|
@ -442,7 +470,21 @@ c all neighbors are of the opposite type
|
|||
t22av = t21
|
||||
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
|
||||
return
|
||||
end
|
||||
|
|
@ -466,6 +508,8 @@ c Number of neighbors for the reference structure
|
|||
Zij = 1
|
||||
else if (latt.eq.'c11') then
|
||||
Zij = 10
|
||||
else if (latt.eq.'l12') then
|
||||
Zij = 12
|
||||
else
|
||||
c call error('Lattice not defined in get_Zij.')
|
||||
endif
|
||||
|
|
@ -504,9 +548,9 @@ c call error('can not do 2NN MEAM for dia')
|
|||
a = sqrt(2.d0)
|
||||
numscr = 4
|
||||
else if (latt.eq.'b1') then
|
||||
Zij2 = 6
|
||||
Zij2 = 12
|
||||
a = sqrt(2.d0)
|
||||
numscr = 4
|
||||
numscr = 2
|
||||
else
|
||||
c call error('Lattice not defined in get_Zij2.')
|
||||
endif
|
||||
|
|
@ -535,7 +579,7 @@ c Calculate density functions, assuming reference configuration
|
|||
integer a,b
|
||||
integer Zij1nn,Zij2nn
|
||||
real*8 rhoa01nn,rhoa02nn
|
||||
real*8 arat,scrn
|
||||
real*8 arat,scrn,denom
|
||||
|
||||
a1 = r/re_meam(a,a) - 1.d0
|
||||
a2 = r/re_meam(b,b) - 1.d0
|
||||
|
|
@ -567,17 +611,6 @@ c Calculate density functions, assuming reference configuration
|
|||
rho01 = 8.d0*rhoa02
|
||||
rho02 = 8.d0*rhoa01
|
||||
else if (lat.eq.'b1') then
|
||||
c Warning: this assumes that atoms of element 1 are not completely
|
||||
c screened from each other, but that atoms of element 2 are. So
|
||||
c Cmin_meam(1,1,2) can be less than 1.0, but Cmin_meam(2,2,1) cannot.
|
||||
c This should be made more general in the future.
|
||||
c Cmin = Cmin_meam(1,1,2)
|
||||
c Cmax = Cmax_meam(1,1,2)
|
||||
c C = (1.0-Cmin)/(Cmax-Cmin)
|
||||
c call fcut(C,fc)
|
||||
c a1 = sqrt(2.0)*r/re_meam(1,1) - 1.d0
|
||||
c rho01 = 6*rhoa02 + 12*fc*fc*rho0_meam(1)*exp(-beta0_meam(1)*a1)
|
||||
c rho02 = 6*rhoa01
|
||||
rho01 = 6*rhoa02
|
||||
rho02 = 6*rhoa01
|
||||
else if (lat.eq.'dia') then
|
||||
|
|
@ -609,6 +642,18 @@ c rho02 = 6*rhoa01
|
|||
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
|
||||
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
|
||||
endif
|
||||
else
|
||||
c call error('Lattice not defined in get_densref.')
|
||||
endif
|
||||
|
|
@ -617,9 +662,9 @@ c call error('Lattice not defined in get_densref.')
|
|||
c Assume that second neighbor is of same type, first neighbor
|
||||
c may be of different type
|
||||
|
||||
if (lat.ne.'bcc') then
|
||||
c if (lat.ne.'bcc') then
|
||||
c call error('Second-neighbor not defined for this lattice')
|
||||
endif
|
||||
c endif
|
||||
|
||||
call get_Zij2(Zij2nn,arat,scrn,lat,
|
||||
$ Cmin_meam(a,a,b),Cmax_meam(a,a,b))
|
||||
|
|
@ -630,8 +675,12 @@ c call error('Second-neighbor not defined for this lattice')
|
|||
rhoa01nn = rho0_meam(a)*exp(-beta0_meam(a)*a1)
|
||||
rhoa02nn = rho0_meam(b)*exp(-beta0_meam(b)*a2)
|
||||
|
||||
rho01 = rho01 + Zij2nn*scrn*rhoa02nn
|
||||
rho02 = rho02 + Zij2nn*scrn*rhoa01nn
|
||||
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
|
||||
|
||||
|
|
@ -731,3 +780,25 @@ c phir interp
|
|||
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
|
||||
|
|
|
|||
|
|
@ -98,6 +98,7 @@ c Set some defaults
|
|||
nn2_meam(:,:) = 0
|
||||
gsmooth_factor = 99.0
|
||||
augt1 = 1
|
||||
ialloy = 0
|
||||
|
||||
return
|
||||
end
|
||||
|
|
|
|||
|
|
@ -27,6 +27,7 @@ c 11 = delr_meam
|
|||
c 12 = augt1
|
||||
c 13 = gsmooth_factor
|
||||
c 14 = re_meam
|
||||
c 15 = ialloy
|
||||
|
||||
subroutine meam_setup_param(which, value, nindex,
|
||||
$ index, errorflag)
|
||||
|
|
@ -71,6 +72,8 @@ c 4 = lattce_meam
|
|||
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'
|
||||
endif
|
||||
|
||||
c 5 = attrac_meam
|
||||
|
|
@ -113,6 +116,10 @@ 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
|
||||
|
||||
else
|
||||
errorflag = 1
|
||||
endif
|
||||
|
|
|
|||
Loading…
Reference in New Issue