/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing authors: Tomas Oppelstrup, LLNL (oppelstrup2@llnl.gov)
   and John Moriarty, LLNL (moriarty2@llnl.gov)

   Fast MGPT algorithm developed by Tomas Oppelstrup (2015) based on the
   matrix MGPT v4.4 FORTRAN routine of John Moriarty (2006) as converted
   to C++ for LAMMPS application by Jamie Marian and Alexander Stukowski
   (2011).  See LLNL copyright notice at bottom of this file.
------------------------------------------------------------------------- */

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cassert>

#include "pair_mgpt.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"


using namespace LAMMPS_NS;

//#define TIMING_ON

#ifdef TIMING_ON
#include <sys/time.h>
#include <time.h>
//#include "rdtsc.h"
#ifdef __bgq__
#include <hwi/include/bqc/A2_inlines.h>
#endif

static double gettime(int x = 0) {
  if(1) {
    /*
      struct timeval tv;
      gettimeofday(&tv,NULL);
      return tv.tv_sec + 1e-6 * tv.tv_usec;
    */
    /*
      const double x = 1.0 / CLOCKS_PER_SEC;
      return clock() * x;
    */

    //const double invfreq = 1.0 / 2394.108e6;
    /*
    const double invfreq = 1.0 / 700e6;
      unsigned long long int x = rdtsc();
      return x*invfreq;
    */

    const double invfreq = 1.0 / 1.6e9;
    unsigned long long int x = GetTimeBase();
    return x*invfreq;


  } else
    return 0.0;
}
#else
static double gettime(int /*x*/ = 0) { return 0.0; }
#endif


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

PairMGPT::PairMGPT(LAMMPS *lmp) : Pair(lmp)
{
	single_enable = 0;
	one_coeff = 1;
	ghostneigh = 1;
}

PairMGPT::~PairMGPT()
{
	if(allocated) {
		memory->destroy(setflag);
		memory->destroy(cutsq);
		memory->destroy(cutghost);
	}
}

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


static double t_make_b2 = 0.0,n_make_b2 = 0.0;

template<typename intype,typename outtype,int ni,int nj> void fmatconv(intype *array) {
  outtype *cast = (outtype *) array;
  for(int i = 0; i<ni; i++)
    for(int j = 0; j<nj; j++)
      cast[i*nj+j] = array[i*nj+j];
}

void PairMGPT::make_bond(const double xx[][3],int i,int j,bond_data *bptr) {
  double rrij[3],rij;
  int p;

  double t0,t1;

  /* Check that alignment requirements for SIMD code are fullfilled */
  assert( (((unsigned long long int) (bptr->H.m )) & 31) == 0 );
  assert( (((unsigned long long int) (bptr->Hx.m)) & 31) == 0 );
  assert( (((unsigned long long int) (bptr->Hy.m)) & 31) == 0 );
  assert( (((unsigned long long int) (bptr->Hz.m)) & 31) == 0 );

  rij = 0.0;
  for(p = 0; p<3; p++) {
    rrij[p] = xx[i][p] - xx[j][p];
    rij = rij + rrij[p]*rrij[p];
  }

  /* Zero all matrix elements */
  for(i = 0; i<8; i++)
    for(j = 0; j<8; j++) {
      bptr->H.m[i][j] = 0.0;
      bptr->Hx.m[i][j] = 0.0;
      bptr->Hy.m[i][j] = 0.0;
      bptr->Hz.m[i][j] = 0.0;
      bptr->Hz.m[j][i] = 0.0;
    }

  if(rij <= rcrit*rcrit) {
    t0 = gettime();
    if(lang == 3) {
      hamltn_5_raw(rrij[0],rrij[1],rrij[2],
		   bptr->H.m ,bptr->Hx.m,
		   bptr->Hy.m,bptr->Hz.m,&bptr->fl_deriv_sum);
    } else {
      hamltn_7_raw(rrij[0],rrij[1],rrij[2],
		   bptr->H.m ,bptr->Hx.m,
		   bptr->Hy.m,bptr->Hz.m,&bptr->fl_deriv_sum);
    }

    t1 = gettime();
    t_make_b2 += t1-t0;
    n_make_b2++;
  } else {
    bptr->fl_deriv_sum = 0.0;
  }

  if(linalg.single) {
    fmatconv<double,float,7,8>(&(bptr->H.m[1][0]));
    fmatconv<double,float,7,8>(&(bptr->Hx.m[1][0]));
    fmatconv<double,float,7,8>(&(bptr->Hy.m[1][0]));
    fmatconv<double,float,7,8>(&(bptr->Hz.m[1][0]));
  }
}

static double t_trace = 0.0,n_trace = 0.0;
/*
static inline double mtrace(int n,double A[8][8],double B[8][8]) {
  double t0,t1;
  double s;

  t0 = gettime();
  if(n == 5) s = mtrace_5(A,B);
  else if(n == 7) s = mtrace_7(A,B);
  else {
    s = 0.0;
    for(int i = 1; i<=n; i++)
      for(int j = 1; j<=n; j++)
	s = s + A[i][j]*B[i][j];
  }
  t1 = gettime();
  t_trace += t1-t0;
  n_trace++;

  return s;
}
*/

void PairMGPT::make_triplet(bond_data *ij_bond,bond_data *ik_bond,
			     triplet_data *triptr) {
  if(1) {
    const trmul_fun tr_mul = linalg.tr_mul;
    tr_mul(&(ij_bond->H.m[1][0]), &(ik_bond->H.m[1][0]) ,&(triptr->H1H2.m[1][0]) );
    tr_mul(&(ij_bond->Hx.m[1][0]),&(ik_bond->H.m[1][0]) ,&(triptr->H1xH2.m[1][0]));
    tr_mul(&(ij_bond->Hy.m[1][0]),&(ik_bond->H.m[1][0]) ,&(triptr->H1yH2.m[1][0]));
    tr_mul(&(ij_bond->Hz.m[1][0]),&(ik_bond->H.m[1][0]) ,&(triptr->H1zH2.m[1][0]));
    tr_mul(&(ij_bond->H.m[1][0]) ,&(ik_bond->Hx.m[1][0]),&(triptr->H1H2x.m[1][0]));
    tr_mul(&(ij_bond->H.m[1][0]) ,&(ik_bond->Hy.m[1][0]),&(triptr->H1H2y.m[1][0]));
    tr_mul(&(ij_bond->H.m[1][0]) ,&(ik_bond->Hz.m[1][0]),&(triptr->H1H2z.m[1][0]));
  } else {
    transprod(ij_bond->H, ik_bond->H ,triptr->H1H2 );
    transprod(ij_bond->Hx,ik_bond->H ,triptr->H1xH2);
    transprod(ij_bond->Hy,ik_bond->H ,triptr->H1yH2);
    transprod(ij_bond->Hz,ik_bond->H ,triptr->H1zH2);
    transprod(ij_bond->H ,ik_bond->Hx,triptr->H1H2x);
    transprod(ij_bond->H ,ik_bond->Hy,triptr->H1H2y);
    transprod(ij_bond->H ,ik_bond->Hz,triptr->H1H2z);
  }
}

static double t_make_t = 0.0,t_make_b = 0.0,n_make = 0.0;

PairMGPT::triplet_data *PairMGPT::get_triplet(const double xx[][3],int i,int j,int k,
						Hash<bond_data,Doublet> *bhash,
						triplet_data *twork,
						double *dvir_ij_p,double *dvir_ik_p) {
  const int recompute = 0;
  static bond_data bij_work,bik_work;

  double t0,t1;

  bond_data *bij = 0,*bik = 0;
  triplet_data *tptr = 0;

  t0 = gettime();
  if(recompute == 0) {
    bij = bhash->Lookup(Doublet(i,j));
    bik = bhash->Lookup(Doublet(i,k));
  }

  if(bij == 0) {
    if(recompute == 0)
      bij = bhash->Insert(Doublet(i,j));
    else
      bij = &bij_work;
    if(i < j)
      make_bond(xx,i,j,bij);
    else
      make_bond(xx,j,i,bij);
  }

  if(bik == 0) {
    if(recompute == 0)
      bik = bhash->Insert(Doublet(i,k));
    else
      bik = &bik_work;
    if(i < k)
      make_bond(xx,i,k,bik);
    else
      make_bond(xx,k,i,bik);
  }
  t1 = gettime();
  t_make_b += t1-t0;

  t0 = gettime();
  if(bij != 0 && bij != 0) {
    tptr = twork;
    make_triplet(bij,bik,tptr);
    *dvir_ij_p = bij->fl_deriv_sum;
    *dvir_ik_p = bik->fl_deriv_sum;
  } else {
    *dvir_ij_p = 0.0;
    *dvir_ik_p = 0.0;
  }
  t1 = gettime();
  t_make_t += t1-t0;
  n_make++;
  return tptr;
}


double PairMGPT::numderiv3t(double xx[][3],int i,int j,int k,int p) {
  static bond_data Bij,Bjk,Bki;
  const double delta = 1e-5;

  const double xsave = xx[i][p];
  double e1,e2;

  const double vc = splinepot.vc;

  xx[i][p] = xsave + delta;
  make_bond(xx,i,j,&Bij);
  make_bond(xx,j,k,&Bjk);
  make_bond(xx,k,i,&Bki);
  e1 = trace(prodmat(Bij.H,Bjk.H),Bki.H) * (vc/anorm3);

  xx[i][p] = xsave - delta;
  make_bond(xx,i,j,&Bij);
  if(0) { /* This bond doesn't change when i is perturbed */
    make_bond(xx,j,k,&Bjk);
  }
  make_bond(xx,k,i,&Bki);
  e2 = trace(prodmat(Bij.H,Bjk.H),Bki.H) * (vc/anorm3);

  xx[i][p] = xsave;
  return (e1 - e2)/(2.0*delta);
}

double PairMGPT::numderiv3v(double xx[][3],int i,int j,int k,int p,int ipert) {
  static bond_data Bij,Bik;
  const double delta = 1e-5;

  const double xsave = xx[ipert][p];
  double e1,e2;

  const double vd = splinepot.vd;

  xx[ipert][p] = xsave + delta;
  make_bond(xx,i,j,&Bij);
  make_bond(xx,i,k,&Bik);
  e1 = trace(prodmat(Bij.H,Bij.H),prodmat(Bik.H,Bik.H)) * (vd/anorm4);

  xx[ipert][p] = xsave - delta;
  make_bond(xx,i,j,&Bij);
  make_bond(xx,i,k,&Bik);
  e2 = trace(prodmat(Bij.H,Bij.H),prodmat(Bik.H,Bik.H)) * (vd/anorm4);

  xx[ipert][p] = xsave;
  return (e1 - e2)/(2.0*delta);
}


double PairMGPT::numderiv4(double xx[][3],int i,int j,int k,int m,int p) {
  static bond_data Bij,Bjk,Bkm,Bmi;
  const double delta = 1e-5;

  const double xsave = xx[i][p];
  double e1,e2;

  const double ve = splinepot.ve;

  xx[i][p] = xsave + delta;
  make_bond(xx,i,j,&Bij);
  make_bond(xx,j,k,&Bjk);
  make_bond(xx,k,m,&Bkm);
  make_bond(xx,m,i,&Bmi);
  e1 = trace(prodmat(Bij.H,Bjk.H),prodmat(Bkm.H,Bmi.H)) * (ve/anorm4);

  xx[i][p] = xsave - delta;
  make_bond(xx,i,j,&Bij);
  if(0) { /* Only the i coordinates changed... */
    make_bond(xx,j,k,&Bjk);
    make_bond(xx,k,m,&Bkm);
  }
  make_bond(xx,m,i,&Bmi);
  e2 = trace(prodmat(Bij.H,Bjk.H),prodmat(Bkm.H,Bmi.H)) * (ve/anorm4);

  xx[i][p] = xsave;
  return (e1 - e2)/(2.0*delta);
}

static double dtol = 1e-6;
void PairMGPT::force_debug_3t(double xx[][3],
		    int i0,int j0,int k0,
		    int i ,int j ,int k ,
		    double dfix,double dfiy,double dfiz,
		    double dfjx,double dfjy,double dfjz,
		    double dfkx,double dfky,double dfkz) {
  double dfi[3],dfj[3],dfk[3];
  dfi[0] = dfix; dfi[1] = dfiy; dfi[2] = dfiz;
  dfj[0] = dfjx; dfj[1] = dfjy; dfj[2] = dfjz;
  dfk[0] = dfkx; dfk[1] = dfky; dfk[2] = dfkz;

  for(int p = 0; p<3; p++) {
    /* Compute numerical derivatives by displacing atoms i,j,k */
    double ndfi,ndfj,ndfk;
    ndfi = -numderiv3t(xx,i,j,k,p);
    ndfj = -numderiv3t(xx,j,k,i,p);
    ndfk = -numderiv3t(xx,k,i,j,p);

    if((fabs(dfi[p] - ndfi) > dtol &&
	fabs(dfi[p] - ndfi) > dtol*fabs(ndfi)) ||
       (fabs(dfj[p] - ndfj) > dtol &&
	fabs(dfj[p] - ndfj) > dtol*fabs(ndfj)) ||
       (fabs(dfk[p] - ndfk) > dtol &&
	fabs(dfk[p] - ndfk) > dtol*fabs(ndfk))) {
      printf("Force error in T12 & T23 & T31 :: i,j,k = %d,%d,%d\n",i0,j0,k0);
      printf("    dE/d%c[i] = %20.10e    %20.10e\n", 'x'+p,ndfi, dfi[p]);
      printf("    dE/d%c[j] = %20.10e    %20.10e\n", 'x'+p,ndfj, dfj[p]);
      printf("    dE/d%c[k] = %20.10e    %20.10e\n", 'x'+p,ndfk, dfk[p]);
      printf("\n");
    }
  }
}

