--- forces-edip.1.8.f	Fri Apr 19 10:43:54 2002
+++ forces-edip.1.9.f	Fri Apr 19 10:44:07 2002
@@ -2,7 +2,7 @@
 C   forces-edip.f
 C   -------------
 
-C   Version 1.8f
+C   Version 1.9f
 
 C   Force and Energy Calculation with the
 C   Environment-Dependent Interatomic Potential
@@ -109,511 +109,510 @@
 C    u4 = 0.66;
 
 C  ******************************************
-        
-        subroutine init_EDIP()
-        implicit none
+ 
+      subroutine init_EDIP()
+      implicit none
 
 C     NOTE: you do not need these header files, since you will have to
 C     customize the way that atoms are passed to the force subroutine 
 C     for your particular MD program.
-          include "edip_pot_include.h"
-          include "edip_neighbors_include.h"
-        
-          call input_edip_params()
-          par_bg=par_a
-          par_eta = par_delta/par_Qo
-          pot_cutoff = par_a
-          delta_safe = 0.2d0
-        
-        end subroutine
-        
+      include "edip_pot_include.h"
+      include "edip_neighbors_include.h"
+
+      call input_edip_params()
+      par_bg=par_a
+      par_eta = par_delta/par_Qo
+      pot_cutoff = par_a
+      delta_safe = 0.2d0
+
+      end subroutine
+
 C ******************************************
 
-        subroutine input_EDIP_params()
-        implicit none
-        
-          include "edip_pot_include.h"
-        
-             print *, "EDIP-Si parameters:"
-C         taken from Justo et al., Phys. Rev. B 58, 2539 (1998).
-          par_cap_A = 5.6714030d0
-          par_cap_B = 2.0002804d0
-            par_rh = 1.2085196d0
-          par_a = 3.1213820d0
-          par_sig = 0.5774108d0
-          par_lam = 1.4533108d0
-          par_gam = 1.1247945d0
-          par_b = 3.1213820d0
-          par_c = 2.5609104d0
-          par_delta = 78.7590539d0
-          par_mu = -0.6966326d0
-          par_Qo = 312.1341346d0
-          par_palp = 1.4074424d0
-          par_bet = -0.0070975d0
-          par_alp = 3.1083847d0
-          print *, par_cap_A,par_cap_B,par_rh,par_a,par_sig
-          print *, par_lam,par_gam,par_b,par_c,par_delta
-          print *, par_mu,par_Qo,par_palp,par_bet,par_alp
+      subroutine input_EDIP_params()
+      implicit none
+ 
+      include "edip_pot_include.h"
+ 
+      print *, "EDIP-Si parameters:"
+C   taken from Justo et al., Phys. Rev. B 58, 2539 (1998).
+      par_cap_A = 5.6714030d0
+      par_cap_B = 2.0002804d0
+      par_rh = 1.2085196d0
+      par_a = 3.1213820d0
+      par_sig = 0.5774108d0
+      par_lam = 1.4533108d0
+      par_gam = 1.1247945d0
+      par_b = 3.1213820d0
+      par_c = 2.5609104d0
+      par_delta = 78.7590539d0
+      par_mu = -0.6966326d0
+      par_Qo = 312.1341346d0
+      par_palp = 1.4074424d0
+      par_bet = -0.0070975d0
+      par_alp = 3.1083847d0
+      print *, par_cap_A,par_cap_B,par_rh,par_a,par_sig
+      print *, par_lam,par_gam,par_b,par_c,par_delta
+      print *, par_mu,par_Qo,par_palp,par_bet,par_alp
+
+      end subroutine
 
-        end subroutine
-        
 C  ******************************************
-        
-        subroutine compute_forces_EDIP(N_own, pos, L_x, L_y, L_z, 
-     &            E_potential, f)
-        implicit none
-          
+
+      subroutine compute_forces_EDIP(N_own, pos, L_x, L_y, L_z, 
+     &           E_potential, f)
+      implicit none
+
 C  ------------------------- VARIABLE DECLARATIONS -------------------------
