--- forces-edip.1.9.f	Fri Apr 19 10:44:07 2002
+++ forces-edip.1.10.f	Fri Apr 19 10:44:21 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.9f
+C   Version 1.10f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -163,6 +163,12 @@
      &           E_potential, f)
       implicit none
 
+
+c#define FIXZ
+#undef FIXZ
+c   switch on/off option using fixed Z
+
+
 C  ------------------------- VARIABLE DECLARATIONS -------------------------
 
 
@@ -177,14 +183,13 @@
 
       integer i,j,k,l,n
       double precision dx,dy,dz,r,rsqr,asqr
-      double precision rinv,rmainv,xinv,xinv3,den,Z,fZ
+      double precision rinv,rmainv,Z
       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
       double precision winv_2lam,TdW_WdT_over_W
       double precision dV3rij,dV3rik,dV3l
-      double precision dVdZ_sum
       double precision dEdrl,dEdrlx,dEdrly,dEdrlz
       double precision bmc,cmb3inv
       double precision fjx,fjy,fjz,fkx,fky,fkz
@@ -215,6 +220,9 @@
       integer num3(MAX_NBRS_1)
 C   atom ID numbers for s3[]
 
+#ifndef FIXZ
+      double precision xinv,xinv3,den,fZ
+      double precision dVdZ_sum
       double precision sz_df(MAX_NBRS_1)
       double precision sz_dx(MAX_NBRS_1)
       double precision sz_dy(MAX_NBRS_1)
@@ -224,6 +232,10 @@
 C   size of sz[]
       integer numz(MAX_NBRS_1)
 C   atom ID numbers for sz[]
+#else
+      integer tricks_Zfix
+      parameter (tricks_Zfix = 5)
+#endif
 
       integer nj,nk,nl
 C   indices for the store arrays
@@ -233,8 +245,6 @@
 
       double precision coord_total
 
-      integer fixZ, tricks_Zfix
-      parameter (fixZ = 0, tricks_Zfix = 5)
 
       L_x_div_2 = L_x/2.0D0
       L_y_div_2 = L_y/2.0D0
@@ -349,8 +359,7 @@
               s3_rinv(n3) = rinv
                s3_r(n3) = r
 
-              if(fixZ .eq. 0) then
-
+#ifndef FIXZ
 
 C   COORDINATION AND NEIGHBOR FUNCTION par_c<r<par_b
 
@@ -372,9 +381,7 @@
                 sz_dz(nz) = dz
                  sz_r(nz) = r
               end if
-C  r < par_c < par_b
-              end if
-C  fixZ .eq. 0
+#endif
             end if
 C  r < par_bg
           end if
@@ -387,31 +394,31 @@
 C  dz < par_a
         end do
 
-        if(fixZ .ne. 0) then
-
-          Z = tricks_Zfix
-          coord_total = coord_total + Z
-          pZ = par_palp*dexp(par_bet*Z*Z)
-          dp = 0.0d0
 
-        else
+#ifndef FIXZ
 
-          coord_total = coord_total + Z
+        coord_total = coord_total + Z
 
 
 C   ZERO ACCUMULATION VARIABLE FOR ENVIRONMENT FORCES
 
-          dVdZ_sum = 0.d0
+        dVdZ_sum = 0.d0
 
 C   ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION
 
-          pZ = par_palp*dexp(par_bet*Z**2)
+        pZ = par_palp*dexp(par_bet*Z**2)
 C   bond order
-          dp = (2.0d0*par_bet)*Z*pZ
+        dp = (2.0d0*par_bet)*Z*pZ
 C   derivative of bond order
 
-        end if
+#else
 
+        Z = tricks_Zfix
+        coord_total = coord_total + Z
+        pZ = par_palp*dexp(par_bet*Z*Z)
+        dp = 0.0d0
+
+#endif
 
 C  --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---
 
@@ -450,40 +457,44 @@
           virial_xyz(3) = virial_xyz(3) - temp2
           virial = virial - (temp0 + temp1 + temp2)
 
-          if(fixZ .eq. 0) then
+#ifndef FIXZ
 
 C  accumulation of pair coordination forces
 
             dVdZ_sum = dVdZ_sum - s2_t0(nj)
 
-          end if
-C  fixZ
+#endif
+
         end do
 
+
+#ifndef FIXZ
+
 C  multiply prefactor to pair coordination forces
 
         dVdZ_sum = dp*dVdZ_sum
 
