jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6544 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2011-07-19 18:31:58 +00:00
parent 2402fbd086
commit 5f6811f4e0
5 changed files with 74 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
lib/meam/meam_data.F
View File

@ -31,6 +31,8 @@ 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 emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0
c bkgd_dyn = 1 if reference densities follows Dynamo, 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
@ -73,6 +75,7 @@ c nrar,rdrar = spline coeff array parameters
real*8 Cmax_meam(maxelt,maxelt,maxelt)
real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt)
integer augt1, ialloy, mix_ref_t, erose_form
integer emb_lin_neg, bkgd_dyn
real*8 gsmooth_factor
integer vind2D(3,3),vind3D(3,3,3)
integer v2D(6),v3D(10)

46
lib/meam/meam_dens_final.F
View File

@ -127,7 +127,11 @@ c Complete the calculation of density
endif
rho_bkgd = rho0_meam(elti)*Z*Gbar
else
rho_bkgd = rho_ref_meam(elti)
if (bkgd_dyn.eq.1) then
rho_bkgd = rho0_meam(elti)*Z
else
rho_bkgd = rho_ref_meam(elti)
endif
endif
rhob = rho(i)/rho_bkgd
denom = 1.d0/rho_bkgd
@ -154,17 +158,33 @@ c included for backward compatibility
B = A_meam(elti)*Ec_meam(elti,elti)
if( rhob .ne. 0.d0 ) then
fp(i) = B*(log(rhob)+1.d0)
if (emb_lin_neg.eq.1 .and. rhob.le.0) then
fp(i) = -B
else
fp(i) = B*(log(rhob)+1.d0)
endif
if (eflag_either.ne.0) then
if (eflag_global.ne.0) then
eng_vdwl = eng_vdwl + B*rhob*log(rhob)
if (emb_lin_neg.eq.1 .and. rhob.le.0) then
eng_vdwl = eng_vdwl - B*rhob
else
eng_vdwl = eng_vdwl + B*rhob*log(rhob)
endif
endif
if (eflag_atom.ne.0) then
eatom(i) = eatom(i) + B*rhob*log(rhob)
if (emb_lin_neg.eq.1 .and. rhob.le.0) then
eatom(i) = eatom(i) - B*rhob
else
eatom(i) = eatom(i) + B*rhob*log(rhob)
endif
endif
endif
else
fp(i) = B
if (emb_lin_neg.eq.1) then
fp(i) = -B
else
fp(i) = B
endif
endif
endif
enddo
@ -181,6 +201,7 @@ c 1 => G = exp(Gamma/2)
c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma))
c 4 => G = sqrt(1+Gamma)
c -5 => G = +-sqrt(abs(1+Gamma))
implicit none
real*8 Gamma,G
real*8 gsmooth_factor, gsmooth_switchpoint
@ -201,6 +222,12 @@ c G = 0.01*(-0.99/Gamma)**99
G = exp(Gamma/2.d0)
else if (ibar.eq.3) then
G = 2.d0/(1.d0+exp(-Gamma))
else if (ibar.eq.-5) then
if ((1.d0+Gamma).ge.0) then
G = sqrt(1.d0+Gamma)
else
G = -sqrt(-1.d0-Gamma)
endif
else
errorflag = 1
endif
@ -216,6 +243,7 @@ c 1 => G = exp(Gamma/2)
c 2 => not implemented
c 3 => G = 2/(1+exp(-Gamma))
c 4 => G = sqrt(1+Gamma)
c -5 => G = +-sqrt(abs(1+Gamma))
real*8 Gamma,G,dG
real*8 gsmooth_factor, gsmooth_switchpoint
integer ibar
@ -239,6 +267,14 @@ c G = 0.01*(-0.99/Gamma)**99
else if (ibar.eq.3) then
G = 2.d0/(1.d0+exp(-Gamma))
dG = G*(2.d0-G)/2
else if (ibar.eq.-5) then
if ((1.d0+Gamma).ge.0) then
G = sqrt(1.d0+Gamma)
dG = 1.d0/(2.d0*G)
else
G = -sqrt(-1.d0-Gamma)
dG = -1.d0/(2.d0*G)
endif
endif
return
end

22
lib/meam/meam_setup_done.F
View File

@ -105,6 +105,7 @@ c here or in the input file)
Ec_meam(i,j) = Ec_meam(j,i)
alpha_meam(i,j) = alpha_meam(j,i)
lattce_meam(i,j) = lattce_meam(j,i)
nn2_meam(i,j) = nn2_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
@ -446,8 +447,13 @@ c use precomputed values
rho_bkgd1 = rho0_1
rho_bkgd2 = rho0_2
else
rho_bkgd1 = rho_ref_meam(a)
rho_bkgd2 = rho_ref_meam(b)
if (bkgd_dyn.eq.1) then
rho_bkgd1 = rho0_meam(a)*Z_meam(a)
rho_bkgd2 = rho0_meam(b)*Z_meam(b)
else
rho_bkgd1 = rho_ref_meam(a)
rho_bkgd2 = rho_ref_meam(b)
endif
endif
rhobar1 = rho01/rho_bkgd1*G1
rhobar2 = rho02/rho_bkgd2*G2
@ -458,12 +464,20 @@ 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)
if (emb_lin_neg.eq.1 .and. rhobar1.le.0) then
F1 = -A_meam(a)*Ec_meam(a,a)*rhobar1
else
F1 = A_meam(a)*Ec_meam(a,a)*rhobar1*log(rhobar1)
endif
endif
if (rhobar2.eq.0.d0) then
F2 = 0.d0
else
F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2)
if (emb_lin_neg.eq.1 .and. rhobar2.le.0) then
F2 = -A_meam(b)*Ec_meam(b,b)*rhobar2
else
F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2)
endif
endif
c compute Rose function, I.16

2
lib/meam/meam_setup_global.F
View File

@ -101,6 +101,8 @@ c Set some defaults
augt1 = 1
ialloy = 0
mix_ref_t = 0
emb_lin_neg = 0
bkgd_dyn = 0
erose_form = 0
return

10
lib/meam/meam_setup_param.F
View File

@ -31,6 +31,8 @@ c 15 = ialloy
c 16 = mixture_ref_t
c 17 = erose_form
c 18 = zbl_meam
c 19 = emb_lin_neg
c 20 = bkgd_dyn
subroutine meam_setup_param(which, value, nindex,
$ index, errorflag)
@ -142,6 +144,14 @@ c 18 = zbl_meam
i2 = max(index(1),index(2))
zbl_meam(i1,i2) = value
c 19 = emb_lin_neg
else if (which.eq.19) then
emb_lin_neg = value
c 20 = bkgd_dyn
else if (which.eq.20) then
bkgd_dyn = value
else
errorflag = 1
endif