C   forces-edip.f
C   -------------

C   Version 1.0f

C   Force and Energy Calculation with the
C   Environment-Dependent Interatomic Potential

C   written by Martin Z. Bazant, 
C   Department of Physics, Harvard University
C   April - October 1997
C   (based on forces.c by Martin Z. Bazant, June 1994)

c   New address (2000):
c   Professor Martin Z. Bazant
c   Department of Mathematics 2-363B
c   Massachusetts Institute of Technology
c   Cambridge, MA 02139-4307

c   E-mail:
c   bazant@math.mit.edu

C   translated from C to FORTRAN 
C   by Noam Bernstein, noamb@cmt.harvard.edu, December 1997


C   COPYRIGHT NOTICE
C   ----------------

C   forces-edip, copyright 1997 by Martin Z. Bazant and Harvard University.
C   Permission is granted to use forces-edip for academic use only, at
C   no cost. Unauthorized sale or commerical use of this software
C   is prohibited by United States copyright law. Any publication describing
C   research involving this software should contain the following citations,
C   at once and in this order, to give proper credit to the theoretical work and
C   fitting that produced EDIP and this subroutine:

C     1.  M. Z. Bazant and E. Kaxiras, Phys. Rev. Lett. 77, 4370 (1996).
C     2.  M. Z. Bazant, E. Kaxiras, J. F. Justo, Phys. Rev. B 56, 8542 (1997).
C     3.  J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, and S. Yip, 
C	    Phys. Rev. B 58, 2539 (1998).
C
C   This software has been extensively tested for molecular dynamics simulations
C   on Sun, SGI and IBM architectures, but no guarantees are made.


C   WEBSITE
C   -------

C   Updated versions of this software are available at the EDIP distribution site,
C   http://pelion.eas.harvard.edu/software/EDIP/
C   Postscript files of related papers are available at the Kaxiras group web site
C   in the Department of Physics at Harvard University, http://pelion.eas.harvard.edu,
C   under 'Empirical Methods'.
C   A description of the algorithm used in this subroutine can be found
C   in the Ph.D. Thesis of M. Z. Bazant (1997), chapter 6, on the web at
C   http://pelion.eas.harvard.edu/~bazant/thesis/


C   INTERFACE
C   ---------

C     compute_forces_EDIP(N, p, lx, ly, lz, E, f)
C     N : number of particles
C     p : array (3,N) of positions in Angstroms
C     E : returned energy in eV
C     f : returned forces array (3,N) in eV/Angstroms
C
C     neighbors(p_nbrs(i)),...,neighbors(p_nbrs(i+1)) are the 
C     atoms which are "neighbors" of atom i, using a standard Verlet 
C     neighbor list. These are a global arrays, not passed, that are declared
C     in "edip_neighbors_include.h". This way of storing atomic positions
C     is not unique, and will require a customized patch to the main MD program.
C
C     The parameters of the potential initialized in input_EDIP_params() are global
C     variables declared in "edip_pot_include.h".
C

C   PARAMETERS
C   ---------- 

C    par_cap_A,par_cap_B,par_rh,par_a,par_sig
C    par_lam,par_gam,par_b,par_c,par_delta
C    par_mu,par_Qo,par_palp,par_bet,par_alp

C    5.6714030     2.0002804     1.2085196     3.1213820     0.5774108
C    1.4533108     1.1247945     3.1213820     2.5609104    78.7590539
C    0.6966326   312.1341346     1.4074424     0.0070975     3.1083847

C    Connection between these parameters and those given in the paper,
C    Justo et al., Phys. Rev. B 58, 2539 (1998):

C    A((B/r)**rh-palp*exp(-bet*Z*Z)) = A'((B'/r)**rh-exp(-bet*Z*Z))