void PairMGPT::force_debug_3v(double xx[][3],
		    int i0,int j0,int k0,
		    int i ,int j ,int k ,
		    double dfix,double dfiy,double dfiz,
		    double dfjx,double dfjy,double dfjz,
		    double dfkx,double dfky,double dfkz) {
  double dfi[3],dfj[3],dfk[3];
  dfi[0] = dfix; dfi[1] = dfiy; dfi[2] = dfiz;
  dfj[0] = dfjx; dfj[1] = dfjy; dfj[2] = dfjz;
  dfk[0] = dfkx; dfk[1] = dfky; dfk[2] = dfkz;

  for(int p = 0; p<3; p++) {
    /* Compute numerical derivatives by displacing atoms i,j,k */
    double ndfi,ndfj,ndfk;
    ndfi = -numderiv3v(xx,i,j,k,p,i0);
    ndfj = -numderiv3v(xx,i,j,k,p,j0);
    ndfk = -numderiv3v(xx,i,j,k,p,k0);

    if((fabs(dfi[p] - ndfi) > dtol &&
	fabs(dfi[p] - ndfi) > dtol*fabs(ndfi)) ||
       (fabs(dfj[p] - ndfj) > dtol &&
	fabs(dfj[p] - ndfj) > dtol*fabs(ndfj)) ||
       (fabs(dfk[p] - ndfk) > dtol &&
	fabs(dfk[p] - ndfk) > dtol*fabs(ndfk))) {
      printf("Force error in T12 :: i,j,k = %d,%d,%d\n",i0,j0,k0);
      printf("    dE/d%c[i] = %20.10e    %20.10e\n", 'x'+p,ndfi, dfi[p]);
      printf("    dE/d%c[j] = %20.10e    %20.10e\n", 'x'+p,ndfj, dfj[p]);
      printf("    dE/d%c[k] = %20.10e    %20.10e\n", 'x'+p,ndfk, dfk[p]);
      printf("\n");
    }
  }
}

void PairMGPT::force_debug_4(double xx[][3],
		   int i0,int j0,int k0,int m0,
		   int i ,int j ,int k ,int m ,
		   double dfix,double dfiy,double dfiz,
		   double dfjx,double dfjy,double dfjz,
		   double dfkx,double dfky,double dfkz,
		   double dfmx,double dfmy,double dfmz) {

  double dfi[3],dfj[3],dfk[3],dfm[3];
  dfi[0] = dfix; dfi[1] = dfiy; dfi[2] = dfiz;
  dfj[0] = dfjx; dfj[1] = dfjy; dfj[2] = dfjz;
  dfk[0] = dfkx; dfk[1] = dfky; dfk[2] = dfkz;
  dfm[0] = dfmx; dfm[1] = dfmy; dfm[2] = dfmz;

  const int ii0[] = {i0,j0,k0,m0},ii[] = {i,j,k,m,i,j,k};

  for(int p = 0; p<3; p++) {
    /* Compute numerical derivatives by displacing atoms i,j,k,m */
    double ndfi,ndfj,ndfk,ndfm;
    if(1) {
      double ndf[] = {0.0,0.0,0.0,0.0};
      for(int s = 0; s<4; s++)
	for(int t = 0; t<4; t++)
	  if(ii[s] == ii0[t])
	    ndf[t] = -numderiv4(xx,ii[s],ii[s+1],ii[s+2],ii[s+3],p);
      ndfi = ndf[0]; ndfj = ndf[1];
      ndfk = ndf[2]; ndfm = ndf[3];
    } else {
      ndfi = -numderiv4(xx,i,j,k,m,p);
      ndfj = -numderiv4(xx,j,k,m,i,p);
      ndfk = -numderiv4(xx,k,m,i,j,p);
      ndfm = -numderiv4(xx,m,i,j,k,p);
    }

    if((fabs(dfi[p] - ndfi) > dtol &&
	fabs(dfi[p] - ndfi) > dtol*fabs(ndfi)) ||
       (fabs(dfj[p] - ndfj) > dtol &&
	fabs(dfj[p] - ndfj) > dtol*fabs(ndfj)) ||
       (fabs(dfk[p] - ndfk) > dtol &&
	fabs(dfk[p] - ndfk) > dtol*fabs(ndfk)) ||
       (fabs(dfm[p] - ndfm) > dtol &&
	fabs(dfm[p] - ndfm) > dtol*fabs(ndfm))) {
      printf("Force error in T31 & T64 :: i,j,k,m = %d,%d,%d,%d\n",i0,j0,k0,m0);
      printf("    dE/d%c[i] = %20.10e    %20.10e\n", 'x'+p,ndfi, dfi[p]);
      printf("    dE/d%c[j] = %20.10e    %20.10e\n", 'x'+p,ndfj, dfj[p]);
      printf("    dE/d%c[k] = %20.10e    %20.10e\n", 'x'+p,ndfk, dfk[p]);
      printf("    dE/d%c[m] = %20.10e    %20.10e\n", 'x'+p,ndfm, dfm[p]);
      printf("\n");
    }
  }
}



/*
#define trd_update_4(T12,T45,coord)                   \
  do {                                                \
    trd1 = transtrace(T12->H1##coord##H2,T45->H1H2 ); \
    trd2 = transtrace(T12->H1H2##coord,T45->H1H2   ); \
    trd3 = transtrace(T12->H1H2 ,T45->H1##coord##H2); \
    trd4 = transtrace(T12->H1H2 ,T45->H1H2##coord  ); \
  } while(0)
*/
#define trd_update_4(T12,T45) \
  do {                                         \
    tr_trace3(&(T45->H1H2.m[1][0]),            \
	      &(T12->H1xH2.m[1][0]),&utr1x.d,  \
	      &(T12->H1yH2.m[1][0]),&utr1y.d,  \
	      &(T12->H1zH2.m[1][0]),&utr1z.d); \
    tr_trace3(&(T45->H1H2.m[1][0]),            \
	      &(T12->H1H2x.m[1][0]),&utr2x.d,  \
	      &(T12->H1H2y.m[1][0]),&utr2y.d,  \
	      &(T12->H1H2z.m[1][0]),&utr2z.d); \
    tr_trace3(&(T12->H1H2.m[1][0]),            \
	      &(T45->H1xH2.m[1][0]),&utr3x.d,  \
	      &(T45->H1yH2.m[1][0]),&utr3y.d,  \
	      &(T45->H1zH2.m[1][0]),&utr3z.d); \
    tr_trace3(&(T12->H1H2.m[1][0]),            \
	      &(T45->H1H2x.m[1][0]),&utr4x.d,  \
	      &(T45->H1H2y.m[1][0]),&utr4y.d,  \
	      &(T45->H1H2z.m[1][0]),&utr4z.d); \
    if(linalg.single) {                        \
      trd1x = utr1x.f; trd2x = utr2x.f; trd3x = utr3x.f; trd4x = utr4x.f; \
      trd1y = utr1y.f; trd2y = utr2y.f; trd3y = utr3y.f; trd4y = utr4y.f; \
      trd1z = utr1z.f; trd2z = utr2z.f; trd3z = utr3z.f; trd4z = utr4z.f; \
    } else {                                                              \
      trd1x = utr1x.d; trd2x = utr2x.d; trd3x = utr3x.d; trd4x = utr4x.d; \
      trd1y = utr1y.d; trd2y = utr2y.d; trd3y = utr3y.d; trd4y = utr4y.d; \
      trd1z = utr1z.d; trd2z = utr2z.d; trd3z = utr3z.d; trd4z = utr4z.d; \
    }                                          \
  } while(0)

#define dfix_update_4a(coord) \
  do { \
    dfi##coord = ( (-sij)*trd1##coord + (-sim)*trd3##coord ) * (ve / anorm4); \
    dfj##coord = ( ( sij)*trd1##coord + (-sjk)*trd2##coord ) * (ve / anorm4); \
    dfk##coord = ( ( sjk)*trd2##coord + (-skm)*trd4##coord ) * (ve / anorm4); \
    dfm##coord = ( ( sim)*trd3##coord + ( skm)*trd4##coord ) * (ve / anorm4); \
  } while(0)


#define dfix_update_4b(coord) \
  do { \
    dfi##coord = ( ( ski)*trd1##coord + (-sim)*trd3##coord ) * (ve / anorm4); \
    dfj##coord = ( (-sjk)*trd2##coord + (-sjm)*trd4##coord ) * (ve / anorm4); \
    dfk##coord = ( (-ski)*trd1##coord + ( sjk)*trd2##coord ) * (ve / anorm4); \
    dfm##coord = ( ( sim)*trd3##coord + ( sjm)*trd4##coord ) * (ve / anorm4); \
  } while(0);

#define dfix_update_4c(coord) \
  do { \
    dfi##coord = ( (-sij)*trd1##coord + ( ski)*trd2##coord ) * (ve / anorm4); \
    dfj##coord = ( ( sij)*trd1##coord + (-sjm)*trd3##coord ) * (ve / anorm4); \
    dfk##coord = ( (-ski)*trd2##coord + (-skm)*trd4##coord ) * (ve / anorm4); \
    dfm##coord = ( ( sjm)*trd3##coord + ( skm)*trd4##coord ) * (ve / anorm4); \
  } while(0);

#define accumulate_forces_2(w) \
  do { \
    fix = fix + dfix*(w);    \
    fiy = fiy + dfiy*(w);    \
    fiz = fiz + dfiz*(w);    \
                             \
    fjx = fjx + dfjx*(w);    \
    fjy = fjy + dfjy*(w);    \
    fjz = fjz + dfjz*(w);    \
  } while(0)

#define accumulate_forces_3(w) \
  do { \
    accumulate_forces_2(w); \
    fkx = fkx + dfkx*(w);    \
    fky = fky + dfky*(w);    \
    fkz = fkz + dfkz*(w);    \
  } while(0)

#define accumulate_forces_4(w) \
  do { \
    accumulate_forces_3(w); \
    fmx = fmx + dfmx*(w);    \
    fmy = fmy + dfmy*(w);    \
    fmz = fmz + dfmz*(w);    \
  } while(0)



#define restrict __restrict__
#ifdef __bg__
#define const
#endif
static int ntr_calls = 0;
static trtrace3_fun tr_internal;
static void tr_count(const double * restrict A,
		     const double * restrict B1,double * restrict t1,
		     const double * restrict B2,double * restrict t2,
		     const double * restrict B3,double * restrict t3) {
  tr_internal(A,B1,t1,B2,t2,B3,t3);
  ntr_calls++;
}
#ifdef __bg__
#undef const
#endif
#undef restrict