-        
-          
-          integer N_own
-          double precision pos(3,N_own)
-          double precision E_potential
-          double precision f(3,N_own)
-          double precision L_x, L_y, L_z
-        
-          include "edip_pot_include.h"
-
-          include "edip_neighbors_include.h"
-        
-          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 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
-        
-          double precision s2_t0(MAX_NBRS_1)
-          double precision s2_t1(MAX_NBRS_1)
-          double precision s2_t2(MAX_NBRS_1)
-          double precision s2_t3(MAX_NBRS_1)
-          double precision s2_dx(MAX_NBRS_1)
-          double precision s2_dy(MAX_NBRS_1)
-          double precision s2_dz(MAX_NBRS_1)
-          double precision s2_r(MAX_NBRS_1)
-          integer n2                
+
+
+      integer N_own
+      double precision pos(3,N_own)
+      double precision E_potential
+      double precision f(3,N_own)
+      double precision L_x, L_y, L_z
+
+      include "edip_pot_include.h"
+      include "edip_neighbors_include.h"
+
+      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 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
+
+      double precision s2_t0(MAX_NBRS_1)
+      double precision s2_t1(MAX_NBRS_1)
+      double precision s2_t2(MAX_NBRS_1)
+      double precision s2_t3(MAX_NBRS_1)
+      double precision s2_dx(MAX_NBRS_1)
+      double precision s2_dy(MAX_NBRS_1)
+      double precision s2_dz(MAX_NBRS_1)
+      double precision s2_r(MAX_NBRS_1)
+      integer n2
 C   size of s2[]
-          integer num2(MAX_NBRS_1)  
+      integer num2(MAX_NBRS_1)
 C   atom ID numbers for s2[]
-        
-          double precision s3_g(MAX_NBRS_1)
-          double precision s3_dg(MAX_NBRS_1)
-          double precision s3_rinv(MAX_NBRS_1)
-          double precision s3_dx(MAX_NBRS_1)
-          double precision s3_dy(MAX_NBRS_1)
-          double precision s3_dz(MAX_NBRS_1)
-          double precision s3_r(MAX_NBRS_1)
-        
-          integer n3                
+
+      double precision s3_g(MAX_NBRS_1)
+      double precision s3_dg(MAX_NBRS_1)
+      double precision s3_rinv(MAX_NBRS_1)
+      double precision s3_dx(MAX_NBRS_1)
+      double precision s3_dy(MAX_NBRS_1)
+      double precision s3_dz(MAX_NBRS_1)
+      double precision s3_r(MAX_NBRS_1)
+
+      integer n3
 C   size of s3[]
-          integer num3(MAX_NBRS_1)  
+      integer num3(MAX_NBRS_1)
 C   atom ID numbers for s3[]
-        
-          double precision sz_df(MAX_NBRS_1)
-          double precision sz_dx(MAX_NBRS_1)
-          double precision sz_dy(MAX_NBRS_1)
-          double precision sz_dz(MAX_NBRS_1)
-          double precision sz_r(MAX_NBRS_1)
-          integer nz                
+
+      double precision sz_df(MAX_NBRS_1)
+      double precision sz_dx(MAX_NBRS_1)
+      double precision sz_dy(MAX_NBRS_1)
+      double precision sz_dz(MAX_NBRS_1)
+      double precision sz_r(MAX_NBRS_1)
+      integer nz
 C   size of sz[]
-          integer numz(MAX_NBRS_1)  
+      integer numz(MAX_NBRS_1)
 C   atom ID numbers for sz[]
-        
-          integer nj,nk,nl         
+
+      integer nj,nk,nl
 C   indices for the store arrays