C    so in the paper (')
C    A' = A*palp
C    B' = B * palp**(-1/rh)
C    eta = detla/Qo

C    Non-adjustable parameters for tau(Z) from Ismail & Kaxiras, 1993,
C    also in Bazant, Kaxiras, Justo, PRB (1997):

C    u1 = -0.165799;
C    u2 = 32.557;
C    u3 = 0.286198;
C    u4 = 0.66;

C  ******************************************
        
        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.2
        
        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.6714030
	  par_cap_B = 2.0002804
  	  par_rh = 1.2085196 
          par_a = 3.1213820
          par_sig = 0.5774108
	  par_lam = 1.4533108
	  par_gam = 1.1247945
	  par_b = 3.1213820
	  par_c = 2.5609104
          par_delta = 78.7590539
	  par_mu = 0.6966326
          par_Qo = 312.1341346
          par_palp = 1.4074424
	  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
          print *, par_mu,par_Qo,par_palp,par_bet,par_alp

        end subroutine
        
C  ******************************************
        
        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,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 dV2dZ,dxdZ,dV3dZ
          double precision dEdrl,dEdrlx,dEdrly,dEdrlz
          double precision bmc,cmbinv
          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)  
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                
C   size of s3[]
          integer num3(MAX_NBRS_1)  
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)
          double precision sz_r(MAX_NBRS_1)
          integer nz                
C   size of sz[]
          integer numz(MAX_NBRS_1)  
C   atom ID numbers for sz[]
        
          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.0
            f(2,i) = 0.0
            f(3,i) = 0.0
          end do
        
          coord_total=0.0
          virial=0.0
          virial_xyz(1)= 0.0
          virial_xyz(2)= 0.0
          virial_xyz(3)= 0.0
          V2=0.0
          V3=0.0
        
          
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
          cmbinv = 1.0D0/(par_c-par_b)
        
        
          
C  --- LEVEL 1: OUTER LOOP OVER ATOMS ---
        
          do i=1, N_own
            
C   RESET COORDINATION AND NEIGHBOR NUMBERS
        
            Z = 0.0
            n2 = 1
            n3 = 1
            nz = 1
        
            
C  --- LEVEL 2: LOOP PREPASS OVER PAIRS ---
        
            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 
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 
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 
C  dy periodic
              if(dabs(dz) < par_a) then
              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
                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
                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
        
                
