ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
sjplimp
2010-10-06 14:56:50 +00:00
parent d2442ba5df
commit ecbd64b062
1 changed files with 12 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
lib/meam/meam_dens_init.F
View File

@ -154,12 +154,17 @@ c Now compute derivatives
xik = rik2/rij2 xik = rik2/rij2
xjk = rjk2/rij2 xjk = rjk2/rij2
a = 1 - (xik-xjk)*(xik-xjk) a = 1 - (xik-xjk)*(xik-xjk)
if (a.eq.0.d0) goto 10 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 cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
Cmax = Cmax_meam(elti,eltj,eltk) Cmax = Cmax_meam(elti,eltj,eltk)
Cmin = Cmin_meam(elti,eltj,eltk) Cmin = Cmin_meam(elti,eltj,eltk)
if (cikj.ge.Cmax) then if (cikj.ge.Cmax) then
goto 10 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) c Note that we never have 0<cikj<Cmin here, else sij=0 (rejected above)
else else
delc = Cmax - Cmin delc = Cmax - Cmin
@ -364,12 +369,17 @@ c if rjksq > ebound*rijsq, atom k is definitely outside the ellipse
xik = riksq/rijsq xik = riksq/rijsq
xjk = rjksq/rijsq xjk = rjksq/rijsq
a = 1 - (xik-xjk)*(xik-xjk) a = 1 - (xik-xjk)*(xik-xjk)
if (a.eq.0.d0) goto 10 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 cikj = (2.d0*(xik+xjk) + a - 2.d0)/a
Cmax = Cmax_meam(elti,eltj,eltk) Cmax = Cmax_meam(elti,eltj,eltk)
Cmin = Cmin_meam(elti,eltj,eltk) Cmin = Cmin_meam(elti,eltj,eltk)
if (cikj.ge.Cmax) then if (cikj.ge.Cmax) then
goto 10 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 else if (cikj.le.Cmin) then
sij = 0.d0 sij = 0.d0
goto 20 goto 20