-        
-          double precision V2, V3, virial, virial_xyz(3)
-          double precision L_x_div_2, L_y_div_2, L_z_div_2
-
-          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
-          L_z_div_2 = L_z/2.0D0
-        
-          do i=1, N_own
-            f(1,i) = 0.0d0
-            f(2,i) = 0.0d0
-            f(3,i) = 0.0d0
-          end do
-        
-          coord_total=0.0d0
-          virial=0.0d0
-          virial_xyz(1)= 0.0d0
-          virial_xyz(2)= 0.0d0
-          virial_xyz(3)= 0.0d0
-          V2=0.0d0
-          V3=0.0d0
-        
-          
+
+      double precision V2, V3, virial, virial_xyz(3)
+      double precision L_x_div_2, L_y_div_2, L_z_div_2
+
+      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
+      L_z_div_2 = L_z/2.0D0
+
+      do i=1, N_own
+        f(1,i) = 0.0d0
+        f(2,i) = 0.0d0
+        f(3,i) = 0.0d0
+      end do
+
+      coord_total=0.0d0
+      virial=0.0d0
+      virial_xyz(1)= 0.0d0
+      virial_xyz(2)= 0.0d0
+      virial_xyz(3)= 0.0d0
+      V2=0.0d0
+      V3=0.0d0
+
+
 C   COMBINE COEFFICIENTS
-        
-          asqr = par_a*par_a
-          Qort = sqrt(par_Qo)
-          muhalf = par_mu*0.5D0
-          u5 = u2*u4
-          bmc = par_b-par_c
-          cmb3inv = 3.0D0/(par_c-par_b)
-        
-        
-          
+
+      asqr = par_a*par_a
+      Qort = sqrt(par_Qo)
+      muhalf = par_mu*0.5D0
+      u5 = u2*u4
+      bmc = par_b-par_c
+      cmb3inv = 3.0D0/(par_c-par_b)
+
+
 C  --- LEVEL 1: OUTER LOOP OVER ATOMS ---
-        
-          do i=1, N_own
-            
+
+      do i=1, N_own
+ 
 C   RESET COORDINATION AND NEIGHBOR NUMBERS
-        
-            Z = 0.0d0
-            n2 = 0
-            n3 = 0
-            nz = 0
-        
-            
+
+        Z = 0.0d0
+        n2 = 0
+        n3 = 0
+        nz = 0
+
+
 C  --- LEVEL 2: LOOP PREPASS OVER PAIRS ---
-        
-            do n=p_nbrs(i), p_nbrs(i+1)-1
-              j = neighbors(n)
-        
-              
+
+        do n=p_nbrs(i), p_nbrs(i+1)-1
+          j = neighbors(n)
+
+
 C   TEST IF WITHIN OUTER CUTOFF
-        
-              dx = pos(1,j)-pos(1,i)
-              if (dx .gt. L_x_div_2) then
-                dx = dx - L_x
-              else if (dx .lt. -L_x_div_2) then
-                dx = dx + L_x
-              end if 
+
+          dx = pos(1,j)-pos(1,i)
+          if (dx .gt. L_x_div_2) then
+            dx = dx - L_x
+          else if (dx .lt. -L_x_div_2) then
+            dx = dx + L_x
+          end if
 C  dx periodic
-              if(dabs(dx) < par_a) then
-              dy = pos(2,j)-pos(2,i)
-              if (dy .gt. L_y_div_2) then
-                dy = dy - L_y
-              else if (dy .lt. -L_y_div_2) then
-                dy = dy + L_y
-              end if 
+          if(dabs(dx) < par_a) then
+          dy = pos(2,j)-pos(2,i)
+          if (dy .gt. L_y_div_2) then
+            dy = dy - L_y
+          else if (dy .lt. -L_y_div_2) then
+            dy = dy + L_y
+          end if
 C  dy periodic
-              if(dabs(dy) < par_a) then
-              dz = pos(3,j)-pos(3,i)
-              if (dz .gt. L_z_div_2) then
-                dz = dz - L_z
-              else if (dz .lt. -L_z_div_2) then
-                dz = dz + L_z
-              end if 
+          if(dabs(dy) < par_a) then
+          dz = pos(3,j)-pos(3,i)
+          if (dz .gt. L_z_div_2) then
+            dz = dz - L_z
+          else if (dz .lt. -L_z_div_2) then
+            dz = dz + L_z
+          end if
 C  dy periodic
