Diffusion lammps (lammm-master-deprecated) Browse

pair_mgpt.cpp
38356949ba42lammm-master-deprecated

pair_mgpt.cpp
View Options

			/* ----------------------------------------------------------------------
			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[inj+j] = array[inj+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 && bik != 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,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 = rmaxrmax,rcut2_bond = rcritrcrit,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] + rhoinve_scalevolvir2;
			vatom[i][pp] -= 0.5 * rhoinve_scaledvir*w2;
			vatom[j][pp] -= 0.5 * rhoinve_scaledvir*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)tr0w3/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)tr1w3/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)tr2w3/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)tr3w3/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.vcv33 + splinepot.vdv43) * e_scale * w3;
			e_triplet = e_triplet + de_triplet;
			e_triplet_c = e_triplet_c + splinepot.vcv33 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) * rhoinve_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 )tr1w4/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 )tr2w4/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 )tr3w4/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 * rhoinve_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] - rhoinvsplinepot.devol0e_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",
			tmem1e3,ntmem,memadj1e3,memadj/ntmem*1e9);
			printf(" sort = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			tsort1e3,ntsort,sortadj1e3,sortadj/ntsort*1e9);
			printf(" pair = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			tpair1e3,ntpair,pairadj1e3,pairadj/ntpair*1e9);
			printf(" lookup = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			tlookup1e3,ntlookup,lookupadj1e3,lookupadj/ntlookup*1e9);
			printf(" triplet = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			ttriplet1e3,nttriplet,tripletadj1e3,tripletadj/nttriplet*1e9);
			printf(" quad = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			tquad1e3,ntquaditer,quadadj1e3,quadadj/ntquaditer*1e9);
			printf(" sum = %10.3fms adj = %10.3fms\n",
			tsum1e3,(tsum - adjnsum)*1e3);
			printf("\n make_b = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			t_make_b1e3,n_make,make_b_adj1e3,make_b_adj/n_make*1e9);
			printf(" make_b2 = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n",
			t_make_b21e3,n_make_b2,make_b2_adj1e3,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_t1e3,n_make,make_t_adj1e3,make_t_adj/n_make*1e9);
			printf(" trace = %10.3fms n = %8.0f adj = %10.3fms one = %10.3fns\n\n",
			t_trace1e3,n_trace,trace_adj1e3,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-adjnbc)1e3,(tbm-adjnbc)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] + nlocpot_input_volpvol0*e_scale;
			virial[i] = virial[i] - nlocrhoinvsplinepot.devol0*e_scale;
			}

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

			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_cpc + (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_vole_quadpe;
			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)
			{
			if(eflag \|\| vflag) ev_setup(eflag, vflag);
			else evflag = vflag_fdotr = eflag_global = vflag_global = eflag_atom = vflag_atom = 0;

			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[3nmax];
			double ffloc = new double[3listfull->inum];
			double ffloc2 = new double[3listfull->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[3ii+p] = -(e1-e2)/(2delta);

			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[3ii+0],ffloc[3ii+1],ffloc[3*ii+2]);
			printf(" Boundary %20.10e %20.10e %20.10e\n",
			ffloc2[3ii+0],ffloc2[3ii+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(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->requests[irequest_half]->half = 0;
			neighbor->requests[irequest_half]->half_from_full = 1;
			neighbor->requests[irequest_half]->otherlist = irequest_full;

			}

			/* ----------------------------------------------------------------------
			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 = -rp1_rpirpi;
			t_p1_ti = 1.0 - r*rpi;
			} else {
			double pni = 1.0/pn;
			double rprpn = pow(r*rpi,pn);
			t = exp(-p1pni(rprpn - 1.0));
			s = -p1rprpnri * t;
			t_rp_ti = p1rprpnrpi;
			t_p1_ti = pni - pnirprpn;// -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 = aluu;
			s = 2alu*r0i;
			qx = s * xhat;
			qy = s * yhat;
			qz = s * zhat;
			q_r0 = -2.0alurr0i*r0i;
			q_al = u*u;

			exp_q = exp(-q);
			if(mode <= 2) {
			fl_0 = exp_q * t;
			fl_x = exp_q(tx - tqx);
			fl_y = exp_q(ty - tqy);
			fl_z = exp_q(tz - tqz);

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

			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 = w4dl4 + w2dl2 + w0*dl0; \
			const double tmpx = w4dl4x + w2dl2x; \
			const double tmpy = w4dl4y + w2dl2y; \
			const double tmpz = w4dl4z + w2dl2z; \
			const double tmpsum = tmpxx + tmpyy + tmpz*z; \
			M [i][j] = fl *tmp; \
			Mx[i][j] = fl_xtmp + fl_ri(tmpx - x*tmpsum); \
			My[i][j] = fl_ytmp + fl_ri(tmpy - y*tmpsum); \
			Mz[i][j] = fl_ztmp + fl_ri(tmpz - z*tmpsum); \
			if (i != j) { \
			M[j][i] = M[i][j]; Mx[j][i] = Mx[i][j]; \
			My[j][i] = My[i][j]; Mz[j][i] = Mz[i][j]; \
			} \
			} 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 = w6dl6 + w4dl4 + w2dl2 + w0dl0; \
			const double tmpx = w6dl6x + w4dl4x + w2*dl2x; \
			const double tmpy = w6dl6y + w4dl4y + w2*dl2y; \
			const double tmpz = w6dl6z + w4dl4z + w2*dl2z; \
			const double tmpsum = tmpxx + tmpyy + tmpz*z; \
			M [i][j] = fl *tmp; \
			Mx[i][j] = fl_xtmp + fl_ri(tmpx - x*tmpsum); \
			My[i][j] = fl_ytmp + fl_ri(tmpy - y*tmpsum); \
			Mz[i][j] = fl_ztmp + fl_ri(tmpz - z*tmpsum); \
			if (i != j) { \
			M[j][i] = M[i][j]; Mx[j][i] = Mx[i][j]; \
			My[j][i] = My[i][j]; Mz[j][i] = Mz[i][j]; \
			} \
			} 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(xinxin + yinyin + zin*zin),ri = 1.0/r;
			const double x = xinri,y = yinri,z = zin*ri;

			// d-d
			// call delndd(x,y,z)
			const double x2 = xx,y2 = yy,z2 = z*z;
			const double xy = xy,xz = xz,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_3ddl_2 + frac_1_3ddl_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.0z2 - 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_rpsplinepot.drp + fl_p1splinepot.dp1 +
			fl_r0splinepot.dr00 + fl_alsplinepot.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.0sr3ix;
			dl2y = 2.0sr3iy;
			dl2z = 0.0;
			MAKE_ELEMENT_5(1,4);

			//del2.m[1][5] = -2.0xysr3i;
			dl2 = -2.0xysr3i;
			dl2x = -2.0ysr3i;
			dl2y = -2.0xsr3i;
			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(xinxin + yinyin + zin*zin),ri = 1.0/r;
			const double x = xinri,y = yinri,z = zin*ri;

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

			//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.5ddl_2 + 0.6ddl_3 - 0.1*ddl_4;
			const double w4 = 0.625ddl_2 - ddl_3 + 0.375ddl_4;
			const double w2 = 0.625*(ddl_2 - ddl_4);
			const double w0 = 0.625ddl_2 + 0.375ddl_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.0n2 - 3.0)*n;
			gsl_1 = 0.5(5.0z2 - 3.0)*z;
			gsl_1x = 0.0;
			gsl_1y = 0.0;
			gsl_1z = 7.5*z2 - 1.5;

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

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

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

			//gslf[5] = sr150lmn;
			gsl_5 = sr150xyz;
			gsl_5x = sr150*yz;
			gsl_5y = sr150*xz;
			gsl_5z = sr150*xy;

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

			//gslf[7] = sr0625(3.0l2 - m2)*m;
			gsl_7 = sr0625(3.0x2 - y2)*y;
			gsl_7x = 6.0sr0625xy;
			gsl_7y = 3.0sr0625(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_rpsplinepot.drp + fl_p1splinepot.dp1 +
			fl_r0splinepot.dr00 + fl_alsplinepot.dal;

			// del2f

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

			//del4f.m[1][1] = 0.60(5.0n2 - 4.0)*n2;
			dl4 = 0.60(5.0z2 - 4.0)*z2;
			dl4x = 0.0;
			dl4y = 0.0;
			dl4z = 0.60(20.0z2 - 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.0n2 - 2.0)*ln;
			dl4 = sr024(5.0z2 - 2.0)*xz;
			dl4x = sr024(5.0z2 - 2.0)*z;
			dl4y = 0.0;
			dl4z = sr024(15.0z2 - 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.0n2 - 2.0)*mn;
			dl4 = sr024(5.0z2 - 2.0)*yz;
			dl4x = 0.0;
			dl4y = sr024(5.0z2 - 2.0)*z;
			dl4z = sr024(15.0z2 - 2.0)*y;

			MAKE_ELEMENT_7(1,3);


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

			//del4f.m[1][4] = -sr06(l2 - m2)n2;
			dl4 = -sr06(x2 - y2)z2;
			dl4x = -2.0sr06x*z2;
			dl4y = 2.0sr06y*z2;
			dl4z = -2.0sr06(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] = -sr24lmn2;
			dl4 = -sr24xyz2;
			dl4x = -sr24yz2;
			dl4y = -sr24xz2;
			dl4z = -2.0sr24xy*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.0m2 - l2)*ln;
			dl4 = sr36(3.0y2 - x2)*xz;
			dl4x = 3.0sr36(y2 - x2)*z;
			dl4y = 6.0sr36y*xz;
			dl4z = sr36(3.0y2 - 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.0l2 - m2)*mn;
			dl4 = -sr36(3.0x2 - y2)*yz;
			dl4x = -6.0sr36x*yz;
			dl4y = -3.0sr36(x2 - y2)*z;
			dl4z = -sr36(3.0x2 - y2)*y;

			MAKE_ELEMENT_7(1,7);


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

			//del4f.m[2][2] = -0.4ll - 2.5(m2 - 0.6l2)*n2;
			dl4 = -0.4x2 - 2.5(y2 - 0.6x2)z2;
			dl4x = -0.8x + 3.0x*z2;
			dl4y = -5.0yz2;
			dl4z = -5.0(y2 - 0.6x2)*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.0n2 - 1.0)*lm;
			dl4 = (4.0z2 - 0.4)xy;
			dl4x = (4.0z2 - 0.4)y;
			dl4y = (4.0z2 - 0.4)x;
			dl4z = 8.0zxy;

			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.0n2 - 8.0m2 - 1.0)ln;
			dl4 = sr01(6.0z2 - 8.0y2 - 1.0)xz;
			dl4x = sr01(6.0z2 - 8.0y2 - 1.0)z;
			dl4y = -16.0sr01y*xz;
			dl4z = sr01(18.0z2 - 8.0y2 - 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.0n2 - 8.0m2 + 3.0)mn;
			dl4 = sr01(2.0z2 - 8.0y2 + 3.0)yz;
			dl4x = 0.0;
			dl4y = sr01(2.0z2 - 24.0y2 + 3.0)z;
			dl4z = sr01(6.0z2 - 8.0y2 + 3.0)y;

			MAKE_ELEMENT_7(2,5);


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

			//del4f.m[2][6] = sr375(l2 - m2 - 1.4l4 + 1.2l2m2 + m4);
			dl4 = sr375(x2 - y2 - 1.4x4 + 1.2x2y2 + y4);
			dl4x = sr375(2.0 - 5.6x2 + 2.4y2)x;
			dl4y = sr375(-2.0 + 2.4x2 + 4.0y2)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.0x2 + 0.25)*y;
			dl4y = sr96(z2 - x2 + 0.25)x;
			dl4z = 2.0sr96z*xy;

			MAKE_ELEMENT_7(2,7);


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

			//del4f.m[3][3] = -0.4m2 - 2.5(l2 - 0.6m2)n2;
			dl4 = -0.4y2 - 2.5(x2 - 0.6y2)z2;
			dl4x = -5.0xz2;
			dl4y = y(3.0z2 - 0.8);
			dl4z = -5.0(x2 - 0.6y2)*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.0n2 - 8.0l2 - 1.0)mn;
			dl4 = -sr01(6.0z2 - 8.0x2 - 1.0)yz;
			dl4x = 16.0sr01x*yz;
			dl4y = -sr01(6.0z2 - 8.0x2 - 1.0)z;
			dl4z = -sr01(18.0z2 - 8.0x2 - 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.0n2 - 8.0l2 + 3.0)ln;
			dl4 = sr01(2.0z2 - 8.0x2 + 3.0)xz;
			dl4x = sr01(2.0z2 - 24.0x2 + 3.0)z;
			dl4y = 0.0;
			dl4z = sr01(6.0z2 - 8.0x2 + 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.0y2 - z2 - 0.25)*x;
			dl4z = -2.0sr96z*xy;

			MAKE_ELEMENT_7(3,6);


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

			//del4f.m[3][7] = sr375(l2 - m2 + 1.4m4 - 1.2l2m2 - l4);
			dl4 = sr375(x2 - y2 + 1.4y4 - 1.2x2y2 - x4);
			dl4x = sr375(2.0 - 2.4y2 - 4.0x2)x;
			dl4y = sr375(-2.0 + 5.6y2 - 2.4x2)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.0n2)n2 - 4.0l2m2;
			dl4 = (2.0 - 3.0z2)z2 - 4.0x2y2;
			dl4x = -8.0xy2;
			dl4y = -8.0x2y;
			dl4z = (4.0 - 12.0z2)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.0x2 - y2)*y;
			dl4y = 2.0(x2 - 3.0y2)*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.0n2 - 1.0)*ln;
			dl4 = -sr15(2.0z2 - 1.0)*xz;
			dl4x = -sr15(2.0z2 - 1.0)*z;
			dl4y = 0.0;
			dl4z = -sr15(6.0z2 - 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.0n2 - 1.0)*mn;
			dl4 = -sr15(2.0z2 - 1.0)*yz;
			dl4x = 0.0;
			dl4y = -sr15(2.0z2 - 1.0)*z;
			dl4z = -sr15(6.0z2 - 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.0n2 - 1.0),2) + 4.0l2*m2;
			t1 = 2.0*z2 - 1.0;
			dl4 = -t1t1 + 4.0x2*y2;
			dl4x = 8.0xy2;
			dl4y = 8.0x2y;
			dl4z = -8.0t1z;

			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.0n2 - 1.0)*mn;
			dl4 = sr15(2.0z2 - 1.0)*yz;
			dl4x = 0.0;
			dl4y = sr15(2.0z2 - 1.0)*z;
			dl4z = sr15(6.0z2 - 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.0n2 - 1.0)*ln;
			dl4 = -sr15(2.0z2 - 1.0)*xz;
			dl4x = -sr15(2.0z2 - 1.0)*z;
			dl4y = 0.0;
			dl4z = -sr15(6.0z2 - 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.5z2 - 1.5)z2;
			dl4x = 0.0;
			dl4y = 0.0;
			dl4z = (6.0z2 - 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.5z2 - 1.5)z2;
			dl4x = 0.0;
			dl4y = 0.0;
			dl4z = (6.0z2 - 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.
			------------------------------------------------------------------------- */

lammps/src/USER-MGPT/pair_mgpt.cpp38356949ba42lammm-master-deprecated

pair_mgpt.cppView Options

lammps/src/USER-MGPT/pair_mgpt.cpp
38356949ba42lammm-master-deprecated

pair_mgpt.cpp
View Options