#if 0 

   forces-edip.c
   -------------

   Version 1.1c

   Force and Energy Calculation with the  
   Environment-Dependent Interatomic Potential 

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

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

   E-mail:
   bazant@math.mit.edu

   optimized by Xianglong Yuan, yuanx@mit.edu, April 2002


   COPYRIGHT NOTICE
   ---------------- 

   forces-edip, copyright 1997 by Martin Z. Bazant and Harvard University.
   Permission is granted to use forces-edip.c for academic use only, at
   no cost. Unauthorized sale or commerical use of this software 
   is prohibited by United States copyright law. Any publication describing 
   research involving this software should contain the following citations,
   at once and in this order, to give proper credit to the theoretical work and 
   fitting that produced EDIP and this subroutine:  
 
     1.  M. Z. Bazant and E. Kaxiras, Phys. Rev. Lett. 77, 4370 (1996).
     2.  M. Z. Bazant, E. Kaxiras, J. F. Justo, Phys. Rev. B 56, 8542 (1997). 
     3.  J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, and S. Yip, 
	   Phys. Rev. B 58, 2539 (1998).

   This software has been extensively tested for molecular dynamics simulations
   on Sun, SGI and IBM architectures, but no guarantees are made. 


   WEBSITE
   -------

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


   INTERFACE - PASSED VARIABLES
   ----------------------------

   void compute_forces_EDIP(measure,version,local,pos,f,box)

   int measure;
     If this flag is nonzero, then the average energy and coordination are measured. 

   int version;
     Not presently used. 

The next three arrays of structures contain the particle information.
The division reflects different communication roles in SPMD parallelism.
The following type declarations should appear elsewhere in the code,
altough for the force routine, only export and import are really needed:

typedef double SCALAR_t;
typedef struct{
  SCALAR_t x,y,z;
} VECTOR_t;
typedef VECTOR_t export_t;    /* r */
typedef VECTOR_t import_t;    /* f */
typedef struct{
        /* derivatives of atomic positions */
   VECTOR_t v;          /* r' */
   VECTOR_t a;          /* r'' */
#if HIGH_ORDER > 0
   VECTOR_t b;          /* r''' */
#endif
#if HIGH_ORDER > 1
   VECTOR_t c;          /* r'''' */
#endif
        /* other local atomic information - coordination and ID numbers */
   double z;
#if OLD
   int global_cell;
#endif
   int global_id;
   int layer;
} local_t;
#define MAX_PART 4097   /* maximum number of particles (no dynamic mem.alloc.) */

   local_t local[MAX_PART];
     This structure is used for "extra" information, not needed for
     force computation. On output, it contains the average EDIP coordination 
     of each atom.  
    
   export_t pos[MAX_PART];
     Contains position vectors of each particle.

   import_t f[MAX_PART];
     On output, contains the total force on each particle.

   DOMAIN_t box;
     Dimensions of the periodic box. Note that the routine takes possible
     interactions from the neighbor list. So, one can "turn off" periodic
     boundary conditions in one or more directions by changing the
     neighbor list. No change would need to be made here.

typedef struct{
  double L_x,L_y,L_z;
  double L_x_div2,L_y_div2,L_z_div2;
} DOMAIN_t;
 

   INTERFACE - GLOBAL VARIABLES  
   ----------------------------

   These are all defined below. In addition to the parameters for the potential,
   which are defined by init_EDIP(), one needs N_own, the total number of particles
   and the Verlet neighbor list. The array neighbors[] contains a list of 
   particle numbers (indices in the pos[], local[], f[] arrays). The neighbors
   of atom i are given by the block of neighbors[] between p_nbrs[i] and 
   p_nbrs[i+1].  Measurements of energy, virial,... are also returned as global
   variables. 

   The atomic arrays and neighbor list are specific to our particular MD program.
   You will need to rewrite the outermost part of the loop over atoms (loop "i") 
   to conform to your MD program. 

#endif

#include <stdio.h>
#include <math.h>

/*---- These are only needed for the "mbmd" program. Can delete. -------*/
#include "parameters.h"
#include "functions.h"
extern PARAMETERS_t parameters;
extern double cct,bct;
/*--------*/