-              if(dabs(dz) < par_a) then
-              rsqr = dx*dx + dy*dy + dz*dz
-              if(rsqr < asqr) then
-                r = sqrt(rsqr)
-                rinv = 1.0d0/r
-                dx = dx * rinv
-                dy = dy * rinv
-                dz = dz * rinv
-                
+          if(dabs(dz) < par_a) then
+          rsqr = dx*dx + dy*dy + dz*dz
+          if(rsqr < asqr) then
+            r = sqrt(rsqr)
+            rinv = 1.0d0/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.0d0/(r-par_a)
-                s2_t0(n2) = par_cap_A*dexp(par_sig*rmainv)
-                s2_t1(n2) = (par_cap_B*rinv)**par_rh
-                s2_t2(n2) = par_rh*rinv
-                s2_t3(n2) = par_sig*rmainv*rmainv
-                s2_dx(n2) = dx
-                s2_dy(n2) = dy
-                s2_dz(n2) = dz
-                 s2_r(n2) = r
-        
-                
+
+            n2 = n2 + 1
+            num2(n2) = j
+            rmainv = 1.0d0/(r-par_a)
+            s2_t0(n2) = par_cap_A*dexp(par_sig*rmainv)
+            s2_t1(n2) = (par_cap_B*rinv)**par_rh
+            s2_t2(n2) = par_rh*rinv
+            s2_t3(n2) = par_sig*rmainv*rmainv
+            s2_dx(n2) = dx
+            s2_dy(n2) = dy
+            s2_dz(n2) = dz
+             s2_r(n2) = r
+
+
 C   RADIAL PARTS OF THREE-BODY INTERACTION r<par_b
-        
-                if(r < par_bg)  then
-        
-                  n3 = n3 + 1
-                  num3(n3) = j
-                  rmbinv = 1.0d0/(r-par_bg)
-                  temp1 = par_gam*rmbinv
-                  temp0 = dexp(temp1)
-                  s3_g(n3) = temp0
-                  s3_dg(n3) = -rmbinv*temp1*temp0
-                  s3_dx(n3) = dx
-                  s3_dy(n3) = dy
-                  s3_dz(n3) = dz
-                  s3_rinv(n3) = rinv
-                   s3_r(n3) = r
-        
-                  if(fixZ .eq. 0) then
-        
-                  
+
+            if(r < par_bg)  then
+
+              n3 = n3 + 1
+              num3(n3) = j
+              rmbinv = 1.0d0/(r-par_bg)
+              temp1 = par_gam*rmbinv
+              temp0 = dexp(temp1)
+              s3_g(n3) = temp0
+              s3_dg(n3) = -rmbinv*temp1*temp0
+              s3_dx(n3) = dx
+              s3_dy(n3) = dy
+              s3_dz(n3) = dz
+              s3_rinv(n3) = rinv
+               s3_r(n3) = r
+
+              if(fixZ .eq. 0) then
+
+
 C   COORDINATION AND NEIGHBOR FUNCTION par_c<r<par_b
-        
-                  if(r .lt. par_c) then
-                    Z = Z + 1.0d0
-                  else if(r .lt. par_b) then
-                    nz = nz + 1
-                    numz(nz) = j
-                    xinv = bmc/(r-par_c)
-                    xinv3 = xinv*xinv*xinv
-                    den = 1.0d0/(1.d0 - xinv3)
-                    temp1 = par_alp*den
-                    fZ = dexp(temp1)
-                    Z = Z + fZ
-                    sz_df(nz) = fZ*temp1*den*xinv3*xinv*cmb3inv   
+
+              if(r .lt. par_c) then
+                Z = Z + 1.0d0
+              else if(r .lt. par_b) then
+                nz = nz + 1
+                numz(nz) = j
+                xinv = bmc/(r-par_c)
+                xinv3 = xinv*xinv*xinv
+                den = 1.0d0/(1.d0 - xinv3)
+                temp1 = par_alp*den
+                fZ = dexp(temp1)
+                Z = Z + fZ
+                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
-                  end if 
+                sz_dx(nz) = dx
+                sz_dy(nz) = dy
+                sz_dz(nz) = dz
+                 sz_r(nz) = r
+              end if
 C  r < par_c < par_b
-                  end if 
+              end if
 C  fixZ .eq. 0
-                end if 
+            end if
 C  r < par_bg
-              end if 
+          end if
 C  rsqr < asqr