-        if(fixZ .ne. 0) then
-          winv = Qort*dexp(muhalf*Z)
-          dwinv = 0.0d0
-          temp0 = dexp(-u4*Z)
-          tau = u1+u2*temp0*(u3-temp0)
-          dtau = 0.0d0
-        else
-
 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
+        dwinv = muhalf
 C   its derivative divided by itself
-          temp0 = exp(-u4*Z)
-          tau = u1+u2*temp0*(u3-temp0)
+        temp0 = exp(-u4*Z)
+        tau = u1+u2*temp0*(u3-temp0)
 C   -cosine of angular minimum
-          dtau = u5*temp0*(2.d0*temp0-u3)
+        dtau = u5*temp0*(2.d0*temp0-u3)
 C   its derivative
-        end if
+
+#else
+
+        winv = Qort*dexp(muhalf*Z)
+        dwinv = 0.0d0
+        temp0 = dexp(-u4*Z)
+        tau = u1+u2*temp0*(u3-temp0)
+        dtau = 0.0d0
+
+#endif
 
         winv_2lam = 2.d0*par_lam*winv
         TdW_WdT_over_W = tau*dwinv + dtau
@@ -568,48 +579,49 @@
             virial_xyz(3) = virial_xyz(3) - temp2
             virial = virial - (temp0 + temp1 + temp2)
 
-            if(fixZ .eq. 0) then
+#ifndef FIXZ
 
 C   accumulation of 3-body coordination forces
 
-              dVdZ_sum = dVdZ_sum+dV3l*(TdW_WdT_over_W+muhalf*lcos)
+            dVdZ_sum = dVdZ_sum+dV3l*(TdW_WdT_over_W+muhalf*lcos)
+
+#endif
 
-            end if
           end do
         end do
 
-        if(fixZ .eq. 0) then
 
+#ifndef FIXZ
 
 C  --- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---
 
-          do nl=1, nz
+        do nl=1, nz
 
-            dEdrl = dVdZ_sum * sz_df(nl)
-            dEdrlx = dEdrl * sz_dx(nl)
-            dEdrly = dEdrl * sz_dy(nl)
-            dEdrlz = dEdrl * sz_dz(nl)
-            f(1,i) = f(1,i) + dEdrlx
-            f(2,i) = f(2,i) + dEdrly
-            f(3,i) = f(3,i) + dEdrlz
-            l = numz(nl)
-            f(1,l) = f(1,l) - dEdrlx
-            f(2,l) = f(2,l) - dEdrly
-            f(3,l) = f(3,l) - dEdrlz
+          dEdrl = dVdZ_sum * sz_df(nl)
+          dEdrlx = dEdrl * sz_dx(nl)
+          dEdrly = dEdrl * sz_dy(nl)
+          dEdrlz = dEdrl * sz_dz(nl)
+          f(1,i) = f(1,i) + dEdrlx
+          f(2,i) = f(2,i) + dEdrly
+          f(3,i) = f(3,i) + dEdrlz
+          l = numz(nl)
+          f(1,l) = f(1,l) - dEdrlx
+          f(2,l) = f(2,l) - dEdrly
+          f(3,l) = f(3,l) - dEdrlz
 
 
 C   dE/dZ*dZ/dr contribution to virial
 
-            temp0 = dEdrlx*sz_dx(nl)*sz_r(nl)
-            temp1 = dEdrly*sz_dy(nl)*sz_r(nl)
-            temp2 = dEdrlz*sz_dz(nl)*sz_r(nl)
-            virial_xyz(1) = virial_xyz(1) - temp0
-            virial_xyz(2) = virial_xyz(2) - temp1
-            virial_xyz(3) = virial_xyz(3) - temp2
-            virial = virial - (temp0 + temp1 + temp2)
-          end do
+          temp0 = dEdrlx*sz_dx(nl)*sz_r(nl)
+          temp1 = dEdrly*sz_dy(nl)*sz_r(nl)
+          temp2 = dEdrlz*sz_dz(nl)*sz_r(nl)
+          virial_xyz(1) = virial_xyz(1) - temp0
+          virial_xyz(2) = virial_xyz(2) - temp1
+          virial_xyz(3) = virial_xyz(3) - temp2
+          virial = virial - (temp0 + temp1 + temp2)
+        end do
 
-        end if
+#endif
 
       end do