C   RADIAL PARTS OF THREE-BODY INTERACTION r<par_b
        
                if(r < par_bg)  then
        
        	  num3(n3) = j
        	  rmbinv = 1.0/(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
        	  n3 = n3 + 1
        
        	  if(fixZ .eq. 0) then
        
                  
C   COORDINATION AND NEIGHBOR FUNCTION par_c<r<par_b
        
                  if(r .lt. par_b) then
         	   if(r .lt. par_c) then
        	    Z = Z + 1.0
        	   else
        	    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   
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 
C  r < par_b
        	  end if 
C  fixZ .eq. 0
                end if 
C  r < par_bg
              end if 
C  rsqr < asqr
              end if 
C  dz < par_a
              end if 
C  dy < par_a
              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.0
        
            else
        
              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   ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION
        
              temp0 = par_bet*Z
              pZ = par_palp*dexp(-temp0*Z)         
C   bond order
              dp = -2.0*temp0*pZ         
C   derivative of bond order
        
            end if
        
            
C  --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---
        
            do nj=1, n2-1
        
              temp0 = s2_t1(nj) - pZ
        
              
C   two-body energy V2(rij,Z)
        
              V2 = V2 + temp0*s2_t0(nj)
              
C   two-body forces
        
              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
        
              
C   dV2/dr contribution to virial
        
              virial_xyz(1) = virial_xyz(1) - 
     &				s2_r(nj)*(dV2ijx*s2_dx(nj))
              virial_xyz(2) = virial_xyz(2) - 
     &				s2_r(nj)*(dV2ijy*s2_dy(nj))
              virial_xyz(3) = virial_xyz(3) - 
     &				s2_r(nj)*(dV2ijz*s2_dz(nj))
              virial = virial - s2_r(nj) * (dV2ijx*s2_dx(nj)+ 
     &        	dV2ijy*s2_dy(nj) + dV2ijz*s2_dz(nj))
        
              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
        
              end if 
C  fixZ
            end do
        
            if(fixZ .ne. 0) then
              winv = Qort*dexp(-muhalf*Z)
              dwinv = 0.0
              temp0 = dexp(-u4*Z)
              tau = u1+u2*temp0*(u3-temp0)
              dtau = 0.0
            else
              
C   COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION
        
              winv = Qort*exp(-muhalf*Z) 
C   inverse width of angular function
              dwinv = -muhalf*winv       
C   its derivative
              temp0 = exp(-u4*Z)
              tau = u1+u2*temp0*(u3-temp0) 
C   -cosine of angular minimum
              dtau = u5*temp0*(2*temp0-u3) 
C   its derivative
            end if
        
            
C  --- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---
        
            do nj=1, n3-2
        
              j=num3(nj)
        
              
C  --- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---
        
              do nk=nj+1, n3-1
        
        	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
                temp0 = exp(-x*x)
        
        	H = par_lam*(1 - temp0 + par_eta*x*x)
        	dHdx = 2*par_lam*x*(temp0 + par_eta)
        
        	dhdl = dHdx*winv
        
        	
C   three-body energy
        
        	temp1 = s3_g(nj) * s3_g(nk)
        	V3 = V3 + temp1*H
        
        	
C   (-) radial force on atom j
        
        	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
        
        	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
        
        	
C   (-) angular force on j
        
        	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
        
        	
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
        	
C   dV3/dR contributions to virial
        
        	 virial = virial - s3_r(nj) * (fjx*s3_dx(nj)  
     &        		   + fjy*s3_dy(nj) + fjz*s3_dz(nj))
        	 virial = virial - s3_r(nk) * (fkx*s3_dx(nk)  
     &        		   + fky*s3_dy(nk) + fkz*s3_dz(nk))
        	 virial_xyz(1) = virial_xyz(1) - 
     &		    s3_r(nj)*(fjx*s3_dx(nj))
        	 virial_xyz(2) = virial_xyz(2) - 
     &		         s3_r(nj)*(fjy*s3_dy(nj))
        	 virial_xyz(3) = virial_xyz(3) - 
     &		         s3_r(nj)*(fjz*s3_dz(nj))
        	 virial_xyz(1) = virial_xyz(1) - 
     &		    s3_r(nk)*(fkx*s3_dx(nk))
        	 virial_xyz(2) = virial_xyz(2) - 
     &		    s3_r(nk)*(fky*s3_dy(nk))
        	 virial_xyz(3) = virial_xyz(3) - 
     &		    s3_r(nk)*(fkz*s3_dz(nk))
        
                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 ---
        
        	  do nl=1, nz-1
        	    sz_sum(nl) = sz_sum(nl) + dV3dZ
        	  end do
                end if
              end do
            end do
        
            if(fixZ .eq. 0) then
        
            
C  --- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---
        
            do nl=1, nz-1
        
                dEdrl = sz_sum(nl) * 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
        
        	virial = virial - sz_r(nl) * (dEdrlx*sz_dx(nl)+ 
     &				dEdrly*sz_dy(nl) + dEdrlz*sz_dz(nl))
        	virial_xyz(1) = virial_xyz(1) - 
     &				sz_r(nl)*(dEdrlx*sz_dx(nl))
        	virial_xyz(2) = virial_xyz(2) - 
     &				sz_r(nl)*(dEdrly*sz_dy(nl))
        	virial_xyz(3) = virial_xyz(3) - 
     &				sz_r(nl)*(dEdrlz*sz_dz(nl))
            end do
        
           end if
          end do
        
          E_potential = V2+V3
          virial = virial / 3.0
        end subroutine