-              end if 
+          end if
 C  dz < par_a
-              end if 
+          end if
 C  dy < par_a
-              end if 
+          end if
 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
-        
-              coord_total = coord_total + Z
-        
-              
+        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
+
+          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
-        
-            
+
+        end if
+
+
 C  --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---
-        
-            do nj=1, n2
-        
-              temp0 = (s2_t1(nj) - pZ) * s2_t0(nj)
-        
-              
+
+        do nj=1, n2
+
+          temp0 = (s2_t1(nj) - pZ) * s2_t0(nj)
+
+
 C   two-body energy V2(rij,Z)
-        
-              V2 = V2 + temp0
-              
+
+          V2 = V2 + temp0
+
 C   two-body forces
-        
-              dV2j = - s2_t0(nj)*s2_t1(nj)*s2_t2(nj) - temp0*s2_t3(nj)
+
+          dV2j = - s2_t0(nj)*s2_t1(nj)*s2_t2(nj) - temp0*s2_t3(nj)
 C   dV2/dr
-              dV2ijx = dV2j * s2_dx(nj)
-              dV2ijy = dV2j * s2_dy(nj)
-              dV2ijz = dV2j * s2_dz(nj)
-              f(1,i) = f(1,i) + dV2ijx
-              f(2,i) = f(2,i) + dV2ijy
-              f(3,i) = f(3,i) + dV2ijz
-              j = num2(nj)
-              f(1,j) = f(1,j) - dV2ijx
-              f(2,j) = f(2,j) - dV2ijy
-              f(3,j) = f(3,j) - dV2ijz
-        
-              
+          dV2ijx = dV2j * s2_dx(nj)
+          dV2ijy = dV2j * s2_dy(nj)
+          dV2ijz = dV2j * s2_dz(nj)
+          f(1,i) = f(1,i) + dV2ijx
+          f(2,i) = f(2,i) + dV2ijy
+          f(3,i) = f(3,i) + dV2ijz
+          j = num2(nj)
+          f(1,j) = f(1,j) - dV2ijx
+          f(2,j) = f(2,j) - dV2ijy
+          f(3,j) = f(3,j) - dV2ijz
+
+
 C   dV2/dr contribution to virial
-        
-              temp0 = dV2ijx*s2_dx(nj)*s2_r(nj)
-              temp1 = dV2ijy*s2_dy(nj)*s2_r(nj)
-              temp2 = dV2ijz*s2_dz(nj)*s2_r(nj)
-              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)
-        
-              if(fixZ .eq. 0) then
-        
+
+          temp0 = dV2ijx*s2_dx(nj)*s2_r(nj)
+          temp1 = dV2ijy*s2_dy(nj)*s2_r(nj)
+          temp2 = dV2ijz*s2_dz(nj)*s2_r(nj)
+          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)
+
+          if(fixZ .eq. 0) then
+
 C  accumulation of pair coordination forces
 
-              dVdZ_sum = dVdZ_sum - s2_t0(nj)
-        
-              end if 
+            dVdZ_sum = dVdZ_sum - s2_t0(nj)
+
+          end if
 C  fixZ
-            end do
+        end do
 
 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
-              
+        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
-        
-            winv_2lam = 2.d0*par_lam*winv
-            TdW_WdT_over_W = tau*dwinv + dtau
-            
+        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 ---
-        
-            do nj=1, n3-1
-        
-              j=num3(nj)
-        
-              
+
+        do nj=1, n3-1
+
+          j=num3(nj)
+
+
 C  --- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---
-        
-              do nk=nj+1, n3
-        
-                k=num3(nk)
-        
-                
+
+          do nk=nj+1, n3
+
+            k=num3(nk)
+
+
 C   angular function h(l,Z)
-        
-                lcos = s3_dx(nj) * s3_dx(nk) + s3_dy(nj) * s3_dy(nk)
-     &                       + s3_dz(nj) * s3_dz(nk)
-                x = (lcos + tau)*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)
-        
-                
+
+            lcos = s3_dx(nj)*s3_dx(nk) + s3_dy(nj)*s3_dy(nk)
+     &           + s3_dz(nj)*s3_dz(nk)
+            x = (lcos + tau)*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
-        
-                V3 = V3 + temp1*H
-        
-                
+
+            V3 = V3 + temp1*H
+
+
 C   (-) radial force on atom j and k