int PairMGPT::Matrix::sz;
void PairMGPT::compute_x(const int *nnei,const int * const *nlist,
			  double *e_s,double *e_p,double *e_t,double *e_q,
			  int evflag,int newton_pair) {
  Hash<bond_data,Doublet> bond_hash(100000);
  int i,j,k,m,ix,jx,kx,mx,itag,jtag,p;

  double e_single,e_pair,e_triplet,e_triplet_c,e_quad;
  double volvir2;

  double nbc = 0.0,tbl = 0.0,tbm = 0.0;
  const int lmax_local = lmax;

  //if(evflag) printf("##### ev flag is set... wasting cycles...\n");

  *e_s = -99.0;
  *e_p = -99.0;
  *e_t = -99.0;
  *e_q = -99.0;

  double t0,t1;

  t0 = gettime(1);
  e_single = e_pair = e_triplet = e_triplet_c = e_quad = 0.0;
  volvir2 = 0.0;

  t_make_t = t_make_b = t_make_b2 = t_trace = 0.0;
  n_make = n_make_b2 = n_trace =  0.0;

  double tx0,tx1,tsort = 0.0,tpair = 0.0,tlookup = 0.0;
  double ttriplet = 0.0,tquad = 0.0,tmem = 0.0;
  double ntsort = 0.0,ntpair = 0.0,ntlookup = 0.0;
  double nttriplet = 0.0,ntquad = 0.0,ntmem = 0.0,ntquaditer = 0.0;
  double mcount = 0.0,mcount2 = 0.0, qcount = 0.0;

  double fix,fjx,fkx,fmx,dfix,dfjx,dfkx,dfmx;
  double fiy,fjy,fky,fmy,dfiy,dfjy,dfky,dfmy;
  double fiz,fjz,fkz,fmz,dfiz,dfjz,dfkz,dfmz;

  double fsave[4][3] = { {0.0} } /* {{0.0}} is to get rid of uninitialized use warning */;

  //const int numerical_pair_forces = (nbody_flag/16)%2;
  const int pair_forces   = (nbody_flag/2)%2,three_body_forces   = (nbody_flag/4)%2,four_body_forces   = (nbody_flag/8)%2;
  const int pair_energies = (nbody_flag/2)%2,three_body_energies = (nbody_flag/4)%2,four_body_energies = (nbody_flag/8)%2;
  const int single_energies = nbody_flag%2;

  const int triplet_debug = 0,quad_debug = 0;

  /* Energy and force scale factor for unit conversion. */
  const double e_scale = 0.5;

#ifdef NEIGHMASK
#define NIDX(x) (x)
#else
#define NIDX(x) ((x) & NEIGHMASK)
#endif

  int nneitot,*first,*nlist_short;
  double w2,w3,w4;
  triplet_data T12work,T23work,T31work,T45work,T56work,T64work;
  triplet_data *T12,*T23,*T31,*T45,*T56,*T64;
  int c_ij,c_jk,c_ki,c_im,c_jm,c_km;
  int mi,mj,mk;
  double tr0,tr1,tr2,tr3;
  double v33,v43;
  double rcut2_pair = rmax*rmax,rcut2_bond = rcrit*rcrit,rij2;
  int ntot,nloc;

  double dvir_ij,dvir_jk,dvir_ki,dvir_im,dvir_jm,dvir_km;
  double vir3t = 0.0,vir3v = 0.0,vir4 = 0.0;

  double (*xx)[3],(*ff)[3],(*ss)[3];

#ifdef TIMING_ON
  tr_internal = linalg.tr_trace; ntr_calls = 0;
  const trtrace3_fun tr_trace3 = tr_count;
#else
  const trtrace3_fun tr_trace3 = linalg.tr_trace;
#endif
  union {
    double d;
    float f;
  } utr1x,utr2x,utr3x,utr4x,utr1y,utr2y,utr3y,utr4y,utr1z,utr2z,utr3z,utr4z;
  double trd1x,trd2x,trd3x,trd4x;
  double trd1y,trd2y,trd3y,trd4y;
  double trd1z,trd2z,trd3z,trd4z;


  tx0 = gettime();

  double rhoinv;
  {
    double vtot = 1.0;
    double ntot = atom->natoms;
    for(i = 0; i<3; i++)
      vtot = vtot * (domain->boxhi[i] - domain->boxlo[i]);
    rhoinv = vtot / ntot;
  }

  /* Make sure triplet data work area is aligned and zeroed out. */ {
    assert(T12work.align_check() == 0);
    assert(T23work.align_check() == 0);
    assert(T31work.align_check() == 0);
    assert(T45work.align_check() == 0);
    assert(T56work.align_check() == 0);
    assert(T64work.align_check() == 0);

    T12work.zero();  T23work.zero();  T31work.zero();
    T45work.zero();  T56work.zero();  T64work.zero();
  }

  ntot = atom->nlocal + atom->nghost;
  nloc = atom->nlocal;
  //printf("[%3d] Allocating local array, size is %d atoms...\n",comm->me,j);
  xx = (double (*)[3]) memory->smalloc(sizeof(double [3]) * ntot,"mgpt: local position vector.");
  ff = (double (*)[3]) memory->smalloc(sizeof(double [3]) * ntot,"mgpt: local force vector.");
  //printf("[%3d] Initializing arrays...\n",comm->me);

  const int triclinic = domain->triclinic;
  double alpha[3] = {0.0,0.0,0.0};
  if(triclinic) {
    double E[3][3],EX[3][3];
    int cyc[] = {0,1,2,0,1};

    ss = (double (*)[3]) memory->smalloc(sizeof(double [3]) * ntot,
					 "mgpt: local reduced coordinate vector.");

    for(i = 0; i<3; i++) {
      for(j = 0; j<3; j++)
	E[i][j] = 0.0;
      E[i][i] = domain->subhi_lamda[i] - domain->sublo_lamda[i];
      domain->lamda2x(E[i],EX[i]);
    }
    for(i = 0; i<3; i++) {
      int i1 = cyc[i+1],i2 = cyc[i+2];
      double dot = 0.0,ns2 = 0.0;
      for(j = 0; j<3; j++) {
	int j1 = cyc[j+1],j2 = cyc[j+2];
	double cj = EX[i1][j1]*EX[i2][j2] - EX[i1][j2]*EX[i2][j1];
	ns2 = ns2 + cj*cj;
	dot = dot + EX[i][j]*cj;
      }
      alpha[i] = E[i][i] / (dot/sqrt(ns2));
      if(comm->me == 0) {
	static int count = 0;
	if(count < 3)
	  printf("@@@ alpha(%d) = %15.5e\n",i+1,alpha[i]);
	count++;
      }
      if(alpha[i] < 0.0) alpha[i] = -alpha[i];
    }
  } else
    ss = xx;

  nneitot = 0;
  for(ix = 0; ix<ntot; ix++) {
    for(p = 0; p<3; p++) {
      xx[ix][p] = atom->x[ix][p];
      ff[ix][p] = 0.0;
    }
    if(triclinic)
      domain->x2lamda(xx[ix],ss[ix]);
    nneitot = nneitot + nnei[ix];
  }

  first = (int *) memory->smalloc(sizeof(int) * (ntot+1),"mgpt: first");
  nlist_short = (int *) memory->smalloc(sizeof(int) * nneitot,"mgpt: nlist_short");

  tx1 = gettime();
  tmem += tx1-tx0;
  ntmem++;

  //printf("[%3d] Starting calculation...\n",comm->me);


  fix = fjx = fkx = fmx = 0.0;
  fiy = fjy = fky = fmy = 0.0;
  fiz = fjz = fkz = fmz = 0.0;

  int c_p = 0, c_t = 0, c_q = 0;

  if(0)
    if(domain->triclinic) {
      if(comm->me == 0)
	printf("Can not handle triclinic box yet\n");
      error->all(__FILE__,__LINE__,"Can not handle triclinic cell with mgpt yet.");
    }

  /*
  for(i = 0; i<nloc; i++) {
    printf("Atom %3d:: %10.3f  %10.3f  %10.3f\n",
	   i,xx[i][0],xx[i][1],xx[i][2]);
  }
  */

  first[0] = 0;
  for(i = 0; i<ntot; i++) {
    fix = fiy = fiz = 0.0;

    first[i+1] = first[i];

    const int c1 = c1_outside(ss[i],triclinic,alpha);

    tx0 = gettime();
    for(jx = 0; jx<nnei[i]; jx++) {
      fjx = fjy = fjz = 0.0;

      j = NIDX( nlist[i][jx] );

      rij2 = 0.0;
      for(p = 0; p<3; p++) {
	double t = xx[i][p] - xx[j][p];
	rij2 = rij2 + t*t;
      }

      if(c1 == 0 && rij2 < rcut2_pair) {
	if(j < i) {
	  w2 = get_weight(triclinic,ss[i],ss[j]);

	  if(w2 > 0.0) {
	    /*
	      Compute pair energy/force
	    */
	    double de_pair,df,rij = sqrt(rij2);
	    splinepot.eval_pot(rij,&de_pair,&df);
	    de_pair = de_pair * e_scale * w2;
	    df = df / rij * w2;


	    if(pair_energies == 0) de_pair = 0.0;
	    e_pair = e_pair + de_pair;
	    c_p++;

	    if(pair_forces == 0) df = 0.0;

	    if(volpres_flag && pair_energies) {
	      double dvir;
	      splinepot.eval_vir(rij,&dvir);
	      volvir2 = volvir2 - dvir * w2;

	      /* Per-atom virial contribution of volumetric energy term */
	      if(vflag_atom)
		for(int pp = 0; pp<3; pp++) {
		  //virial[i] = virial[i] + rhoinv*e_scale*volvir2;
		  vatom[i][pp] -= 0.5 * rhoinv*e_scale*dvir*w2;
		  vatom[j][pp] -= 0.5 * rhoinv*e_scale*dvir*w2;
		}
	    }

	    double drijx = xx[j][0] - xx[i][0];
	    double drijy = xx[j][1] - xx[i][1];
	    double drijz = xx[j][2] - xx[i][2];

	    fix = fix + df*drijx;
	    fjx = fjx - df*drijx;

	    fiy = fiy + df*drijy;
	    fjy = fjy - df*drijy;

	    fiz = fiz + df*drijz;
	    fjz = fjz - df*drijz;

	    if(evflag) {
	      //ev_tally(i,j,nloc,newton_pair,de_pair,0.0,df,-drijx,-drijy,-drijz);
	      /* To fix stress-per-atom scaling, and sign */
	      ev_tally(i,j,nloc,newton_pair,de_pair,0.0,-df * e_scale,-drijx,-drijy,-drijz);
	    }

	    ff[j][0] += fjx * e_scale;
	    ff[j][1] += fjy * e_scale;
	    ff[j][2] += fjz * e_scale;

	  }
	}
      }

      if(rij2 < rcut2_bond && c2_outside(ss[i],ss[j],triclinic,alpha) == 0) {
	/*
	  Add j to short neighbor list for i.
	  Insert j to keep list sorted.
	*/

	p = first[i+1]-1;
	while(p >= first[i] && nlist_short[p] > j) {
	  nlist_short[p+1] = nlist_short[p];
	  p = p - 1;
	}
	nlist_short[p+1] = j;
	first[i+1] = first[i+1] + 1;
	if(first[i+1] > nneitot) {
	  printf("nneitot = %d, short list full. i=%d\n",
		 nneitot,i);
	  error->one(__FILE__,__LINE__,"Shit! Short list full\n");
	}

      }
    }

    ff[i][0] += fix * e_scale;
    ff[i][1] += fiy * e_scale;
    ff[i][2] += fiz * e_scale;

    tx1 = gettime();
    tpair += tx1-tx0;
    ntpair += nnei[i];
  }

  for(i = 0; i<ntot; i++) {
    fix = fiy = fiz = 0.0;

    /*
      Use short lists for triplets and quadruplets.
      For open (2-bonded) triplets, can only use k<j, but not k<i.
      For closed (3-bonded) triplets, we can assume k<j<i.

      Quadruplets:
      Always use k<j<i, and require m<i.
      If 5-bonded with im bond, ignore the quadruplet.
      If 6-bonded, require m<k.

      For 4-bonded quadruplets, we can still use k<j, but also
      assume max(m,j)<i
    */

    if(three_body_energies || three_body_forces ||
       four_body_energies || four_body_forces)
      for(jx = first[i]; jx<first[i+1]; jx++) {
	fjx = fjy = fjz = 0.0;

	j = nlist_short[jx];

	for(kx = first[i]; kx<jx; kx++) {
	  fkx = fky = fkz = 0.0;

	  k = nlist_short[kx];

	  /*
	    Search lists of j and k, and see if
	    1) j is in k-list (closed triplet)
	    2) j and k have a common neighbor (closed quadruplet)
	  */

	  c_ij = c_ki = 1;

	  const int sij = (i < j) ? 1 : -1;
	  const int sjk = (j < k) ? 1 : -1;
	  const int ski = (k < i) ? 1 : -1;


	  T12 = T23 = T31 = 0;

	  mj = first[j];
	  /*
	    Since i is in the j-list, and i > k and the list
	    is sorted, the loop below terminates:-)
	  */
	  while(mj < first[j+1] && nlist_short[mj] < k) mj = mj + 1;
	  if(mj < first[j+1] && nlist_short[mj] == k) {
	    /* Closed triplet */
	    c_jk = 1;

	    if(j > i) continue; /* Require k<j<i for closed triplets */
	  } else {
	    /* Open triplet */
	    c_jk = 0;
	  }

	  tx0 = gettime();

	  w3 = get_weight(triclinic,ss[i],ss[j],ss[k]);

	  int triplet_defer;
	  if(w3 > 0.0) {
	    triplet_defer = 0;

	    dvir_ij = dvir_jk = dvir_ki = 0.0;
	    if(c_ij && c_jk)
	      T12 = get_triplet(xx,j,i,k,&bond_hash,&T12work,&dvir_ij,&dvir_jk);
	    if(c_ki && c_jk)
	      T23 = get_triplet(xx,k,i,j,&bond_hash,&T23work,&dvir_ki,&dvir_jk);
	    if(c_ij && c_ki)
	      T31 = get_triplet(xx,i,j,k,&bond_hash,&T31work,&dvir_ij,&dvir_ki);

	    if(evflag) {
	      fsave[0][0] = fix; fsave[0][1] = fiy; fsave[0][2] = fiz;
	      fsave[1][0] = fjx; fsave[1][1] = fjy; fsave[1][2] = fjz;
	      fsave[2][0] = fkx; fsave[2][1] = fky; fsave[2][2] = fkz;
	      fix = fiy = fiz = 0.0;
	      fjx = fjy = fjz = 0.0;
	      fkx = fky = fkz = 0.0;
	    }

	    tr0 = tr1 = tr2 = tr3 = 0.0;
	    double xvir3t,xvir3v;
	    xvir3t = xvir3v = 0.0;

	    if(T12 && T23) {
	      bond_data *bki = bond_hash.Lookup(Doublet(k,i));

	      if(three_body_energies && evflag) {
		tr0 = transtrace(T12->H1H2,bki->H);
		double dvir = ((dvir_ij + dvir_jk + bki->fl_deriv_sum)*splinepot.vc +
				 splinepot.dvc)*tr0*w3/anorm3;
		vir3t = vir3t + dvir;
		xvir3t = xvir3t + dvir;
	      }
	      mcount2++;

	      {
		const double vc = splinepot.vc;
		tr_trace3(&(bki->H.m[1][0]),
			    &(T12->H1xH2.m[1][0]),&utr1x.d,
			    &(T12->H1yH2.m[1][0]),&utr1y.d,
			    &(T12->H1zH2.m[1][0]),&utr1z.d);

		tr_trace3(&(bki->H.m[1][0]),
			    &(T12->H1H2x.m[1][0]),&utr2x.d,
			    &(T12->H1H2y.m[1][0]),&utr2y.d,
			    &(T12->H1H2z.m[1][0]),&utr2z.d);

		tr_trace3(&(T12->H1H2.m[1][0]),
			  &(bki->Hx.m[1][0]),&utr3x.d,
			  &(bki->Hy.m[1][0]),&utr3y.d,
			  &(bki->Hz.m[1][0]),&utr3z.d);

		if(linalg.single) {
		  trd1x = utr1x.f; trd2x = utr2x.f; trd3x = utr3x.f;
		  trd1y = utr1y.f; trd2y = utr2y.f; trd3y = utr3y.f;
		  trd1z = utr1z.f; trd2z = utr2z.f; trd3z = utr3z.f;
		} else {
		  trd1x = utr1x.d; trd2x = utr2x.d; trd3x = utr3x.d;
		  trd1y = utr1y.d; trd2y = utr2y.d; trd3y = utr3y.d;
		  trd1z = utr1z.d; trd2z = utr2z.d; trd3z = utr3z.d;
		}

		dfix = ( (-sij)*trd1x + ( ski)*trd3x ) * (vc / anorm3);
		dfjx = ( ( sij)*trd1x + (-sjk)*trd2x ) * (vc / anorm3);
		dfkx = ( ( sjk)*trd2x + (-ski)*trd3x ) * (vc / anorm3);

		dfiy = ( (-sij)*trd1y + ( ski)*trd3y ) * (vc / anorm3);
		dfjy = ( ( sij)*trd1y + (-sjk)*trd2y ) * (vc / anorm3);
		dfky = ( ( sjk)*trd2y + (-ski)*trd3y ) * (vc / anorm3);

		dfiz = ( (-sij)*trd1z + ( ski)*trd3z ) * (vc / anorm3);
		dfjz = ( ( sij)*trd1z + (-sjk)*trd2z ) * (vc / anorm3);
		dfkz = ( ( sjk)*trd2z + (-ski)*trd3z ) * (vc / anorm3);
	      }

	      if(triplet_debug)
		force_debug_3t(xx,i,j,k, i,j,k,
			       dfix,dfiy,dfiz,
			       dfjx,dfjy,dfjz,
			       dfkx,dfky,dfkz);

	      if(three_body_forces)
		accumulate_forces_3(w3);
	    }

	    if(T12 != 0) {
	      //printf("T12 i,j,k = %d,%d,%d\n",i,j,k);
	      mcount++;
	      if(three_body_energies && evflag) {
		tr1 = transtrace(T12->H1H2,T12->H1H2);
		double dvir = (2.0*(dvir_ij + dvir_jk)*splinepot.vd +
				 splinepot.dvd)*tr1*w3/anorm4;
		vir3v = vir3v + dvir;
		xvir3v = xvir3v + dvir;
	      }

	      {
		const double vd = splinepot.vd;

		tr_trace3(&(T12->H1H2.m[1][0]),
			  &(T12->H1xH2.m[1][0]),&utr1x.d,
			  &(T12->H1yH2.m[1][0]),&utr1y.d,
			  &(T12->H1zH2.m[1][0]),&utr1z.d);
		tr_trace3(&(T12->H1H2.m[1][0]),
			  &(T12->H1H2x.m[1][0]),&utr2x.d,
			  &(T12->H1H2y.m[1][0]),&utr2y.d,
			  &(T12->H1H2z.m[1][0]),&utr2z.d);
		if(linalg.single) {
		  trd1x = utr1x.f; trd2x = utr2x.f;
		  trd1y = utr1y.f; trd2y = utr2y.f;
		  trd1z = utr1z.f; trd2z = utr2z.f;
		} else {
		  trd1x = utr1x.d; trd2x = utr2x.d;
		  trd1y = utr1y.d; trd2y = utr2y.d;
		  trd1z = utr1z.d; trd2z = utr2z.d;
		}

		dfix = 2.0*(-sij)*trd1x * (vd / anorm4);
		dfkx = 2.0*( sjk)*trd2x * (vd / anorm4);
		dfjx = -(dfix + dfkx);

		dfiy = 2.0*(-sij)*trd1y * (vd / anorm4);
		dfky = 2.0*( sjk)*trd2y * (vd / anorm4);
		dfjy = -(dfiy + dfky);

		dfiz = 2.0*(-sij)*trd1z * (vd / anorm4);
		dfkz = 2.0*( sjk)*trd2z * (vd / anorm4);
		dfjz = -(dfiz + dfkz);
	      }

	      if(triplet_debug) /* Compare forces to numerical derivatives */
		force_debug_3v(xx,i,j,k, j,i,k,
			       dfix,dfiy,dfiz,
			       dfjx,dfjy,dfjz,
			       dfkx,dfky,dfkz);


	      if(three_body_forces)
		accumulate_forces_3(w3);
	    }

	    if(T23 != 0) {
	      //printf("T23 i,j,k = %d,%d,%d\n",i,j,k);
	      mcount++;
	      if(three_body_energies && evflag) {
		tr2 = transtrace(T23->H1H2,T23->H1H2);
		double dvir = (2.0*(dvir_jk + dvir_ki)*splinepot.vd +
				 splinepot.dvd)*tr2*w3/anorm4;
		vir3v = vir3v + dvir;
		xvir3v = xvir3v + dvir;
	      }

	      {
		const double vd = splinepot.vd;

		tr_trace3(&(T23->H1H2.m[1][0]),
			  &(T23->H1xH2.m[1][0]),&utr1x.d,
			  &(T23->H1yH2.m[1][0]),&utr1y.d,
			  &(T23->H1zH2.m[1][0]),&utr1z.d);
		tr_trace3(&(T23->H1H2.m[1][0]),
			  &(T23->H1H2x.m[1][0]),&utr2x.d,
			  &(T23->H1H2y.m[1][0]),&utr2y.d,
			  &(T23->H1H2z.m[1][0]),&utr2z.d);
		if(linalg.single) {
		  trd1x = utr1x.f; trd2x = utr2x.f;
		  trd1y = utr1y.f; trd2y = utr2y.f;
		  trd1z = utr1z.f; trd2z = utr2z.f;
		} else {
		  trd1x = utr1x.d; trd2x = utr2x.d;
		  trd1y = utr1y.d; trd2y = utr2y.d;
		  trd1z = utr1z.d; trd2z = utr2z.d;
		}

		dfix = 2.0*( ski)*trd1x * (vd / anorm4);
		dfjx = 2.0*(-sjk)*trd2x * (vd / anorm4);
		dfkx = -(dfix + dfjx);

		dfiy = 2.0*( ski)*trd1y * (vd / anorm4);
		dfjy = 2.0*(-sjk)*trd2y * (vd / anorm4);
		dfky = -(dfiy + dfjy);

		dfiz = 2.0*( ski)*trd1z * (vd / anorm4);
		dfjz = 2.0*(-sjk)*trd2z * (vd / anorm4);
		dfkz = -(dfiz + dfjz);
	      }

	      if(triplet_debug) /* Compare forces to numerical derivatives */
		force_debug_3v(xx,i,j,k, k,i,j,
			       dfix,dfiy,dfiz,
			       dfjx,dfjy,dfjz,
			       dfkx,dfky,dfkz);

	      if(three_body_forces)
		accumulate_forces_3(w3);

	    }

	    if(T31 != 0) {
	      //printf("T31 i,j,k = %d,%d,%d\n",i,j,k);
	      mcount++;
	      if(three_body_energies && evflag) {
		tr3 = transtrace(T31->H1H2,T31->H1H2);
		double dvir = (2.0*(dvir_ki + dvir_ij)*splinepot.vd +
				 splinepot.dvd)*tr3*w3/anorm4;
		vir3v = vir3v + dvir;
		xvir3v = xvir3v + dvir;
	      }

	      {
		const double vd = splinepot.vd;

		tr_trace3(&(T31->H1H2.m[1][0]),
			  &(T31->H1xH2.m[1][0]),&utr1x.d,
			  &(T31->H1yH2.m[1][0]),&utr1y.d,
			  &(T31->H1zH2.m[1][0]),&utr1z.d);
		tr_trace3(&(T31->H1H2.m[1][0]),
			  &(T31->H1H2x.m[1][0]),&utr2x.d,
			  &(T31->H1H2y.m[1][0]),&utr2y.d,
			  &(T31->H1H2z.m[1][0]),&utr2z.d);
		if(linalg.single) {
		  trd1x = utr1x.f; trd2x = utr2x.f;
		  trd1y = utr1y.f; trd2y = utr2y.f;
		  trd1z = utr1z.f; trd2z = utr2z.f;
		} else {
		  trd1x = utr1x.d; trd2x = utr2x.d;
		  trd1y = utr1y.d; trd2y = utr2y.d;
		  trd1z = utr1z.d; trd2z = utr2z.d;
		}

		dfjx = 2.0*( sij)*trd1x * (vd / anorm4);
		dfkx = 2.0*(-ski)*trd2x * (vd / anorm4);
		dfix = -(dfjx + dfkx);

		dfjy = 2.0*( sij)*trd1y * (vd / anorm4);
		dfky = 2.0*(-ski)*trd2y * (vd / anorm4);
		dfiy = -(dfjy + dfky);

		dfjz = 2.0*( sij)*trd1z * (vd / anorm4);
		dfkz = 2.0*(-ski)*trd2z * (vd / anorm4);
		dfiz = -(dfjz + dfkz);

	      }

	      if(triplet_debug) /* Compare forces to numerical derivatives */
		force_debug_3v(xx,i,j,k, i,j,k,
			       dfix,dfiy,dfiz,
			       dfjx,dfjy,dfjz,
			       dfkx,dfky,dfkz);

	      if(three_body_forces)
		accumulate_forces_3(w3);
	    }

	    v33 = tr0 / anorm3;
	    v43 = (tr1 + tr2 + tr3) / anorm4;
	    double de_triplet = (splinepot.vc*v33 + splinepot.vd*v43) * e_scale * w3;
	    e_triplet = e_triplet + de_triplet;
	    e_triplet_c = e_triplet_c + splinepot.vc*v33 * e_scale * w3;
	    c_t++;

	    //printf("xxxx %6d %6d %6d :: %20.10e\n",1,2,3,de_triplet);

	    if(evflag) {
	      double drji[3],drki[3];
	      double fj[3] = {fjx,fjy,fjz},fk[3] = {fkx,fky,fkz};
	      for(int p = 0; p<3; p++) {
		drji[p] = xx[j][p] - xx[i][p];
		drki[p] = xx[k][p] - xx[i][p];
		/* To fix stress-per-atom scaling. */
		fj[p] *= e_scale;
		fk[p] *= e_scale;
	      }

	      ev_tally3(i,j,k,de_triplet,0.0,fj,fk,drji,drki);

	      if(volpres_flag && vflag_atom) {
		//virial[i] = virial[i] - (vir3v + vir3t) * rhoinv*e_scale;
		double dvir = -(xvir3v + xvir3t) * rhoinv*e_scale * (1.0/3.0);
		for(int pp = 0; pp<3; pp++) {
		  vatom[i][pp] += dvir;
		  vatom[j][pp] += dvir;
		  vatom[k][pp] += dvir;
		}
	      }

	      fix = fix+fsave[0][0]; fiy = fiy+fsave[0][1]; fiz = fiz+fsave[0][2];
	      fjx = fjx+fsave[1][0]; fjy = fjy+fsave[1][1]; fjz = fjz+fsave[1][2];
	      fkx = fkx+fsave[2][0]; fky = fky+fsave[2][1]; fkz = fkz+fsave[2][2];
	    }

	    tx1 = gettime();
	    ttriplet += tx1 - tx0;
	    nttriplet++;
	  } else {
	    triplet_defer = 1;
	  }

	  if(four_body_energies || four_body_forces)
	    if(j < i) { /* Search for quadruplet */
	      tx0 = gettime();

	      mj = first[j];
	      mk = first[k];
	      /*
		i is in both the j-list and the k-list, and i > k,
		and lists are sorted, so the loop terminates.
	      */
	      while(nlist_short[mj] < i && nlist_short[mk] < i) {

		if(mj >= first[j+1] || mk >= first[k+1]) {
		  printf("Illegal quad...\n"
			 "  j=%d  first[j]=%d  first[j+1]=%d  mj=%d\n"
			 "  k=%d  first[k]=%d  first[k+1]=%d  mk=%d\n",
			 j,first[j],first[j+1],mj,
			 k,first[k],first[k+1],mk);
		  error->one(__FILE__,__LINE__,"Shit, brkoen quad loop");
		}

		if(nlist_short[mj] == nlist_short[mk]) {
		  /* Closed quadruplet */
		  m = nlist_short[mj];
		  c_jm = c_km = 1;

		  const int sim = (i < m) ? 1 : -1;
		  const int sjm = (j < m) ? 1 : -1;
		  const int skm = (k < m) ? 1 : -1;

		  w4 = get_weight(triclinic,ss[i],ss[j],ss[k],ss[m]);

		  if(w4 > 0.0) {

		    /* Alrady know ij,jk,ki,jm,km bonds. Look for im bond. */
		    mi = first[i];
		    while(mi < first[i+1] && nlist_short[mi] < m) mi = mi + 1;
		    if(mi < first[i+1] && nlist_short[mi] == m)
		      c_im = 1;
		    else
		      c_im = 0;

		    if(c_im == 0 || c_jk == 0 || (c_jk && c_im && m < k)) {

		      if(triplet_defer) {
			dvir_ij = dvir_jk = dvir_ki = 0.0;
			if(c_ij && c_jk)
			  T12 = get_triplet(xx,j,i,k,&bond_hash,&T12work,&dvir_ij,&dvir_jk);
			if(c_ki && c_jk)
			  T23 = get_triplet(xx,k,i,j,&bond_hash,&T23work,&dvir_ki,&dvir_jk);
			if(c_ij && c_ki)
			  T31 = get_triplet(xx,i,j,k,&bond_hash,&T31work,&dvir_ij,&dvir_ki);
			triplet_defer = 0;
		      }


		      fmx = fmy = fmz = 0.0;
		      double xvir4 = 0.0;

		      if(evflag) {
			fsave[0][0] = fix; fsave[0][1] = fiy; fsave[0][2] = fiz;
			fsave[1][0] = fjx; fsave[1][1] = fjy; fsave[1][2] = fjz;
			fsave[2][0] = fkx; fsave[2][1] = fky; fsave[2][2] = fkz;
			fsave[3][0] = fmx; fsave[3][1] = fmy; fsave[3][2] = fmz;
			fix = fiy = fiz = 0.0;
			fjx = fjy = fjz = 0.0;
			fkx = fky = fkz = 0.0;
			fmx = fmy = fmz = 0.0;
		      }

		      tr1 = tr2 = tr3 = 0.0;

		      dvir_im = dvir_jm = dvir_km = 0.0;
		      T45 = T56 = T64 = 0;
		      if(T12 != 0 && c_km && c_im)
			T45 = get_triplet(xx,m,i,k,&bond_hash,&T45work,&dvir_im,&dvir_km);
		      if(T23 != 0 && c_im && c_jm)
			T56 = get_triplet(xx,m,i,j,&bond_hash,&T56work,&dvir_im,&dvir_jm);
		      if(T31 != 0 && c_jm && c_km)
			T64 = get_triplet(xx,m,j,k,&bond_hash,&T64work,&dvir_jm,&dvir_km);

		      if(T12 != 0 && T45 != 0) {
			if(four_body_energies && evflag) {
			  tr1 = transtrace(T12->H1H2,T45->H1H2);
			  double dvir = ( (dvir_ij + dvir_jk + dvir_im + dvir_km)*splinepot.ve +
					  splinepot.dve )*tr1*w4/anorm4;
			  vir4 = vir4 + dvir;
			  xvir4 = xvir4 + dvir;
			}
			qcount++;

			{
			  const double ve = splinepot.ve;

			  trd_update_4(T12,T45);

			  dfix_update_4a(x);
			  dfix_update_4a(y);
			  dfix_update_4a(z);
			}

			if(quad_debug) /* Compare forces to numerical derivatives */
			  force_debug_4(xx,i,j,k,m, i,j,k,m,
					dfix,dfiy,dfiz , dfjx,dfjy,dfjz,
					dfkx,dfky,dfkz , dfmx,dfmy,dfmz);

			if(four_body_forces)
			  accumulate_forces_4(w4);
		      }

		      if(T23 != 0 && T56 != 0) {
			if(four_body_energies && evflag) {
			  tr2 = transtrace(T23->H1H2,T56->H1H2);
			  double dvir = ( (dvir_ki + dvir_jk + dvir_im + dvir_jm)*splinepot.ve +
					  splinepot.dve )*tr2*w4/anorm4;
			  vir4 = vir4 + dvir;
			  xvir4 = xvir4 + dvir;
			}
			qcount++;

			{
			  const double ve = splinepot.ve;

			  trd_update_4(T23,T56);

			  dfix_update_4b(x);
			  dfix_update_4b(y);
			  dfix_update_4b(z);
			}

			if(quad_debug) /* Compare forces to numerical derivatives */
			  force_debug_4(xx,i,j,k,m, i,m,j,k,
					dfix,dfiy,dfiz , dfjx,dfjy,dfjz,
					dfkx,dfky,dfkz , dfmx,dfmy,dfmz);

			if(four_body_forces)
			  accumulate_forces_4(w4);

		      }

		      if(T31 != 0 && T64 != 0) {
			if(four_body_energies && evflag) {
			  tr3 = transtrace(T31->H1H2,T64->H1H2);
			  double dvir = ( (dvir_ki + dvir_ij + dvir_jm + dvir_km)*splinepot.ve +
					  splinepot.dve )*tr3*w4/anorm4;
			  vir4 = vir4 + dvir;
			  xvir4 = xvir4 + dvir;
			}
			qcount++;

			{
			  const double ve = splinepot.ve;

			  /* X */
			  trd_update_4(T31,T64);

			  dfix_update_4c(x);
			  dfix_update_4c(y);
			  dfix_update_4c(z);
			}

			if(quad_debug) /* Compare forces to numerical derivatives */
			  force_debug_4(xx,i,j,k,m, i,j,m,k,
					dfix,dfiy,dfiz , dfjx,dfjy,dfjz,
					dfkx,dfky,dfkz , dfmx,dfmy,dfmz);

			if(four_body_forces)
			  accumulate_forces_4(w4);
		      }

		      double de_quad = splinepot.ve*(tr1 + tr2 + tr3)/anorm4 * e_scale * w4;
		      e_quad = e_quad + de_quad;
		      if((T12 && T45) ||
			 (T23 && T56) ||
			 (T31 && T64)) {
			c_q++;
		      }

		      if(evflag) {
			double drim[3],drjm[3],drkm[3];
			double fi[3] = {fix,fiy,fiz};
			double fj[3] = {fjx,fjy,fjz};
			double fk[3] = {fkx,fky,fkz};
			for(int p = 0; p<3; p++) {
			  drim[p] = xx[i][p] - xx[m][p];
			  drjm[p] = xx[j][p] - xx[m][p];
			  drkm[p] = xx[k][p] - xx[m][p];
			  fi[p] *= e_scale;
			  fj[p] *= e_scale;
			  fk[p] *= e_scale;
			}

			ev_tally4(i,j,k,m,de_quad,fi,fj,fk,drim,drjm,drkm);

			if(volpres_flag && vflag_atom) {
			  //virial[i] = virial[i] - vir4 * rhoinv*e_scale;
			  double dvir = -xvir4 * rhoinv*e_scale * (1.0/4.0);
			  for(int pp = 0; pp<3; pp++) {
			    vatom[i][pp] += dvir;
			    vatom[j][pp] += dvir;
			    vatom[k][pp] += dvir;
			    vatom[m][pp] += dvir;
			  }
			}

			fix = fix+fsave[0][0]; fiy = fiy+fsave[0][1]; fiz = fiz+fsave[0][2];
			fjx = fjx+fsave[1][0]; fjy = fjy+fsave[1][1]; fjz = fjz+fsave[1][2];
			fkx = fkx+fsave[2][0]; fky = fky+fsave[2][1]; fkz = fkz+fsave[2][2];
			fmx = fmx+fsave[3][0]; fmy = fmy+fsave[3][1]; fmz = fmz+fsave[3][2];
		      }

		      ff[m][0] += fmx * e_scale;
		      ff[m][1] += fmy * e_scale;
		      ff[m][2] += fmz * e_scale;

		    }
		  }
		  mj = mj + 1;
		  mk = mk + 1;
		} else if(nlist_short[mj] < nlist_short[mk]) {
		  mj = mj + 1;
		} else {
		  mk = mk + 1;
		}

	      }
	      tx1 = gettime();
	      tquad += tx1 - tx0;
	      ntquad++;
	      ntquaditer++;
	    }


	  ff[k][0] += fkx * e_scale;
	  ff[k][1] += fky * e_scale;
	  ff[k][2] += fkz * e_scale;

	}
#undef transtrace

	ff[j][0] += fjx * e_scale;
	ff[j][1] += fjy * e_scale;
	ff[j][2] += fjz * e_scale;

      }

    ff[i][0] += fix * e_scale;
    ff[i][1] += fiy * e_scale;
    ff[i][2] += fiz * e_scale;

    if(single_energies == 1 && i < nloc) {
      const double evol0 = splinepot.evol0;
      if(eflag_global) {
	e_single = e_single + evol0 * e_scale;
	eng_vdwl = eng_vdwl + evol0 * e_scale;
      }
      if(eflag_atom) eatom[i] = eatom[i] + evol0 * e_scale;
      if(volpres_flag && vflag_atom) {
	for(int pp = 0; pp<3; pp++)
	  vatom[i][pp] = vatom[i][pp] - rhoinv*splinepot.devol0*e_scale;
      }

    }

  }

  tx0 = gettime();
  for(i = 0; i<ntot; i++)
    for(p = 0; p<3; p++)
      atom->f[i][p] = atom->f[i][p] + ff[i][p];

  memory->sfree(nlist_short);
  memory->sfree(first);
  if(ss != xx) memory->sfree(ss);
  memory->sfree(ff);
  memory->sfree(xx);
  tx1 = gettime();
  tmem += tx1-tx0;
  ntmem++;

  t1 = gettime(1);

  //printf("compute_x: c_p = %d    c_t = %d    c_q = %d\n",c_p,c_t,c_q);


#ifdef TIMING_ON
  if(comm->me == 0) {
    double tsum = (tmem+tsort+tpair+tlookup+ttriplet+tquad);
    double nsum = (ntmem+ntsort+ntpair+ntlookup+nttriplet+ntquad);
    //double adj = ((t1-t0)-tsum)/nsum;
    /* Use adj = 6ns for RDTSC, and 58ns for gettimeofday,
       on monkfish.llnl.gov, 2.4GHz Intel

       Use adj = 35.945ns for RDTSC on uBGL (assumed rate set to 700MHz)
    */
    double adj = 35.945e-9;

    double
      memadj     = tmem     - adj*ntmem    ,
      sortadj    = tsort    - adj*ntsort   ,
      pairadj    = tpair    - adj*ntpair   ,
      lookupadj  = tlookup  - adj*ntlookup ,
      tripletadj = ttriplet - adj*nttriplet,
      quadadj    = tquad    - adj*ntquad   ,

      make_b_adj = t_make_b - adj*n_make,
      make_t_adj = t_make_t - adj*n_make,
      make_b2_adj = t_make_b2 - adj*n_make_b2,
      trace_adj = t_trace - adj*n_trace;

    printf("mgpt engy = %10.3fms\n",(t1-t0)*1e3);
    printf("       mem = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   tmem*1e3,ntmem,memadj*1e3,memadj/ntmem*1e9);
    printf("      sort = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   tsort*1e3,ntsort,sortadj*1e3,sortadj/ntsort*1e9);
    printf("      pair = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   tpair*1e3,ntpair,pairadj*1e3,pairadj/ntpair*1e9);
    printf("    lookup = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   tlookup*1e3,ntlookup,lookupadj*1e3,lookupadj/ntlookup*1e9);
    printf("   triplet = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   ttriplet*1e3,nttriplet,tripletadj*1e3,tripletadj/nttriplet*1e9);
    printf("      quad = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   tquad*1e3,ntquaditer,quadadj*1e3,quadadj/ntquaditer*1e9);
    printf("       sum = %10.3fms                adj = %10.3fms\n",
	   tsum*1e3,(tsum - adj*nsum)*1e3);
    printf("\n    make_b = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   t_make_b*1e3,n_make,make_b_adj*1e3,make_b_adj/n_make*1e9);
    printf("   make_b2 = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n",
	   t_make_b2*1e3,n_make_b2,make_b2_adj*1e3,make_b2_adj/n_make_b2*1e9);
    printf("    make_t = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n\n",
	   t_make_t*1e3,n_make,make_t_adj*1e3,make_t_adj/n_make*1e9);
    printf("     trace = %10.3fms  n = %8.0f  adj = %10.3fms  one = %10.3fns\n\n",
	   t_trace*1e3,n_trace,trace_adj*1e3,trace_adj/n_trace*1e9);

    printf("mcount (transpose + trace for triplet) = %.0f , %.0f  qcount = %.0f  lmax = %d\n",
	   mcount,mcount2,qcount,lmax);

    printf("nbc=%.0f  tbl=%.3fms  tbm=%.3fms  one tbl=%.3fns  one tbm=%.3fns\n",
	   nbc,(tbl-adj*nbc)*1e3,(tbm-adj*nbc)*1e3,(tbl/nbc-adj)*1e9,
	   (tbm/nbc-adj)*1e9);
    printf("\n\nForces:\n");
    printf("fix = %.3f  fiy=%.3f  fiz=%.3f\n",fix,fiy,fiz);
    printf("fjx = %.3f  fjy=%.3f  fjz=%.3f\n",fjx,fjy,fjz);
    printf("fkx = %.3f  fky=%.3f  fkz=%.3f\n",fkx,fky,fkz);
    printf("\n");

    printf("Bonds   : nsearch=%d maxlen=%d avg.len=%.3f\n",
	   bond_hash.NSearch(),bond_hash.MaxLength(),
	   bond_hash.NStep()/(double) bond_hash.NSearch());

    printf("compute_x: c_p = %d    c_t = %d    c_q = %d\n",c_p,c_t,c_q);

    printf("@@ Total number of trace3 calls is %d, total number of make_triplet is %.1f\n",
	   ntr_calls,n_make);

    {
      Hash<bond_data,Doublet>::Iterator iter = bond_hash.begin();
      int nitem = 0,nhit = 0;
      while(iter != bond_hash.end()) {
	nitem++;
	nhit += iter.link()->hits;
	iter.next();
      }
      printf("bond_hash    hits: nitems=%d   nhits=%d  hits/item = %.3f\n",
	     nitem,nhit,nhit/(double) nitem);
    }
  }
#endif

  if(volpres_flag) {
    /*
      Include contributions to the pressure due to derivatines
      of the energy with respect to the potential input volume.
    */

    /* The following lines have moved to beginning of functions,
       since they are used in calculating per-atom virial contributions */
    /*
      double vtot = 1.0;
      double ntot = atom->natoms;
      for(i = 0; i<3; i++)
      vtot = vtot * (domain->boxhi[i] - domain->boxlo[i]);
      double rhoinv = vtot / ntot;
    */

    if(single_energies) // Virial correction for self energy
      for(i = 0; i<3; i++) {
	//virial[i] = virial[i] + nloc*pot_input_vol*pvol0*e_scale;
	virial[i] = virial[i] - nloc*rhoinv*splinepot.devol0*e_scale;
      }

    if(pair_energies) // Virial correction for pair energy
      for(i = 0; i<3; i++)
	virial[i] = virial[i] + rhoinv*e_scale*volvir2;

    if(three_body_energies) // Virial correction for three body enegries
      for(i = 0; i<3; i++) {
	//virial[i] = virial[i] - pot_input_vol*(e_triplet_c*pc + (e_triplet-e_triplet_c)*pd);
	virial[i] = virial[i] - (vir3v + vir3t) * rhoinv*e_scale;
      }

    if(four_body_energies) // Virial correction for four body enegries
      for(i = 0; i<3; i++) {
	//virial[i] = virial[i] - pot_input_vol*e_quad*pe;
	virial[i] = virial[i] - vir4 * rhoinv*e_scale;
      }

  }

  *e_s = e_single;
  *e_p = e_pair;
  *e_t = e_triplet;
  *e_q = e_quad;
}

void PairMGPT::compute(int eflag, int vflag)
{
  ev_init(eflag, vflag);

  int newton_pair = force->newton_pair;
  double e_s,e_p,e_t,e_q;

  //printf("newton_pair = %d, newton = %d, tag_enable = %d\n",force->newton_pair,force->newton,atom->tag_enable);

  if(newton_pair == 0) {
    printf("This is a problem. MGPT requires newton_pair flag to be on. Exiting...\n");
    exit(1);
  }

  if(atom->tag_enable == 0) {
    printf("This is a problem. MGPT requires tag_enable flag to be on. Exiting...\n");
    exit(1);
  }

  compute_x(listfull->numneigh,listfull->firstneigh,&e_s,&e_p,&e_t,&e_q,evflag,newton_pair);

  if(0) { // Stupid force calculation / verification
    int ii,nmax=-1;
    for(ii = 0; ii<listfull->inum + listfull->gnum; ii++) {
      int i = listfull->ilist[ii];
      if(i > nmax) nmax = i;
    }
    nmax++;
    double *ffwork = new double[3*nmax];
    double *ffloc = new double[3*listfull->inum];
    double *ffloc2 = new double[3*listfull->inum];
    double **ffptr = new double *[nmax];
    for(ii = 0; ii<listfull->inum + listfull->gnum; ii++)
      ffptr[ii] = &ffwork[3*ii];

    printf("Computing boundary forces\n");
    for(ii = 0; ii<listfull->inum; ii++) {
      ffloc2[3*ii] = 0.0;
      ffloc2[3*ii+1] = 0.0;
      ffloc2[3*ii+2] = 0.0;
      int i = listfull->ilist[ii];
      for(int jj = 0; jj<listfull->inum+listfull->gnum; jj++) {
	int j = listfull->ilist[jj];
	if(atom->tag[i] == atom->tag[j])
	  for(int p = 0; p<3; p++)
	    ffloc2[3*ii+p] += atom->f[j][p];
      }
    }

    printf("Starting main displacement force calculation\n");
    for(ii = 0; ii<listfull->inum; ii++) {
      int i = listfull->ilist[ii];

      double **atom_f_save = atom->f;
      atom->f = ffptr;

      for(int p = 0; p<3; p++) {
	double xsave = atom->x[i][p];
	const double delta = 1e-3;

	atom->x[i][p] = xsave + delta;
	for(int jj = 0; jj<3*nmax; jj++) ffwork[jj] = 0.0;
	compute_x(listfull->numneigh,
		  listfull->firstneigh,
		  &e_s,&e_p,&e_t,&e_q,evflag,newton_pair);
	double e1 = e_s + e_p + e_t + e_q;

	atom->x[i][p] = xsave - delta;
	for(int jj = 0; jj<3*nmax; jj++) ffwork[jj] = 0.0;
	compute_x(listfull->numneigh,
		  listfull->firstneigh,
		  &e_s,&e_p,&e_t,&e_q,evflag,newton_pair);
	double e2 = e_s + e_p + e_t + e_q;

	ffloc[3*ii+p] = -(e1-e2)/(2*delta);

	atom->x[i][p] = xsave;
      }

      atom->f = atom_f_save;
      printf("Force on i=%4d:\n",i);
      printf("  Position   %20.10e  %20.10e  %20.10e\n",
	     atom->x[i][0],atom->x[i][1],atom->x[i][2]);
      printf("  Exact      %20.10e  %20.10e  %20.10e\n",
	     atom->f[i][0],atom->f[i][1],atom->f[i][2]);
      printf("  Numerical  %20.10e  %20.10e  %20.10e\n",
	     ffloc[3*ii+0],ffloc[3*ii+1],ffloc[3*ii+2]);
      printf("  Boundary   %20.10e  %20.10e  %20.10e\n",
	     ffloc2[3*ii+0],ffloc2[3*ii+1],ffloc2[3*ii+2]);
    }


    delete[] ffloc2;
    delete[] ffloc;
    delete[] ffptr;
    delete[] ffwork;
  }


  if(0) {
    printf("\nForces MGPT:\n");
    const int iimax = (listfull->inum < 10) ? listfull->inum : 10;
    for(int ii = 0; ii<iimax; ii++) {
      int i = listfull->ilist[ii];
      printf("%4d  =  %20.10e  %20.10e  %20.10e\n",
             i,atom->f[i][0],atom->f[i][1],atom->f[i][2]);
    }
    printf("\n\n");
  }

  if(vflag_fdotr) {
    //printf("##### Using virial_compute!!!\n");
    virial_fdotr_compute();
  }
}

void PairMGPT::allocate()
{
	allocated = 1;
	int n = atom->ntypes;

	memory->create(setflag,n+1,n+1,"pair:setflag");
	for(int i = 0; i <= n; i++)
		for(int j = 0; j <= n; j++)
			setflag[i][j] = 0;

	memory->create(cutsq,n+1,n+1,"pair:cutsq");
	memory->create(cutghost,n+1,n+1,"pair:cutsq");
}

/* ----------------------------------------------------------------------
   global settings
------------------------------------------------------------------------- */
void PairMGPT::settings(int narg, char **/*arg*/)
{
  if(narg != 0) error->all(__FILE__,__LINE__,"Illegal pair_style command");
}

/* ----------------------------------------------------------------------
   set coeffs for one or more type pairs
------------------------------------------------------------------------- */

void PairMGPT::coeff(int narg, char **arg)
{
  int single_precision = 0;

  if(narg < 5)
    error->all(__FILE__,__LINE__,
	       "Not enough arguments for mgpt (MGPT) pair coefficients.");

  if(!allocated) allocate();

  // Make sure I,J args are * *
  if(strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
    error->all(__FILE__,__LINE__,"Incorrect args for pair coefficients");

  double vol;
  if(sscanf(arg[4], "%lg", &vol) != 1 || vol <= 0.0)
    error->all(__FILE__,__LINE__,"Invalid volume in mgpt (MGPT) pair coefficients.");

  volpres_flag = 1;
  single_precision = 0;

  /* Parse arguments */ {
    int volpres_tag = 0,precision_tag = 0,nbody_tag = 0;

    int iarg = 5;
    while (iarg < narg) {
      if(strcmp(arg[iarg],"volpress") == 0) { /* Volumetric pressure flag */
        if (iarg+2 > narg)
          error->all(FLERR,"Incorrect args for pair coefficients");
	if(strcmp(arg[iarg+1],"yes") == 0) volpres_flag = 1;
	else if(strcmp(arg[iarg+1],"no") == 0) volpres_flag = 0;
	else {
	  char line[1024];
	  sprintf(line,"(In %s:%d) Invalid value for volumetric pressure argument.\n"
		  "It should be \"volpress yes\" or \"volpress no\".\n"
		  "The value is \"%s\".\n",__FILE__,__LINE__,arg[iarg+1]);
	  error->all(__FILE__,__LINE__,line);
	}
	volpres_tag = 1;
        iarg += 2;
	if(comm->me == 0) printf("* volpress: volpres_flag = %d [%s %s]\n",volpres_flag,arg[iarg-2],arg[iarg-1]);
      } else if(strcmp(arg[iarg],"nbody") == 0) {
        if (iarg+2 > narg)
          error->all(FLERR,"Incorrect args for pair coefficients");
	if(strspn(arg[iarg+1],"1234") == strlen(arg[iarg+1])) {
	  nbody_flag = 0;
	  for(int i = 0; i<4; i++)
	    if(strchr(arg[iarg+1],'1'+i) != NULL) {
	      nbody_flag = nbody_flag + (1<<i);
	      if(comm->me == 0) printf("Explicitly adding %d-tuple forces.\n",i+1);
	    }
	} else {
	  char line[1024];
	  sprintf(line,"(In %s:%d) Invalid value for nbody flag.\n"
		  "It should be e.g. \"nbody=1234\" (for single, pair, triple, and quad forces/energiers)\n"
		  "For e.g. only pair and triple forces/energies, use \"nbody=23\".\n"
		  "The default is \"nbody=1234\".\n"
		  "The current value is \"%s\".\n",__FILE__,__LINE__,arg[iarg+1]);
	  error->all(__FILE__,__LINE__,line);
	}
	nbody_tag = 1;
        iarg += 2;
      } else if(strcmp(arg[iarg],"precision") == 0) {
        if (iarg+2 > narg)
          error->all(FLERR,"Incorrect args for pair coefficients");
	if(strcmp(arg[iarg+1],"single") == 0) single_precision = 1;
	else if(strcmp(arg[iarg+1],"double") == 0) single_precision = 0;
	else {
	  char line[1024];
	  sprintf(line,"(In %s:%d) Invalid value for precision argument.\n"
		  "It should be \"precision single\" or \"precision double\".\n"
		  "The value is \"%s\".\n",__FILE__,__LINE__,arg[iarg+1]);
	  error->all(__FILE__,__LINE__,line);
	}
	precision_tag = 1;
        iarg += 2;
	if(comm->me == 0) printf("* precision: single_flag = %d [%s %s]\n",single_precision,arg[iarg-2],arg[iarg-1]);
      } else {
	char line[1024];
	sprintf(line,"(In %s:%d) Invalid argument. Allowed arguments are:\n"
		"    volpress {yes|no} , default = yes\n"
		"    precision {single|double} , default = double\n"
		"    nbody {[1234,]*} , default = whichever terms potential require\n"
		"The invalid argument is \"%s\".\n",__FILE__,__LINE__,arg[iarg]);
	error->all(__FILE__,__LINE__,line);
      }
    }

    if(comm->me == 0)
      printf("Volumetric pressure is %s.\n",volpres_flag ? "on" : "off");

    if(comm->me == 0) {
      FILE *parmin_fp = force->open_potential(arg[2]);
      FILE *potin_fp = force->open_potential(arg[3]);
      if (parmin_fp == NULL || potin_fp == NULL) {
        char str[128];
        sprintf(str,"Cannot open MGPT potential files %s %s",arg[2],arg[3]);
        error->one(FLERR,str);
      }
      fclose(parmin_fp);
      fclose(potin_fp);

      splinepot.readpot(arg[2],arg[3],vol);
      printf("evol0 = %.10e\n",splinepot.evol0);

      /* Set up default and requested nbody forces to include */ {
	int nbody_default = (1<<0) + (1<<1) + (1<<2) + (1<<3);

	if(splinepot.vd == 0.0 && splinepot.dvd == 0.0)
	  nbody_default -= (1<<2); // No 3-body contributions
	if(splinepot.ve == 0.0 && splinepot.dve == 0.0)
	  nbody_default -= (1<<3); // No 4-body contributions

	if(nbody_tag == 0) nbody_flag = nbody_default;

	if(nbody_flag != nbody_default) {
	  printf("Warning: nbody=%d (suggested=%d) set to disregard multibody-forces in potential.\n",
		 nbody_flag,nbody_default);
	}
      }
    }
  }

  MPI_Bcast(&nbody_flag,sizeof(nbody_flag),MPI_BYTE,0,world);

  /*
    Broadcast structure to all processes. In receiving
    processes, pointes will be screwed up. We allocate
    memory, and then broadcast contents of arrays.
  */
  MPI_Bcast(&splinepot,sizeof(splinepot),MPI_BYTE,0,world);
  if(comm->me != 0) {
    splinepot.vpair_spline = new double[splinepot.nr-1][4];
    splinepot.dvpair_spline = new double[splinepot.nr-1][4];
  }
  MPI_Bcast(splinepot.vpair_spline,4*(splinepot.nr-1),MPI_DOUBLE,0,world);
  MPI_Bcast(splinepot.dvpair_spline,4*(splinepot.nr-1),MPI_DOUBLE,0,world);
  anorm3 = splinepot.anorm3;
  anorm4 = splinepot.anorm4;
  lmax = splinepot.lmax;
  lang = splinepot.lang;
  //ipot = splinepot.ipot;
  for(int i = 0; i<(int) (sizeof(ddl)/sizeof(double)); i++)
    ddl[i] = splinepot.ddl[i];
  for(int i = 0; i<lmax; i++) {
    for(int j = 0; j<lmax; j++)
      del0.m[i+1][j+1] = 0.0;
    del0.m[i+1][i+1] = 1.0;
  }

  /* Set matrix param, cutoff, LAMMPS param */
  Matrix::sz = lmax;

  rcrit = splinepot.rcrit;
  rmax = splinepot.rmax;
  cutoff = rmax;
  if(rcrit > rmax) cutoff = rcrit;

  // Set LAMMPS pair interaction flags.
  for(int i = 1; i <= atom->ntypes; i++) {
    for(int j = 1; j <= atom->ntypes; j++) {
      setflag[i][j] = 1;
      cutsq[i][j] = cutoff;
      cutghost[i][j] = cutoff;
    }
  }

  // Set atomic mass.
  for(int i = 1; i <= atom->ntypes; i++)
    atom->set_mass(FLERR,i, splinepot.mass);

  // Initialize linear algebra routines.
  linalg = mgpt_linalg(lmax,single_precision);
  if(comm->me == 0)
    printf("%s",linalg.msg);
}


/* ----------------------------------------------------------------------
   init specific to this pair style
------------------------------------------------------------------------- */
void PairMGPT::init_style()
{
	if(force->newton_pair == 0)
	  error->all(__FILE__,__LINE__,"Pair style mgpt requires newton pair on.");

	// Need full neighbor list.
	int irequest_full = neighbor->request(this);
	neighbor->requests[irequest_full]->id = 1;
	neighbor->requests[irequest_full]->half = 0;
	neighbor->requests[irequest_full]->full = 1;
	neighbor->requests[irequest_full]->ghost = 1;

	// Also need half neighbor list.
	int irequest_half = neighbor->request(this);
	neighbor->requests[irequest_half]->id = 2;
}

/* ----------------------------------------------------------------------
   neighbor callback to inform pair style of neighbor list to use
   half or full
------------------------------------------------------------------------- */
void PairMGPT::init_list(int id, NeighList *ptr)
{
	if(id == 1) listfull = ptr;
	else if(id == 2) listhalf = ptr;
}

/* ----------------------------------------------------------------------
   init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairMGPT::init_one(int /*i*/, int /*j*/)
{
	return cutoff;
}

/************************************************************************
 ****  REIMPLEMENTATION OF FL AND HAMLTN WITH ANALYTICAL DERIVATIVES ****
 ************************************************************************/
/*
  Reimplementation of bond length potential, including
  derivatives with respect to x,y, and z.
*/
void PairMGPT::fl_deriv_new(double r,double ri,double xhat,double yhat,double zhat,
			    double &fl_0,double &fl_x,double &fl_y,double &fl_z,
			    double &fl_rp,double &fl_p1,double &fl_r0,double &fl_al) {
  const double rp = splinepot.rp,p1 = splinepot.p1,r0 = splinepot.r00,al = splinepot.al;
  const int mode = splinepot.mode;
  const double pn = splinepot.pn;

  double t,tx,ty,tz,t_rp_ti,t_p1_ti;
  double s;

  /*
  // Original code
    double term;
    double pn=1.0;
    if (mode <= 4)
      term = pow(rp/r, p1);
    else
      term = exp(-p1*(pow(r/rp, pn) - 1.0)/pn);
  */

  double rpi = 1.0/rp;
  if(mode <= 4) {
    t = pow(rp*ri,p1);
    s = -p1 * t * ri;
    t_rp_ti = p1*rpi;
    t_p1_ti = log(rp*ri);
  } else {
    if(pn == 1.0) {
      double p1_rpi = -p1*rpi;
      t = exp(p1 + r*p1_rpi);
      s = p1_rpi * t;
      t_rp_ti = -r*p1_rpi*rpi;
      t_p1_ti = 1.0 - r*rpi;
    } else {
      double pni = 1.0/pn;
      double rprpn = pow(r*rpi,pn);
      t = exp(-p1*pni*(rprpn - 1.0));
      s = -p1*rprpn*ri * t;
      t_rp_ti = p1*rprpn*rpi;
      t_p1_ti = pni - pni*rprpn;// -pni*(rprpn - 1.0);
    }
  }
  tx = s * xhat;
  ty = s * yhat;
  tz = s * zhat;

  fl_rp = t_rp_ti;
  fl_p1 = t_p1_ti;

  if (r <= r0) {
    fl_0 = t;
    fl_x = tx;
    fl_y = ty;
    fl_z = tz;

    fl_r0 = 0.0;
    fl_al = 0.0;
  } else {
    double q,qx,qy,qz,exp_q,q_r0,q_al;
    double r0i,u;

    r0i = 1.0/r0;
    u = r*r0i - 1.0;
    q = al*u*u;
    s = 2*al*u*r0i;
    qx = s * xhat;
    qy = s * yhat;
    qz = s * zhat;
    q_r0 = -2.0*al*u*r*r0i*r0i;
    q_al = u*u;

    exp_q = exp(-q);
    if(mode <= 2) {
      fl_0 = exp_q * t;
      fl_x = exp_q*(tx - t*qx);
      fl_y = exp_q*(ty - t*qy);
      fl_z = exp_q*(tz - t*qz);

      fl_r0 = -q_r0;
      fl_al = -q_al;
    } else {
      fl_0 = exp_q * (1.0 + q) * t;
      fl_x = exp_q * (tx + q*(tx - t*qx));
      fl_y = exp_q * (ty + q*(ty - t*qy));
      fl_z = exp_q * (tz + q*(tz - t*qz));

      fl_r0 = -q_r0 * q/(1.0 + q);
      fl_al = -q_al * q/(1.0 + q);
    }
  }
}


/*
  Macros to build elements of the bond matrix, and also
  its derivatives with repsect to x,y, and z.
*/
#define MAKE_ELEMENT_5(i,j)                                          \
  do {                                                               \
    const double dl0  = del0.m[i][j];                                \
    const double dl4  = gsl_##i    * gsl_##j;                        \
    const double dl4x = gsl_##i##x * gsl_##j + gsl_##i * gsl_##j##x; \
    const double dl4y = gsl_##i##y * gsl_##j + gsl_##i * gsl_##j##y; \
    const double dl4z = gsl_##i##z * gsl_##j + gsl_##i * gsl_##j##z; \
                                                                     \
    const double tmp  = w4*dl4  + w2*dl2 + w0*dl0;                   \
    const double tmpx = w4*dl4x + w2*dl2x;                           \
    const double tmpy = w4*dl4y + w2*dl2y;                           \
    const double tmpz = w4*dl4z + w2*dl2z;                           \
    const double tmpsum = tmpx*x + tmpy*y + tmpz*z;                  \
    M [j][i] = M[i][j] =  fl  *tmp;                                  \
    Mx[j][i] = Mx[i][j] = fl_x*tmp + fl_ri*(tmpx - x*tmpsum);        \
    My[j][i] = My[i][j] = fl_y*tmp + fl_ri*(tmpy - y*tmpsum);        \
    Mz[j][i] = Mz[i][j] = fl_z*tmp + fl_ri*(tmpz - z*tmpsum);        \
  } while(0)

#define MAKE_ELEMENT_7(i,j)                                          \
  do {                                                               \
    const double dl0  = del0.m[i][j];                                \
    const double dl6  = gsl_##i    * gsl_##j;                        \
    const double dl6x = gsl_##i##x * gsl_##j + gsl_##i * gsl_##j##x; \
    const double dl6y = gsl_##i##y * gsl_##j + gsl_##i * gsl_##j##y; \
    const double dl6z = gsl_##i##z * gsl_##j + gsl_##i * gsl_##j##z; \
                                                                     \
    const double tmp  = w6*dl6  + w4*dl4  + w2*dl2 + w0*dl0;         \
    const double tmpx = w6*dl6x + w4*dl4x + w2*dl2x;                 \
    const double tmpy = w6*dl6y + w4*dl4y + w2*dl2y;                 \
    const double tmpz = w6*dl6z + w4*dl4z + w2*dl2z;                 \
    const double tmpsum = tmpx*x + tmpy*y + tmpz*z;                  \
    M [j][i] = M[i][j] =  fl  *tmp;                                  \
    Mx[j][i] = Mx[i][j] = fl_x*tmp + fl_ri*(tmpx - x*tmpsum);        \
    My[j][i] = My[i][j] = fl_y*tmp + fl_ri*(tmpy - y*tmpsum);        \
    Mz[j][i] = Mz[i][j] = fl_z*tmp + fl_ri*(tmpz - z*tmpsum);        \
  } while(0)
/* End of bond matrix macros */


/*
  Construction of bond matrix, and its derivatives
  with respect to the coordinates
*/
void PairMGPT::hamltn_5_raw(const double xin,const double yin,const double zin,
			     double M [8][8],double Mx[8][8],
			     double My[8][8],double Mz[8][8],
			     double *fl_deriv_sum_p) {

  const double r = sqrt(xin*xin + yin*yin + zin*zin),ri = 1.0/r;
  const double x = xin*ri,y = yin*ri,z = zin*ri;

  // d-d
  // call delndd(x,y,z)
  const double x2 = x*x,y2 = y*y,z2 = z*z;
  const double xy = x*y,xz = x*z,yz = y*z;

  const double sr3 = sqrt(3.0),sr3i = 1.0/sr3;
  const double frac_1_3 = 1.0/3.0,frac_2_3 = 2.0/3.0,frac_4_3 = 4.0/3.0;

  const double ddl_1 = ddl[1],ddl_2 = ddl[2],ddl_3 = ddl[3];
  const double w4 = ddl_1 - frac_4_3*ddl_2 + frac_1_3*ddl_3;
  const double w2 = ddl_2 - ddl_3;
  const double w0 = ddl_2;


  //del4
  double gsl_1 ,gsl_2 ,gsl_3 ,gsl_4 ,gsl_5;
  double gsl_1x,gsl_2x,gsl_3x,gsl_4x,gsl_5x;
  double gsl_1y,gsl_2y,gsl_3y,gsl_4y,gsl_5y;
  double gsl_1z,gsl_2z,gsl_3z,gsl_4z,gsl_5z;

  double dl2,dl2x,dl2y,dl2z;

  double fl,fl_x,fl_y,fl_z,fl_ri;
  double fl_rp,fl_p1,fl_r0,fl_al;

  gsl_1  = 0.5*(3.0*z2 - 1.0);
  gsl_1x = 0.0;
  gsl_1y = 0.0;
  gsl_1z = 3.0*z;

  gsl_2  = sr3*xz;
  gsl_2x = sr3*z;
  gsl_2y = 0.0;
  gsl_2z = sr3*x;

  gsl_3  = sr3*yz;
  gsl_3x = 0.0;
  gsl_3y = sr3*z;
  gsl_3z = sr3*y;

  gsl_4  = sr3*(x2 - y2)*0.5;
  gsl_4x = sr3*x;
  gsl_4y = -sr3*y;
  gsl_4z = 0.0;

  gsl_5  = sr3*xy;
  gsl_5x = sr3*y;
  gsl_5y = sr3*x;
  gsl_5z = 0.0;

  // Compute bond length potential
  fl_deriv_new(r,ri,x,y,z,fl,fl_x,fl_y,fl_z , fl_rp,fl_p1,fl_r0,fl_al);
  fl_ri = fl*ri;

  *fl_deriv_sum_p =
    fl_rp*splinepot.drp  + fl_p1*splinepot.dp1 +
    fl_r0*splinepot.dr00 + fl_al*splinepot.dal;

  // del2

  //del2.m[1][1] = z2 - 2.0/3.0;
  dl2 = z2 - frac_2_3;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 2*z;
  MAKE_ELEMENT_5(1,1);

  //del2.m[1][2] = xz/sr3;
  dl2  = xz*sr3i;
  dl2x = z*sr3i;
  dl2y = 0.0;
  dl2z = x*sr3i;
  MAKE_ELEMENT_5(1,2);

  //del2.m[1][3] = yz/sr3;
  dl2  = yz*sr3i;
  dl2x = 0.0;
  dl2y = z*sr3i;
  dl2z = y*sr3i;
  MAKE_ELEMENT_5(1,3);

  //del2.m[1][4] = -(x2 - y2)*sr3i;
  dl2  = -(x2 - y2)*sr3i;
  dl2x = -2.0*sr3i*x;
  dl2y =  2.0*sr3i*y;
  dl2z = 0.0;
  MAKE_ELEMENT_5(1,4);

  //del2.m[1][5] = -2.0*xy*sr3i;
  dl2 = -2.0*xy*sr3i;
  dl2x = -2.0*y*sr3i;
  dl2y = -2.0*x*sr3i;
  dl2z =  0.0;
  MAKE_ELEMENT_5(1,5);

  //del2.m[2][2] = -y2;
  dl2  = -y2;
  dl2x =  0.0;
  dl2y = -2.0*y;
  dl2z =  0.0;
  MAKE_ELEMENT_5(2,2);

  //del2.m[2][3] = xy;
  dl2  = xy;
  dl2x = y;
  dl2y = x;
  dl2z = 0.0;
  MAKE_ELEMENT_5(2,3);

  //del2.m[2][4] = xz;
  dl2  = xz;
  dl2x = z;
  dl2y = 0.0;
  dl2z = x;
  MAKE_ELEMENT_5(2,4);

  //del2.m[2][5] = yz;
  dl2  = yz;
  dl2x = 0.0;
  dl2y = z;
  dl2z = y;
  MAKE_ELEMENT_5(2,5);

  //del2.m[3][3] = -x2;
  dl2  = -x2;
  dl2x = -2.0*x;
  dl2y =  0.0;
  dl2z =  0.0;
  MAKE_ELEMENT_5(3,3);

  //del2.m[3][4] = -yz;
  dl2  = -yz;
  dl2x =  0.0;
  dl2y = -z;
  dl2z = -y;
  MAKE_ELEMENT_5(3,4);

  //del2.m[3][5] = xz;
  dl2  = xz;
  dl2x = z;
  dl2y = 0.0;
  dl2z = x;
  MAKE_ELEMENT_5(3,5);

  //del2.m[4][4] = -z2;
  dl2  = -z2;
  dl2x =  0.0;
  dl2y =  0.0;
  dl2z = -2.0*z;
  MAKE_ELEMENT_5(4,4);

  //del2.m[4][5] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;
  MAKE_ELEMENT_5(4,5);

  //del2.m[5][5] = -z2;
  dl2  = -z2;
  dl2x =  0.0;
  dl2y =  0.0;
  dl2z = -2.0*z;
  MAKE_ELEMENT_5(5,5);
}


void PairMGPT::hamltn_7_raw(const double xin,const double yin,const double zin,
		  double M [8][8],double Mx[8][8],
		  double My[8][8],double Mz[8][8],
			     double *fl_deriv_sum_p) {

  const double r = sqrt(xin*xin + yin*yin + zin*zin),ri = 1.0/r;
  const double x = xin*ri,y = yin*ri,z = zin*ri;

  // d-d
  // call delndd(x,y,z)
  const double x2 = x*x,y2 = y*y,z2 = z*z;
  const double xy = x*y,xz = x*z,yz = y*z;
  const double x4 = x2*x2,y4 = y2*y2;

  //const double sr3 = sqrt(3.0);//,sr3i = 1.0/sr3;
  //const double frac_1_3 = 1.0/3.0,frac_2_3 = 2.0/3.0,frac_4_3 = 4.0/3.0;

  const double sr01   = sqrt(0.1);
  const double sr015  = sqrt(0.15);
  const double sr024  = sqrt(0.24);
  const double sr0375 = sqrt(0.375);
  const double sr06   = sqrt(0.6);
  const double sr0625 = sqrt(0.625);
  const double sr09   = sqrt(0.9);
  const double sr15   = sqrt(1.5);
  const double sr24   = sqrt(2.4);
  const double sr36   = sqrt(3.6);
  const double sr375  = sqrt(3.75);
  const double sr96   = sqrt(9.6);
  const double sr150  = sqrt(15.0);



  const double ddl_1 = ddl[1],ddl_2 = ddl[2],ddl_3 = ddl[3],ddl_4 = ddl[4];
  const double w6 = ddl_1 - 1.5*ddl_2 + 0.6*ddl_3 - 0.1*ddl_4;
  const double w4 = 0.625*ddl_2 - ddl_3 + 0.375*ddl_4;
  const double w2 = 0.625*(ddl_2 - ddl_4);
  const double w0 = 0.625*ddl_2 + 0.375*ddl_4;


  //del6
  double gsl_1 ,gsl_2 ,gsl_3 ,gsl_4 ,gsl_5 ,gsl_6, gsl_7;
  double gsl_1x,gsl_2x,gsl_3x,gsl_4x,gsl_5x,gsl_6x,gsl_7x;
  double gsl_1y,gsl_2y,gsl_3y,gsl_4y,gsl_5y,gsl_6y,gsl_7y;
  double gsl_1z,gsl_2z,gsl_3z,gsl_4z,gsl_5z,gsl_6z,gsl_7z;

  double dl2,dl2x,dl2y,dl2z;
  double dl4,dl4x,dl4y,dl4z;
  double t1;

  double fl,fl_x,fl_y,fl_z,fl_ri;
  double fl_rp,fl_p1,fl_r0,fl_al;

  //gslf[1] = 0.5*(5.0*n2 - 3.0)*n;
  gsl_1  = 0.5*(5.0*z2 - 3.0)*z;
  gsl_1x = 0.0;
  gsl_1y = 0.0;
  gsl_1z = 7.5*z2 - 1.5;

  //gslf[2] = sr0375*(5.0*n2 - 1.0)*l;
  gsl_2  = sr0375*(5.0*z2 - 1.0)*x;
  gsl_2x = sr0375*(5.0*z2 - 1.0);
  gsl_2y = 0.0;
  gsl_2z = sr0375*10.0*xz;

  //gslf[3] = sr0375*(5.0*n2 - 1.0)*m;
  gsl_3  = sr0375*(5.0*z2 - 1.0)*y;
  gsl_3x = 0.0;
  gsl_3y = sr0375*(5.0*z2 - 1.0);
  gsl_3z = sr0375*10.0*yz;

  //gslf[4] = sr375*(l2 - m2)*n;
  gsl_4  =  sr375*(x2 - y2)*z;
  gsl_4x =  2.0*sr375*xz;
  gsl_4y = -2.0*sr375*yz;
  gsl_4z =  sr375*(x2 - y2);

  //gslf[5] = sr150*lm*n;
  gsl_5  = sr150*xy*z;
  gsl_5x = sr150*yz;
  gsl_5y = sr150*xz;
  gsl_5z = sr150*xy;

  //gslf[6] = sr0625*(l2 - 3.0*m2)*l;
  gsl_6  =  sr0625*(x2 - 3.0*y2)*x;
  gsl_6x =  3.0*sr0625*(x2 - y2);
  gsl_6y = -6.0*sr0625*xy;
  gsl_6z =  0.0;

  //gslf[7] = sr0625*(3.0*l2 - m2)*m;
  gsl_7  = sr0625*(3.0*x2 - y2)*y;
  gsl_7x = 6.0*sr0625*xy;
  gsl_7y = 3.0*sr0625*(x2 - y2);
  gsl_7z = 0.0;


  // Compute bond length potential
  fl_deriv_new(r,ri,x,y,z,fl,fl_x,fl_y,fl_z , fl_rp,fl_p1,fl_r0,fl_al);
  fl_ri = fl*ri;

  *fl_deriv_sum_p =
    fl_rp*splinepot.drp  + fl_p1*splinepot.dp1 +
    fl_r0*splinepot.dr00 + fl_al*splinepot.dal;

  // del2f

  //del2f.m[1][1] = 0.4*(3.0*n2 - 1.0);
  dl2  = 0.4*(3.0*z2 - 1.0);
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 2.4*z;

  //del4f.m[1][1] = 0.60*(5.0*n2 - 4.0)*n2;
  dl4  = 0.60*(5.0*z2 - 4.0)*z2;
  dl4x = 0.0;
  dl4y = 0.0;
  dl4z = 0.60*(20.0*z2 - 8.0)*z;

  MAKE_ELEMENT_7(1,1);


  //del2f.m[1][2] = sr024*ln;
  dl2  = sr024*xz;
  dl2x = sr024*z;
  dl2y = 0.0;
  dl2z = sr024*x;

  //del4f.m[1][2] = sr024*(5.0*n2 - 2.0)*ln;
  dl4  = sr024*(5.0*z2 - 2.0)*xz;
  dl4x = sr024*(5.0*z2 - 2.0)*z;
  dl4y = 0.0;
  dl4z = sr024*(15.0*z2 - 2.0)*x;

  MAKE_ELEMENT_7(1,2);


  //del2f.m[1][3] = sr024*mn;
  dl2  = sr024*yz;
  dl2x = 0.0;
  dl2y = sr024*z;
  dl2z = sr024*y;

  //del4f.m[1][3] = sr024*(5.0*n2 - 2.0)*mn;
  dl4  = sr024*(5.0*z2 - 2.0)*yz;
  dl4x = 0.0;
  dl4y = sr024*(5.0*z2 - 2.0)*z;
  dl4z = sr024*(15.0*z2 - 2.0)*y;

  MAKE_ELEMENT_7(1,3);


  //del2f.m[1][4] = -sr06*(l2 - m2);
  dl2  = -sr06*(x2 - y2);
  dl2x = -2.0*sr06*x;
  dl2y =  2.0*sr06*y;
  dl2z = 0.0;

  //del4f.m[1][4] = -sr06*(l2 - m2)*n2;
  dl4  = -sr06*(x2 - y2)*z2;
  dl4x = -2.0*sr06*x*z2;
  dl4y =  2.0*sr06*y*z2;
  dl4z = -2.0*sr06*(x2 - y2)*z;

  MAKE_ELEMENT_7(1,4);


  //del2f.m[1][5] = -sr24*lm;
  dl2  = -sr24*xy;
  dl2x = -sr24*y;
  dl2y = -sr24*x;
  dl2z =  0.0;

  //del4f.m[1][5] = -sr24*lm*n2;
  dl4  = -sr24*xy*z2;
  dl4x = -sr24*y*z2;
  dl4y = -sr24*x*z2;
  dl4z = -2.0*sr24*xy*z;

  MAKE_ELEMENT_7(1,5);


  //del2f.m[1][6] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;

  //del4f.m[1][6] = sr36*(3.0*m2 - l2)*ln;
  dl4  = sr36*(3.0*y2 - x2)*xz;
  dl4x = 3.0*sr36*(y2 - x2)*z;
  dl4y = 6.0*sr36*y*xz;
  dl4z = sr36*(3.0*y2 - x2)*x;

  MAKE_ELEMENT_7(1,6);


  //del2f.m[1][7] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;

  //del4f.m[1][7] = -sr36*(3.0*l2 - m2)*mn;
  dl4  = -sr36*(3.0*x2 - y2)*yz;
  dl4x = -6.0*sr36*x*yz;
  dl4y = -3.0*sr36*(x2 - y2)*z;
  dl4z = -sr36*(3.0*x2 - y2)*y;

  MAKE_ELEMENT_7(1,7);


  //del2f.m[2][2] = 0.3*(1.0 - 4.0*m2 + n2);
  dl2  =  0.3 - 1.2*y2 + 0.3*z2;
  dl2x =  0.0;
  dl2y = -2.4*y;
  dl2z =  0.6*z;

  //del4f.m[2][2] = -0.4*l*l - 2.5*(m2 - 0.6*l2)*n2;
  dl4  = -0.4*x2 - 2.5*(y2 - 0.6*x2)*z2;
  dl4x = -0.8*x + 3.0*x*z2;
  dl4y = -5.0*y*z2;
  dl4z = -5.0*(y2 - 0.6*x2)*z;

  MAKE_ELEMENT_7(2,2);


  //del2f.m[2][3] = 1.2*lm;
  dl2  = 1.2*xy;
  dl2x = 1.2*y;
  dl2y = 1.2*x;
  dl2z = 0.0;

  //del4f.m[2][3] = 0.4*(10.0*n2 - 1.0)*lm;
  dl4  = (4.0*z2 - 0.4)*xy;
  dl4x = (4.0*z2 - 0.4)*y;
  dl4y = (4.0*z2 - 0.4)*x;
  dl4z = 8.0*z*xy;

  MAKE_ELEMENT_7(2,3);


  //del2f.m[2][4] = sr09*ln;
  dl2  = sr09*xz;
  dl2x = sr09*z;
  dl2y = 0.0;
  dl2z = sr09*x;

  //del4f.m[2][4] = sr01*(6.0*n2 - 8.0*m2 - 1.0)*ln;
  dl4  =  sr01*(6.0*z2 - 8.0*y2 - 1.0)*xz;
  dl4x =  sr01*(6.0*z2 - 8.0*y2 - 1.0)*z;
  dl4y = -16.0*sr01*y*xz;
  dl4z =  sr01*(18.0*z2 - 8.0*y2 - 1.0)*x;

  MAKE_ELEMENT_7(2,4);


  //del2f.m[2][5] = sr09*mn;
  dl2  = sr09*yz;
  dl2x = 0.0;
  dl2y = sr09*z;
  dl2z = sr09*y;

  //del4f.m[2][5] = sr01*(2.0*n2 - 8.0*m2 + 3.0)*mn;
  dl4  = sr01*(2.0*z2 - 8.0*y2 + 3.0)*yz;
  dl4x = 0.0;
  dl4y = sr01*(2.0*z2 - 24.0*y2 + 3.0)*z;
  dl4z = sr01*(6.0*z2 - 8.0*y2 + 3.0)*y;

  MAKE_ELEMENT_7(2,5);


  //del2f.m[2][6] = -sr015*(l2 - m2);
  dl2  = -sr015*(x2 - y2);
  dl2x = -2.0*sr015*x;
  dl2y =  2.0*sr015*y;
  dl2z =  0.0;

  //del4f.m[2][6] = sr375*(l2 - m2 - 1.4*l4 + 1.2*l2*m2 + m4);
  dl4  = sr375*(x2 - y2 - 1.4*x4 + 1.2*x2*y2 + y4);
  dl4x = sr375*(2.0 - 5.6*x2 + 2.4*y2)*x;
  dl4y = sr375*(-2.0 + 2.4*x2 + 4.0*y2)*y;
  dl4z = 0.0;

  MAKE_ELEMENT_7(2,6);


  //del2f.m[2][7] = -sr06*lm;
  dl2  = -sr06*xy;
  dl2x = -sr06*y;
  dl2y = -sr06*x;
  dl2z =  0.0;

  //del4f.m[2][7] = sr96*(n2 - l2 + 0.25)*lm;
  dl4  = sr96*(z2 - x2 + 0.25)*xy;
  dl4x = sr96*(z2 - 3.0*x2 + 0.25)*y;
  dl4y = sr96*(z2 - x2 + 0.25)*x;
  dl4z = 2.0*sr96*z*xy;

  MAKE_ELEMENT_7(2,7);


  //del2f.m[3][3] = 0.30*(1.0 - 4.0*l2 + n2);
  dl2  =  0.3 - 1.2*x2 + 0.3*z2;
  dl2x = -2.4*x;
  dl2y =  0.0;
  dl2z =  0.6*z;

  //del4f.m[3][3] = -0.4*m2 - 2.5*(l2 - 0.6*m2)*n2;
  dl4  = -0.4*y2 - 2.5*(x2 - 0.6*y2)*z2;
  dl4x = -5.0*x*z2;
  dl4y = y*(3.0*z2 - 0.8);
  dl4z = -5.0*(x2 - 0.6*y2)*z;

  MAKE_ELEMENT_7(3,3);


  //del2f.m[3][4] = -sr09*mn;
  dl2  = -sr09*yz;
  dl2x =  0.0;
  dl2y = -sr09*z;
  dl2z = -sr09*y;

  //del4f.m[3][4] = -sr01*(6.0*n2 - 8.0*l2 - 1.0)*mn;
  dl4  = -sr01*(6.0*z2 - 8.0*x2 - 1.0)*yz;
  dl4x =  16.0*sr01*x*yz;
  dl4y = -sr01*(6.0*z2 - 8.0*x2 - 1.0)*z;
  dl4z = -sr01*(18.0*z2 - 8.0*x2 - 1.0)*y;

  MAKE_ELEMENT_7(3,4);


  //del2f.m[3][5] = sr09*ln;
  dl2  = sr09*xz;
  dl2x = sr09*z;
  dl2y = 0.0;
  dl2z = sr09*x;

  //del4f.m[3][5] = sr01*(2.0*n2 - 8.0*l2 + 3.0)*ln;
  dl4  = sr01*(2.0*z2 - 8.0*x2 + 3.0)*xz;
  dl4x = sr01*(2.0*z2 - 24.0*x2 + 3.0)*z;
  dl4y = 0.0;
  dl4z = sr01*(6.0*z2 - 8.0*x2 + 3.0)*x;

  MAKE_ELEMENT_7(3,5);


  //del2f.m[3][6] = sr06*lm;
  dl2  = sr06*xy;
  dl2x = sr06*y;
  dl2y = sr06*x;
  dl2z = 0.0;

  //del4f.m[3][6] = sr96*(m2 - n2 - 0.25)*lm;
  dl4  =  sr96*(y2 - z2 - 0.25)*xy;
  dl4x =  sr96*(y2 - z2 - 0.25)*y;
  dl4y =  sr96*(3.0*y2 - z2 - 0.25)*x;
  dl4z = -2.0*sr96*z*xy;

  MAKE_ELEMENT_7(3,6);


  //del2f.m[3][7] = -sr015*(l2 - m2);
  dl2  = -sr015*(x2 - y2);
  dl2x = -2.0*sr015*x;
  dl2y =  2.0*sr015*y;
  dl2z =  0.0;

  //del4f.m[3][7] = sr375*(l2 - m2 + 1.4*m4 - 1.2*l2*m2 - l4);
  dl4  = sr375*(x2 - y2 + 1.4*y4 - 1.2*x2*y2 - x4);
  dl4x = sr375*(2.0 - 2.4*y2 - 4.0*x2)*x;
  dl4y = sr375*(-2.0 + 5.6*y2 - 2.4*x2)*y;
  dl4z = 0.0;

  MAKE_ELEMENT_7(3,7);


  //del2f.m[4][4] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;

  //del4f.m[4][4] = (2.0 - 3.0*n2)*n2 - 4.0*l2*m2;
  dl4  = (2.0 - 3.0*z2)*z2 - 4.0*x2*y2;
  dl4x = -8.0*x*y2;
  dl4y = -8.0*x2*y;
  dl4z =  (4.0 - 12.0*z2)*z;

  MAKE_ELEMENT_7(4,4);

  //del2f.m[4][5] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;

  //del4f.m[4][5] = 2.0*(l2 - m2)*lm;
  dl4  = 2.0*(x2 - y2)*xy;
  dl4x = 2.0*(3.0*x2 - y2)*y;
  dl4y = 2.0*(x2 - 3.0*y2)*x;
  dl4z = 0.0;

  MAKE_ELEMENT_7(4,5);


  //del2f.m[4][6] = sr15*ln;
  dl2  = sr15*xz;
  dl2x = sr15*z;
  dl2y = 0.0;
  dl2z = sr15*x;

  //del4f.m[4][6] = -sr15*(2.0*n2 - 1.0)*ln;
  dl4  = -sr15*(2.0*z2 - 1.0)*xz;
  dl4x = -sr15*(2.0*z2 - 1.0)*z;
  dl4y =  0.0;
  dl4z = -sr15*(6.0*z2 - 1.0)*x;

  MAKE_ELEMENT_7(4,6);


  //del2f.m[4][7] = sr15*mn;
  dl2  = sr15*yz;
  dl2x = 0.0;
  dl2y = sr15*z;
  dl2z = sr15*y;

  //del4f.m[4][7] = -sr15*(2.0*n2 - 1.0)*mn;
  dl4  = -sr15*(2.0*z2 - 1.0)*yz;
  dl4x =  0.0;
  dl4y = -sr15*(2.0*z2 - 1.0)*z;
  dl4z = -sr15*(6.0*z2 - 1.0)*y;

  MAKE_ELEMENT_7(4,7);


  //del2f.m[5][5] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;

  //del4f.m[5][5] = -pow((2.0*n2 - 1.0),2) + 4.0*l2*m2;
  t1 = 2.0*z2 - 1.0;
  dl4  = -t1*t1 + 4.0*x2*y2;
  dl4x =  8.0*x*y2;
  dl4y =  8.0*x2*y;
  dl4z = -8.0*t1*z;

  MAKE_ELEMENT_7(5,5);


  //del2f.m[5][6] = -sr15*mn;
  dl2  = -sr15*yz;
  dl2x =  0.0;
  dl2y = -sr15*z;
  dl2z = -sr15*y;

  //del4f.m[5][6] = sr15*(2.0*n2 - 1.0)*mn;
  dl4  = sr15*(2.0*z2 - 1.0)*yz;
  dl4x = 0.0;
  dl4y = sr15*(2.0*z2 - 1.0)*z;
  dl4z = sr15*(6.0*z2 - 1.0)*y;

  MAKE_ELEMENT_7(5,6);

  //del2f.m[5][7] = sr15*ln;
  dl2  = sr15*xz;
  dl2x = sr15*z;
  dl2y = 0.0;
  dl2z = sr15*x;

  //del4f.m[5][7] = -sr15*(2.0*n2 - 1.0)*ln;
  dl4  = -sr15*(2.0*z2 - 1.0)*xz;
  dl4x = -sr15*(2.0*z2 - 1.0)*z;
  dl4y =  0.0;
  dl4z = -sr15*(6.0*z2 - 1.0)*x;

  MAKE_ELEMENT_7(5,7);


  //del2f.m[6][6] = -(3.0*n2 - 1.0)/2.0;
  dl2  =  0.5 - 1.5*z2;
  dl2x =  0.0;
  dl2y =  0.0;
  dl2z = -3.0*z;

  //del4f.m[6][6] = 1.5*(n2 - 1.0)*n2;
  dl4  = (1.5*z2 - 1.5)*z2;
  dl4x = 0.0;
  dl4y = 0.0;
  dl4z = (6.0*z2 - 3.0)*z;

  MAKE_ELEMENT_7(6,6);


  //del2f.m[6][7] = 0.0;
  dl2  = 0.0;
  dl2x = 0.0;
  dl2y = 0.0;
  dl2z = 0.0;

  //del4f.m[6][7] = 0.0;
  dl4  = 0.0;
  dl4x = 0.0;
  dl4y = 0.0;
  dl4z = 0.0;

  MAKE_ELEMENT_7(6,7);


  //del2f.m[7][7] = -(3.0*n2 - 1.0)/2.0;
  dl2  =  0.5 - 1.5*z2;
  dl2x =  0.0;
  dl2y =  0.0;
  dl2z = -3.0*z;

  //del4f.m[7][7] = 1.5*(n2 - 1.0)*n2;
  dl4  = (1.5*z2 - 1.5)*z2;
  dl4x = 0.0;
  dl4y = 0.0;
  dl4z = (6.0*z2 - 3.0)*z;

  MAKE_ELEMENT_7(7,7);
}

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

/* ----------------------------------------------------------------------
 * Fast Model Generalized Pseudopotential Theory (MGPT) interatomic
 * potential routine.
 *
 * Copyright (2015) Lawrence Livermore National Security, LLC.
 * Produced at the Lawrence Livermore National Laboratory.
 * Written by Tomas Oppelstrup (oppelstrup2@llnl.gov) and John Moriarty
 * (moriarty2@llnl.gov)
 * LLNL-CODE-674031  All rights reserved.
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License (as published by the
 * Free Software Foundation) version 2, dated June 1991.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE. See the terms and conditions of the
 * GNU General Public License for more details.
 *
 * LLNL Preamble Notice
 * A. This notice is required to be provided under our contract with the
 * U.S. Department of Energy (DOE). This work was performed under the auspices
 * of the DOE by Lawrence Livermore National Laboratory under Contract No.
 * DE-AC52-07NA27344.
 *
 * B. Neither the United States Government nor Lawrence Livermore National
 * Security, LLC nor any of their employees, makes any warranty, express or
 * implied, or assumes any liability or responsibility for the accuracy,
 * completeness, or usefulness of any information, apparatus, product, or
 * process disclosed, or represents that its use would not infringe
 * privately-owned rights.
 *
 * C. Also, reference herein to any specific commercial products, process,
 * or services by trade name, trademark, manufacturer or otherwise does not
 * necessarily constitute or imply its endorsement, recommendation, or
 * favoring by the United States Government or Lawrence Livermore National
 * Security, LLC. The views and opinions of authors expressed herein do not
 * necessarily state or reflect those of the United States Government or
 * Lawrence Livermore National Security, LLC, and shall not be used for
 * advertising or product endorsement purposes.
------------------------------------------------------------------------- */