--- forces-edip.1.f	Wed Apr 17 21:49:57 2002
+++ forces-edip.1.1.f	Fri Apr 19 10:42:17 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.0f
+C   Version 1.1f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -24,6 +24,8 @@
 C   translated from C to FORTRAN 
 C   by Noam Bernstein, noamb@cmt.harvard.edu, December 1997
 
+C   optimized by Xianglong Yuan, yuanx@mit.edu, April 2002
+
 
 C   COPYRIGHT NOTICE
 C   ----------------
@@ -185,7 +187,7 @@
           double precision dV3rik,dV3rikx,dV3riky,dV3rikz
           double precision dV3l,dV3ljx,dV3ljy,dV3ljz,dV3lkx,
      &                            dV3lky,dV3lkz
-          double precision dV2dZ,dxdZ,dV3dZ
+          double precision dVdZ_sum
           double precision dEdrl,dEdrlx,dEdrly,dEdrlz
           double precision bmc,cmbinv
           double precision fjx,fjy,fjz,fkx,fky,fkz
@@ -217,7 +219,6 @@
 C   atom ID numbers for s3[]
         
           double precision sz_df(MAX_NBRS_1)
-          double precision sz_sum(MAX_NBRS_1)
           double precision sz_dx(MAX_NBRS_1)
           double precision sz_dy(MAX_NBRS_1)
           double precision sz_dz(MAX_NBRS_1)
@@ -406,12 +407,9 @@
               coord_total = coord_total + Z
         
               
-C   ZERO ACCUMULATION ARRAY FOR ENVIRONMENT FORCES
-        
-              do nl=1, nz-1
-                sz_sum(nl)=0.0
-              end do
-        
+C   ZERO ACCUMULATION VARIABLE FOR ENVIRONMENT FORCES
+
+              dVdZ_sum = 0.d0
               
 C   ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION
         
@@ -465,17 +463,17 @@
         
               if(fixZ .eq. 0) then
         
-              
-C  --- LEVEL 3: LOOP FOR PAIR COORDINATION FORCES ---
-        
-              dV2dZ = - dp * s2_t0(nj)
-              do nl=1, nz-1
-                 sz_sum(nl) =  sz_sum(nl) + dV2dZ
-              end do
+C  accumulation of pair coordination forces
+
+              dVdZ_sum = dVdZ_sum - s2_t0(nj)
         
               end if 
 C  fixZ
             end do
+
+C  multiply prefactor to pair coordination forces
+
+            dVdZ_sum = dp*dVdZ_sum
         
             if(fixZ .ne. 0) then
               winv = Qort*dexp(-muhalf*Z)
@@ -614,17 +612,11 @@
         
                 if(fixZ .eq. 0) then
         
-                
-C   prefactor for 4-body forces from coordination
-                  dxdZ = dwinv*(lcos + tau) + winv*dtau
-                  dV3dZ = temp1*dHdx*dxdZ
-        
-                
-C  --- LEVEL 4: LOOP FOR THREE-BODY COORDINATION FORCES ---
+C   accumulation of 3-body coordination forces
+
+                  dVdZ_sum = dVdZ_sum
+     &                     + temp1*dHdx*(dwinv*(lcos+tau) + winv*dtau)
         
-                  do nl=1, nz-1
-                    sz_sum(nl) = sz_sum(nl) + dV3dZ
-                  end do
                 end if
               end do
             end do
@@ -636,7 +628,7 @@
         
             do nl=1, nz-1
         
-                dEdrl = sz_sum(nl) * sz_df(nl)
+                dEdrl = dVdZ_sum * sz_df(nl)
                 dEdrlx = dEdrl * sz_dx(nl)
                 dEdrly = dEdrl * sz_dy(nl)
                 dEdrlz = dEdrl * sz_dz(nl)