--- forces-edip.1.1.f	Fri Apr 19 10:42:17 2002
+++ forces-edip.1.2.f	Fri Apr 19 10:42:33 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.1f
+C   Version 1.2f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -146,10 +146,10 @@
           par_b = 3.1213820
           par_c = 2.5609104
           par_delta = 78.7590539
-          par_mu = 0.6966326
+          par_mu = -0.6966326
           par_Qo = 312.1341346
           par_palp = 1.4074424
-          par_bet = 0.0070975
+          par_bet = -0.0070975
           par_alp = 3.1083847 
           print *, par_cap_A,par_cap_B,par_rh,par_a,par_sig
           print *, par_lam,par_gam,par_b,par_c,par_delta
@@ -189,7 +189,7 @@
      &                            dV3lky,dV3lkz
           double precision dVdZ_sum
           double precision dEdrl,dEdrlx,dEdrly,dEdrlz
-          double precision bmc,cmbinv
+          double precision bmc,cmb3inv
           double precision fjx,fjy,fjz,fkx,fky,fkz
         
           double precision s2_t0(MAX_NBRS_1)
@@ -265,7 +265,7 @@
           muhalf = par_mu*0.5D0
           u5 = u2*u4
           bmc = par_b-par_c
-          cmbinv = 1.0D0/(par_c-par_b)
+          cmb3inv = 3.0D0/(par_c-par_b)
         
         
           
@@ -276,9 +276,9 @@
 C   RESET COORDINATION AND NEIGHBOR NUMBERS
         
             Z = 0.0
-            n2 = 1
-            n3 = 1
-            nz = 1
+            n2 = 0
+            n3 = 0
+            nz = 0
         
             
 C  --- LEVEL 2: LOOP PREPASS OVER PAIRS ---
@@ -316,15 +316,15 @@
               rsqr = dx*dx + dy*dy + dz*dz
               if(rsqr < asqr) then
                 r = sqrt(rsqr)
-        
-                
-C   PARTS OF TWO-BODY INTERACTION r<par_a
-        
-                num2(n2) = j
                 rinv = 1.0/r
                 dx = dx * rinv
                 dy = dy * rinv
                 dz = dz * rinv
+                
+C   PARTS OF TWO-BODY INTERACTION r<par_a
+        
+                n2 = n2 + 1
+                num2(n2) = j
                 rmainv = 1.0/(r-par_a)
                 s2_t0(n2) = par_cap_A*dexp(par_sig*rmainv)
                 s2_t1(n2) = (par_cap_B*rinv)**par_rh
@@ -334,13 +334,13 @@
                 s2_dy(n2) = dy
                 s2_dz(n2) = dz
                  s2_r(n2) = r
-                n2 = n2 + 1
         
                 
 C   RADIAL PARTS OF THREE-BODY INTERACTION r<par_b
         
                 if(r < par_bg)  then
         
+                  n3 = n3 + 1
                   num3(n3) = j
                   rmbinv = 1.0/(r-par_bg)
                   temp1 = par_gam*rmbinv
@@ -352,7 +352,6 @@
                   s3_dz(n3) = dz
                   s3_rinv(n3) = rinv
                    s3_r(n3) = r
-                  n3 = n3 + 1
         
                   if(fixZ .eq. 0) then
         
@@ -363,20 +362,20 @@
                     if(r .lt. par_c) then
                     Z = Z + 1.0
                    else
+                    nz = nz + 1
+                    numz(nz) = j
                     xinv = bmc/(r-par_c)
                     xinv3 = xinv*xinv*xinv
                     den = 1.0/(1 - xinv3)
                     temp1 = par_alp*den
                     fZ = dexp(temp1)
                     Z = Z + fZ
-                    numz(nz) = j
-                    sz_df(nz) = fZ*temp1*den*3.0*xinv3*xinv*cmbinv   
+                    sz_df(nz) = fZ*temp1*den*xinv3*xinv*cmb3inv   
 C   df/dr
                     sz_dx(nz) = dx
                     sz_dy(nz) = dy
                     sz_dz(nz) = dz
                      sz_r(nz) = r
-                    nz = nz + 1
                    end if 
 C  r < par_C
                   end if 
@@ -399,7 +398,7 @@
         
               Z = tricks_Zfix
               coord_total = coord_total + Z
-              pZ = par_palp*dexp(-par_bet*Z*Z)
+              pZ = par_palp*dexp(par_bet*Z*Z)
               dp = 0.0
         
             else
@@ -413,10 +412,9 @@
               
 C   ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION
         
-              temp0 = par_bet*Z
-              pZ = par_palp*dexp(-temp0*Z)         
+              pZ = par_palp*dexp(par_bet*Z**2)         
 C   bond order
-              dp = -2.0*temp0*pZ         
+              dp = (2.0d0*par_bet)*Z*pZ         
 C   derivative of bond order
         
             end if
@@ -424,7 +422,7 @@
             
 C  --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---
         
-            do nj=1, n2-1
+            do nj=1, n2
         
               temp0 = s2_t1(nj) - pZ
         
@@ -476,7 +474,7 @@
             dVdZ_sum = dp*dVdZ_sum
         
             if(fixZ .ne. 0) then
-              winv = Qort*dexp(-muhalf*Z)
+              winv = Qort*dexp(muhalf*Z)
               dwinv = 0.0
               temp0 = dexp(-u4*Z)
               tau = u1+u2*temp0*(u3-temp0)
@@ -485,9 +483,9 @@
               
 C   COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION
         
-              winv = Qort*exp(-muhalf*Z) 
+              winv = Qort*exp(muhalf*Z) 
 C   inverse width of angular function
-              dwinv = -muhalf*winv       
+              dwinv = muhalf*winv       
 C   its derivative
               temp0 = exp(-u4*Z)
               tau = u1+u2*temp0*(u3-temp0) 
@@ -499,14 +497,14 @@
             
 C  --- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---
         
-            do nj=1, n3-2
+            do nj=1, n3-1
         
               j=num3(nj)
         
               
 C  --- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---
         
-              do nk=nj+1, n3-1
+              do nk=nj+1, n3
         
                 k=num3(nk)
         
@@ -626,7 +624,7 @@
             
 C  --- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---
         
-            do nl=1, nz-1
+            do nl=1, nz
         
                 dEdrl = dVdZ_sum * sz_df(nl)
                 dEdrlx = dEdrl * sz_dx(nl)