pair_airebo_omp.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, Jun 27, 05:18

pair_airebo_omp.cpp
View Options

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

	This software is distributed under the GNU General Public License.

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

	/* ----------------------------------------------------------------------
	Contributing author: Axel Kohlmeyer (Temple U)
	------------------------------------------------------------------------- */

	#include <math.h>
	#include "pair_airebo_omp.h"
	#include "atom.h"
	#include "comm.h"
	#include "error.h"
	#include "force.h"
	#include "memory.h"
	#include "my_page.h"
	#include "math_special.h"
	#include "neighbor.h"
	#include "neigh_list.h"

	#if defined(_OPENMP)
	#include <omp.h>
	#endif

	#include "suffix.h"
	using namespace LAMMPS_NS;
	using namespace MathSpecial;

	#define TOL 1.0e-9

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

	PairAIREBOOMP::PairAIREBOOMP(LAMMPS *lmp) :
	PairAIREBO(lmp), ThrOMP(lmp, THR_PAIR)
	{
	suffix_flag \|= Suffix::OMP;
	respa_enable = 0;
	}

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

	void PairAIREBOOMP::compute(int eflag, int vflag)
	{
	if (eflag \|\| vflag) {
	ev_setup(eflag,vflag);
	} else evflag = vflag_fdotr = vflag_atom = 0;

	REBO_neigh_thr();

	const int nall = atom->nlocal + atom->nghost;
	const int nthreads = comm->nthreads;
	const int inum = list->inum;

	#if defined(_OPENMP)
	#pragma omp parallel default(none) shared(eflag,vflag)
	#endif
	{
	int ifrom, ito, tid;

	loop_setup_thr(ifrom, ito, tid, inum, nthreads);
	ThrData *thr = fix->get_thr(tid);
	thr->timer(Timer::START);
	ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr);

	FREBO_thr(ifrom,ito,evflag,eflag,vflag_atom,thr);
	if (ljflag) FLJ_thr(ifrom,ito,evflag,eflag,vflag_atom,thr);
	if (torflag) TORSION_thr(ifrom,ito,evflag,eflag,thr);

	thr->timer(Timer::PAIR);
	reduce_thr(this, eflag, vflag, thr);
	} // end of omp parallel region
	}

	/* ----------------------------------------------------------------------
	Bij function
	------------------------------------------------------------------------- */

	double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
	double VA, int vflag_atom, ThrData * const thr)
	{
	int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
	int itype,jtype,ktype,ltype,ntype;
	double rik[3],rjl[3],rkn[3],rji[3],rki[3],rlj[3],rknmag,dNki,dwjl,bij;
	double NijC,NijH,NjiC,NjiH,wik,dwik,dwkn,wjl;
	double rikmag,rjlmag,cosjik,cosijl,g,tmp2,tmp3;
	double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ,Nki,Nlj,dS;
	double lamdajik,lamdaijl,dgdc,dgdN,pji,Nijconj,piRC;
	double dcosjikdri[3],dcosijldri[3],dcosjikdrk[3];
	double dN2[2],dN3[3];
	double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3];
	double Tij;
	double r32[3],r32mag,cos321,r43[3],r13[3];
	double dNlj;
	double om1234,rln[3];
	double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag;
	double w21,dw21,r34[3],r34mag,cos234,w34,dw34;
	double cross321[3],cross234[3],prefactor,SpN;
	double fcijpc,fcikpc,fcjlpc,fcjkpc,fcilpc;
	double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa1,aaa2,at2,cw,cwnum,cwnom;
	double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2;
	double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i;
	double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij;
	double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3];
	double f1[3],f2[3],f3[3],f4[4];
	double dcut321,PijS,PjiS;
	double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp;
	int REBO_neighs,REBO_neighs_i,REBO_neighs_j,REBO_neighs_k,*REBO_neighs_l;

	const double * const * const x = atom->x;
	double * const * const f = thr->get_f();
	int *type = atom->type;

	atomi = i;
	atomj = j;
	itype = map[type[i]];
	jtype = map[type[j]];
	wij = Sp(rijmag,rcmin[itype][jtype],rcmax[itype][jtype],dwij);
	NijC = nC[i]-(wij*kronecker(jtype,0));
	NijH = nH[i]-(wij*kronecker(jtype,1));
	NjiC = nC[j]-(wij*kronecker(itype,0));
	NjiH = nH[j]-(wij*kronecker(itype,1));
	bij = 0.0;
	tmp = 0.0;
	tmp2 = 0.0;
	tmp3 = 0.0;
	dgdc = 0.0;
	dgdN = 0.0;
	NconjtmpI = 0.0;
	NconjtmpJ = 0.0;
	Etmp = 0.0;

	REBO_neighs = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs[k];
	if (atomk != atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	lamdajik = 4.0kronecker(itype,1)
	((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dS);
	Nki = nC[atomk]-(wik*kronecker(itype,0))+nH[atomk] -
	(wik*kronecker(itype,1));
	cosjik = ((rij[0]rik[0])+(rij[1]rik[1])+(rij[2]*rik[2])) /
	(rijmag*rikmag);
	cosjik = MIN(cosjik,1.0);
	cosjik = MAX(cosjik,-1.0);

	// evaluate splines g and derivatives dg

	g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN);
	Etmp = Etmp+(wikgexp(lamdajik));
	tmp3 = tmp3+(wikdgdNexp(lamdajik));
	NconjtmpI = NconjtmpI+(kronecker(ktype,0)wikSp(Nki,Nmin,Nmax,dS));
	}
	}

	PijS = 0.0;
	dN2[0] = 0.0;
	dN2[1] = 0.0;
	PijS = PijSpline(NijC,NijH,itype,jtype,dN2);
	pij = 1.0/sqrt(1.0+Etmp+PijS);
	tmp = -0.5pijpij*pij;

	const double invrijm = 1.0/rijmag;
	const double invrijm2 = invrijm*invrijm;

	// pij forces

	REBO_neighs = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs[k];
	if (atomk != atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	lamdajik = 4.0kronecker(itype,1)
	((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);

	const double invrikm = 1.0/rikmag;
	const double invrijkm = invrijm*invrikm;
	const double invrikm2 = invrikm*invrikm;

	cosjik = (rij[0]rik[0] + rij[1]rik[1] + rij[2]rik[2]) invrijkm;
	cosjik = MIN(cosjik,1.0);
	cosjik = MAX(cosjik,-1.0);

	dcosjikdri[0] = ((rij[0]+rik[0])*invrijkm) -
	(cosjik((rij[0]invrijm2)+(rik[0]*invrikm2)));
	dcosjikdri[1] = ((rij[1]+rik[1])*invrijkm) -
	(cosjik((rij[1]invrijm2)+(rik[1]*invrikm2)));
	dcosjikdri[2] = ((rij[2]+rik[2])*invrijkm) -
	(cosjik((rij[2]invrijm2)+(rik[2]*invrikm2)));
	dcosjikdrk[0] = (-rij[0]invrijkm) + (cosjik(rik[0]*invrikm2));
	dcosjikdrk[1] = (-rij[1]invrijkm) + (cosjik(rik[1]*invrikm2));
	dcosjikdrk[2] = (-rij[2]invrijkm) + (cosjik(rik[2]*invrikm2));
	dcosjikdrj[0] = (-rik[0]invrijkm) + (cosjik(rij[0]*invrijm2));
	dcosjikdrj[1] = (-rik[1]invrijkm) + (cosjik(rij[1]*invrijm2));
	dcosjikdrj[2] = (-rik[2]invrijkm) + (cosjik(rij[2]*invrijm2));

	g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN);
	tmp2 = VA.5(tmpwikdgdc*exp(lamdajik));
	fj[0] = -tmp2*dcosjikdrj[0];
	fj[1] = -tmp2*dcosjikdrj[1];
	fj[2] = -tmp2*dcosjikdrj[2];
	fi[0] = -tmp2*dcosjikdri[0];
	fi[1] = -tmp2*dcosjikdri[1];
	fi[2] = -tmp2*dcosjikdri[2];
	fk[0] = -tmp2*dcosjikdrk[0];
	fk[1] = -tmp2*dcosjikdrk[1];
	fk[2] = -tmp2*dcosjikdrk[2];

	tmp2 = VA.5(tmpwikgexp(lamdajik)4.0*kronecker(itype,1));
	fj[0] -= tmp2(-rij[0]invrijm);
	fj[1] -= tmp2(-rij[1]invrijm);
	fj[2] -= tmp2(-rij[2]invrijm);
	fi[0] -= tmp2((-rik[0]invrikm)+(rij[0]*invrijm));
	fi[1] -= tmp2((-rik[1]invrikm)+(rij[1]*invrijm));
	fi[2] -= tmp2((-rik[2]invrikm)+(rij[2]*invrijm));
	fk[0] -= tmp2(rik[0]invrikm);
	fk[1] -= tmp2(rik[1]invrikm);
	fk[2] -= tmp2(rik[2]invrikm);

	// coordination forces

	// dwik forces

	tmp2 = VA.5(tmpdwikgexp(lamdajik))invrikm;
	fi[0] -= tmp2*rik[0];
	fi[1] -= tmp2*rik[1];
	fi[2] -= tmp2*rik[2];
	fk[0] += tmp2*rik[0];
	fk[1] += tmp2*rik[1];
	fk[2] += tmp2*rik[2];

	// PIJ forces

	tmp2 = VA.5(tmpdN2[ktype]dwik)*invrikm;
	fi[0] -= tmp2*rik[0];
	fi[1] -= tmp2*rik[1];
	fi[2] -= tmp2*rik[2];
	fk[0] += tmp2*rik[0];
	fk[1] += tmp2*rik[1];
	fk[2] += tmp2*rik[2];

	// dgdN forces

	tmp2 = VA.5(tmptmp3dwik)*invrikm;
	fi[0] -= tmp2*rik[0];
	fi[1] -= tmp2*rik[1];
	fi[2] -= tmp2*rik[2];
	fk[0] += tmp2*rik[0];
	fk[1] += tmp2*rik[1];
	fk[2] += tmp2*rik[2];

	f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2];
	f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
	f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2];

	if (vflag_atom) {
	rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
	rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2];
	v_tally3_thr(atomi,atomj,atomk,fj,fk,rji,rki,thr);
	}
	}
	}

	tmp = 0.0;
	tmp2 = 0.0;
	tmp3 = 0.0;
	Etmp = 0.0;

	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml != atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	lamdaijl = 4.0kronecker(jtype,1)
	((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dS);
	Nlj = nC[atoml]-(wjl*kronecker(jtype,0)) +
	nH[atoml]-(wjl*kronecker(jtype,1));
	cosijl = -1.0((rij[0]rjl[0])+(rij[1]rjl[1])+(rij[2]rjl[2])) /
	(rijmag*rjlmag);
	cosijl = MIN(cosijl,1.0);
	cosijl = MAX(cosijl,-1.0);

	// evaluate splines g and derivatives dg

	g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN);
	Etmp = Etmp+(wjlgexp(lamdaijl));
	tmp3 = tmp3+(wjldgdNexp(lamdaijl));
	NconjtmpJ = NconjtmpJ+(kronecker(ltype,0)wjlSp(Nlj,Nmin,Nmax,dS));
	}
	}

	PjiS = 0.0;
	dN2[0] = 0.0;
	dN2[1] = 0.0;
	PjiS = PijSpline(NjiC,NjiH,jtype,itype,dN2);
	pji = 1.0/sqrt(1.0+Etmp+PjiS);
	tmp = -0.5pjipji*pji;

	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml != atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	lamdaijl = 4.0kronecker(jtype,1)
	((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);

	const double invrjlm = 1.0/rjlmag;
	const double invrijlm = invrijm*invrjlm;
	const double invrjlm2 = invrjlm*invrjlm;

	cosijl = (-1.0((rij[0]rjl[0])+(rij[1]rjl[1])+(rij[2]rjl[2])))
	* invrijlm;

	cosijl = MIN(cosijl,1.0);
	cosijl = MAX(cosijl,-1.0);

	dcosijldri[0] = (-rjl[0]invrijlm) - (cosijlrij[0]*invrijm2);
	dcosijldri[1] = (-rjl[1]invrijlm) - (cosijlrij[1]*invrijm2);
	dcosijldri[2] = (-rjl[2]invrijlm) - (cosijlrij[2]*invrijm2);
	dcosijldrj[0] = ((-rij[0]+rjl[0])*invrijlm) +
	(cosijl((rij[0]invrijm2)-(rjl[0]*invrjlm2)));
	dcosijldrj[1] = ((-rij[1]+rjl[1])*invrijlm) +
	(cosijl((rij[1]invrijm2)-(rjl[1]*invrjlm2)));
	dcosijldrj[2] = ((-rij[2]+rjl[2])*invrijlm) +
	(cosijl((rij[2]invrijm2)-(rjl[2]*invrjlm2)));
	dcosijldrl[0] = (rij[0]invrijlm)+(cosijlrjl[0]*invrjlm2);
	dcosijldrl[1] = (rij[1]invrijlm)+(cosijlrjl[1]*invrjlm2);
	dcosijldrl[2] = (rij[2]invrijlm)+(cosijlrjl[2]*invrjlm2);

	// evaluate splines g and derivatives dg

	g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN);
	tmp2 = VA.5(tmpwjldgdc*exp(lamdaijl));
	fi[0] = -tmp2*dcosijldri[0];
	fi[1] = -tmp2*dcosijldri[1];
	fi[2] = -tmp2*dcosijldri[2];
	fj[0] = -tmp2*dcosijldrj[0];
	fj[1] = -tmp2*dcosijldrj[1];
	fj[2] = -tmp2*dcosijldrj[2];
	fl[0] = -tmp2*dcosijldrl[0];
	fl[1] = -tmp2*dcosijldrl[1];
	fl[2] = -tmp2*dcosijldrl[2];

	tmp2 = VA.5(tmpwjlgexp(lamdaijl)4.0*kronecker(jtype,1));
	fi[0] -= tmp2(rij[0]invrijm);
	fi[1] -= tmp2(rij[1]invrijm);
	fi[2] -= tmp2(rij[2]invrijm);
	fj[0] -= tmp2((-rjl[0]invrjlm)-(rij[0]*invrijm));
	fj[1] -= tmp2((-rjl[1]invrjlm)-(rij[1]*invrijm));
	fj[2] -= tmp2((-rjl[2]invrjlm)-(rij[2]*invrijm));
	fl[0] -= tmp2(rjl[0]invrjlm);
	fl[1] -= tmp2(rjl[1]invrjlm);
	fl[2] -= tmp2(rjl[2]invrjlm);

	// coordination forces

	// dwik forces

	tmp2 = VA.5(tmpdwjlgexp(lamdaijl))invrjlm;
	fj[0] -= tmp2*rjl[0];
	fj[1] -= tmp2*rjl[1];
	fj[2] -= tmp2*rjl[2];
	fl[0] += tmp2*rjl[0];
	fl[1] += tmp2*rjl[1];
	fl[2] += tmp2*rjl[2];

	// PIJ forces

	tmp2 = VA.5(tmpdN2[ltype]dwjl)*invrjlm;
	fj[0] -= tmp2*rjl[0];
	fj[1] -= tmp2*rjl[1];
	fj[2] -= tmp2*rjl[2];
	fl[0] += tmp2*rjl[0];
	fl[1] += tmp2*rjl[1];
	fl[2] += tmp2*rjl[2];

	// dgdN forces

	tmp2 = VA.5(tmptmp3dwjl)*invrjlm;
	fj[0] -= tmp2*rjl[0];
	fj[1] -= tmp2*rjl[1];
	fj[2] -= tmp2*rjl[2];
	fl[0] += tmp2*rjl[0];
	fl[1] += tmp2*rjl[1];
	fl[2] += tmp2*rjl[2];

	f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2];
	f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
	f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2];

	if (vflag_atom) {
	rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2];
	v_tally3_thr(atomi,atomj,atoml,fi,fl,rij,rlj,thr);
	}
	}
	}

	// evaluate Nij conj

	Nijconj = 1.0+(NconjtmpINconjtmpI)+(NconjtmpJNconjtmpJ);
	piRC = piRCSpline(NijC+NijH,NjiC+NjiH,Nijconj,itype,jtype,dN3);

	// piRC forces

	REBO_neighs_i = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs_i[k];
	if (atomk !=atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
	Nki = nC[atomk]-(wik*kronecker(itype,0))+nH[atomk] -
	(wik*kronecker(itype,1));
	SpN = Sp(Nki,Nmin,Nmax,dNki);

	tmp2 = VAdN3[0]dwik/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	tmp2 = VAdN3[2](2.0NconjtmpIdwik*SpN)/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	if (fabs(dNki) > TOL) {
	REBO_neighs_k = REBO_firstneigh[atomk];
	for (n = 0; n < REBO_numneigh[atomk]; n++) {
	atomn = REBO_neighs_k[n];
	if (atomn != atomi) {
	ntype = map[type[atomn]];
	rkn[0] = x[atomk][0]-x[atomn][0];
	rkn[1] = x[atomk][1]-x[atomn][1];
	rkn[2] = x[atomk][2]-x[atomn][2];
	rknmag = sqrt((rkn[0]rkn[0])+(rkn[1]rkn[1])+(rkn[2]*rkn[2]));
	Sp(rknmag,rcmin[ktype][ntype],rcmax[ktype][ntype],dwkn);

	tmp2 = VAdN3[2](2.0NconjtmpIwikdNkidwkn)/rknmag;
	f[atomk][0] -= tmp2*rkn[0];
	f[atomk][1] -= tmp2*rkn[1];
	f[atomk][2] -= tmp2*rkn[2];
	f[atomn][0] += tmp2*rkn[0];
	f[atomn][1] += tmp2*rkn[1];
	f[atomn][2] += tmp2*rkn[2];

	if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
	}
	}
	}
	}
	}

	// piRC forces

	REBO_neighs = REBO_firstneigh[atomj];
	for (l = 0; l < REBO_numneigh[atomj]; l++) {
	atoml = REBO_neighs[l];
	if (atoml !=atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);
	Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] -
	(wjl*kronecker(jtype,1));
	SpN = Sp(Nlj,Nmin,Nmax,dNlj);

	tmp2 = VAdN3[1]dwjl/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	tmp2 = VAdN3[2](2.0NconjtmpJdwjl*SpN)/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	if (fabs(dNlj) > TOL) {
	REBO_neighs_l = REBO_firstneigh[atoml];
	for (n = 0; n < REBO_numneigh[atoml]; n++) {
	atomn = REBO_neighs_l[n];
	if (atomn != atomj) {
	ntype = map[type[atomn]];
	rln[0] = x[atoml][0]-x[atomn][0];
	rln[1] = x[atoml][1]-x[atomn][1];
	rln[2] = x[atoml][2]-x[atomn][2];
	rlnmag = sqrt((rln[0]rln[0])+(rln[1]rln[1])+(rln[2]*rln[2]));
	Sp(rlnmag,rcmin[ltype][ntype],rcmax[ltype][ntype],dwln);

	tmp2 = VAdN3[2](2.0NconjtmpJwjldNljdwln)/rlnmag;
	f[atoml][0] -= tmp2*rln[0];
	f[atoml][1] -= tmp2*rln[1];
	f[atoml][2] -= tmp2*rln[2];
	f[atomn][0] += tmp2*rln[0];
	f[atomn][1] += tmp2*rln[1];
	f[atomn][2] += tmp2*rln[2];

	if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
	}
	}
	}
	}
	}

	Tij = 0.0;
	dN3[0] = 0.0;
	dN3[1] = 0.0;
	dN3[2] = 0.0;
	if (itype == 0 && jtype == 0)
	Tij=TijSpline((NijC+NijH),(NjiC+NjiH),Nijconj,dN3);
	Etmp = 0.0;

	if (fabs(Tij) > TOL) {
	atom2 = atomi;
	atom3 = atomj;
	r32[0] = x[atom3][0]-x[atom2][0];
	r32[1] = x[atom3][1]-x[atom2][1];
	r32[2] = x[atom3][2]-x[atom2][2];
	r32mag = sqrt((r32[0]r32[0])+(r32[1]r32[1])+(r32[2]*r32[2]));
	r23[0] = -r32[0];
	r23[1] = -r32[1];
	r23[2] = -r32[2];
	r23mag = r32mag;
	REBO_neighs_i = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs_i[k];
	atom1 = atomk;
	ktype = map[type[atomk]];
	if (atomk != atomj) {
	r21[0] = x[atom2][0]-x[atom1][0];
	r21[1] = x[atom2][1]-x[atom1][1];
	r21[2] = x[atom2][2]-x[atom1][2];
	r21mag = sqrt(r21[0]r21[0] + r21[1]r21[1] + r21[2]*r21[2]);
	cos321 = -1.0((r21[0]r32[0])+(r21[1]r32[1])+(r21[2]r32[2])) /
	(r21mag*r32mag);
	cos321 = MIN(cos321,1.0);
	cos321 = MAX(cos321,-1.0);
	Sp2(cos321,thmin,thmax,dcut321);
	sin321 = sqrt(1.0 - cos321*cos321);
	sink2i = 1.0/(sin321*sin321);
	rik2i = 1.0/(r21mag*r21mag);
	if (sin321 != 0.0) {
	rr = (r23magr23mag)-(r21magr21mag);
	rjk[0] = r21[0]-r23[0];
	rjk[1] = r21[1]-r23[1];
	rjk[2] = r21[2]-r23[2];
	rjk2 = (rjk[0]rjk[0])+(rjk[1]rjk[1])+(rjk[2]*rjk[2]);
	rijrik = 2.0r23magr21mag;
	rik2 = r21mag*r21mag;
	dctik = (-rr+rjk2)/(rijrik*rik2);
	dctij = (rr+rjk2)/(rijrikr23magr23mag);
	dctjk = -2.0/rijrik;
	w21 = Sp(r21mag,rcmin[itype][ktype],rcmaxp[itype][ktype],dw21);
	rijmag = r32mag;
	rikmag = r21mag;
	rij2 = r32mag*r32mag;
	rik2 = r21mag*r21mag;
	costmp = 0.5*(rij2+rik2-rjk2)/rijmag/rikmag;
	tspjik = Sp2(costmp,thmin,thmax,dtsjik);
	dtsjik = -dtsjik;

	REBO_neighs_j = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs_j[l];
	atom4 = atoml;
	ltype = map[type[atoml]];
	if (!(atoml == atomi \|\| atoml == atomk)) {
	r34[0] = x[atom3][0]-x[atom4][0];
	r34[1] = x[atom3][1]-x[atom4][1];
	r34[2] = x[atom3][2]-x[atom4][2];
	r34mag = sqrt((r34[0]r34[0])+(r34[1]r34[1])+(r34[2]*r34[2]));
	cos234 = (r32[0]r34[0] + r32[1]r34[1] + r32[2]*r34[2]) /
	(r32mag*r34mag);
	cos234 = MIN(cos234,1.0);
	cos234 = MAX(cos234,-1.0);
	sin234 = sqrt(1.0 - cos234*cos234);
	sinl2i = 1.0/(sin234*sin234);
	rjl2i = 1.0/(r34mag*r34mag);

	if (sin234 != 0.0) {
	w34 = Sp(r34mag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dw34);
	rr = (r23magr23mag)-(r34magr34mag);
	ril[0] = r23[0]+r34[0];
	ril[1] = r23[1]+r34[1];
	ril[2] = r23[2]+r34[2];
	ril2 = (ril[0]ril[0])+(ril[1]ril[1])+(ril[2]*ril[2]);
	rijrjl = 2.0r23magr34mag;
	rjl2 = r34mag*r34mag;
	dctjl = (-rr+ril2)/(rijrjl*rjl2);
	dctji = (rr+ril2)/(rijrjlr23magr23mag);
	dctil = -2.0/rijrjl;
	rjlmag = r34mag;
	rjl2 = r34mag*r34mag;
	costmp = 0.5*(rij2+rjl2-ril2)/rijmag/rjlmag;
	tspijl = Sp2(costmp,thmin,thmax,dtsijl);
	dtsijl = -dtsijl;
	prefactor = VA*Tij;

	cross321[0] = (r32[1]r21[2])-(r32[2]r21[1]);
	cross321[1] = (r32[2]r21[0])-(r32[0]r21[2]);
	cross321[2] = (r32[0]r21[1])-(r32[1]r21[0]);
	cross234[0] = (r23[1]r34[2])-(r23[2]r34[1]);
	cross234[1] = (r23[2]r34[0])-(r23[0]r34[2]);
	cross234[2] = (r23[0]r34[1])-(r23[1]r34[0]);

	cwnum = (cross321[0]*cross234[0]) +
	(cross321[1]cross234[1]) + (cross321[2]cross234[2]);
	cwnom = r21magr34magr23magr23magsin321*sin234;
	om1234 = cwnum/cwnom;
	cw = om1234;
	Etmp += ((1.0-(om1234om1234))w21w34)
	(1.0-tspjik)*(1.0-tspijl);

	dt1dik = (rik2i)-(dctiksink2icos321);
	dt1djk = (-dctjksink2icos321);
	dt1djl = (rjl2i)-(dctjlsinl2icos234);
	dt1dil = (-dctilsinl2icos234);
	dt1dij = (2.0/(r23magr23mag))-(dctijsink2i*cos321) -
	(dctjisinl2icos234);

	dt2dik[0] = (-r23[2]cross234[1])+(r23[1]cross234[2]);
	dt2dik[1] = (-r23[0]cross234[2])+(r23[2]cross234[0]);
	dt2dik[2] = (-r23[1]cross234[0])+(r23[0]cross234[1]);

	dt2djl[0] = (-r23[1]cross321[2])+(r23[2]cross321[1]);
	dt2djl[1] = (-r23[2]cross321[0])+(r23[0]cross321[2]);
	dt2djl[2] = (-r23[0]cross321[1])+(r23[1]cross321[0]);

	dt2dij[0] = (r21[2]cross234[1])-(r34[2]cross321[1]) -
	(r21[1]cross234[2])+(r34[1]cross321[2]);
	dt2dij[1] = (r21[0]cross234[2])-(r34[0]cross321[2]) -
	(r21[2]cross234[0])+(r34[2]cross321[0]);
	dt2dij[2] = (r21[1]cross234[0])-(r34[1]cross321[0]) -
	(r21[0]cross234[1])+(r34[0]cross321[1]);

	aa = (prefactor2.0cw/cwnom)w21w34 *
	(1.0-tspjik)*(1.0-tspijl);
	aaa1 = -prefactor(1.0-(om1234om1234)) *
	(1.0-tspjik)*(1.0-tspijl);
	aaa2 = aaa1w21w34;
	at2 = aa*cwnum;

	fcijpc = (-dt1dijat2)+(aaa2dtsjikdctij(1.0-tspijl)) +
	(aaa2dtsijldctji*(1.0-tspjik));
	fcikpc = (-dt1dikat2)+(aaa2dtsjikdctik(1.0-tspijl));
	fcjlpc = (-dt1djlat2)+(aaa2dtsijldctjl(1.0-tspjik));
	fcjkpc = (-dt1djkat2)+(aaa2dtsjikdctjk(1.0-tspijl));
	fcilpc = (-dt1dilat2)+(aaa2dtsijldctil(1.0-tspjik));

	F23[0] = (fcijpcr23[0])+(aadt2dij[0]);
	F23[1] = (fcijpcr23[1])+(aadt2dij[1]);
	F23[2] = (fcijpcr23[2])+(aadt2dij[2]);

	F12[0] = (fcikpcr21[0])+(aadt2dik[0]);
	F12[1] = (fcikpcr21[1])+(aadt2dik[1]);
	F12[2] = (fcikpcr21[2])+(aadt2dik[2]);

	F34[0] = (fcjlpcr34[0])+(aadt2djl[0]);
	F34[1] = (fcjlpcr34[1])+(aadt2djl[1]);
	F34[2] = (fcjlpcr34[2])+(aadt2djl[2]);

	F31[0] = (fcjkpc*rjk[0]);
	F31[1] = (fcjkpc*rjk[1]);
	F31[2] = (fcjkpc*rjk[2]);

	F24[0] = (fcilpc*ril[0]);
	F24[1] = (fcilpc*ril[1]);
	F24[2] = (fcilpc*ril[2]);

	f1[0] = -F12[0]-F31[0];
	f1[1] = -F12[1]-F31[1];
	f1[2] = -F12[2]-F31[2];
	f2[0] = F23[0]+F12[0]+F24[0];
	f2[1] = F23[1]+F12[1]+F24[1];
	f2[2] = F23[2]+F12[2]+F24[2];
	f3[0] = -F23[0]+F34[0]+F31[0];
	f3[1] = -F23[1]+F34[1]+F31[1];
	f3[2] = -F23[2]+F34[2]+F31[2];
	f4[0] = -F34[0]-F24[0];
	f4[1] = -F34[1]-F24[1];
	f4[2] = -F34[2]-F24[2];

	// coordination forces

	tmp2 = VATij((1.0-(om1234om1234)))
	(1.0-tspjik)(1.0-tspijl)dw21*w34/r21mag;
	f2[0] -= tmp2*r21[0];
	f2[1] -= tmp2*r21[1];
	f2[2] -= tmp2*r21[2];
	f1[0] += tmp2*r21[0];
	f1[1] += tmp2*r21[1];
	f1[2] += tmp2*r21[2];

	tmp2 = VATij((1.0-(om1234om1234)))
	(1.0-tspjik)(1.0-tspijl)w21*dw34/r34mag;
	f3[0] -= tmp2*r34[0];
	f3[1] -= tmp2*r34[1];
	f3[2] -= tmp2*r34[2];
	f4[0] += tmp2*r34[0];
	f4[1] += tmp2*r34[1];
	f4[2] += tmp2*r34[2];

	f[atom1][0] += f1[0]; f[atom1][1] += f1[1];
	f[atom1][2] += f1[2];
	f[atom2][0] += f2[0]; f[atom2][1] += f2[1];
	f[atom2][2] += f2[2];
	f[atom3][0] += f3[0]; f[atom3][1] += f3[1];
	f[atom3][2] += f3[2];
	f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
	f[atom4][2] += f4[2];

	if (vflag_atom) {
	r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2];
	r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2];
	v_tally4_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
	}
	}
	}
	}
	}
	}
	}

	// Tij forces now that we have Etmp

	REBO_neighs = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs[k];
	if (atomk != atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
	Nki = nC[atomk]-(wik*kronecker(itype,0))+nH[atomk] -
	(wik*kronecker(itype,1));
	SpN = Sp(Nki,Nmin,Nmax,dNki);

	tmp2 = VAdN3[0]dwik*Etmp/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	tmp2 = VAdN3[2](2.0NconjtmpIdwikSpN)Etmp/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	if (fabs(dNki) > TOL) {
	REBO_neighs_k = REBO_firstneigh[atomk];
	for (n = 0; n < REBO_numneigh[atomk]; n++) {
	atomn = REBO_neighs_k[n];
	ntype = map[type[atomn]];
	if (atomn != atomi) {
	rkn[0] = x[atomk][0]-x[atomn][0];
	rkn[1] = x[atomk][1]-x[atomn][1];
	rkn[2] = x[atomk][2]-x[atomn][2];
	rknmag = sqrt((rkn[0]rkn[0])+(rkn[1]rkn[1])+(rkn[2]*rkn[2]));
	Sp(rknmag,rcmin[ktype][ntype],rcmax[ktype][ntype],dwkn);

	tmp2 = VAdN3[2](2.0NconjtmpIwikdNkidwkn)*Etmp/rknmag;
	f[atomk][0] -= tmp2*rkn[0];
	f[atomk][1] -= tmp2*rkn[1];
	f[atomk][2] -= tmp2*rkn[2];
	f[atomn][0] += tmp2*rkn[0];
	f[atomn][1] += tmp2*rkn[1];
	f[atomn][2] += tmp2*rkn[2];

	if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
	}
	}
	}
	}
	}

	// Tij forces

	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml != atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);
	Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] -
	(wjl*kronecker(jtype,1));
	SpN = Sp(Nlj,Nmin,Nmax,dNlj);

	tmp2 = VAdN3[1]dwjl*Etmp/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	tmp2 = VAdN3[2](2.0NconjtmpJdwjlSpN)Etmp/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	if (fabs(dNlj) > TOL) {
	REBO_neighs_l = REBO_firstneigh[atoml];
	for (n = 0; n < REBO_numneigh[atoml]; n++) {
	atomn = REBO_neighs_l[n];
	ntype = map[type[atomn]];
	if (atomn !=atomj) {
	rln[0] = x[atoml][0]-x[atomn][0];
	rln[1] = x[atoml][1]-x[atomn][1];
	rln[2] = x[atoml][2]-x[atomn][2];
	rlnmag = sqrt((rln[0]rln[0])+(rln[1]rln[1])+(rln[2]*rln[2]));
	Sp(rlnmag,rcmin[ltype][ntype],rcmax[ltype][ntype],dwln);

	tmp2 = VAdN3[2](2.0NconjtmpJwjldNljdwln)*Etmp/rlnmag;
	f[atoml][0] -= tmp2*rln[0];
	f[atoml][1] -= tmp2*rln[1];
	f[atoml][2] -= tmp2*rln[2];
	f[atomn][0] += tmp2*rln[0];
	f[atomn][1] += tmp2*rln[1];
	f[atomn][2] += tmp2*rln[2];

	if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
	}
	}
	}
	}
	}
	}

	bij = (0.5(pij+pji))+piRC+(TijEtmp);
	return bij;
	}

	/* ----------------------------------------------------------------------
	Bij* function
	------------------------------------------------------------------------- */

	double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag,
	double VA, double rij0[3], double rij0mag,
	int vflag_atom, ThrData * const thr)
	{
	int k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
	int atomi,atomj,itype,jtype,ktype,ltype,ntype;
	double rik[3], rjl[3], rkn[3],rknmag,dNki;
	double NijC,NijH,NjiC,NjiH,wik,dwik,dwkn,wjl;
	double rikmag,rjlmag,cosjik,cosijl,g,tmp2,tmp3;
	double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ;
	double Nki,Nlj,dS,lamdajik,lamdaijl,dgdc,dgdN,pji,Nijconj,piRC;
	double dcosjikdri[3],dcosijldri[3],dcosjikdrk[3];
	double dN2[2],dN3[3];
	double dcosijldrj[3],dcosijldrl[3],dcosjikdrj[3],dwjl;
	double Tij,crosskij[3],crosskijmag;
	double crossijl[3],crossijlmag,omkijl;
	double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3];
	double bij,tmp3pij,tmp3pji,Stb,dStb;
	double r32[3],r32mag,cos321;
	double om1234,rln[3];
	double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag;
	double w21,dw21,r34[3],r34mag,cos234,w34,dw34;
	double cross321[3],cross234[3],prefactor,SpN;
	double fcijpc,fcikpc,fcjlpc,fcjkpc,fcilpc;
	double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa1,aaa2,at2,cw,cwnum,cwnom;
	double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2;
	double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i;
	double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij;
	double dNlj;
	double PijS,PjiS;
	double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp;
	int REBO_neighs,REBO_neighs_i,REBO_neighs_j,REBO_neighs_k,*REBO_neighs_l;
	double F12[3],F23[3],F34[3],F31[3],F24[3];
	double fi[3],fj[3],fk[3],fl[3],f1[3],f2[3],f3[3],f4[4];
	double rji[3],rki[3],rlj[3],r13[3],r43[3];

	const double * const * const x = atom->x;
	double * const * const f = thr->get_f();
	const int * const type = atom->type;

	atomi = i;
	atomj = j;
	itype = map[type[atomi]];
	jtype = map[type[atomj]];
	wij = Sp(rij0mag,rcmin[itype][jtype],rcmax[itype][jtype],dwij);
	NijC = nC[atomi]-(wij*kronecker(jtype,0));
	NijH = nH[atomi]-(wij*kronecker(jtype,1));
	NjiC = nC[atomj]-(wij*kronecker(itype,0));
	NjiH = nH[atomj]-(wij*kronecker(itype,1));

	bij = 0.0;
	tmp = 0.0;
	tmp2 = 0.0;
	tmp3 = 0.0;
	dgdc = 0.0;
	dgdN = 0.0;
	NconjtmpI = 0.0;
	NconjtmpJ = 0.0;
	Etmp = 0.0;

	REBO_neighs = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs[k];
	if (atomk != atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	lamdajik = 4.0kronecker(itype,1)
	((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dS);
	Nki = nC[atomk]-(wik*kronecker(itype,0)) +
	nH[atomk]-(wik*kronecker(itype,1));
	cosjik = ((rij[0]rik[0])+(rij[1]rik[1])+(rij[2]*rik[2])) /
	(rijmag*rikmag);
	cosjik = MIN(cosjik,1.0);
	cosjik = MAX(cosjik,-1.0);

	// evaluate splines g and derivatives dg

	g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN);
	Etmp += (wikgexp(lamdajik));
	tmp3 += (wikdgdNexp(lamdajik));
	NconjtmpI = NconjtmpI+(kronecker(ktype,0)wikSp(Nki,Nmin,Nmax,dS));
	}
	}

	PijS = 0.0;
	dN2PIJ[0] = 0.0;
	dN2PIJ[1] = 0.0;
	PijS = PijSpline(NijC,NijH,itype,jtype,dN2PIJ);
	pij = 1.0/sqrt(1.0+Etmp+PijS);
	tmppij = -.5pijpij*pij;
	tmp3pij = tmp3;
	tmp = 0.0;
	tmp2 = 0.0;
	tmp3 = 0.0;
	Etmp = 0.0;

	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml != atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	lamdaijl = 4.0kronecker(jtype,1)
	((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dS);
	Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] -
	(wjl*kronecker(jtype,1));
	cosijl = -1.0((rij[0]rjl[0])+(rij[1]rjl[1])+(rij[2]rjl[2])) /
	(rijmag*rjlmag);
	cosijl = MIN(cosijl,1.0);
	cosijl = MAX(cosijl,-1.0);

	// evaluate splines g and derivatives dg

	g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN);
	Etmp = Etmp+(wjlgexp(lamdaijl));
	tmp3 = tmp3+(wjldgdNexp(lamdaijl));
	NconjtmpJ = NconjtmpJ+(kronecker(ltype,0)wjlSp(Nlj,Nmin,Nmax,dS));
	}
	}

	PjiS = 0.0;
	dN2PJI[0] = 0.0;
	dN2PJI[1] = 0.0;
	PjiS = PijSpline(NjiC,NjiH,jtype,itype,dN2PJI);
	pji = 1.0/sqrt(1.0+Etmp+PjiS);
	tmppji = -.5pjipji*pji;
	tmp3pji = tmp3;

	// evaluate Nij conj

	Nijconj = 1.0+(NconjtmpINconjtmpI)+(NconjtmpJNconjtmpJ);
	piRC = piRCSpline(NijC+NijH,NjiC+NjiH,Nijconj,itype,jtype,dN3piRC);
	Tij = 0.0;
	dN3Tij[0] = 0.0;
	dN3Tij[1] = 0.0;
	dN3Tij[2] = 0.0;
	if (itype == 0 && jtype == 0)
	Tij=TijSpline((NijC+NijH),(NjiC+NjiH),Nijconj,dN3Tij);

	Etmp = 0.0;
	if (fabs(Tij) > TOL) {
	REBO_neighs_i = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs_i[k];
	ktype = map[type[atomk]];
	if (atomk != atomj) {
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	cos321 = ((rij[0]rik[0])+(rij[1]rik[1])+(rij[2]*rik[2])) /
	(rijmag*rikmag);
	cos321 = MIN(cos321,1.0);
	cos321 = MAX(cos321,-1.0);

	rjk[0] = rik[0]-rij[0];
	rjk[1] = rik[1]-rij[1];
	rjk[2] = rik[2]-rij[2];
	rjk2 = (rjk[0]rjk[0])+(rjk[1]rjk[1])+(rjk[2]*rjk[2]);
	rij2 = rijmag*rijmag;
	rik2 = rikmag*rikmag;
	costmp = 0.5*(rij2+rik2-rjk2)/rijmag/rikmag;
	tspjik = Sp2(costmp,thmin,thmax,dtsjik);

	if (sqrt(1.0 - cos321*cos321) > sqrt(TOL)) {
	wik = Sp(rikmag,rcmin[itype][ktype],rcmaxp[itype][ktype],dwik);
	REBO_neighs_j = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs_j[l];
	ltype = map[type[atoml]];
	if (!(atoml == atomi \|\| atoml == atomk)) {
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt(rjl[0]rjl[0] + rjl[1]rjl[1] + rjl[2]*rjl[2]);
	cos234 = -((rij[0]rjl[0])+(rij[1]rjl[1])+(rij[2]*rjl[2])) /
	(rijmag*rjlmag);
	cos234 = MIN(cos234,1.0);
	cos234 = MAX(cos234,-1.0);

	ril[0] = rij[0]+rjl[0];
	ril[1] = rij[1]+rjl[1];
	ril[2] = rij[2]+rjl[2];
	ril2 = (ril[0]ril[0])+(ril[1]ril[1])+(ril[2]*ril[2]);
	rijrjl = 2.0rijmagrjlmag;
	rjl2 = rjlmag*rjlmag;
	costmp = 0.5*(rij2+rjl2-ril2)/rijmag/rjlmag;
	tspijl = Sp2(costmp,thmin,thmax,dtsijl);

	if (sqrt(1.0 - cos234*cos234) > sqrt(TOL)) {
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dS);
	crosskij[0] = (rij[1]rik[2]-rij[2]rik[1]);
	crosskij[1] = (rij[2]rik[0]-rij[0]rik[2]);
	crosskij[2] = (rij[0]rik[1]-rij[1]rik[0]);
	crosskijmag = sqrt(crosskij[0]*crosskij[0] +
	crosskij[1]*crosskij[1] +
	crosskij[2]*crosskij[2]);
	crossijl[0] = (rij[1]rjl[2]-rij[2]rjl[1]);
	crossijl[1] = (rij[2]rjl[0]-rij[0]rjl[2]);
	crossijl[2] = (rij[0]rjl[1]-rij[1]rjl[0]);
	crossijlmag = sqrt(crossijl[0]*crossijl[0] +
	crossijl[1]*crossijl[1] +
	crossijl[2]*crossijl[2]);
	omkijl = -1.0(((crosskij[0]crossijl[0]) +
	(crosskij[1]*crossijl[1]) +
	(crosskij[2]*crossijl[2])) /
	(crosskijmag*crossijlmag));
	Etmp += ((1.0-(omkijlomkijl))wikwjl)
	(1.0-tspjik)*(1.0-tspijl);
	}
	}
	}
	}
	}
	}
	}

	bij = (.5(pij+pji))+piRC+(TijEtmp);
	Stb = Sp2(bij,bLJmin[itype][jtype],bLJmax[itype][jtype],dStb);
	VA = VA*dStb;

	if (dStb != 0.0) {
	tmp = tmppij;
	dN2[0] = dN2PIJ[0];
	dN2[1] = dN2PIJ[1];
	tmp3 = tmp3pij;

	const double invrijm = 1.0/rijmag;
	const double invrijm2 = invrijm*invrijm;

	// pij forces

	REBO_neighs_i = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs_i[k];
	ktype = map[type[atomk]];
	if (atomk != atomj) {
	lamdajik = 0.0;
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt(rik[0]rik[0] + rik[1]rik[1] + rik[2]*rik[2]);
	lamdajik = 4.0kronecker(itype,1)
	((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);

	const double invrikm = 1.0/rikmag;
	const double invrijkm = invrijm*invrikm;
	const double invrikm2 = invrikm*invrikm;

	cosjik = ((rij[0]rik[0])+(rij[1]rik[1])+(rij[2]*rik[2]))
	* invrijkm;
	cosjik = MIN(cosjik,1.0);
	cosjik = MAX(cosjik,-1.0);

	dcosjikdri[0] = ((rij[0]+rik[0])*invrijkm) -
	(cosjik((rij[0]invrijm2)+(rik[0]*invrikm2)));
	dcosjikdri[1] = ((rij[1]+rik[1])*invrijkm) -
	(cosjik((rij[1]invrijm2)+(rik[1]*invrikm2)));
	dcosjikdri[2] = ((rij[2]+rik[2])*invrijkm) -
	(cosjik((rij[2]invrijm2)+(rik[2]*invrikm2)));
	dcosjikdrk[0] = (-rij[0]invrijkm) + (cosjik(rik[0]*invrikm2));
	dcosjikdrk[1] = (-rij[1]invrijkm) + (cosjik(rik[1]*invrikm2));
	dcosjikdrk[2] = (-rij[2]invrijkm) + (cosjik(rik[2]*invrikm2));
	dcosjikdrj[0] = (-rik[0]invrijkm) + (cosjik(rij[0]*invrijm2));
	dcosjikdrj[1] = (-rik[1]invrijkm) + (cosjik(rij[1]*invrijm2));
	dcosjikdrj[2] = (-rik[2]invrijkm) + (cosjik(rij[2]*invrijm2));

	g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN);

	tmp2 = VA.5(tmpwikdgdc*exp(lamdajik));
	fj[0] = -tmp2*dcosjikdrj[0];
	fj[1] = -tmp2*dcosjikdrj[1];
	fj[2] = -tmp2*dcosjikdrj[2];
	fi[0] = -tmp2*dcosjikdri[0];
	fi[1] = -tmp2*dcosjikdri[1];
	fi[2] = -tmp2*dcosjikdri[2];
	fk[0] = -tmp2*dcosjikdrk[0];
	fk[1] = -tmp2*dcosjikdrk[1];
	fk[2] = -tmp2*dcosjikdrk[2];

	tmp2 = VA.5(tmpwikgexp(lamdajik)4.0*kronecker(itype,1));
	fj[0] -= tmp2(-rij[0]invrijm);
	fj[1] -= tmp2(-rij[1]invrijm);
	fj[2] -= tmp2(-rij[2]invrijm);
	fi[0] -= tmp2((-rik[0]/rikmag)+(rij[0]invrijm));
	fi[1] -= tmp2((-rik[1]/rikmag)+(rij[1]invrijm));
	fi[2] -= tmp2((-rik[2]/rikmag)+(rij[2]invrijm));
	fk[0] -= tmp2*(rik[0]/rikmag);
	fk[1] -= tmp2*(rik[1]/rikmag);
	fk[2] -= tmp2*(rik[2]/rikmag);

	// coordination forces

	// dwik forces

	tmp2 = VA.5(tmpdwikg*exp(lamdajik))/rikmag;
	fi[0] -= tmp2*rik[0];
	fi[1] -= tmp2*rik[1];
	fi[2] -= tmp2*rik[2];
	fk[0] += tmp2*rik[0];
	fk[1] += tmp2*rik[1];
	fk[2] += tmp2*rik[2];

	// PIJ forces

	tmp2 = VA.5(tmpdN2[ktype]dwik)/rikmag;
	fi[0] -= tmp2*rik[0];
	fi[1] -= tmp2*rik[1];
	fi[2] -= tmp2*rik[2];
	fk[0] += tmp2*rik[0];
	fk[1] += tmp2*rik[1];
	fk[2] += tmp2*rik[2];

	// dgdN forces

	tmp2 = VA.5(tmptmp3dwik)/rikmag;
	fi[0] -= tmp2*rik[0];
	fi[1] -= tmp2*rik[1];
	fi[2] -= tmp2*rik[2];
	fk[0] += tmp2*rik[0];
	fk[1] += tmp2*rik[1];
	fk[2] += tmp2*rik[2];

	f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2];
	f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
	f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2];

	if (vflag_atom) {
	rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
	rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2];
	v_tally3_thr(atomi,atomj,atomk,fj,fk,rji,rki,thr);
	}
	}
	}

	tmp = tmppji;
	tmp3 = tmp3pji;
	dN2[0] = dN2PJI[0];
	dN2[1] = dN2PJI[1];
	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml !=atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	lamdaijl = 4.0kronecker(jtype,1)
	((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);

	const double invrjlm = 1.0/rjlmag;
	const double invrijlm = invrijm*invrjlm;
	const double invrjlm2 = invrjlm*invrjlm;

	cosijl = (-1.0((rij[0]rjl[0])+(rij[1]rjl[1])+(rij[2]rjl[2]))) *
	invrijlm;
	cosijl = MIN(cosijl,1.0);
	cosijl = MAX(cosijl,-1.0);

	dcosijldri[0] = (-rjl[0]invrijlm) - (cosijlrij[0]*invrijm2);
	dcosijldri[1] = (-rjl[1]invrijlm) - (cosijlrij[1]*invrijm2);
	dcosijldri[2] = (-rjl[2]invrijlm) - (cosijlrij[2]*invrijm2);
	dcosijldrj[0] = ((-rij[0]+rjl[0])*invrijlm) +
	(cosijl((rij[0]invrijm2)-(rjl[0]*invrjlm2)));
	dcosijldrj[1] = ((-rij[1]+rjl[1])*invrijlm) +
	(cosijl((rij[1]invrijm2)-(rjl[1]*invrjlm2)));
	dcosijldrj[2] = ((-rij[2]+rjl[2])*invrijlm) +
	(cosijl((rij[2]invrijm2)-(rjl[2]*invrjlm2)));
	dcosijldrl[0] = (rij[0]invrijlm) + (cosijlrjl[0]*invrjlm2);
	dcosijldrl[1] = (rij[1]invrijlm) + (cosijlrjl[1]*invrjlm2);
	dcosijldrl[2] = (rij[2]invrijlm) + (cosijlrjl[2]*invrjlm2);

	// evaluate splines g and derivatives dg

	g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN);
	tmp2 = VA.5(tmpwjldgdc*exp(lamdaijl));
	fi[0] = -tmp2*dcosijldri[0];
	fi[1] = -tmp2*dcosijldri[1];
	fi[2] = -tmp2*dcosijldri[2];
	fj[0] = -tmp2*dcosijldrj[0];
	fj[1] = -tmp2*dcosijldrj[1];
	fj[2] = -tmp2*dcosijldrj[2];
	fl[0] = -tmp2*dcosijldrl[0];
	fl[1] = -tmp2*dcosijldrl[1];
	fl[2] = -tmp2*dcosijldrl[2];

	tmp2 = VA.5(tmpwjlgexp(lamdaijl)4.0*kronecker(jtype,1));
	fi[0] -= tmp2(rij[0]invrijm);
	fi[1] -= tmp2(rij[1]invrijm);
	fi[2] -= tmp2(rij[2]invrijm);
	fj[0] -= tmp2((-rjl[0]invrjlm)-(rij[0]*invrijm));
	fj[1] -= tmp2((-rjl[1]invrjlm)-(rij[1]*invrijm));
	fj[2] -= tmp2((-rjl[2]invrjlm)-(rij[2]*invrijm));
	fl[0] -= tmp2(rjl[0]invrjlm);
	fl[1] -= tmp2(rjl[1]invrjlm);
	fl[2] -= tmp2(rjl[2]invrjlm);

	// coordination forces
	// dwik forces

	tmp2 = VA.5(tmpdwjlgexp(lamdaijl))invrjlm;
	fj[0] -= tmp2*rjl[0];
	fj[1] -= tmp2*rjl[1];
	fj[2] -= tmp2*rjl[2];
	fl[0] += tmp2*rjl[0];
	fl[1] += tmp2*rjl[1];
	fl[2] += tmp2*rjl[2];

	// PIJ forces

	tmp2 = VA.5(tmpdN2[ltype]dwjl)*invrjlm;
	fj[0] -= tmp2*rjl[0];
	fj[1] -= tmp2*rjl[1];
	fj[2] -= tmp2*rjl[2];
	fl[0] += tmp2*rjl[0];
	fl[1] += tmp2*rjl[1];
	fl[2] += tmp2*rjl[2];

	// dgdN forces

	tmp2=VA.5(tmptmp3dwjl)*invrjlm;
	fj[0] -= tmp2*rjl[0];
	fj[1] -= tmp2*rjl[1];
	fj[2] -= tmp2*rjl[2];
	fl[0] += tmp2*rjl[0];
	fl[1] += tmp2*rjl[1];
	fl[2] += tmp2*rjl[2];

	f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2];
	f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2];
	f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2];

	if (vflag_atom) {
	rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2];
	v_tally3_thr(atomi,atomj,atoml,fi,fl,rij,rlj,thr);
	}
	}
	}

	// piRC forces

	dN3[0] = dN3piRC[0];
	dN3[1] = dN3piRC[1];
	dN3[2] = dN3piRC[2];

	REBO_neighs_i = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs_i[k];
	if (atomk != atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
	Nki = nC[atomk]-(wik*kronecker(itype,0))+nH[atomk] -
	(wik*kronecker(itype,1));
	SpN = Sp(Nki,Nmin,Nmax,dNki);

	tmp2 = VAdN3[0]dwik/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	tmp2 = VAdN3[2](2.0NconjtmpIdwik*SpN)/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	if (fabs(dNki) > TOL) {
	REBO_neighs_k = REBO_firstneigh[atomk];
	for (n = 0; n < REBO_numneigh[atomk]; n++) {
	atomn = REBO_neighs_k[n];
	if (atomn != atomi) {
	ntype = map[type[atomn]];
	rkn[0] = x[atomk][0]-x[atomn][0];
	rkn[1] = x[atomk][1]-x[atomn][1];
	rkn[2] = x[atomk][2]-x[atomn][2];
	rknmag = sqrt((rkn[0]rkn[0])+(rkn[1]rkn[1])+(rkn[2]*rkn[2]));
	Sp(rknmag,rcmin[ktype][ntype],rcmax[ktype][ntype],dwkn);

	tmp2 = VAdN3[2](2.0NconjtmpIwikdNkidwkn)/rknmag;
	f[atomk][0] -= tmp2*rkn[0];
	f[atomk][1] -= tmp2*rkn[1];
	f[atomk][2] -= tmp2*rkn[2];
	f[atomn][0] += tmp2*rkn[0];
	f[atomn][1] += tmp2*rkn[1];
	f[atomn][2] += tmp2*rkn[2];

	if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
	}
	}
	}
	}
	}

	// piRC forces to J side

	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml != atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);
	Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] -
	(wjl*kronecker(jtype,1));
	SpN = Sp(Nlj,Nmin,Nmax,dNlj);

	tmp2 = VAdN3[1]dwjl/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	tmp2 = VAdN3[2](2.0NconjtmpJdwjl*SpN)/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	if (fabs(dNlj) > TOL) {
	REBO_neighs_l = REBO_firstneigh[atoml];
	for (n = 0; n < REBO_numneigh[atoml]; n++) {
	atomn = REBO_neighs_l[n];
	if (atomn != atomj) {
	ntype = map[type[atomn]];
	rln[0] = x[atoml][0]-x[atomn][0];
	rln[1] = x[atoml][1]-x[atomn][1];
	rln[2] = x[atoml][2]-x[atomn][2];
	rlnmag = sqrt((rln[0]rln[0])+(rln[1]rln[1])+(rln[2]*rln[2]));
	Sp(rlnmag,rcmin[ltype][ntype],rcmax[ltype][ntype],dwln);

	tmp2 = VAdN3[2](2.0NconjtmpJwjldNljdwln)/rlnmag;
	f[atoml][0] -= tmp2*rln[0];
	f[atoml][1] -= tmp2*rln[1];
	f[atoml][2] -= tmp2*rln[2];
	f[atomn][0] += tmp2*rln[0];
	f[atomn][1] += tmp2*rln[1];
	f[atomn][2] += tmp2*rln[2];

	if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
	}
	}
	}
	}
	}

	if (fabs(Tij) > TOL) {
	dN3[0] = dN3Tij[0];
	dN3[1] = dN3Tij[1];
	dN3[2] = dN3Tij[2];
	atom2 = atomi;
	atom3 = atomj;
	r32[0] = x[atom3][0]-x[atom2][0];
	r32[1] = x[atom3][1]-x[atom2][1];
	r32[2] = x[atom3][2]-x[atom2][2];
	r32mag = sqrt((r32[0]r32[0])+(r32[1]r32[1])+(r32[2]*r32[2]));
	r23[0] = -r32[0];
	r23[1] = -r32[1];
	r23[2] = -r32[2];
	r23mag = r32mag;

	REBO_neighs_i = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs_i[k];
	atom1 = atomk;
	ktype = map[type[atomk]];
	if (atomk != atomj) {
	r21[0] = x[atom2][0]-x[atom1][0];
	r21[1] = x[atom2][1]-x[atom1][1];
	r21[2] = x[atom2][2]-x[atom1][2];
	r21mag = sqrt(r21[0]r21[0] + r21[1]r21[1] + r21[2]*r21[2]);
	cos321 = ((r21[0]rij[0])+(r21[1]rij[1])+(r21[2]*rij[2])) /
	(r21mag*rijmag);
	cos321 = MIN(cos321,1.0);
	cos321 = MAX(cos321,-1.0);
	sin321 = sqrt(1.0 - cos321*cos321);
	sink2i = 1.0/(sin321*sin321);
	rik2i = 1.0/(r21mag*r21mag);

	if (sin321 != 0.0) {
	rr = (rijmagrijmag)-(r21magr21mag);
	rjk[0] = r21[0]-rij[0];
	rjk[1] = r21[1]-rij[1];
	rjk[2] = r21[2]-rij[2];
	rjk2 = (rjk[0]rjk[0])+(rjk[1]rjk[1])+(rjk[2]*rjk[2]);
	rijrik = 2.0rijmagr21mag;
	rik2 = r21mag*r21mag;
	dctik = (-rr+rjk2)/(rijrik*rik2);
	dctij = (rr+rjk2)/(rijrikrijmagrijmag);
	dctjk = -2.0/rijrik;
	w21 = Sp(r21mag,rcmin[itype][ktype],rcmaxp[itype][ktype],dw21);
	rikmag = r21mag;
	rij2 = r32mag*r32mag;
	rik2 = r21mag*r21mag;
	costmp = 0.5*(rij2+rik2-rjk2)/rijmag/rikmag;
	tspjik = Sp2(costmp,thmin,thmax,dtsjik);
	dtsjik = -dtsjik;

	REBO_neighs_j = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs_j[l];
	atom4 = atoml;
	ltype = map[type[atoml]];
	if (!(atoml == atomi \|\| atoml == atomk)) {
	r34[0] = x[atom3][0]-x[atom4][0];
	r34[1] = x[atom3][1]-x[atom4][1];
	r34[2] = x[atom3][2]-x[atom4][2];
	r34mag = sqrt(r34[0]r34[0] + r34[1]r34[1] + r34[2]*r34[2]);
	cos234 = -1.0((rij[0]r34[0])+(rij[1]*r34[1]) +
	(rij[2]r34[2]))/(rijmagr34mag);
	cos234 = MIN(cos234,1.0);
	cos234 = MAX(cos234,-1.0);
	sin234 = sqrt(1.0 - cos234*cos234);
	sinl2i = 1.0/(sin234*sin234);
	rjl2i = 1.0/(r34mag*r34mag);

	if (sin234 != 0.0) {
	w34 = Sp(r34mag,rcmin[jtype][ltype],
	rcmaxp[jtype][ltype],dw34);
	rr = (r23magr23mag)-(r34magr34mag);
	ril[0] = r23[0]+r34[0];
	ril[1] = r23[1]+r34[1];
	ril[2] = r23[2]+r34[2];
	ril2 = (ril[0]ril[0])+(ril[1]ril[1])+(ril[2]*ril[2]);
	rijrjl = 2.0r23magr34mag;
	rjl2 = r34mag*r34mag;
	dctjl = (-rr+ril2)/(rijrjl*rjl2);
	dctji = (rr+ril2)/(rijrjlr23magr23mag);
	dctil = -2.0/rijrjl;
	rjlmag = r34mag;
	rjl2 = r34mag*r34mag;
	costmp = 0.5*(rij2+rjl2-ril2)/rijmag/rjlmag;
	tspijl = Sp2(costmp,thmin,thmax,dtsijl);
	dtsijl = -dtsijl; //need minus sign
	prefactor = VA*Tij;

	cross321[0] = (r32[1]r21[2])-(r32[2]r21[1]);
	cross321[1] = (r32[2]r21[0])-(r32[0]r21[2]);
	cross321[2] = (r32[0]r21[1])-(r32[1]r21[0]);
	cross234[0] = (r23[1]r34[2])-(r23[2]r34[1]);
	cross234[1] = (r23[2]r34[0])-(r23[0]r34[2]);
	cross234[2] = (r23[0]r34[1])-(r23[1]r34[0]);

	cwnum = (cross321[0]*cross234[0]) +
	(cross321[1]cross234[1])+(cross321[2]cross234[2]);
	cwnom = r21magr34magr23magr23magsin321*sin234;
	om1234 = cwnum/cwnom;
	cw = om1234;
	Etmp += ((1.0-(om1234om1234))w21w34)
	(1.0-tspjik)*(1.0-tspijl);

	dt1dik = (rik2i)-(dctiksink2icos321);
	dt1djk = (-dctjksink2icos321);
	dt1djl = (rjl2i)-(dctjlsinl2icos234);
	dt1dil = (-dctilsinl2icos234);
	dt1dij = (2.0/(r23mag*r23mag)) -
	(dctijsink2icos321)-(dctjisinl2icos234);

	dt2dik[0] = (-r23[2]cross234[1])+(r23[1]cross234[2]);
	dt2dik[1] = (-r23[0]cross234[2])+(r23[2]cross234[0]);
	dt2dik[2] = (-r23[1]cross234[0])+(r23[0]cross234[1]);

	dt2djl[0] = (-r23[1]cross321[2])+(r23[2]cross321[1]);
	dt2djl[1] = (-r23[2]cross321[0])+(r23[0]cross321[2]);
	dt2djl[2] = (-r23[0]cross321[1])+(r23[1]cross321[0]);

	dt2dij[0] = (r21[2]*cross234[1]) -
	(r34[2]cross321[1])-(r21[1]cross234[2]) +
	(r34[1]*cross321[2]);
	dt2dij[1] = (r21[0]*cross234[2]) -
	(r34[0]cross321[2])-(r21[2]cross234[0]) +
	(r34[2]*cross321[0]);
	dt2dij[2] = (r21[1]*cross234[0]) -
	(r34[1]cross321[0])-(r21[0]cross234[1]) +
	(r34[0]*cross321[1]);

	aa = (prefactor2.0cw/cwnom)w21w34 *
	(1.0-tspjik)*(1.0-tspijl);
	aaa1 = -prefactor(1.0-(om1234om1234)) *
	(1.0-tspjik)*(1.0-tspijl);
	aaa2 = aaa1w21w34;
	at2 = aa*cwnum;

	fcijpc = (-dt1dijat2)+(aaa2dtsjikdctij(1.0-tspijl)) +
	(aaa2dtsijldctji*(1.0-tspjik));
	fcikpc = (-dt1dikat2)+(aaa2dtsjikdctik(1.0-tspijl));
	fcjlpc = (-dt1djlat2)+(aaa2dtsijldctjl(1.0-tspjik));
	fcjkpc = (-dt1djkat2)+(aaa2dtsjikdctjk(1.0-tspijl));
	fcilpc = (-dt1dilat2)+(aaa2dtsijldctil(1.0-tspjik));

	F23[0] = (fcijpcr23[0])+(aadt2dij[0]);
	F23[1] = (fcijpcr23[1])+(aadt2dij[1]);
	F23[2] = (fcijpcr23[2])+(aadt2dij[2]);

	F12[0] = (fcikpcr21[0])+(aadt2dik[0]);
	F12[1] = (fcikpcr21[1])+(aadt2dik[1]);
	F12[2] = (fcikpcr21[2])+(aadt2dik[2]);

	F34[0] = (fcjlpcr34[0])+(aadt2djl[0]);
	F34[1] = (fcjlpcr34[1])+(aadt2djl[1]);
	F34[2] = (fcjlpcr34[2])+(aadt2djl[2]);

	F31[0] = (fcjkpc*rjk[0]);
	F31[1] = (fcjkpc*rjk[1]);
	F31[2] = (fcjkpc*rjk[2]);

	F24[0] = (fcilpc*ril[0]);
	F24[1] = (fcilpc*ril[1]);
	F24[2] = (fcilpc*ril[2]);

	f1[0] = -F12[0]-F31[0];
	f1[1] = -F12[1]-F31[1];
	f1[2] = -F12[2]-F31[2];
	f2[0] = F23[0]+F12[0]+F24[0];
	f2[1] = F23[1]+F12[1]+F24[1];
	f2[2] = F23[2]+F12[2]+F24[2];
	f3[0] = -F23[0]+F34[0]+F31[0];
	f3[1] = -F23[1]+F34[1]+F31[1];
	f3[2] = -F23[2]+F34[2]+F31[2];
	f4[0] = -F34[0]-F24[0];
	f4[1] = -F34[1]-F24[1];
	f4[2] = -F34[2]-F24[2];

	// coordination forces

	tmp2 = VATij((1.0-(om1234om1234)))
	(1.0-tspjik)(1.0-tspijl)dw21*w34/r21mag;
	f2[0] -= tmp2*r21[0];
	f2[1] -= tmp2*r21[1];
	f2[2] -= tmp2*r21[2];
	f1[0] += tmp2*r21[0];
	f1[1] += tmp2*r21[1];
	f1[2] += tmp2*r21[2];

	tmp2 = VATij((1.0-(om1234om1234)))
	(1.0-tspjik)(1.0-tspijl)w21*dw34/r34mag;
	f3[0] -= tmp2*r34[0];
	f3[1] -= tmp2*r34[1];
	f3[2] -= tmp2*r34[2];
	f4[0] += tmp2*r34[0];
	f4[1] += tmp2*r34[1];
	f4[2] += tmp2*r34[2];

	f[atom1][0] += f1[0]; f[atom1][1] += f1[1];
	f[atom1][2] += f1[2];
	f[atom2][0] += f2[0]; f[atom2][1] += f2[1];
	f[atom2][2] += f2[2];
	f[atom3][0] += f3[0]; f[atom3][1] += f3[1];
	f[atom3][2] += f3[2];
	f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
	f[atom4][2] += f4[2];

	if (vflag_atom) {
	r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2];
	r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2];
	v_tally4_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
	}
	}
	}
	}
	}
	}
	}

	REBO_neighs = REBO_firstneigh[i];
	for (k = 0; k < REBO_numneigh[i]; k++) {
	atomk = REBO_neighs[k];
	if (atomk != atomj) {
	ktype = map[type[atomk]];
	rik[0] = x[atomi][0]-x[atomk][0];
	rik[1] = x[atomi][1]-x[atomk][1];
	rik[2] = x[atomi][2]-x[atomk][2];
	rikmag = sqrt((rik[0]rik[0])+(rik[1]rik[1])+(rik[2]*rik[2]));
	wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
	Nki = nC[atomk]-(wik*kronecker(itype,0))+nH[atomk] -
	(wik*kronecker(itype,1));
	SpN = Sp(Nki,Nmin,Nmax,dNki);

	tmp2 = VAdN3[0]dwik*Etmp/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	tmp2 = VAdN3[2](2.0NconjtmpIdwikSpN)Etmp/rikmag;
	f[atomi][0] -= tmp2*rik[0];
	f[atomi][1] -= tmp2*rik[1];
	f[atomi][2] -= tmp2*rik[2];
	f[atomk][0] += tmp2*rik[0];
	f[atomk][1] += tmp2*rik[1];
	f[atomk][2] += tmp2*rik[2];

	if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);

	if (fabs(dNki) >TOL) {
	REBO_neighs_k = REBO_firstneigh[atomk];
	for (n = 0; n < REBO_numneigh[atomk]; n++) {
	atomn = REBO_neighs_k[n];
	ntype = map[type[atomn]];
	if (atomn !=atomi) {
	rkn[0] = x[atomk][0]-x[atomn][0];
	rkn[1] = x[atomk][1]-x[atomn][1];
	rkn[2] = x[atomk][2]-x[atomn][2];
	rknmag = sqrt((rkn[0]rkn[0])+(rkn[1]rkn[1])+(rkn[2]*rkn[2]));
	Sp(rknmag,rcmin[ktype][ntype],rcmax[ktype][ntype],dwkn);

	tmp2 = VAdN3[2](2.0NconjtmpIwikdNkidwkn)*Etmp/rknmag;
	f[atomk][0] -= tmp2*rkn[0];
	f[atomk][1] -= tmp2*rkn[1];
	f[atomk][2] -= tmp2*rkn[2];
	f[atomn][0] += tmp2*rkn[0];
	f[atomn][1] += tmp2*rkn[1];
	f[atomn][2] += tmp2*rkn[2];

	if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
	}
	}
	}
	}
	}

	// Tij forces

	REBO_neighs = REBO_firstneigh[j];
	for (l = 0; l < REBO_numneigh[j]; l++) {
	atoml = REBO_neighs[l];
	if (atoml != atomi) {
	ltype = map[type[atoml]];
	rjl[0] = x[atomj][0]-x[atoml][0];
	rjl[1] = x[atomj][1]-x[atoml][1];
	rjl[2] = x[atomj][2]-x[atoml][2];
	rjlmag = sqrt((rjl[0]rjl[0])+(rjl[1]rjl[1])+(rjl[2]*rjl[2]));
	wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);
	Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] -
	(wjl*kronecker(jtype,1));
	SpN = Sp(Nlj,Nmin,Nmax,dNlj);

	tmp2 = VAdN3[1]dwjl*Etmp/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	tmp2 = VAdN3[2](2.0NconjtmpJdwjlSpN)Etmp/rjlmag;
	f[atomj][0] -= tmp2*rjl[0];
	f[atomj][1] -= tmp2*rjl[1];
	f[atomj][2] -= tmp2*rjl[2];
	f[atoml][0] += tmp2*rjl[0];
	f[atoml][1] += tmp2*rjl[1];
	f[atoml][2] += tmp2*rjl[2];

	if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);

	if (fabs(dNlj) > TOL) {
	REBO_neighs_l = REBO_firstneigh[atoml];
	for (n = 0; n < REBO_numneigh[atoml]; n++) {
	atomn = REBO_neighs_l[n];
	ntype = map[type[atomn]];
	if (atomn != atomj) {
	rln[0] = x[atoml][0]-x[atomn][0];
	rln[1] = x[atoml][1]-x[atomn][1];
	rln[2] = x[atoml][2]-x[atomn][2];
	rlnmag = sqrt((rln[0]rln[0])+(rln[1]rln[1])+(rln[2]*rln[2]));
	Sp(rlnmag,rcmin[ltype][ntype],rcmax[ltype][ntype],dwln);

	tmp2 = VAdN3[2](2.0NconjtmpJwjldNljdwln)*Etmp/rlnmag;
	f[atoml][0] -= tmp2*rln[0];
	f[atoml][1] -= tmp2*rln[1];
	f[atoml][2] -= tmp2*rln[2];
	f[atomn][0] += tmp2*rln[0];
	f[atomn][1] += tmp2*rln[1];
	f[atomn][2] += tmp2*rln[2];

	if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
	}
	}
	}
	}
	}
	}
	}

	return Stb;
	}

	/* ----------------------------------------------------------------------
	REBO forces and energy
	------------------------------------------------------------------------- */

	void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
	int vflag_atom, ThrData * const thr)
	{
	int i,j,k,m,ii,itype,jtype;
	tagint itag,jtag;
	double delx,dely,delz,evdwl,fpair,xtmp,ytmp,ztmp;
	double rsq,rij,wij;
	double Qij,Aij,alphaij,VR,pre,dVRdi,VA,term,bij,dVAdi,dVA;
	double dwij,del[3];
	int ilist,REBO_neighs;

	evdwl = 0.0;

	const double * const * const x = atom->x;
	double * const * const f = thr->get_f();
	int *type = atom->type;
	tagint *tag = atom->tag;
	int nlocal = atom->nlocal;

	ilist = list->ilist;

	// two-body interactions from REBO neighbor list, skip half of them

	for (ii = ifrom; ii < ito; ii++) {
	i = ilist[ii];
	itag = tag[i];
	itype = map[type[i]];
	xtmp = x[i][0];
	ytmp = x[i][1];
	ztmp = x[i][2];
	REBO_neighs = REBO_firstneigh[i];

	for (k = 0; k < REBO_numneigh[i]; k++) {
	j = REBO_neighs[k];
	jtag = tag[j];

	if (itag > jtag) {
	if ((itag+jtag) % 2 == 0) continue;
	} else if (itag < jtag) {
	if ((itag+jtag) % 2 == 1) continue;
	} else {
	if (x[j][2] < ztmp) continue;
	if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
	if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
	}

	jtype = map[type[j]];

	delx = x[i][0] - x[j][0];
	dely = x[i][1] - x[j][1];
	delz = x[i][2] - x[j][2];
	rsq = delxdelx + delydely + delz*delz;
	rij = sqrt(rsq);
	wij = Sp(rij,rcmin[itype][jtype],rcmax[itype][jtype],dwij);
	if (wij <= TOL) continue;

	Qij = Q[itype][jtype];
	Aij = A[itype][jtype];
	alphaij = alpha[itype][jtype];

	VR = wij(1.0+(Qij/rij)) Aijexp(-alphaijrij);
	pre = wijAij exp(-alphaij*rij);
	dVRdi = pre * ((-alphaij)-(Qij/rsq)-(Qij*alphaij/rij));
	dVRdi += VR/wij * dwij;

	VA = dVA = 0.0;
	for (m = 0; m < 3; m++) {
	term = -wij * BIJc[itype][jtype][m] * exp(-Beta[itype][jtype][m]*rij);
	VA += term;
	dVA += -Beta[itype][jtype][m] * term;
	}
	dVA += VA/wij * dwij;
	del[0] = delx;
	del[1] = dely;
	del[2] = delz;
	bij = bondorder_thr(i,j,del,rij,VA,vflag_atom,thr);
	dVAdi = bij*dVA;

	fpair = -(dVRdi+dVAdi) / rij;
	f[i][0] += delx*fpair;
	f[i][1] += dely*fpair;
	f[i][2] += delz*fpair;
	f[j][0] -= delx*fpair;
	f[j][1] -= dely*fpair;
	f[j][2] -= delz*fpair;

	if (eflag) evdwl = VR + bij*VA;
	if (evflag) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
	evdwl,0.0,fpair,delx,dely,delz,thr);
	}
	}
	}

	/* ----------------------------------------------------------------------
	compute LJ forces and energy
	find 3- and 4-step paths between atoms I,J via REBO neighbor lists
	------------------------------------------------------------------------- */

	void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
	int vflag_atom, ThrData * const thr)
	{
	int i,j,k,m,ii,jj,kk,mm,jnum,itype,jtype,ktype,mtype;
	tagint itag,jtag;
	int atomi,atomj,atomk,atomm;
	int testpath,npath,done;
	double evdwl,fpair,xtmp,ytmp,ztmp;
	double rsq,best,wik,wkm,cij,rij,dwij,dwik,dwkj,dwkm,dwmj;
	double delij[3],rijsq,delik[3],rik,deljk[3];
	double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb;
	double vdw,slw,dvdw,dslw,drij,swidth,tee,tee2;
	double rljmin,rljmax,sigcut,sigmin,sigwid;
	double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3];
	int ilist,jlist,numneigh,*firstneigh;
	int REBO_neighs_i,REBO_neighs_k;
	double delikS[3],deljkS[3],delkmS[3],deljmS[3],delimS[3];
	double rikS,rkjS,rkmS,rmjS,wikS,dwikS;
	double wkjS,dwkjS,wkmS,dwkmS,wmjS,dwmjS;
	double fpair1,fpair2,fpair3;
	double fi[3],fj[3],fk[3],fm[3];

	// I-J interaction from full neighbor list
	// skip 1/2 of interactions since only consider each pair once

	evdwl = 0.0;
	rljmin = 0.0;
	rljmax = 0.0;
	sigcut = 0.0;
	sigmin = 0.0;
	sigwid = 0.0;

	const double * const * const x = atom->x;
	double * const * const f = thr->get_f();
	tagint *tag = atom->tag;
	int *type = atom->type;
	int nlocal = atom->nlocal;

	ilist = list->ilist;
	numneigh = list->numneigh;
	firstneigh = list->firstneigh;

	// loop over neighbors of my atoms

	for (ii = ifrom; ii < ito; ii++) {
	i = ilist[ii];
	itag = tag[i];
	itype = map[type[i]];
	atomi = i;
	xtmp = x[i][0];
	ytmp = x[i][1];
	ztmp = x[i][2];
	jlist = firstneigh[i];
	jnum = numneigh[i];

	for (jj = 0; jj < jnum; jj++) {
	j = jlist[jj];
	j &= NEIGHMASK;
	jtag = tag[j];

	if (itag > jtag) {
	if ((itag+jtag) % 2 == 0) continue;
	} else if (itag < jtag) {
	if ((itag+jtag) % 2 == 1) continue;
	} else {
	if (x[j][2] < ztmp) continue;
	if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
	if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
	}

	jtype = map[type[j]];
	atomj = j;

	delij[0] = xtmp - x[j][0];
	delij[1] = ytmp - x[j][1];
	delij[2] = ztmp - x[j][2];
	rijsq = delij[0]delij[0] + delij[1]delij[1] + delij[2]*delij[2];

	// if outside of LJ cutoff, skip
	// if outside of 4-path cutoff, best = 0.0, no need to test paths
	// if outside of 2-path cutoff but inside 4-path cutoff,
	// best = 0.0, test 3-,4-paths
	// if inside 2-path cutoff, best = wij, only test 3-,4-paths if best < 1

	if (rijsq >= cutljsq[itype][jtype]) continue;
	rij = sqrt(rijsq);
	if (rij >= cut3rebo) {
	best = 0.0;
	testpath = 0;
	} else if (rij >= rcmax[itype][jtype]) {
	best = 0.0;
	testpath = 1;
	} else {
	best = Sp(rij,rcmin[itype][jtype],rcmax[itype][jtype],dwij);
	npath = 2;
	if (best < 1.0) testpath = 1;
	else testpath = 0;
	}

	done = 0;
	if (testpath) {

	// test all 3-body paths = I-K-J
	// I-K interactions come from atom I's REBO neighbors
	// if wik > current best, compute wkj
	// if best = 1.0, done

	REBO_neighs_i = REBO_firstneigh[i];
	for (kk = 0; kk < REBO_numneigh[i] && done==0; kk++) {
	k = REBO_neighs_i[kk];
	if (k == j) continue;
	ktype = map[type[k]];

	delik[0] = x[i][0] - x[k][0];
	delik[1] = x[i][1] - x[k][1];
	delik[2] = x[i][2] - x[k][2];
	rsq = delik[0]delik[0] + delik[1]delik[1] + delik[2]*delik[2];
	if (rsq < rcmaxsq[itype][ktype]) {
	rik = sqrt(rsq);
	wik = Sp(rik,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
	} else wik = 0.0;

	if (wik > best) {
	deljk[0] = x[j][0] - x[k][0];
	deljk[1] = x[j][1] - x[k][1];
	deljk[2] = x[j][2] - x[k][2];
	rsq = deljk[0]deljk[0] + deljk[1]deljk[1] + deljk[2]*deljk[2];
	if (rsq < rcmaxsq[ktype][jtype]) {
	rkj = sqrt(rsq);
	wkj = Sp(rkj,rcmin[ktype][jtype],rcmax[ktype][jtype],dwkj);
	if (wik*wkj > best) {
	best = wik*wkj;
	npath = 3;
	atomk = k;
	delikS[0] = delik[0];
	delikS[1] = delik[1];
	delikS[2] = delik[2];
	rikS = rik;
	wikS = wik;
	dwikS = dwik;
	deljkS[0] = deljk[0];
	deljkS[1] = deljk[1];
	deljkS[2] = deljk[2];
	rkjS = rkj;
	wkjS = wkj;
	dwkjS = dwkj;
	if (best == 1.0) {
	done = 1;
	break;
	}
	}
	}

	// test all 4-body paths = I-K-M-J
	// K-M interactions come from atom K's REBO neighbors
	// if wik*wkm > current best, compute wmj
	// if best = 1.0, done

	REBO_neighs_k = REBO_firstneigh[k];
	for (mm = 0; mm < REBO_numneigh[k] && done==0; mm++) {
	m = REBO_neighs_k[mm];
	if (m == i \|\| m == j) continue;
	mtype = map[type[m]];
	delkm[0] = x[k][0] - x[m][0];
	delkm[1] = x[k][1] - x[m][1];
	delkm[2] = x[k][2] - x[m][2];
	rsq = delkm[0]delkm[0] + delkm[1]delkm[1] + delkm[2]*delkm[2];
	if (rsq < rcmaxsq[ktype][mtype]) {
	rkm = sqrt(rsq);
	wkm = Sp(rkm,rcmin[ktype][mtype],rcmax[ktype][mtype],dwkm);
	} else wkm = 0.0;

	if (wik*wkm > best) {
	deljm[0] = x[j][0] - x[m][0];
	deljm[1] = x[j][1] - x[m][1];
	deljm[2] = x[j][2] - x[m][2];
	rsq = deljm[0]deljm[0] + deljm[1]deljm[1] +
	deljm[2]*deljm[2];
	if (rsq < rcmaxsq[mtype][jtype]) {
	rmj = sqrt(rsq);
	wmj = Sp(rmj,rcmin[mtype][jtype],rcmax[mtype][jtype],dwmj);
	if (wikwkmwmj > best) {
	best = wikwkmwmj;
	npath = 4;
	atomk = k;
	delikS[0] = delik[0];
	delikS[1] = delik[1];
	delikS[2] = delik[2];
	rikS = rik;
	wikS = wik;
	dwikS = dwik;
	atomm = m;
	delkmS[0] = delkm[0];
	delkmS[1] = delkm[1];
	delkmS[2] = delkm[2];
	rkmS = rkm;
	wkmS = wkm;
	dwkmS = dwkm;
	deljmS[0] = deljm[0];
	deljmS[1] = deljm[1];
	deljmS[2] = deljm[2];
	rmjS = rmj;
	wmjS = wmj;
	dwmjS = dwmj;
	if (best == 1.0) {
	done = 1;
	break;
	}
	}
	}
	}
	}
	}
	}
	}

	cij = 1.0 - best;
	if (cij == 0.0) continue;

	// compute LJ forces and energy

	sigwid = 0.84;
	sigcut = 3.0;
	sigmin = sigcut - sigwid;

	rljmin = sigma[itype][jtype];
	rljmax = sigcut * rljmin;
	rljmin = sigmin * rljmin;

	if (rij > rljmax) {
	slw = 0.0;
	dslw = 0.0;
	} else if (rij > rljmin) {
	drij = rij - rljmin;
	swidth = rljmax - rljmin;
	tee = drij / swidth;
	tee2 = tee*tee;
	slw = 1.0 - tee2 * (3.0 - 2.0 * tee);
	dslw = 6.0 * tee * (1.0 - tee) / rij / swidth;
	} else {
	slw = 1.0;
	dslw = 0.0;
	}

	if (morseflag) {

	const double exr = exp(-rij*lj4[itype][jtype]);
	vdw = lj1[itype][jtype]exr(lj2[itype][jtype]*exr - 2);
	dvdw = lj3[itype][jtype]exr(1-lj2[itype][jtype]*exr);

	} else {

	r2inv = 1.0/rijsq;
	r6inv = r2invr2invr2inv;

	vdw = r6inv(lj3[itype][jtype]r6inv-lj4[itype][jtype]);
	dvdw = -r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) / rij;
	}

	// VLJ now becomes vdw * slw, derivaties, etc.

	VLJ = vdw * slw;
	dVLJ = dvdw * slw + vdw * dslw;

	Str = Sp2(rij,rcLJmin[itype][jtype],rcLJmax[itype][jtype],dStr);
	VA = StrcijVLJ;
	if (Str > 0.0) {
	scale = rcmin[itype][jtype] / rij;
	delscale[0] = scale * delij[0];
	delscale[1] = scale * delij[1];
	delscale[2] = scale * delij[2];
	Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA,
	delij,rij,vflag_atom,thr);
	} else Stb = 0.0;

	fpair = -(dStr * (StbcijVLJ - cij*VLJ) +
	dVLJ * (StrStbcij + cij - Str*cij)) / rij;

	f[i][0] += delij[0]*fpair;
	f[i][1] += delij[1]*fpair;
	f[i][2] += delij[2]*fpair;
	f[j][0] -= delij[0]*fpair;
	f[j][1] -= delij[1]*fpair;
	f[j][2] -= delij[2]*fpair;

	if (eflag) evdwl = VAStb + (1.0-Str)cij*VLJ;
	if (evflag) ev_tally_thr(this,i,j,nlocal,/* newton_pair */ 1,
	evdwl,0.0,fpair,delij[0],delij[1],delij[2],thr);

	if (cij < 1.0) {
	dC = StrStbVLJ + (1.0-Str)*VLJ;
	if (npath == 2) {
	fpair = dC*dwij / rij;
	f[atomi][0] += delij[0]*fpair;
	f[atomi][1] += delij[1]*fpair;
	f[atomi][2] += delij[2]*fpair;
	f[atomj][0] -= delij[0]*fpair;
	f[atomj][1] -= delij[1]*fpair;
	f[atomj][2] -= delij[2]*fpair;

	if (vflag_atom) v_tally2_thr(atomi,atomj,fpair,delij,thr);

	} else if (npath == 3) {
	fpair1 = dCdwikSwkjS / rikS;
	fi[0] = delikS[0]*fpair1;
	fi[1] = delikS[1]*fpair1;
	fi[2] = delikS[2]*fpair1;
	fpair2 = dCwikSdwkjS / rkjS;
	fj[0] = deljkS[0]*fpair2;
	fj[1] = deljkS[1]*fpair2;
	fj[2] = deljkS[2]*fpair2;

	f[atomi][0] += fi[0];
	f[atomi][1] += fi[1];
	f[atomi][2] += fi[2];
	f[atomj][0] += fj[0];
	f[atomj][1] += fj[1];
	f[atomj][2] += fj[2];
	f[atomk][0] -= fi[0] + fj[0];
	f[atomk][1] -= fi[1] + fj[1];
	f[atomk][2] -= fi[2] + fj[2];

	if (vflag_atom)
	v_tally3_thr(atomi,atomj,atomk,fi,fj,delikS,deljkS,thr);

	} else if (npath == 4) {
	fpair1 = dCdwikSwkmS*wmjS / rikS;
	fi[0] = delikS[0]*fpair1;
	fi[1] = delikS[1]*fpair1;
	fi[2] = delikS[2]*fpair1;

	fpair2 = dCwikSdwkmS*wmjS / rkmS;
	fk[0] = delkmS[0]*fpair2 - fi[0];
	fk[1] = delkmS[1]*fpair2 - fi[1];
	fk[2] = delkmS[2]*fpair2 - fi[2];

	fpair3 = dCwikSwkmS*dwmjS / rmjS;
	fj[0] = deljmS[0]*fpair3;
	fj[1] = deljmS[1]*fpair3;
	fj[2] = deljmS[2]*fpair3;

	fm[0] = -delkmS[0]*fpair2 - fj[0];
	fm[1] = -delkmS[1]*fpair2 - fj[1];
	fm[2] = -delkmS[2]*fpair2 - fj[2];

	f[atomi][0] += fi[0];
	f[atomi][1] += fi[1];
	f[atomi][2] += fi[2];
	f[atomj][0] += fj[0];
	f[atomj][1] += fj[1];
	f[atomj][2] += fj[2];
	f[atomk][0] += fk[0];
	f[atomk][1] += fk[1];
	f[atomk][2] += fk[2];
	f[atomm][0] += fm[0];
	f[atomm][1] += fm[1];
	f[atomm][2] += fm[2];

	if (vflag_atom) {
	delimS[0] = delikS[0] + delkmS[0];
	delimS[1] = delikS[1] + delkmS[1];
	delimS[2] = delikS[2] + delkmS[2];
	v_tally4_thr(atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr);
	}
	}
	}
	}
	}
	}

	/* ----------------------------------------------------------------------
	torsional forces and energy
	------------------------------------------------------------------------- */

	void PairAIREBOOMP::TORSION_thr(int ifrom, int ito,
	int evflag, int eflag, ThrData * const thr)
	{
	int i,j,k,l,ii;
	tagint itag,jtag;
	double evdwl,fpair,xtmp,ytmp,ztmp;
	double cos321;
	double w21,dw21,cos234,w34,dw34;
	double cross321[3],cross321mag,cross234[3],cross234mag;
	double w23,dw23,cw2,ekijl,Ec;
	double cw,cwnum,cwnom;
	double rij,rij2,rik,rjl,tspjik,dtsjik,tspijl,dtsijl,costmp,fcpc;
	double sin321,sin234,rjk2,rik2,ril2,rjl2;
	double rjk,ril;
	double Vtors;
	double dndij[3],tmpvec[3],dndik[3],dndjl[3];
	double dcidij,dcidik,dcidjk,dcjdji,dcjdjl,dcjdil;
	double dsidij,dsidik,dsidjk,dsjdji,dsjdjl,dsjdil;
	double dxidij,dxidik,dxidjk,dxjdji,dxjdjl,dxjdil;
	double ddndij,ddndik,ddndjk,ddndjl,ddndil,dcwddn,dcwdn,dvpdcw,Ftmp[3];
	double del32[3],rsq,r32,del23[3],del21[3],r21;
	double deljk[3],del34[3],delil[3],delkl[3],r23,r34;
	double fi[3],fj[3],fk[3],fl[3];
	int itype,jtype,ktype,ltype,kk,ll,jj;
	int ilist,REBO_neighs_i,*REBO_neighs_j;

	const double * const * const x = atom->x;
	double * const * const f = thr->get_f();
	int *type = atom->type;
	tagint *tag = atom->tag;

	ilist = list->ilist;

	for (ii = ifrom; ii < ito; ii++) {
	i = ilist[ii];
	itag = tag[i];
	itype = map[type[i]];
	if (itype != 0) continue;
	xtmp = x[i][0];
	ytmp = x[i][1];
	ztmp = x[i][2];
	REBO_neighs_i = REBO_firstneigh[i];

	for (jj = 0; jj < REBO_numneigh[i]; jj++) {
	j = REBO_neighs_i[jj];
	jtag = tag[j];

	if (itag > jtag) {
	if ((itag+jtag) % 2 == 0) continue;
	} else if (itag < jtag) {
	if ((itag+jtag) % 2 == 1) continue;
	} else {
	if (x[j][2] < ztmp) continue;
	if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
	if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
	}

	jtype = map[type[j]];
	if (jtype != 0) continue;

	del32[0] = x[j][0]-x[i][0];
	del32[1] = x[j][1]-x[i][1];
	del32[2] = x[j][2]-x[i][2];
	rsq = del32[0]del32[0] + del32[1]del32[1] + del32[2]*del32[2];
	r32 = sqrt(rsq);
	del23[0] = -del32[0];
	del23[1] = -del32[1];
	del23[2] = -del32[2];
	r23 = r32;
	w23 = Sp(r23,rcmin[itype][jtype],rcmax[itype][jtype],dw23);

	for (kk = 0; kk < REBO_numneigh[i]; kk++) {
	k = REBO_neighs_i[kk];
	ktype = map[type[k]];
	if (k == j) continue;
	del21[0] = x[i][0]-x[k][0];
	del21[1] = x[i][1]-x[k][1];
	del21[2] = x[i][2]-x[k][2];
	rsq = del21[0]del21[0] + del21[1]del21[1] + del21[2]*del21[2];
	r21 = sqrt(rsq);
	cos321 = - ((del21[0]del32[0]) + (del21[1]del32[1]) +
	(del21[2]del32[2])) / (r21r32);
	cos321 = MIN(cos321,1.0);
	cos321 = MAX(cos321,-1.0);
	sin321 = sqrt(1.0 - cos321*cos321);
	if (sin321 < TOL) continue;

	deljk[0] = del21[0]-del23[0];
	deljk[1] = del21[1]-del23[1];
	deljk[2] = del21[2]-del23[2];
	rjk2 = deljk[0]deljk[0] + deljk[1]deljk[1] + deljk[2]*deljk[2];
	rjk=sqrt(rjk2);
	rik2 = r21*r21;
	w21 = Sp(r21,rcmin[itype][ktype],rcmax[itype][ktype],dw21);

	rij = r32;
	rik = r21;
	rij2 = r32*r32;
	rik2 = r21*r21;
	costmp = 0.5*(rij2+rik2-rjk2)/rij/rik;
	tspjik = Sp2(costmp,thmin,thmax,dtsjik);
	dtsjik = -dtsjik;

	REBO_neighs_j = REBO_firstneigh[j];
	for (ll = 0; ll < REBO_numneigh[j]; ll++) {
	l = REBO_neighs_j[ll];
	ltype = map[type[l]];
	if (l == i \|\| l == k) continue;
	del34[0] = x[j][0]-x[l][0];
	del34[1] = x[j][1]-x[l][1];
	del34[2] = x[j][2]-x[l][2];
	rsq = del34[0]del34[0] + del34[1]del34[1] + del34[2]*del34[2];
	r34 = sqrt(rsq);
	cos234 = (del32[0]del34[0] + del32[1]del34[1] +
	del32[2]del34[2]) / (r32r34);
	cos234 = MIN(cos234,1.0);
	cos234 = MAX(cos234,-1.0);
	sin234 = sqrt(1.0 - cos234*cos234);
	if (sin234 < TOL) continue;
	w34 = Sp(r34,rcmin[jtype][ltype],rcmax[jtype][ltype],dw34);
	delil[0] = del23[0] + del34[0];
	delil[1] = del23[1] + del34[1];
	delil[2] = del23[2] + del34[2];
	ril2 = delil[0]delil[0] + delil[1]delil[1] + delil[2]*delil[2];
	ril=sqrt(ril2);
	rjl2 = r34*r34;

	rjl = r34;
	rjl2 = r34*r34;
	costmp = 0.5*(rij2+rjl2-ril2)/rij/rjl;
	tspijl = Sp2(costmp,thmin,thmax,dtsijl);
	dtsijl = -dtsijl; //need minus sign
	cross321[0] = (del32[1]del21[2])-(del32[2]del21[1]);
	cross321[1] = (del32[2]del21[0])-(del32[0]del21[2]);
	cross321[2] = (del32[0]del21[1])-(del32[1]del21[0]);
	cross321mag = sqrt(cross321[0]*cross321[0]+
	cross321[1]*cross321[1]+
	cross321[2]*cross321[2]);
	cross234[0] = (del23[1]del34[2])-(del23[2]del34[1]);
	cross234[1] = (del23[2]del34[0])-(del23[0]del34[2]);
	cross234[2] = (del23[0]del34[1])-(del23[1]del34[0]);
	cross234mag = sqrt(cross234[0]*cross234[0]+
	cross234[1]*cross234[1]+
	cross234[2]*cross234[2]);
	cwnum = (cross321[0]*cross234[0]) +
	(cross321[1]cross234[1])+(cross321[2]cross234[2]);
	cwnom = r21r34r32r32sin321*sin234;
	cw = cwnum/cwnom;

	cw2 = (.5*(1.0-cw));
	ekijl = epsilonT[ktype][ltype];
	Ec = 256.0*ekijl/405.0;
	Vtors = (Ec*(powint(cw2,5)))-(ekijl/10.0);

	if (eflag) evdwl = Vtorsw21w23w34(1.0-tspjik)*(1.0-tspijl);

	dndij[0] = (cross234[1]del21[2])-(cross234[2]del21[1]);
	dndij[1] = (cross234[2]del21[0])-(cross234[0]del21[2]);
	dndij[2] = (cross234[0]del21[1])-(cross234[1]del21[0]);

	tmpvec[0] = (del34[1]cross321[2])-(del34[2]cross321[1]);
	tmpvec[1] = (del34[2]cross321[0])-(del34[0]cross321[2]);
	tmpvec[2] = (del34[0]cross321[1])-(del34[1]cross321[0]);

	dndij[0] = dndij[0]+tmpvec[0];
	dndij[1] = dndij[1]+tmpvec[1];
	dndij[2] = dndij[2]+tmpvec[2];

	dndik[0] = (del23[1]cross234[2])-(del23[2]cross234[1]);
	dndik[1] = (del23[2]cross234[0])-(del23[0]cross234[2]);
	dndik[2] = (del23[0]cross234[1])-(del23[1]cross234[0]);

	dndjl[0] = (cross321[1]del23[2])-(cross321[2]del23[1]);
	dndjl[1] = (cross321[2]del23[0])-(cross321[0]del23[2]);
	dndjl[2] = (cross321[0]del23[1])-(cross321[1]del23[0]);

	dcidij = ((r23r23)-(r21r21)+(rjkrjk))/(2.0r23r23r21);
	dcidik = ((r21r21)-(r23r23)+(rjkrjk))/(2.0r23r21r21);
	dcidjk = (-rjk)/(r23*r21);
	dcjdji = ((r23r23)-(r34r34)+(rilril))/(2.0r23r23r34);
	dcjdjl = ((r34r34)-(r23r23)+(rilril))/(2.0r23r34r34);
	dcjdil = (-ril)/(r23*r34);

	dsidij = (-cos321/sin321)*dcidij;
	dsidik = (-cos321/sin321)*dcidik;
	dsidjk = (-cos321/sin321)*dcidjk;

	dsjdji = (-cos234/sin234)*dcjdji;
	dsjdjl = (-cos234/sin234)*dcjdjl;
	dsjdil = (-cos234/sin234)*dcjdil;

	dxidij = (r21sin321)+(r23r21*dsidij);
	dxidik = (r23sin321)+(r23r21*dsidik);
	dxidjk = (r23r21dsidjk);

	dxjdji = (r34sin234)+(r23r34*dsjdji);
	dxjdjl = (r23sin234)+(r23r34*dsjdjl);
	dxjdil = (r23r34dsjdil);

	ddndij = (dxidijcross234mag)+(cross321magdxjdji);
	ddndik = dxidik*cross234mag;
	ddndjk = dxidjk*cross234mag;
	ddndjl = cross321mag*dxjdjl;
	ddndil = cross321mag*dxjdil;
	dcwddn = -cwnum/(cwnom*cwnom);
	dcwdn = 1.0/cwnom;
	dvpdcw = (-1.0)Ec(-.5)5.0powint(cw2,4) *
	w23w21w34(1.0-tspjik)(1.0-tspijl);

	Ftmp[0] = dvpdcw((dcwdndndij[0])+(dcwddnddndijdel23[0]/r23));
	Ftmp[1] = dvpdcw((dcwdndndij[1])+(dcwddnddndijdel23[1]/r23));
	Ftmp[2] = dvpdcw((dcwdndndij[2])+(dcwddnddndijdel23[2]/r23));
	fi[0] = Ftmp[0];
	fi[1] = Ftmp[1];
	fi[2] = Ftmp[2];
	fj[0] = -Ftmp[0];
	fj[1] = -Ftmp[1];
	fj[2] = -Ftmp[2];

	Ftmp[0] = dvpdcw((dcwdndndik[0])+(dcwddnddndikdel21[0]/r21));
	Ftmp[1] = dvpdcw((dcwdndndik[1])+(dcwddnddndikdel21[1]/r21));
	Ftmp[2] = dvpdcw((dcwdndndik[2])+(dcwddnddndikdel21[2]/r21));
	fi[0] += Ftmp[0];
	fi[1] += Ftmp[1];
	fi[2] += Ftmp[2];
	fk[0] = -Ftmp[0];
	fk[1] = -Ftmp[1];
	fk[2] = -Ftmp[2];

	Ftmp[0] = (dvpdcwdcwddnddndjk*deljk[0])/rjk;
	Ftmp[1] = (dvpdcwdcwddnddndjk*deljk[1])/rjk;
	Ftmp[2] = (dvpdcwdcwddnddndjk*deljk[2])/rjk;
	fj[0] += Ftmp[0];
	fj[1] += Ftmp[1];
	fj[2] += Ftmp[2];
	fk[0] -= Ftmp[0];
	fk[1] -= Ftmp[1];
	fk[2] -= Ftmp[2];

	Ftmp[0] = dvpdcw((dcwdndndjl[0])+(dcwddnddndjldel34[0]/r34));
	Ftmp[1] = dvpdcw((dcwdndndjl[1])+(dcwddnddndjldel34[1]/r34));
	Ftmp[2] = dvpdcw((dcwdndndjl[2])+(dcwddnddndjldel34[2]/r34));
	fj[0] += Ftmp[0];
	fj[1] += Ftmp[1];
	fj[2] += Ftmp[2];
	fl[0] = -Ftmp[0];
	fl[1] = -Ftmp[1];
	fl[2] = -Ftmp[2];

	Ftmp[0] = (dvpdcwdcwddnddndil*delil[0])/ril;
	Ftmp[1] = (dvpdcwdcwddnddndil*delil[1])/ril;
	Ftmp[2] = (dvpdcwdcwddnddndil*delil[2])/ril;
	fi[0] += Ftmp[0];
	fi[1] += Ftmp[1];
	fi[2] += Ftmp[2];
	fl[0] -= Ftmp[0];
	fl[1] -= Ftmp[1];
	fl[2] -= Ftmp[2];

	// coordination forces

	fpair = Vtorsdw21w23w34(1.0-tspjik)*(1.0-tspijl) / r21;
	fi[0] -= del21[0]*fpair;
	fi[1] -= del21[1]*fpair;
	fi[2] -= del21[2]*fpair;
	fk[0] += del21[0]*fpair;
	fk[1] += del21[1]*fpair;
	fk[2] += del21[2]*fpair;

	fpair = Vtorsw21dw23w34(1.0-tspjik)*(1.0-tspijl) / r23;
	fi[0] -= del23[0]*fpair;
	fi[1] -= del23[1]*fpair;
	fi[2] -= del23[2]*fpair;
	fj[0] += del23[0]*fpair;
	fj[1] += del23[1]*fpair;
	fj[2] += del23[2]*fpair;

	fpair = Vtorsw21w23dw34(1.0-tspjik)*(1.0-tspijl) / r34;
	fj[0] -= del34[0]*fpair;
	fj[1] -= del34[1]*fpair;
	fj[2] -= del34[2]*fpair;
	fl[0] += del34[0]*fpair;
	fl[1] += del34[1]*fpair;
	fl[2] += del34[2]*fpair;

	// additional cut off function forces

	fcpc = -Vtorsw21w23w34dtsjik*(1.0-tspijl);
	fpair = fcpc*dcidij/rij;
	fi[0] += fpair*del23[0];
	fi[1] += fpair*del23[1];
	fi[2] += fpair*del23[2];
	fj[0] -= fpair*del23[0];
	fj[1] -= fpair*del23[1];
	fj[2] -= fpair*del23[2];

	fpair = fcpc*dcidik/rik;
	fi[0] += fpair*del21[0];
	fi[1] += fpair*del21[1];
	fi[2] += fpair*del21[2];
	fk[0] -= fpair*del21[0];
	fk[1] -= fpair*del21[1];
	fk[2] -= fpair*del21[2];

	fpair = fcpc*dcidjk/rjk;
	fj[0] += fpair*deljk[0];
	fj[1] += fpair*deljk[1];
	fj[2] += fpair*deljk[2];
	fk[0] -= fpair*deljk[0];
	fk[1] -= fpair*deljk[1];
	fk[2] -= fpair*deljk[2];

	fcpc = -Vtorsw21w23w34(1.0-tspjik)*dtsijl;
	fpair = fcpc*dcjdji/rij;
	fi[0] += fpair*del23[0];
	fi[1] += fpair*del23[1];
	fi[2] += fpair*del23[2];
	fj[0] -= fpair*del23[0];
	fj[1] -= fpair*del23[1];
	fj[2] -= fpair*del23[2];

	fpair = fcpc*dcjdjl/rjl;
	fj[0] += fpair*del34[0];
	fj[1] += fpair*del34[1];
	fj[2] += fpair*del34[2];
	fl[0] -= fpair*del34[0];
	fl[1] -= fpair*del34[1];
	fl[2] -= fpair*del34[2];

	fpair = fcpc*dcjdil/ril;
	fi[0] += fpair*delil[0];
	fi[1] += fpair*delil[1];
	fi[2] += fpair*delil[2];
	fl[0] -= fpair*delil[0];
	fl[1] -= fpair*delil[1];
	fl[2] -= fpair*delil[2];

	// sum per-atom forces into atom force array

	f[i][0] += fi[0]; f[i][1] += fi[1]; f[i][2] += fi[2];
	f[j][0] += fj[0]; f[j][1] += fj[1]; f[j][2] += fj[2];
	f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2];
	f[l][0] += fl[0]; f[l][1] += fl[1]; f[l][2] += fl[2];

	if (evflag) {
	delkl[0] = delil[0] - del21[0];
	delkl[1] = delil[1] - del21[1];
	delkl[2] = delil[2] - del21[2];
	ev_tally4_thr(this,i,j,k,l,evdwl,fi,fj,fk,delil,del34,delkl,thr);
	}
	}
	}
	}
	}
	}

	/* ----------------------------------------------------------------------
	create REBO neighbor list from main neighbor list
	REBO neighbor list stores neighbors of ghost atoms
	------------------------------------------------------------------------- */

	void PairAIREBOOMP::REBO_neigh_thr()
	{
	const int nthreads = comm->nthreads;

	if (atom->nmax > maxlocal) {
	maxlocal = atom->nmax;
	memory->destroy(REBO_numneigh);
	memory->sfree(REBO_firstneigh);
	memory->destroy(nC);
	memory->destroy(nH);
	memory->create(REBO_numneigh,maxlocal,"AIREBO:numneigh");
	REBO_firstneigh = (int *) memory->smalloc(maxlocalsizeof(int *),
	"AIREBO:firstneigh");
	memory->create(nC,maxlocal,"AIREBO:nC");
	memory->create(nH,maxlocal,"AIREBO:nH");
	}

	#if defined(_OPENMP)
	#pragma omp parallel default(none)
	#endif
	{
	int i,j,ii,jj,n,jnum,itype,jtype;
	double xtmp,ytmp,ztmp,delx,dely,delz,rsq,dS;
	int ilist,jlist,numneigh,*firstneigh;
	int *neighptr;

	double **x = atom->x;
	int *type = atom->type;

	const int allnum = list->inum + list->gnum;
	ilist = list->ilist;
	numneigh = list->numneigh;
	firstneigh = list->firstneigh;

	#if defined(_OPENMP)
	const int tid = omp_get_thread_num();
	#else
	const int tid = 0;
	#endif

	const int iidelta = 1 + allnum/nthreads;
	const int iifrom = tid*iidelta;
	const int iito = ((iifrom+iidelta)>allnum) ? allnum : (iifrom+iidelta);

	// store all REBO neighs of owned and ghost atoms
	// scan full neighbor list of I

	// each thread has its own page allocator
	MyPage<int> &ipg = ipage[tid];
	ipg.reset();

	for (ii = iifrom; ii < iito; ii++) {
	i = ilist[ii];

	n = 0;
	neighptr = ipg.vget();

	xtmp = x[i][0];
	ytmp = x[i][1];
	ztmp = x[i][2];
	itype = map[type[i]];
	nC[i] = nH[i] = 0.0;
	jlist = firstneigh[i];
	jnum = numneigh[i];

	for (jj = 0; jj < jnum; jj++) {
	j = jlist[jj];
	j &= NEIGHMASK;
	jtype = map[type[j]];
	delx = xtmp - x[j][0];
	dely = ytmp - x[j][1];
	delz = ztmp - x[j][2];
	rsq = delxdelx + delydely + delz*delz;

	if (rsq < rcmaxsq[itype][jtype]) {
	neighptr[n++] = j;
	if (jtype == 0)
	nC[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
	else
	nH[i] += Sp(sqrt(rsq),rcmin[itype][jtype],rcmax[itype][jtype],dS);
	}
	}

	REBO_firstneigh[i] = neighptr;
	REBO_numneigh[i] = n;
	ipg.vgot(n);
	if (ipg.status())
	error->one(FLERR,"REBO list overflow, boost neigh_modify one");
	}
	}
	}


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

	double PairAIREBOOMP::memory_usage()
	{
	double bytes = memory_usage_thr();
	bytes += PairAIREBO::memory_usage();

	return bytes;
	}

pair_airebo_omp.cppNo OneTemporaryActions

File Metadata

pair_airebo_omp.cppView Options

Event Timeline

pair_airebo_omp.cpp
No OneTemporary
Actions

pair_airebo_omp.cpp
View Options