--- forces-edip.1.6.f	Fri Apr 19 10:43:29 2002
+++ forces-edip.1.7.f	Fri Apr 19 10:43:41 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.6f
+C   Version 1.7f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -182,7 +182,8 @@
           double precision dV2j,dV2ijx,dV2ijy,dV2ijz,pZ,dp
           double precision temp0,temp1,temp2
           double precision Qort,muhalf,u5
-          double precision rmbinv,winv,dwinv,tau,dtau,lcos,x,H,dHdx,dhdl
+          double precision rmbinv,winv,dwinv,tau,dtau,lcos,x,H
+          double precision winv_2lam,TdW_WdT_over_W
           double precision dV3rij,dV3rik,dV3l
           double precision dVdZ_sum
           double precision dEdrl,dEdrlx,dEdrly,dEdrlz
@@ -478,8 +479,8 @@
         
               winv = Qort*exp(muhalf*Z) 
 C   inverse width of angular function
-              dwinv = muhalf*winv       
-C   its derivative
+              dwinv = muhalf
+C   its derivative divided by itself
               temp0 = exp(-u4*Z)
               tau = u1+u2*temp0*(u3-temp0) 
 C   -cosine of angular minimum
@@ -487,6 +488,8 @@
 C   its derivative
             end if
         
+            winv_2lam = 2.d0*par_lam*winv
+            TdW_WdT_over_W = tau*dwinv + dtau
             
 C  --- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---
         
@@ -507,17 +510,16 @@
                 lcos = s3_dx(nj) * s3_dx(nk) + s3_dy(nj) * s3_dy(nk)
      &                       + s3_dz(nj) * s3_dz(nk)
                 x = (lcos + tau)*winv
-                temp0 = exp(-x*x)
-        
-                H = par_lam*(1.d0 - temp0 + par_eta*x*x)
-                dHdx = 2.d0*par_lam*x*(temp0 + par_eta)
-        
-                dhdl = dHdx*winv
+                temp1 = x*x
+                temp0 = exp(-temp1)
+                H = par_lam*(1.d0 - temp0 + par_eta*temp1)
+
+                temp1 = s3_g(nj) * s3_g(nk)
+                dV3l = temp1*winv_2lam*x*(temp0 + par_eta)
         
                 
 C   three-body energy
         
-                temp1 = s3_g(nj) * s3_g(nk)
                 V3 = V3 + temp1*H
         
                 
@@ -536,8 +538,6 @@
                 
 C   (-) angular force on j and k
         
-                dV3l = temp1*dhdl
-
                 temp0 = dV3l * s3_rinv(nj)
                 fjx = fjx + temp0 * (s3_dx(nk) - lcos * s3_dx(nj))
                 fjy = fjy + temp0 * (s3_dy(nk) - lcos * s3_dy(nj))
@@ -584,8 +584,7 @@
         
 C   accumulation of 3-body coordination forces
 
-                  dVdZ_sum = dVdZ_sum
-     &                     + temp1*dHdx*(dwinv*(lcos+tau) + winv*dtau)
+                  dVdZ_sum = dVdZ_sum+dV3l*(TdW_WdT_over_W+muhalf*lcos)
         
                 end if
               end do