From ecbd64b0627f791c40828e8206fd36e234829e8c Mon Sep 17 00:00:00 2001 From: sjplimp Date: Wed, 6 Oct 2010 14:56:50 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4975 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/meam/meam_dens_init.F | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/lib/meam/meam_dens_init.F b/lib/meam/meam_dens_init.F index a64c56298f..53ad0956df 100755 --- a/lib/meam/meam_dens_init.F +++ b/lib/meam/meam_dens_init.F @@ -154,12 +154,17 @@ c Now compute derivatives xik = rik2/rij2 xjk = rjk2/rij2 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 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 ebound*rijsq, atom k is definitely outside the ellipse xik = riksq/rijsq xjk = rjksq/rijsq 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 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