--- forces-edip.1.5.f	Fri Apr 19 10:43:16 2002
+++ forces-edip.1.6.f	Fri Apr 19 10:43:29 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.5f
+C   Version 1.6f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -183,10 +183,7 @@
           double precision temp0,temp1,temp2
           double precision Qort,muhalf,u5
           double precision rmbinv,winv,dwinv,tau,dtau,lcos,x,H,dHdx,dhdl
-          double precision dV3rij,dV3rijx,dV3rijy,dV3rijz
-          double precision dV3rik,dV3rikx,dV3riky,dV3rikz
-          double precision dV3l,dV3ljx,dV3ljy,dV3ljz,dV3lkx,
-     &                            dV3lky,dV3lkz
+          double precision dV3rij,dV3rik,dV3l
           double precision dVdZ_sum
           double precision dEdrl,dEdrlx,dEdrly,dEdrlz
           double precision bmc,cmb3inv
@@ -524,53 +521,32 @@
                 V3 = V3 + temp1*H
         
                 
-C   (-) radial force on atom j
+C   (-) radial force on atom j and k
         
                 dV3rij = s3_dg(nj) * s3_g(nk) * H
-                dV3rijx = dV3rij * s3_dx(nj)
-                dV3rijy = dV3rij * s3_dy(nj)
-                dV3rijz = dV3rij * s3_dz(nj)
-                fjx = dV3rijx
-                fjy = dV3rijy
-                fjz = dV3rijz
-        
-                
-C   (-) radial force on atom k
+                fjx = dV3rij * s3_dx(nj)
+                fjy = dV3rij * s3_dy(nj)
+                fjz = dV3rij * s3_dz(nj)
         
                 dV3rik = s3_g(nj) * s3_dg(nk) * H
-                dV3rikx = dV3rik * s3_dx(nk)
-                dV3riky = dV3rik * s3_dy(nk)
-                dV3rikz = dV3rik * s3_dz(nk)
-                fkx = dV3rikx
-                fky = dV3riky
-                fkz = dV3rikz
+                fkx = dV3rik * s3_dx(nk)
+                fky = dV3rik * s3_dy(nk)
+                fkz = dV3rik * s3_dz(nk)
         
                 
-C   (-) angular force on j
+C   (-) angular force on j and k
         
                 dV3l = temp1*dhdl
-                dV3ljx = dV3l * (s3_dx(nk) - 
-     &                        lcos * s3_dx(nj)) * s3_rinv(nj)
-                dV3ljy = dV3l * (s3_dy(nk) - 
-     &                        lcos * s3_dy(nj)) * s3_rinv(nj)
-                dV3ljz = dV3l * (s3_dz(nk) - 
-     &                        lcos * s3_dz(nj)) * s3_rinv(nj)
-                fjx = fjx + dV3ljx
-                fjy = fjy + dV3ljy
-                fjz = fjz + dV3ljz
-        
-                
-C   (-) angular force on k
-        
-                dV3lkx = dV3l * (s3_dx(nj) - lcos * s3_dx(nk)) 
-     &                            * s3_rinv(nk)
-                dV3lky = dV3l * (s3_dy(nj) - lcos * s3_dy(nk))
-     &                            * s3_rinv(nk)
-                dV3lkz = dV3l * (s3_dz(nj) - lcos * s3_dz(nk))
-     &                            * s3_rinv(nk)
-                fkx = fkx + dV3lkx
-                fky = fky + dV3lky
-                fkz = fkz + dV3lkz
+
+                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))
+                fjz = fjz + temp0 * (s3_dz(nk) - lcos * s3_dz(nj))
+        
+                temp0 = dV3l * s3_rinv(nk)
+                fkx = fkx + temp0 * (s3_dx(nj) - lcos * s3_dx(nk)) 
+                fky = fky + temp0 * (s3_dy(nj) - lcos * s3_dy(nk))
+                fkz = fkz + temp0 * (s3_dz(nj) - lcos * s3_dz(nk))
         
                 
 C   apply radial + angular forces to i, j, k