/* macro to inline minimum image distance checking */
#define MIN_IMAGE_DISTANCE(dx,L_x,L_x_div2) ((dx)+(((dx) > (L_x_div2)) ? -(L_x) : (((dx) < -(L_x_div2)) ? (L_x) : 0)))

/* DEBUGGING FLAGS - systematically turn on various pieces of the potential. */
#define V2_on	1	
#define V2Z_on	1	
#define V3_on 	1	
#define V3g_on	1	
#define V3h_on	1		
#define V3Z_on	1 
#define Zfast	1	

/* Total number of particles */
extern int N_own;

/* Verlet neighbor list with bonds double-counted */
extern int neighbors[MAX_NBRS],p_nbrs[MAX_PART]; 

/* EDIP Si PARAMETERS Justo et al., Phys. Rev. B 58, 2539 (1998).

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

connection between these parameters and Justo et al., Phys. Rev. B 58, 2539 (1998):

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

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

*/
static double A,B,rh,sig,a;
static double lam,gam,b,c;
static double mu,Qo,bet,alp;
static double bg;   /* cutoff for g(r) */
static double palp;    /* justo prefactor for bond order - delete later */
static double delta,eta,zet;

/* tau(Z) (Ismail & Kaxiras, 1993) */
static const double u1 = -0.165799;
static const double u2 = 32.557;
static const double u3 = 0.286198;
static const double u4 = 0.66; 

/* measurements in force routine */
extern double virial,v2sum,E_potential,coord_total,V2,V3;

/********************************************/

void init_EDIP(params)
  int params;   /* Old parameter sets have been removed. This is unnecessary. */
{
    A = 5.6714030;
    B = 2.0002804;
    rh = 1.2085196;
    a = 3.1213820;
    sig = 0.5774108;
    lam = 1.4533108;
    gam = 1.1247945;
    b = 3.1213820;
    c = 2.5609104;
    delta = 78.7590539;
    mu = 0.6966326;
    Qo = 312.1341346;
    palp = 1.4074424;
    bet = 0.0070975;
    alp = 3.1083847;
    printf("\n EDIP-Si Parameters: \n");
    printf("%lf %lf %lf %lf %lf \n",A,B,rh,a,sig);
    printf("%lf %lf %lf %lf %lf \n",lam,gam,b,c,delta);
    printf("%lf %lf %lf %lf %lf \n",mu,Qo,palp,bet,alp);
    bg=a;
    eta = delta/Qo;

    /* needed for mbmd program. can delete. */
    parameters.pot_cutoff = a;
    bct=b; cct=c; 
}

/********************************************/