-        
-                dV3rij = s3_dg(nj) * s3_g(nk) * H
-                fjx = dV3rij * s3_dx(nj)
-                fjy = dV3rij * s3_dy(nj)
-                fjz = dV3rij * s3_dz(nj)
-        
-                dV3rik = s3_g(nj) * s3_dg(nk) * H
-                fkx = dV3rik * s3_dx(nk)
-                fky = dV3rik * s3_dy(nk)
-                fkz = dV3rik * s3_dz(nk)
-        
-                
+
+            dV3rij = s3_dg(nj) * s3_g(nk) * H
+            fjx = dV3rij * s3_dx(nj)
+            fjy = dV3rij * s3_dy(nj)
+            fjz = dV3rij * s3_dz(nj)
+
+            dV3rik = s3_g(nj) * s3_dg(nk) * H
+            fkx = dV3rik * s3_dx(nk)
+            fky = dV3rik * s3_dy(nk)
+            fkz = dV3rik * s3_dz(nk)
+
+
 C   (-) angular force on j and k
-        
-                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))
-        
-                
+
+            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
-        
-                f(1,j) = f(1,j) - fjx
-                f(2,j) = f(2,j) - fjy
-                f(3,j) = f(3,j) - fjz
-                f(1,k) = f(1,k) - fkx
-                f(2,k) = f(2,k) - fky
-                f(3,k) = f(3,k) - fkz
-                f(1,i) = f(1,i) + fjx + fkx
-                f(2,i) = f(2,i) + fjy + fky
-                f(3,i) = f(3,i) + fjz + fkz
-                
+
+            f(1,j) = f(1,j) - fjx
+            f(2,j) = f(2,j) - fjy
+            f(3,j) = f(3,j) - fjz
+            f(1,k) = f(1,k) - fkx
+            f(2,k) = f(2,k) - fky
+            f(3,k) = f(3,k) - fkz
+            f(1,i) = f(1,i) + fjx + fkx
+            f(2,i) = f(2,i) + fjy + fky
+            f(3,i) = f(3,i) + fjz + fkz
+
 C   dV3/dR contributions to virial
 
-                temp0 = fjx*s3_dx(nj)*s3_r(nj) + fkx*s3_dx(nk)*s3_r(nk)
-                temp1 = fjy*s3_dy(nj)*s3_r(nj) + fky*s3_dy(nk)*s3_r(nk)
-                temp2 = fjz*s3_dz(nj)*s3_r(nj) + fkz*s3_dz(nk)*s3_r(nk)
-                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)
-        
-                if(fixZ .eq. 0) then
-        
-C   accumulation of 3-body coordination forces
+            temp0 = fjx*s3_dx(nj)*s3_r(nj) + fkx*s3_dx(nk)*s3_r(nk)
+            temp1 = fjy*s3_dy(nj)*s3_r(nj) + fky*s3_dy(nk)*s3_r(nk)
+            temp2 = fjz*s3_dz(nj)*s3_r(nj) + fkz*s3_dz(nk)*s3_r(nk)
+            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)
 
-                  dVdZ_sum = dVdZ_sum+dV3l*(TdW_WdT_over_W+muhalf*lcos)
-        
-                end if
-              end do
-            end do
-        
             if(fixZ .eq. 0) then
-        
-            
+
+C   accumulation of 3-body coordination forces
+
+              dVdZ_sum = dVdZ_sum+dV3l*(TdW_WdT_over_W+muhalf*lcos)
+
+            end if
+          end do
+        end do
+
+        if(fixZ .eq. 0) then
+
+
 C  --- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---
-        
-            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
-        
-                
+
+          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
+
+
 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
-        
-           end if
+
+            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
-        
-          E_potential = V2+V3
-          virial = virial / 3.0d0
-        end subroutine
+
+        end if
+
+      end do
+
+      E_potential = V2+V3
+      virial = virial / 3.0d0
+      end subroutine