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

C   Version 1.9f

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   optimized by Xianglong Yuan, yuanx@mit.edu, April 2002


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.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

      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
      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)
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_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.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)


C  --- LEVEL 1: OUTER LOOP OVER ATOMS ---

      do i=1, N_own
 
C   RESET COORDINATION AND NEIGHBOR NUMBERS

        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)


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)
            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


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


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
C   df/dr
                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
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.0d0

        else

          coord_total = coord_total + Z


C   ZERO ACCUMULATION VARIABLE FOR ENVIRONMENT FORCES

          dVdZ_sum = 0.d0

C   ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION

          pZ = par_palp*dexp(par_bet*Z**2)
C   bond order
          dp = (2.0d0*par_bet)*Z*pZ
C   derivative of bond order

        end if


C  --- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---

        do nj=1, n2

          temp0 = (s2_t1(nj) - pZ) * s2_t0(nj)


C   two-body energy V2(rij,Z)

          V2 = V2 + temp0

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

          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
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)
          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)
C   inverse width of angular function
          dwinv = muhalf
C   its derivative divided by itself
          temp0 = exp(-u4*Z)
          tau = u1+u2*temp0*(u3-temp0)
C   -cosine of angular minimum
          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

C  --- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---

        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)


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)


C   three-body energy

            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)


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))


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

            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

              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


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

      end do

      E_potential = V2+V3
      virial = virial / 3.0d0
      end subroutine