void compute_forces_EDIP(measure,version,local,pos,f,box)
  int measure;
  int version;
  local_t local[MAX_PART];
  export_t pos[MAX_PART];
  import_t f[MAX_PART];
  DOMAIN_t box;
{
  /*------------------------- VARIABLE DECLARATIONS -------------------------*/

  int i,j,k,l,n;
  double dx,dy,dz,r,rsqr,asqr;
  double rinv,rmainv,xinv,xinv3,den,Z,fZ;
  double dV2j,dV2ijx,dV2ijy,dV2ijz,pZ,dp;
  double temp0,temp1,temp2;
  double Qort,muhalf,u5;
  double rmbinv,winv,dwinv,tau,dtau,lcos,x,H,dHdx,dhdl;
  double dV3rij,dV3rijx,dV3rijy,dV3rijz;
  double dV3rik,dV3rikx,dV3riky,dV3rikz;
  double dV3l,dV3ljx,dV3ljy,dV3ljz,dV3lkx,dV3lky,dV3lkz;
  double dV2dZ,dxdZ,dV3dZ,dVdZ_sum;
  double dEdrl,dEdrlx,dEdrly,dEdrlz;
  double bmc,cmbinv;
  double fjx,fjy,fjz,fkx,fky,fkz;

  typedef struct{      
    double t0,t1,t2,t3;  /* various V2 functions and derivatives */
    double dx,dy,dz;     /* unit separation vector */
    double r;		 /* bond length (only needed for virial) */
  } store2;
  store2 s2[MAX_NBRS_1]; /* two-body interaction quantities, r<a */  
  int n2;                /* size of s2[] */
  int num2[MAX_NBRS_1];  /* atom ID numbers for s2[] */

  typedef struct{
    double g,dg;         /* 3-body radial function and its derivative */
    double rinv;         /* 1/r */
    double dx,dy,dz;     /* unit separation vector */
    double r;		 /* bond length (only needed for virial) */
  } store3;
  store3 s3[MAX_NBRS_1]; /* three-body interaction quantities, r<b */
  int n3;                /* size of s3[] */
  int num3[MAX_NBRS_1];  /* atom ID numbers for s3[] */
    
  typedef struct{
    double df;           /* derivative of neighbor function f'(r) */
    double dx,dy,dz;     /* unit separation vector */
    double r;		 /* bond length (only needed for virial) */
  } storez;
  storez sz[MAX_NBRS_1]; /* coordination number stuff, c<r<b */
  int nz;                /* size of sz[] */
  int numz[MAX_NBRS_1];  /* atom ID numbers for sz[] */

  int nj,nk,nl;         /* indices for the store arrays */

  /*---------------------------------------------------------------------------*/

  /* INITIALIZE FORCES AND GLOBAL SUMS */

  for(i=0;i<N_own;i++) {
    f[i].x=0.0;
    f[i].y=0.0;
    f[i].z=0.0;
  }
  if(measure) {
    coord_total=0.0;
    virial=0.0;
    V2=0.0;            /* E_potential = V2 + V3 */
    V3=0.0;
  }

  /* COMBINE COEFFICIENTS */

  asqr = a*a;
  Qort = sqrt(Qo);
  muhalf = mu*0.5;
  u5 = u2*u4;    
  bmc = b-c;
  cmbinv = 1.0/(c-b);


  /*--- LEVEL 1: OUTER LOOP OVER ATOMS ---*/

  for(i=0; i<N_own; i++) {

    /* RESET COORDINATION AND NEIGHBOR NUMBERS */

    Z = 0.0;
    n2 = 0;  
    n3 = 0;
    nz = 0;

    /*--- LEVEL 2: LOOP PREPASS OVER PAIRS ---*/

    for(n=p_nbrs[i]; n<p_nbrs[i+1]; n++) {
      j = neighbors[n];

      /* TEST IF WITHIN OUTER CUTOFF */

      dx = pos[j].x-pos[i].x;
      dx = MIN_IMAGE_DISTANCE(dx,box.L_x,box.L_x_div2); 
      if(fabs(dx) < a) {
      dy = pos[j].y-pos[i].y;
      dy = MIN_IMAGE_DISTANCE(dy,box.L_y,box.L_y_div2); 
      if(fabs(dy) < a) {
      dz = pos[j].z-pos[i].z;
      dz = MIN_IMAGE_DISTANCE(dz,box.L_z,box.L_z_div2);
      if(fabs(dz) < a) {
      rsqr = dx*dx + dy*dy + dz*dz;
      if(rsqr < asqr) {
        r = sqrt(rsqr);

        /* PARTS OF TWO-BODY INTERACTION r<a */
 
        num2[n2] = j;
        rinv = 1.0/r;
        dx *= rinv;
        dy *= rinv;
        dz *= rinv;
        rmainv = 1.0/(r-a);
        s2[n2].t0 = A*exp(sig*rmainv);
        s2[n2].t1 = pow(B*rinv,rh);
        s2[n2].t2 = rh*rinv;
        s2[n2].t3 = sig*rmainv*rmainv;
        s2[n2].dx = dx;
        s2[n2].dy = dy;
        s2[n2].dz = dz;
 	s2[n2].r = r;
        n2++;

        /* RADIAL PARTS OF THREE-BODY INTERACTION r<b */

        if(r < bg) {

	  num3[n3] = j;
	  rmbinv = 1.0/(r-bg);
	  temp1 = gam*rmbinv;
	  temp0 = exp(temp1);
#if V3g_on
	  s3[n3].g = temp0;
	  s3[n3].dg = -rmbinv*temp1*temp0;
#else
	  s3[n3].g = 1;
	  s3[n3].dg = 0;
#endif
	  s3[n3].dx = dx;
	  s3[n3].dy = dy;
	  s3[n3].dz = dz;
	  s3[n3].rinv = rinv;
 	  s3[n3].r = r;
	  n3++;

          /* COORDINATION AND NEIGHBOR FUNCTION c<r<b */

          if(r<b) {
 	   if(r < c) 
	    Z += 1.0;
	   else {
	    xinv = bmc/(r-c);
            xinv3 = xinv*xinv*xinv;
            den = 1.0/(1 - xinv3);
	    temp1 = alp*den;
	    fZ = exp(temp1);
            Z += fZ;
	    numz[nz] = j;
            sz[nz].df = fZ*temp1*den*3.0*xinv3*xinv*cmbinv;   /* df/dr */
	    sz[nz].dx = dx;
	    sz[nz].dy = dy;
	    sz[nz].dz = dz;
 	    sz[nz].r = r;
	    nz++;
	   }      
	  }
        }
      }
      }
      }
      }
    }
    local[i].z = Z;

    if(measure) coord_total += Z;

    /* ZERO ACCUMULATION VARIABLE FOR ENVIRONMENT FORCES */

    dVdZ_sum = 0.0;

    /* ENVIRONMENT-DEPENDENCE OF PAIR INTERACTION */

#if V2_on

#if V2Z_on
    temp0 = bet*Z;
    pZ = palp*exp(-temp0*Z);         /* bond order */
    dp = -2.0*temp0*pZ;         /* derivative of bond order */
#else
    pZ = palp*exp(-bet*16);
    dp=0.0;
#endif

    /*--- LEVEL 2: LOOP FOR PAIR INTERACTIONS ---*/

    for(nj=0; nj<n2; nj++) {

      temp0 = s2[nj].t1 - pZ;

      /* two-body energy V2(rij,Z) */

      if(measure)  V2 += temp0*s2[nj].t0; 

      /* two-body forces */

      dV2j = - (s2[nj].t0) * ((s2[nj].t1)*(s2[nj].t2) 
			    + temp0 * (s2[nj].t3));   /* dV2/dr */
      dV2ijx = dV2j * s2[nj].dx;
      dV2ijy = dV2j * s2[nj].dy;
      dV2ijz = dV2j * s2[nj].dz;
      f[i].x += dV2ijx;
      f[i].y += dV2ijy;
      f[i].z += dV2ijz;
      j = num2[nj];
      f[j].x -= dV2ijx;
      f[j].y -= dV2ijy;
      f[j].z -= dV2ijz;

      /* dV2/dr contribution to virial */

      if(measure) virial -= s2[nj].r * (dV2ijx*s2[nj].dx 
	+ dV2ijy*s2[nj].dy + dV2ijz*s2[nj].dz);

      /* accumulation of pair coordination forces */

      dV2dZ = - dp * s2[nj].t0;
      dVdZ_sum += dV2dZ;

    }
#endif

    /* COORDINATION-DEPENDENCE OF THREE-BODY INTERACTION */

#if V3_on

#if V3Z_on
    winv = Qort*exp(-muhalf*Z); /* inverse width of angular function */    
    dwinv = -muhalf*winv;       /* its derivative */
    temp0 = exp(-u4*Z);
    tau = u1+u2*temp0*(u3-temp0); /* -cosine of angular minimum */
    dtau = u5*temp0*(2*temp0-u3); /* its derivative */
#else
    winv =  Qort*exp(-muhalf*4);
    dwinv = 0.0;
    tau = 1.0/3.0;
    dtau=0.0;
#endif

    /*--- LEVEL 2: FIRST LOOP FOR THREE-BODY INTERACTIONS ---*/

    for(nj=0; nj<(n3-1); nj++) {

      j=num3[nj];
      
      /*--- LEVEL 3: SECOND LOOP FOR THREE-BODY INTERACTIONS ---*/

      for(nk=nj+1; nk<n3; nk++) {

	k=num3[nk];

	/* angular function h(l,Z) */

	lcos = s3[nj].dx * s3[nk].dx + s3[nj].dy * s3[nk].dy
	  + s3[nj].dz * s3[nk].dz;
	x = (lcos + tau)*winv;
        temp0 = exp(-x*x);
#if V3h_on
        H = lam*(1 - temp0 + eta*x*x);
        dHdx = 2*lam*x*(temp0 + eta);
	dhdl = dHdx*winv;
#else
 	H=1.0;
	dhdl=0.0;
#endif

	/* three-body energy */

	temp1 = s3[nj].g * s3[nk].g;
	if(measure) {
	  V3 += temp1*H; 
	}

	/* (-) radial force on atom j */

	dV3rij = s3[nj].dg * s3[nk].g * H;
	dV3rijx = dV3rij * s3[nj].dx;
	dV3rijy = dV3rij * s3[nj].dy;
	dV3rijz = dV3rij * s3[nj].dz;
	fjx = dV3rijx;
	fjy = dV3rijy;
	fjz = dV3rijz;

	/* (-) radial force on atom k */

	dV3rik = s3[nj].g * s3[nk].dg * H;
	dV3rikx = dV3rik * s3[nk].dx;
	dV3riky = dV3rik * s3[nk].dy;
	dV3rikz = dV3rik * s3[nk].dz;
	fkx = dV3rikx;
	fky = dV3riky;
	fkz = dV3rikz;

	/* (-) angular force on j */

	dV3l = temp1*dhdl;
	dV3ljx = dV3l * (s3[nk].dx - lcos * s3[nj].dx) * s3[nj].rinv;
	dV3ljy = dV3l * (s3[nk].dy - lcos * s3[nj].dy) * s3[nj].rinv;
	dV3ljz = dV3l * (s3[nk].dz - lcos * s3[nj].dz) * s3[nj].rinv;
	fjx += dV3ljx;
	fjy += dV3ljy;
	fjz += dV3ljz;

	/* (-) angular force on k */

	dV3lkx = dV3l * (s3[nj].dx - lcos * s3[nk].dx) * s3[nk].rinv;
	dV3lky = dV3l * (s3[nj].dy - lcos * s3[nk].dy) * s3[nk].rinv;
	dV3lkz = dV3l * (s3[nj].dz - lcos * s3[nk].dz) * s3[nk].rinv;
	fkx += dV3lkx;
	fky += dV3lky;
	fkz += dV3lkz;

	/* apply radial + angular forces to i, j, k */

        f[j].x -= fjx;
        f[j].y -= fjy;
        f[j].z -= fjz;
        f[k].x -= fkx;
        f[k].y -= fky;
        f[k].z -= fkz;
	f[i].x += fjx + fkx;
	f[i].y += fjy + fky;
	f[i].z += fjz + fkz;

	/* dV3/dR contributions to virial */

     	if(measure) {
	   virial -= s3[nj].r * (fjx*s3[nj].dx 
	             + fjy*s3[nj].dy + fjz*s3[nj].dz); 
	   virial -= s3[nk].r * (fkx*s3[nk].dx 
	             + fky*s3[nk].dy + fkz*s3[nk].dz); 
	}

	/* prefactor for 4-body forces from coordination */
#if V3Z_on
	dxdZ = dwinv*(lcos + tau) + winv*dtau;
	dV3dZ = temp1*dHdx*dxdZ;

	/* accumulation of 3-body coordination forces */

        dVdZ_sum += dV3dZ;
#endif
      }
    }
#endif

    /*--- LEVEL 2: LOOP TO APPLY COORDINATION FORCES ---*/

    for(nl=0; nl<nz; nl++) {

        dEdrl = dVdZ_sum * sz[nl].df;
        dEdrlx = dEdrl * sz[nl].dx;
        dEdrly = dEdrl * sz[nl].dy;
        dEdrlz = dEdrl * sz[nl].dz;
        f[i].x += dEdrlx;
        f[i].y += dEdrly;
        f[i].z += dEdrlz;
        l = numz[nl];
        f[l].x -= dEdrlx;
        f[l].y -= dEdrly;
        f[l].z -= dEdrlz;

        /* dE/dZ*dZ/dr contribution to virial */

        if(measure) virial -= sz[nl].r * (dEdrlx*sz[nl].dx 
                              + dEdrly*sz[nl].dy + dEdrlz*sz[nl].dz); 
   }
  }

  E_potential = V2+V3;
  virial /= 3.0;
}