pair_airebo_intel.cpp
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Apr 26, 11:36

pair_airebo_intel.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 author: Markus Hohnerbach (RWTH)
	------------------------------------------------------------------------- */

	#ifdef __INTEL_OFFLOAD
	#pragma offload_attribute(push, target(mic))
	#endif
	#include <unistd.h>
	#include <stdlib.h>
	#include <math.h>
	#include <stdio.h>
	#include <stdint.h>
	#include <assert.h>
	#include <stddef.h>
	#include "lmptype.h"
	#include "intel_preprocess.h"
	#include "intel_intrinsics_airebo.h"
	#ifdef __INTEL_OFFLOAD
	#pragma offload_attribute(pop)
	#endif

	#include <omp.h>
	#include <string.h>
	#include "pair_airebo_intel.h"
	#include "atom.h"
	#include "neighbor.h"
	#include "neigh_list.h"
	#include "neigh_request.h"
	#include "force.h"
	#include "comm.h"
	#include "memory.h"
	#include "error.h"
	#include "group.h"
	#include "kspace.h"
	#include "modify.h"
	#include "suffix.h"

	using namespace LAMMPS_NS;

	#ifdef __INTEL_OFFLOAD
	#pragma offload_attribute(push, target(mic))
	#endif

	template<typename flt_t, typename acc_t>
	struct LAMMPS_NS::PairAIREBOIntelParam {
	flt_t cutlj, cutljrebosq, cut3rebo;
	flt_t sigmin, sigcut;
	flt_t cutljsq[2][2];
	flt_t lj1[2][2], lj2[2][2], lj3[2][2], lj4[2][2];

	flt_t smin, Nmin, Nmax, NCmin, NCmax, thmin, thmax;
	flt_t rcmin[2][2], rcmax[2][2], rcmaxsq[2][2], rcmaxp[2][2];
	flt_t Q[2][2], alpha[2][2], A[2][2], rho[2][2], BIJc[2][2][3],
	Beta[2][2][3];
	flt_t rcLJmin[2][2], rcLJmax[2][2], rcLJmaxsq[2][2], bLJmin[2][2],
	bLJmax[2][2];
	flt_t epsilon[2][2], sigma[2][2], epsilonT[2][2];

	// spline coefficients

	flt_t gCdom[5], gC1[4][6], gC2[4][6], gHdom[4], gH[3][6];
	flt_t gDom[5+4];
	flt_t gVal[(4+4+3)*6];
	flt_t pCCdom[2][2], pCHdom[2][2], pCC[4][4][16], pCH[4][4][16];
	flt_t piCCdom[3][2], piCHdom[3][2], piHHdom[3][2];
	acc_t piCC[4][4][9][64], piCH[4][4][9][64], piHH[4][4][9][64];
	flt_t Tijdom[3][2];
	acc_t Tijc[4][4][9][64];

	// spline knot values

	flt_t PCCf[5][5], PCCdfdx[5][5], PCCdfdy[5][5], PCHf[5][5];
	flt_t PCHdfdx[5][5], PCHdfdy[5][5];
	flt_t piCCf[5][5][11], piCCdfdx[5][5][11];
	flt_t piCCdfdy[5][5][11], piCCdfdz[5][5][11];
	flt_t piCHf[5][5][11], piCHdfdx[5][5][11];
	flt_t piCHdfdy[5][5][11], piCHdfdz[5][5][11];
	flt_t piHHf[5][5][11], piHHdfdx[5][5][11];
	flt_t piHHdfdy[5][5][11], piHHdfdz[5][5][11];
	flt_t Tf[5][5][10], Tdfdx[5][5][10], Tdfdy[5][5][10], Tdfdz[5][5][10];
	};

	namespace {

	struct NeighListAIREBO {
	int * num; /* num_all */
	int * num_half; /* num_all */
	int * offset; /* num_all */
	int * entries; /* num_all * num_neighs_per_atom */
	};

	template<typename flt_t>
	struct AtomAIREBOT {
	flt_t x, y, z;
	int w;
	};

	template<typename acc_t>
	struct ResultForceT {
	acc_t x, y, z, w;
	};

	template<typename flt_t, typename acc_t>
	struct KernelArgsAIREBOT {
	int num_local;
	int num_all;
	int num_neighs_per_atom;
	int num_types;
	int frebo_from_atom, frebo_to_atom;
	int neigh_from_atom, neigh_to_atom;
	int rebuild_flag;
	flt_t skin;
	struct NeighListAIREBO neigh_lmp;
	struct NeighListAIREBO neigh_rebo;
	PairAIREBOIntelParam<flt_t,acc_t> params;
	struct AtomAIREBOT<flt_t> * x; /* num_all */
	int * tag; /* num_all */
	flt_t * nC, * nH; /* num_all */
	int * map; /* num_types+1 */
	struct ResultForceT<acc_t> * result_f; /* num_all */
	acc_t result_eng;
	};

	template<typename flt_t, typename acc_t>
	void aut_lennard_jones(KernelArgsAIREBOT<flt_t,acc_t> * ka, int morseflag);
	template<typename flt_t, typename acc_t>
	void aut_rebo_neigh(KernelArgsAIREBOT<flt_t,acc_t> * ka);
	template<typename flt_t, typename acc_t>
	void aut_frebo(KernelArgsAIREBOT<flt_t,acc_t> * ka, int torsion_flag);

	}

	#ifdef __INTEL_OFFLOAD
	#pragma offload_attribute(pop)
	#endif

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

	PairAIREBOIntel::PairAIREBOIntel(LAMMPS *lmp) : PairAIREBO(lmp)
	{
	suffix_flag \|= Suffix::INTEL;
	REBO_cnumneigh = NULL;
	REBO_num_skin = NULL;
	REBO_list_data = NULL;
	fix = NULL;
	}

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

	PairAIREBOIntel::~PairAIREBOIntel()
	{
	memory->destroy(REBO_cnumneigh);
	memory->destroy(REBO_num_skin);
	memory->destroy(REBO_list_data);
	}

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

	void PairAIREBOIntel::init_style()
	{
	PairAIREBO::init_style();
	neighbor->requests[neighbor->nrequest-1]->intel = 1;

	int ifix = modify->find_fix("package_intel");
	if (ifix < 0)
	error->all(FLERR,
	"The 'package intel' command is required for /intel styles");
	fix = static_cast<FixIntel *>(modify->fix[ifix]);

	fix->pair_init_check();
	#ifdef _LMP_INTEL_OFFLOAD
	_cop = fix->coprocessor_number();
	#endif

	if (fix->precision() == FixIntel::PREC_MODE_MIXED) {
	pack_force_const(fix->get_mixed_buffers());
	fix->get_mixed_buffers()->need_tag(1);
	} else if (fix->precision() == FixIntel::PREC_MODE_DOUBLE) {
	pack_force_const(fix->get_double_buffers());
	fix->get_double_buffers()->need_tag(1);
	} else {
	pack_force_const(fix->get_single_buffers());
	fix->get_single_buffers()->need_tag(1);
	}

	#ifdef _LMP_INTEL_OFFLOAD
	if (fix->offload_noghost())
	error->all(FLERR,"The 'ghost no' option cannot be used with airebo/intel.");
	#endif
	}

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

	template<typename T>
	T * calloc_it(size_t size) {
	return static_cast<T*>(calloc(size, sizeof(T)));
	}

	void PairAIREBOIntel::compute(int eflag, int vflag)
	{
	if (fix->precision()==FixIntel::PREC_MODE_MIXED)
	compute<float,double>(eflag, vflag, fix->get_mixed_buffers());
	else if (fix->precision()==FixIntel::PREC_MODE_DOUBLE)
	compute<double,double>(eflag, vflag, fix->get_double_buffers());
	else
	compute<float,float>(eflag, vflag, fix->get_single_buffers());

	fix->balance_stamp();
	vflag_fdotr = 0;
	}

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

	template<class flt_t, class acc_t>
	PairAIREBOIntelParam<flt_t,acc_t> PairAIREBOIntel::get_param()
	{
	PairAIREBOIntelParam<flt_t,acc_t> fc;

	#define A(a) \
	for (int i = 0; i < sizeof(this->a)/sizeof(double); i++) { \
	reinterpret_cast<flt_t*>(&fc.a)[i] = \
	reinterpret_cast<double*>(&this->a)[i]; \
	}
	#define A0(a) \
	for (int i = 0; i < sizeof(fc.a)/sizeof(flt_t); i++) { \
	reinterpret_cast<flt_t*>(&fc.a)[i] = \
	reinterpret_cast<double*>(this->a[0])[i]; \
	}
	#define B(a) \
	for (int i = 0; i < sizeof(this->a)/sizeof(double); i++) { \
	reinterpret_cast<acc_t*>(&fc.a)[i] = \
	reinterpret_cast<double*>(&this->a)[i]; \
	}

	A(cutlj) A(cutljrebosq) A(cut3rebo) A(sigmin);
	A(sigcut) A0(cutljsq) A0(lj1) A0(lj2) A0(lj3);
	A0(lj4) A(smin) A(Nmin) A(Nmax) A(NCmin) A(NCmax) A(thmin) A(thmax);
	A(rcmin) A(rcmax) A(rcmaxsq) A(rcmaxp) A(Q) A(alpha) A(A) A(rho) A(BIJc);
	A(Beta) A(rcLJmin) A(rcLJmax) A(rcLJmaxsq) A(bLJmin) A(bLJmax) A(epsilon);
	A(sigma) A(epsilonT) A(gCdom) A(gC1) A(gC2) A(gHdom) A(gH) A(pCCdom);
	A(pCHdom) A(pCC) A(pCH) A(piCCdom) A(piCHdom) A(piHHdom) B(piCC);
	B(piCH) B(piHH) A(Tijdom) B(Tijc) A(PCCf) A(PCCdfdx) A(PCCdfdy) A(PCHf);
	A(PCHdfdx) A(PCHdfdy) A(piCCf) A(piCCdfdx) A(piCCdfdy) A(piCCdfdz);
	A(piCHf) A(piCHdfdx) A(piCHdfdy) A(piCHdfdz) A(piHHf) A(piHHdfdx);
	A(piHHdfdy) A(piHHdfdz) A(Tf) A(Tdfdx) A(Tdfdy) A(Tdfdz);

	#undef A
	#undef A0
	#undef B
	for (int i = 0; i < 5; i++) fc.gDom[i] = fc.gCdom[i];
	for (int i = 0; i < 4; i++) fc.gDom[5+i] = fc.gHdom[i];
	for (int i = 0; i < 4; i++) for (int j = 0; j < 6; j++)
	fc.gVal[6*i+j] = fc.gC1[i][j];
	for (int i = 0; i < 4; i++) for (int j = 0; j < 6; j++)
	fc.gVal[46+6i+j] = fc.gC2[i][j];
	for (int i = 0; i < 3; i++) for (int j = 0; j < 6; j++)
	fc.gVal[86+6i+j] = fc.gH[i][j];

	return fc;
	}

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

	template<class flt_t, class acc_t>
	void PairAIREBOIntel::compute(
	int eflag, int vflag, IntelBuffers<flt_t,acc_t> * buffers
	) {
	if (eflag \|\| vflag) ev_setup(eflag,vflag);
	else evflag = vflag_fdotr = vflag_atom = 0;
	pvector[0] = pvector[1] = pvector[2] = 0.0;

	const int inum = list->inum;
	const int nthreads = comm->nthreads;
	const int host_start = fix->host_start_pair();
	const int offload_end = fix->offload_end_pair();
	const int ago = neighbor->ago;

	if (ago != 0 && fix->separate_buffers() == 0) {
	fix->start_watch(TIME_PACK);
	int packthreads;
	if (nthreads > INTEL_HTHREADS) packthreads = nthreads;
	else packthreads = 1;
	#if defined(_OPENMP)
	#pragma omp parallel if(packthreads > 1)
	#endif
	{
	int ifrom, ito, tid;
	IP_PRE_omp_range_id_align(ifrom, ito, tid, atom->nlocal + atom->nghost,
	packthreads, sizeof(ATOM_T));
	buffers->thr_pack(ifrom,ito,ago);
	}
	fix->stop_watch(TIME_PACK);
	}

	if (atom->nmax > maxlocal) {
	#ifdef LMP_INTEL_OFFLOAD
	if (maxlocal > 0 && _cop >= 0) {
	int * const REBO_numneigh = this->REBO_numneigh;
	int * const REBO_num_skin = this->REBO_num_skin;
	int * const REBO_cnumneigh = this->REBO_cnumneigh;
	int * const REBO_list_data = this->REBO_list_data;
	double * const nC = this->nC;
	double * const nH = this->nH;
	#pragma offload_transfer target(mic:_cop) \
	nocopy(REBO_numneigh: alloc_if(0) free_if(1)) \
	nocopy(REBO_cnumneigh: alloc_if(0) free_if(1)) \
	nocopy(REBO_num_skin: alloc_if(0) free_if(1)) \
	nocopy(REBO_list_data: alloc_if(0) free_if(1)) \
	nocopy(nH: alloc_if(0) free_if(1)) \
	nocopy(nC: alloc_if(0) free_if(1))
	}
	#endif
	maxlocal = atom->nmax;
	memory->destroy(REBO_numneigh);
	memory->destroy(REBO_cnumneigh);
	memory->destroy(REBO_list_data);
	memory->sfree(REBO_firstneigh);
	memory->destroy(nC);
	memory->destroy(nH);
	memory->create(REBO_numneigh,maxlocal,"AIREBO:numneigh");
	memory->create(REBO_cnumneigh,maxlocal,"AIREBO:cnumneigh");
	memory->create(REBO_num_skin,maxlocal,"AIREBO:cnumneigh");
	int max_nbors = buffers->get_max_nbors();
	memory->create(REBO_list_data,maxlocal * max_nbors,"AIREBO:list_data");
	REBO_firstneigh = (int *) memory->smalloc(maxlocalsizeof(int *),
	"AIREBO:firstneigh");
	memory->create(nC,maxlocal,"AIREBO:nC");
	memory->create(nH,maxlocal,"AIREBO:nH");
	#ifdef _LMP_INTEL_OFFLOAD
	if (_cop >= 0) {
	int * const REBO_numneigh = this->REBO_numneigh;
	int * const REBO_num_skin = this->REBO_num_skin;
	int * const REBO_cnumneigh = this->REBO_cnumneigh;
	int * const REBO_list_data = this->REBO_list_data;
	double * const nC = this->nC;
	double * const nH = this->nH;
	const int mnml = max_nbors * maxlocal;
	#pragma offload_transfer target(mic:_cop) \
	nocopy(REBO_numneigh: length(maxlocal) alloc_if(1) free_if(0)) \
	nocopy(REBO_cnumneigh:length(maxlocal) alloc_if(1) free_if(0)) \
	nocopy(REBO_num_skin: length(maxlocal) alloc_if(1) free_if(0)) \
	nocopy(REBO_list_data:length(mnml) alloc_if(1) free_if(0)) \
	nocopy(nH: length(maxlocal) alloc_if(1) free_if(0)) \
	nocopy(nC: length(maxlocal) alloc_if(1) free_if(0))
	}
	#endif
	}

	if (evflag \|\| vflag_fdotr) {
	int ovflag = 0;
	if (vflag_fdotr) ovflag = 2;
	else if (vflag) ovflag = 1;
	if (eflag) {
	eval<1,1>(1, ovflag, buffers, 0, offload_end);
	eval<1,1>(0, ovflag, buffers, host_start, inum);
	} else {
	eval<1,0>(1, ovflag, buffers, 0, offload_end);
	eval<1,0>(0, ovflag, buffers, host_start, inum);
	}
	} else {
	eval<0,0>(1, 0, buffers, 0, offload_end);
	eval<0,0>(0, 0, buffers, host_start, inum);
	}
	}

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

	template<int EVFLAG, int EFLAG, class flt_t, class acc_t>
	void PairAIREBOIntel::eval(
	const int offload, const int vflag,
	IntelBuffers<flt_t,acc_t> * buffers,
	const int astart, const int aend
	) {
	const int inum = aend - astart;
	if (inum == 0) {
	return;
	}
	int nlocal, nall, minlocal;
	fix->get_buffern(offload, nlocal, nall, minlocal);

	const int ago = neighbor->ago;
	IP_PRE_pack_separate_buffers(fix, buffers, ago, offload, nlocal, nall);

	ATOM_T * _noalias const x = buffers->get_x(offload);
	const int * _noalias const numneighhalf = buffers->get_atombin();
	const int * _noalias const numneigh = list->numneigh;
	const int * _noalias const cnumneigh = buffers->cnumneigh(list);
	const int * _noalias const firstneigh = buffers->firstneigh(list);
	int * const tag = atom->tag;

	const int ntypes = atom->ntypes + 1;
	const int eatom = this->eflag_atom;

	int x_size, q_size, f_stride, ev_size, separate_flag;
	IP_PRE_get_transfern(ago, 1 /NEWTON_PAIR/, EFLAG, vflag,
	buffers, offload, fix, separate_flag,
	x_size, q_size, ev_size, f_stride);

	int tc;
	FORCE_T * _noalias f_start;
	acc_t * _noalias ev_global;
	IP_PRE_get_buffers(offload, buffers, fix, tc, f_start, ev_global);

	const int nthreads = tc;
	const double skin = neighbor->skin;
	const int max_nbor = buffers->get_max_nbors();
	const PairAIREBOIntelParam<flt_t,acc_t> param = get_param<flt_t,acc_t>();

	// offload here
	#ifdef _LMP_INTEL_OFFLOAD
	int *overflow = fix->get_off_overflow_flag();
	double *timer_compute = fix->off_watch_pair();

	int * const REBO_numneigh = this->REBO_numneigh;
	int * const REBO_num_skin = this->REBO_num_skin;
	int * const REBO_cnumneigh = this->REBO_cnumneigh;
	int * const REBO_list_data = this->REBO_list_data;
	double * const nC = this->nC;
	double * const nH = this->nH;
	const int torflag = this->torflag;
	const int ljflag = this->ljflag;
	const int morseflag = this->morseflag;
	int * const map = this->map;

	if (offload) fix->start_watch(TIME_OFFLOAD_LATENCY);

	#pragma offload target(mic:_cop) if(offload) \
	in(firstneigh:length(0) alloc_if(0) free_if(0)) \
	in(cnumneigh:length(0) alloc_if(0) free_if(0)) \
	in(numneigh:length(0) alloc_if(0) free_if(0)) \
	in(numneighhalf:length(0) alloc_if(0) free_if(0)) \
	in(x:length(x_size) alloc_if(0) free_if(0)) \
	in(overflow:length(0) alloc_if(0) free_if(0)) \
	in(astart,nthreads,inum,nall,ntypes,vflag,eatom) \
	in(f_stride,nlocal,minlocal,separate_flag,offload) \
	out(f_start:length(f_stride) alloc_if(0) free_if(0)) \
	out(ev_global:length(ev_size) alloc_if(0) free_if(0)) \
	out(timer_compute:length(1) alloc_if(0) free_if(0)) \
	in(param,skin,max_nbor) \
	in(tag: length(0) alloc_if(0) free_if(0)) \
	in(torflag, ljflag, morseflag, ago) \
	in(nC: length(0) alloc_if(0) free_if(0)) \
	in(nH: length(0) alloc_if(0) free_if(0)) \
	in(REBO_numneigh: length(0) alloc_if(0) free_if(0)) \
	in(REBO_cnumneigh: length(0) alloc_if(0) free_if(0)) \
	in(REBO_num_skin: length(0) alloc_if(0) free_if(0)) \
	in(REBO_list_data: length(0) alloc_if(0) free_if(0)) \
	in(map: length(0) alloc_if(0) free_if(0)) \
	signal(f_start)
	#endif
	{
	#if defined(__MIC__) && defined(_LMP_INTEL_OFFLOAD)
	*timer_compute = MIC_Wtime();
	#endif

	IP_PRE_repack_for_offload(1 /NEWTON_PAIR/, separate_flag, nlocal, nall,
	f_stride, x, 0/q/);

	acc_t oevdwl, oecoul, ov0, ov1, ov2, ov3, ov4, ov5;
	if (EVFLAG) {
	oevdwl = oecoul = (acc_t)0;
	if (vflag) ov0 = ov1 = ov2 = ov3 = ov4 = ov5 = (acc_t)0;
	}

	// loop over neighbors of my atoms
	#if defined(_OPENMP)
	#pragma omp parallel \
	shared(f_start,f_stride,nlocal,nall,minlocal) \
	reduction(+:oevdwl,oecoul,ov0,ov1,ov2,ov3,ov4,ov5)
	#endif
	{
	int iifrom, iito, tid;
	IP_PRE_omp_range_id(iifrom, iito, tid, inum, nthreads);
	iifrom += astart;
	iito += astart;

	int neigh_iifrom, neigh_iito;
	IP_PRE_omp_range(neigh_iifrom, neigh_iito, tid, nall, nthreads);

	FORCE_T * _noalias const f = f_start - minlocal + (tid * f_stride);
	memset(f + minlocal, 0, f_stride * sizeof(FORCE_T));

	KernelArgsAIREBOT<flt_t,acc_t> args;
	args.num_local = nlocal;
	args.num_all = nall;
	args.num_neighs_per_atom = max_nbor;
	args.num_types = ntypes;
	args.frebo_from_atom = 0;
	args.frebo_to_atom = args.num_local;
	args.neigh_from_atom = 0;
	args.neigh_to_atom = args.num_all;
	args.rebuild_flag = ago == 0;
	args.skin = skin;
	args.neigh_lmp.num = const_cast<int*>(numneigh);
	args.neigh_lmp.num_half = const_cast<int*>(numneighhalf);
	args.neigh_lmp.offset = const_cast<int*>(cnumneigh);
	args.neigh_lmp.entries = const_cast<int*>(firstneigh);
	args.neigh_rebo.num = REBO_numneigh;
	args.neigh_rebo.num_half = REBO_num_skin;
	args.neigh_rebo.offset = REBO_cnumneigh;
	args.neigh_rebo.entries = REBO_list_data;
	args.params = param;
	args.tag = tag;
	args.nC = reinterpret_cast<flt_t*>(nC);
	args.nH = reinterpret_cast<flt_t*>(nH);
	args.map = map;
	args.result_eng = 0;
	args.x = (AtomAIREBOT<flt_t>*) x;

	args.result_f = (ResultForceT<acc_t> *) f;
	args.neigh_from_atom = neigh_iifrom;
	args.neigh_to_atom = neigh_iito;
	args.frebo_from_atom = iifrom;
	args.frebo_to_atom = iito;

	aut_rebo_neigh(&args);
	#if defined(_OPENMP)
	#pragma omp barrier
	#endif
	aut_frebo(&args, torflag);
	if (ljflag) aut_lennard_jones(&args, morseflag);

	oevdwl += args.result_eng;

	IP_PRE_fdotr_reduce_omp(1, nall, minlocal, nthreads, f_start, f_stride, x,
	offload, vflag, ov0, ov1, ov2, ov3, ov4, ov5);
	} // end of omp parallel region
	IP_PRE_fdotr_reduce(1, nall, nthreads, f_stride, vflag,
	ov0, ov1, ov2, ov3, ov4, ov5);
	if (EVFLAG) {
	if (EFLAG) {
	ev_global[0] = oevdwl;
	ev_global[1] = oecoul;
	}
	if (vflag) {
	ev_global[2] = ov0;
	ev_global[3] = ov1;
	ev_global[4] = ov2;
	ev_global[5] = ov3;
	ev_global[6] = ov4;
	ev_global[7] = ov5;
	}
	}
	#if defined(__MIC__) && defined(_LMP_INTEL_OFFLOAD)
	timer_compute = MIC_Wtime() - timer_compute;
	#endif
	} // end of offload region

	if (offload)
	fix->stop_watch(TIME_OFFLOAD_LATENCY);
	else
	fix->stop_watch(TIME_HOST_PAIR);

	if (EVFLAG)
	fix->add_result_array(f_start, ev_global, offload, eatom, 0, vflag);
	else
	fix->add_result_array(f_start, 0, offload);
	}

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

	template<class flt_t, class acc_t>
	void PairAIREBOIntel::pack_force_const(IntelBuffers<flt_t,acc_t> * buffers) {
	int tp1 = atom->ntypes + 1;

	buffers->set_ntypes(tp1,1);
	flt_t **cutneighsq = buffers->get_cutneighsq();
	flt_t **cutneighghostsq = buffers->get_cutneighghostsq();

	// Repeat cutsq calculation because done after call to init_style
	double cut, cutneigh;
	for (int i = 1; i <= atom->ntypes; i++) {
	for (int j = i; j <= atom->ntypes; j++) {
	if (setflag[i][j] != 0 \|\| (setflag[i][i] != 0 && setflag[j][j] != 0)) {
	cut = init_one(i, j);
	cutneigh = cut + neighbor->skin;
	cutsq[i][j] = cutsq[j][i] = cut*cut;
	cutneighsq[i][j] = cutneighsq[j][i] = cutneigh * cutneigh;
	cut = cutghost[i][j] + neighbor->skin;
	cutneighghostsq[i][j] = cutneighghostsq[j][i] = cut*cut;
	}
	}
	}

	#ifdef _LMP_INTEL_OFFLOAD
	if (_cop < 0) return;
	flt_t * ocutneighsq = cutneighsq[0];
	size_t VL = 512 / 8 / sizeof(flt_t);
	int ntypes = tp1;
	int tp1sq = tp1 * tp1;
	// TODO the lifecycle of "map" is currently not 100% correct
	// it might not be freed if this method is called more than once
	int * map = this->map;
	#pragma offload_transfer target(mic:_cop) \
	in(ocutneighsq: length(tp1sq) alloc_if(0) free_if(0)) \
	in(map: length(tp1) alloc_if(1) free_if(0))
	#endif

	}

	/* ----------------------------------------------------------------------
	Implementation
	---------------------------------------------------------------------- */

	namespace {

	#ifdef __INTEL_OFFLOAD
	#pragma offload_attribute(push, target(mic))
	#endif

	namespace overloaded {
	double sqrt(double a) { return ::sqrt(a); }
	float sqrt(float a) { return ::sqrtf(a); }
	double sin(double a) { return ::sin(a); }
	float sin(float a) { return ::sinf(a); }
	double cos(double a) { return ::cos(a); }
	float cos(float a) { return ::cosf(a); }
	double exp(double a) { return ::exp(a); }
	float exp(float a) { return ::expf(a); }
	double pow(double a, double b) { return ::pow(a, b); }
	float pow(float a, float b) { return ::powf(a, b); }
	}

	/* ----------------------------------------------------------------------
	Scalar AIREBO implementation, standalone, with massive code reuse
	compared to original code.
	---------------------------------------------------------------------- */

	#define M_PI 3.14159265358979323846 /* pi */

	#define CARBON 0
	#define HYDROGEN 1
	#define TOL 1.0e-9

	template<typename T>
	inline T fmin_nonan(T a, T b) {
	return a < b ? a : b;
	}
	template<typename T>
	inline T fmax_nonan(T a, T b) {
	return a > b ? a : b;
	}

	template<typename flt_t>
	inline flt_t Sp(flt_t r, flt_t lo, flt_t hi, flt_t * del) {
	flt_t t = (r - lo) / (hi - lo);
	if (t <= 0) {
	if (del) *del = 0;
	return 1;
	} else if (t >= 1) {
	if (del) *del = 0;
	return 0;
	} else {
	t *= static_cast<flt_t>(M_PI);
	if (del) del = static_cast<flt_t>(-0.5 M_PI)
	* overloaded::sin(t) / (hi - lo);
	return static_cast<flt_t>(0.5) * (1 + overloaded::cos(t));
	}
	}

	template<typename flt_t>
	inline flt_t Sp2(flt_t r, flt_t lo, flt_t hi, flt_t * del) {
	flt_t t = (r - lo) / (hi - lo);
	if (t <= 0) {
	if (del) *del = 0;
	return 1;
	} else if (t >= 1) {
	if (del) *del = 0;
	return 0;
	} else {
	if (del) del = 6 (t * t - t) / (hi - lo);
	return 1 - t * t * (3 - 2 * t);
	}
	}

	template<typename flt_t>
	inline flt_t eval_poly_lin(int n, flt_t * coeffs, flt_t x, flt_t * deriv) {
	flt_t result = coeffs[n - 1];
	deriv = coeffs[n - 1] (n - 1);
	for (int i = n - 2; i > 0; i--) {
	result = coeffs[i] + x * result;
	deriv = coeffs[i] i + x * (*deriv);
	}
	result = coeffs[0] + x * result;
	return result;
	}

	template<typename flt_t, typename acc_t>
	inline flt_t gSpline(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype, flt_t cos, flt_t N, flt_t * dgdc, flt_t * dgdN) {
	flt_t NCmin = ka->params.NCmin;
	flt_t NCmax = ka->params.NCmax;
	int index = 0;
	flt_t * gDom = NULL;
	int nDom = 0;
	int offs = 0;
	if (itype == 0) {
	nDom = 4;
	gDom = &ka->params.gCdom[0];
	if (N > NCmin) offs = 4 * 6;
	} else {
	nDom = 3;
	gDom = &ka->params.gHdom[0];
	offs = 8 * 6;
	}
	cos = fmax_nonan(gDom[0], fmin_nonan(gDom[nDom], cos));
	int i;
	for (i = 0; i < nDom; i++) {
	if (cos >= gDom[i] && cos <= gDom[i + 1]) {
	index = i;
	}
	}
	flt_t g = eval_poly_lin(6, &ka->params.gVal[offs+index*6], cos, dgdc);
	*dgdN = 0;
	if (itype == 0 && N > NCmin && N < NCmax) {
	flt_t dg1;
	flt_t g1 = eval_poly_lin(6, &ka->params.gVal[index*6], cos, &dg1);
	flt_t dS;
	flt_t cut = Sp(N, NCmin, NCmax, &dS);
	dgdN = dS (g1 - g);
	g = g + cut * (g1 - g);
	dgdc = dgdc + cut * (dg1 - *dgdc);
	}
	return g;
	}

	template<typename flt_t>
	inline flt_t eval_poly_bi(int n, flt_t * coeffs, flt_t x, flt_t y,
	flt_t * deriv) {
	flt_t dy;
	flt_t vy = eval_poly_lin(n, &coeffs[n * (n - 1)], y, &dy);
	flt_t result = vy;
	deriv[0] = vy * (n - 1);
	deriv[1] = dy;
	for (int i = n - 2; i > 0; i--) {
	vy = eval_poly_lin(n, &coeffs[n * i], y, &dy);
	result = vy + x * result;
	deriv[0] = vy * i + x * deriv[0];
	deriv[1] = dy + x * deriv[1];
	}
	result = eval_poly_lin(n, &coeffs[0], y, &dy) + x * result;
	deriv[1] = dy + x * deriv[1];
	return result;
	}

	template<typename flt_t>
	inline flt_t eval_poly_tri(int n, flt_t * coeffs, flt_t x, flt_t y, flt_t z,
	flt_t * deriv) {
	flt_t dyz[2];
	flt_t vyz = eval_poly_bi(n, &coeffs[n * n * (n - 1)], y, z, &dyz[0]);
	flt_t result = vyz;
	deriv[0] = vyz * (n - 1);
	deriv[1] = dyz[0];
	deriv[2] = dyz[1];
	for (int i = n - 2; i > 0; i--) {
	vyz = eval_poly_bi(n, &coeffs[n * n * i], y, z, &dyz[0]);
	result = vyz + x * result;
	deriv[0] = vyz * i + x * deriv[0];
	deriv[1] = dyz[0] + x * deriv[1];
	deriv[2] = dyz[1] + x * deriv[2];
	}
	result = eval_poly_bi(n, &coeffs[0], y, z, &dyz[0]) + x * result;
	deriv[1] = dyz[0] + x * deriv[1];
	deriv[2] = dyz[1] + x * deriv[2];
	return result;
	}

	template<typename flt_t, typename acc_t>
	inline flt_t PijSpline(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, flt_t NC, flt_t NH, flt_t * dN) {
	dN[0] = 0.0;
	dN[1] = 0.0;
	if (itype == HYDROGEN) return 0;
	flt_t *pCJdom = jtype == CARBON ? &ka->params.pCCdom[0][0] :
	&ka->params.pCHdom[0][0];
	NC = fmax_nonan(pCJdom[0], fmin_nonan(pCJdom[1], NC));
	NH = fmax_nonan(pCJdom[2], fmin_nonan(pCJdom[3], NH));
	int nC = floor(NC);
	int nH = floor(NH);
	#define PijSelect(a, b) (jtype == CARBON ? ka->params.a : ka->params.b)
	if (fabs(NC - nC) < TOL && fabs(NH - nH) < TOL) {
	dN[0] = PijSelect(PCCdfdx, PCHdfdx)[nC][nH];
	dN[1] = PijSelect(PCCdfdy, PCHdfdy)[nC][nH];
	return PijSelect(PCCf, PCHf)[nC][nH];
	}
	if (NC == pCJdom[1]) nC -= 1;
	if (NH == pCJdom[3]) nH -= 1;
	return eval_poly_bi(4, &PijSelect(pCC, pCH)[nC][nH][0], NC, NH, dN);
	#undef PijSelect
	}

	template<typename flt_t, typename acc_t>
	inline flt_t TijSpline(KernelArgsAIREBOT<flt_t,acc_t> * ka, flt_t Nij,
	flt_t Nji, flt_t Nijconj, acc_t * dN3) {
	flt_t * Tijdom = &ka->params.Tijdom[0][0];
	Nij = fmax_nonan(Tijdom[0], fmin_nonan(Tijdom[1], Nij));
	Nji = fmax_nonan(Tijdom[2], fmin_nonan(Tijdom[3], Nji));
	Nijconj = fmax_nonan(Tijdom[4], fmin_nonan(Tijdom[5], Nijconj));
	int nij = floor(Nij);
	int nji = floor(Nji);
	int nijconj = floor(Nijconj);
	if (fabs(Nij - nij) < TOL && fabs(Nji - nji) <
	TOL && fabs(Nijconj - nijconj) < TOL) {
	dN3[0] = ka->params.Tdfdx[nij][nji][nijconj];
	dN3[1] = ka->params.Tdfdy[nij][nji][nijconj];
	dN3[2] = ka->params.Tdfdz[nij][nji][nijconj];
	return ka->params.Tf[nij][nji][nijconj];
	}
	if (Nij == Tijdom[1]) nij -= 1;
	if (Nji == Tijdom[3]) nji -= 1;
	if (Nijconj == Tijdom[5]) nijconj -= 1;
	return eval_poly_tri<acc_t>(4, &ka->params.Tijc[nij][nji][nijconj][0], Nij,
	Nji, Nijconj, dN3);
	}

	template<typename flt_t, typename acc_t>
	inline flt_t piRCSpline(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, flt_t Nij, flt_t Nji, flt_t Nijconj, acc_t * dN3) {
	const int HH = 2;
	const int CH = 1;
	/* const int CC = 0; */
	int select = itype + jtype;
	#define piRCSelect(a, b, c) (select == HH ? ka->params.a : select == CH ? \
	ka->params.b : ka->params.c)
	flt_t * piIJdom = &piRCSelect(piHHdom, piCHdom, piCCdom)[0][0];
	if (select == HH) {
	if (Nij < piIJdom[0] \|\| Nij > piIJdom[1] \|\| Nji < piIJdom[2] \|\|
	Nji > piIJdom[3] \|\| Nijconj < piIJdom[4] \|\| Nijconj > piIJdom[5]) {
	Nij = 0;
	Nji = 0;
	Nijconj = 0;
	}
	}
	Nij = fmax_nonan(piIJdom[0], fmin_nonan(piIJdom[1], Nij));
	Nji = fmax_nonan(piIJdom[2], fmin_nonan(piIJdom[3], Nji));
	Nijconj = fmax_nonan(piIJdom[4], fmin_nonan(piIJdom[5], Nijconj));
	int nij = floor(Nij);
	int nji = floor(Nji);
	int nijconj = floor(Nijconj);
	if (fabs(Nij - nij) < TOL && fabs(Nji - nji) <
	TOL && fabs(Nijconj - nijconj) < TOL) {
	dN3[0] = piRCSelect(piHHdfdx, piCHdfdx, piCCdfdx)[nij][nji][nijconj];
	dN3[1] = piRCSelect(piHHdfdy, piCHdfdy, piCCdfdy)[nij][nji][nijconj];
	dN3[2] = piRCSelect(piHHdfdz, piCHdfdz, piCCdfdz)[nij][nji][nijconj];
	return piRCSelect(piHHf, piCHf, piCCf)[nij][nji][nijconj];
	}
	if (Nij == piIJdom[1]) nij -= 1;
	if (Nji == piIJdom[3]) nji -= 1;
	if (Nijconj == piIJdom[5]) nijconj -= 1;
	return eval_poly_tri<acc_t>(4,
	&piRCSelect(piHH, piCH, piCC)[nij][nji][nijconj][0], Nij, Nji, Nijconj,
	dN3);
	#undef piRCSelect
	}

	/*
	* Implements the p_ij term in airebo, which occurs on 4 different occasions
	* in the original lammps code.
	*/
	template<typename flt_t, typename acc_t>
	inline flt_t frebo_pij(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i, int j,
	flt_t rijx, flt_t rijy, flt_t rijz, flt_t rijmag, flt_t wij, flt_t VA,
	flt_t * sum_N, acc_t fij[3]) {
	ResultForceT<acc_t> * result_f = ka->result_f;
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	flt_t * nC = ka->nC;
	flt_t * nH = ka->nH;
	flt_t x_i = x[i].x;
	flt_t y_i = x[i].y;
	flt_t z_i = x[i].z;
	int itype = map[x[i].w];
	int jtype = map[x[j].w];
	flt_t invrijm = 1 / rijmag;
	flt_t invrijm2 = invrijm * invrijm;
	flt_t rcminij = ka->params.rcmin[itype][jtype];
	flt_t rcmaxij = ka->params.rcmax[itype][jtype];
	flt_t Nmin = ka->params.Nmin;
	flt_t Nmax = ka->params.Nmax;
	flt_t Nij = nC[i] + nH[i] - wij;
	flt_t NijC = nC[i] - wij * (1 - jtype);
	flt_t NijH = nH[i] - wij * jtype;
	flt_t sum_pij = 0;
	flt_t sum_dpij_dN = 0;
	flt_t dN2[2] = {0};
	flt_t pij = 0;
	*sum_N = 0;
	int * neighs = ka->neigh_rebo.entries + ka->neigh_rebo.offset[i];
	int pass;
	for (pass = 0; pass < 2; pass++) {
	int kk;
	int knum = ka->neigh_rebo.num[i];
	for (kk = 0; kk < knum; kk++) {
	int k = neighs[kk];
	if (k == j) continue;
	flt_t rikx = x_i - x[k].x;
	flt_t riky = y_i - x[k].y;
	flt_t rikz = z_i - x[k].z;
	int ktype = map[x[k].w];
	flt_t rikmag = overloaded::sqrt(rikx * rikx + riky * riky + rikz * rikz);
	flt_t rho_k = ka->params.rho[ktype][1];
	flt_t rho_j = ka->params.rho[jtype][1];
	flt_t lamdajik = 4 * itype * ((rho_k - rikmag) - (rho_j - rijmag));
	flt_t ex_lam = exp(lamdajik);
	flt_t rcminik = ka->params.rcmin[itype][ktype];
	flt_t rcmaxik = ka->params.rcmax[itype][ktype];
	flt_t dwik;
	flt_t wik = Sp(rikmag, rcminik, rcmaxik, &dwik);
	flt_t Nki = nC[k] + nH[k] - wik;
	flt_t cosjik = (rijx * rikx + rijy * riky + rijz * rikz) /
	(rijmag * rikmag);
	cosjik = fmin_nonan<flt_t>(1, fmax_nonan<flt_t>(-1, cosjik));
	flt_t dgdc, dgdN;
	flt_t g = gSpline(ka, itype, cosjik, Nij, &dgdc, &dgdN);
	if (pass == 0) {
	sum_pij += wik * g * ex_lam;
	sum_dpij_dN += wik * dgdN * ex_lam;
	flt_t cutN = Sp<flt_t>(Nki, Nmin, Nmax, NULL);
	sum_N += (1 - ktype) wik * cutN;
	} else {
	flt_t tmp = -0.5 * pij * pij * pij;
	flt_t invrikm = 1 / rikmag;
	flt_t rjkx = rikx - rijx;
	flt_t rjky = riky - rijy;
	flt_t rjkz = rikz - rijz;
	flt_t rjkmag = sqrt(rjkx * rjkx + rjky * rjky + rjkz * rjkz);
	flt_t rijrik = 2 * rijmag * rikmag;
	flt_t rr = rijmag * rijmag - rikmag * rikmag;
	flt_t dctdjk = -2 / rijrik;
	flt_t dctdik = (-rr + rjkmag * rjkmag) / (rijrik * rikmag * rikmag);
	flt_t dctdij = (rr + rjkmag * rjkmag) / (rijrik * rijmag * rijmag);

	acc_t fi[3], fj[3], fk[3];
	flt_t pref = 0.5 * VA * tmp;
	flt_t tmp20 = pref * wik * dgdc * ex_lam;
	fj[0] = fj[1] = fj[2] = 0;
	fi[0] = -tmp20 * dctdik * rikx;
	fi[1] = -tmp20 * dctdik * riky;
	fi[2] = -tmp20 * dctdik * rikz;
	fk[0] = tmp20 * dctdik * rikx;
	fk[1] = tmp20 * dctdik * riky;
	fk[2] = tmp20 * dctdik * rikz;

	fij[0] += -tmp20 * dctdij * rijx;
	fij[1] += -tmp20 * dctdij * rijy;
	fij[2] += -tmp20 * dctdij * rijz;

	fi[0] += -tmp20 * dctdjk * rjkx;
	fi[1] += -tmp20 * dctdjk * rjky;
	fi[2] += -tmp20 * dctdjk * rjkz;
	fk[0] += tmp20 * dctdjk * rjkx;
	fk[1] += tmp20 * dctdjk * rjky;
	fk[2] += tmp20 * dctdjk * rjkz;
	fij[0] -= -tmp20 * dctdjk * rjkx;
	fij[1] -= -tmp20 * dctdjk * rjky;
	fij[2] -= -tmp20 * dctdjk * rjkz;

	flt_t tmp21 = pref * (wik * g * ex_lam * 4 * itype);
	fij[0] -= 1 * tmp21 * rijx * invrijm;
	fij[1] -= 1 * tmp21 * rijy * invrijm;
	fij[2] -= 1 * tmp21 * rijz * invrijm;
	fi[0] -= tmp21 * (-rikx * invrikm);
	fi[1] -= tmp21 * (-riky * invrikm);
	fi[2] -= tmp21 * (-rikz * invrikm);
	fk[0] -= tmp21 * (rikx * invrikm);
	fk[1] -= tmp21 * (riky * invrikm);
	fk[2] -= tmp21 * (rikz * invrikm);

	// coordination forces

	// dwik forces
	flt_t tmp22 = pref * dwik * g * ex_lam * invrikm;
	fi[0] -= tmp22 * rikx;
	fi[1] -= tmp22 * riky;
	fi[2] -= tmp22 * rikz;
	fk[0] += tmp22 * rikx;
	fk[1] += tmp22 * riky;
	fk[2] += tmp22 * rikz;

	// PIJ forces
	flt_t tmp23 = pref * dN2[ktype] * dwik * invrikm;
	fi[0] -= tmp23 * rikx;
	fi[1] -= tmp23 * riky;
	fi[2] -= tmp23 * rikz;
	fk[0] += tmp23 * rikx;
	fk[1] += tmp23 * riky;
	fk[2] += tmp23 * rikz;

	// dgdN forces
	flt_t tmp24 = pref * sum_dpij_dN * dwik * invrikm;
	fi[0] -= tmp24 * rikx;
	fi[1] -= tmp24 * riky;
	fi[2] -= tmp24 * rikz;
	fk[0] += tmp24 * rikx;
	fk[1] += tmp24 * riky;
	fk[2] += tmp24 * rikz;

	result_f[i].x += fi[0];
	result_f[i].y += fi[1];
	result_f[i].z += fi[2];
	result_f[j].x += fj[0];
	result_f[j].y += fj[1];
	result_f[j].z += fj[2];
	result_f[k].x += fk[0];
	result_f[k].y += fk[1];
	result_f[k].z += fk[2];
	}
	}
	if (pass == 0) {
	flt_t PijS = PijSpline(ka, itype, jtype, NijC, NijH, dN2);
	pij = 1 / overloaded::sqrt(1 + sum_pij + PijS);
	}
	}
	return pij;
	}

	template<typename flt_t, typename acc_t>
	inline flt_t frebo_pi_rc(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, flt_t Nij, flt_t Nji, flt_t Nijconj, flt_t * dN3) {
	acc_t dN3tmp[3] = {0};
	flt_t ret = piRCSpline(ka, itype, jtype, Nij, Nji, Nijconj, dN3tmp);
	dN3[0] = dN3tmp[0];
	dN3[1] = dN3tmp[1];
	dN3[2] = dN3tmp[2];
	return ret;
	}

	template<typename flt_t, typename acc_t>
	inline flt_t frebo_Tij(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, flt_t Nij, flt_t Nji, flt_t Nijconj, flt_t * dN3) {
	dN3[0] = 0;
	dN3[1] = 0;
	dN3[2] = 0;
	if (itype == HYDROGEN \|\| jtype == HYDROGEN) return 0;
	acc_t dN3tmp[3] = {0};
	flt_t ret = TijSpline(ka, Nij, Nji, Nijconj, dN3tmp);
	dN3[0] = dN3tmp[0];
	dN3[1] = dN3tmp[1];
	dN3[2] = dN3tmp[2];
	return ret;
	}

	/*
	* Implements a scalar version of the sum cos^1(omega) term used in pi^dh_ij.
	* Occurs in both bondorder and bondorderLJ.
	*/
	template<typename flt_t, typename acc_t>
	inline flt_t frebo_sum_omega(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i, int j,
	flt_t r23x, flt_t r23y, flt_t r23z, flt_t r23mag, flt_t VA, acc_t fij[3]) {
	ResultForceT<acc_t> * result_f = ka->result_f;
	acc_t sum_omega = 0;
	int a2 = i;
	int a3 = j;
	flt_t r32x = - r23x;
	flt_t r32y = - r23y;
	flt_t r32z = - r23z;
	int * map = ka->map;
	AtomAIREBOT<flt_t> * x = ka->x;
	flt_t thmin = ka->params.thmin;
	flt_t thmax = ka->params.thmax;
	int itype = map[x[i].w];
	int jtype = map[x[j].w];
	int * neighs_i = ka->neigh_rebo.entries + ka->neigh_rebo.offset[i];
	int * neighs_j = ka->neigh_rebo.entries + ka->neigh_rebo.offset[j];
	int num_i = ka->neigh_rebo.num[i];
	int num_j = ka->neigh_rebo.num[j];
	int kk;
	for (kk = 0; kk < num_i; kk++) {
	int k = neighs_i[kk];
	if (k == j) continue;
	int a1 = k;
	int ktype = map[x[k].w];
	flt_t r21x = x[a2].x - x[a1].x;
	flt_t r21y = x[a2].y - x[a1].y;
	flt_t r21z = x[a2].z - x[a1].z;
	flt_t r21mag = overloaded::sqrt(r21x * r21x + r21y * r21y + r21z * r21z);
	flt_t cos321 = (r23x * r21x + r23y * r21y + r23z * r21z) /
	(r23mag * r21mag);
	cos321 = fmin_nonan<flt_t>(1, fmax_nonan<flt_t>(-1, cos321));
	flt_t sin321 = overloaded::sqrt(1 - cos321 * cos321);
	if (sin321 == 0) continue;
	flt_t sink2i = 1 / (sin321 * sin321);
	flt_t rik2i = 1 / (r21mag * r21mag);
	flt_t rr = r23mag * r23mag - r21mag * r21mag;
	flt_t r31x = r21x - r23x;
	flt_t r31y = r21y - r23y;
	flt_t r31z = r21z - r23z;
	flt_t r31mag2 = r31x * r31x + r31y * r31y + r31z * r31z;
	flt_t rijrik = 2 * r23mag * r21mag;
	flt_t r21mag2 = r21mag * r21mag;
	flt_t dctik = (-rr + r31mag2) / (rijrik * r21mag2);
	flt_t dctij = (rr + r31mag2) / (rijrik * r23mag * r23mag);
	flt_t dctjk = -2 / rijrik;
	flt_t rcmin21 = ka->params.rcmin [itype][ktype];
	flt_t rcmaxp21 = ka->params.rcmaxp[itype][ktype];
	flt_t dw21;
	flt_t w21 = Sp(r21mag, rcmin21, rcmaxp21, &dw21);
	// why does this additional cutoff in the cosine exist?
	// the original code by stuart answers this:
	// it avoid issues when bonds in the dihedral are linear
	// by switching the dihedral off beforehand.
	// This is the reason for both the sin == 0 checks and the
	// tspjik = Sp2(..) calls.
	// Unfortunately, this is not exactly stated in the original paper.
	// It might be similar in purpose to the H(sin - s^min) term that
	// appears in that paper, but can not be found in original REBO papers.
	flt_t dtsjik;
	flt_t tspjik = Sp2(cos321, thmin, thmax, &dtsjik);
	dtsjik = - dtsjik;
	int ll;
	for (ll = 0; ll < num_j; ll++) {
	int l = neighs_j[ll];
	if (l == i \|\| l == k) continue;
	int ltype = map[x[l].w];
	int a4 = l;
	flt_t r34x = x[a3].x - x[a4].x;
	flt_t r34y = x[a3].y - x[a4].y;
	flt_t r34z = x[a3].z - x[a4].z;
	flt_t r34mag = overloaded::sqrt(r34x * r34x + r34y * r34y + r34z * r34z);
	flt_t cos234 = (r32x * r34x + r32y * r34y + r32z * r34z) /
	(r23mag * r34mag);
	cos234 = fmin_nonan<flt_t>(1, fmax_nonan<flt_t>(-1, cos234));
	flt_t sin234 = overloaded::sqrt(1 - cos234 * cos234);
	if (sin234 == 0) continue;
	flt_t sinl2i = 1 / (sin234 * sin234);
	flt_t rjl2i = 1 / (r34mag * r34mag);

	flt_t rcminjl = ka->params.rcmin[jtype][ltype];
	flt_t rcmaxpjl = ka->params.rcmaxp[jtype][ltype];
	flt_t dw34;
	flt_t w34 = Sp(r34mag, rcminjl, rcmaxpjl, &dw34);
	flt_t rr = (r23mag * r23mag) - (r34mag * r34mag);
	flt_t r24x = r23x + r34x;
	flt_t r24y = r23y + r34y;
	flt_t r24z = r23z + r34z;
	flt_t r242 =
	(r24x * r24x) + (r24y * r24y) + (r24z * r24z);
	flt_t rijrjl = 2 * r23mag * r34mag;
	flt_t rjl2 = r34mag * r34mag;
	flt_t dctjl = (-rr + r242) / (rijrjl * rjl2);
	flt_t dctji = (rr + r242) / (rijrjl * r23mag * r23mag);
	flt_t dctil = -2 / rijrjl;
	flt_t dtsijl;
	flt_t tspijl = Sp2(cos234, thmin, thmax, &dtsijl);
	dtsijl = -dtsijl; // need minus sign
	flt_t prefactor = VA;

	flt_t cross321x = (r32y * r21z) - (r32z * r21y);
	flt_t cross321y = (r32z * r21x) - (r32x * r21z);
	flt_t cross321z = (r32x * r21y) - (r32y * r21x);
	flt_t cross234x = (r23y * r34z) - (r23z * r34y);
	flt_t cross234y = (r23z * r34x) - (r23x * r34z);
	flt_t cross234z = (r23x * r34y) - (r23y * r34x);

	flt_t cwnum = (cross321x * cross234x) +
	(cross321y * cross234y) +
	(cross321z * cross234z);
	flt_t cwnom = r21mag * r34mag * r23mag * r23mag * sin321 * sin234;
	flt_t om1234 = cwnum / cwnom;
	flt_t cw = om1234;
	sum_omega += ((1 - (om1234 * om1234)) * w21 * w34) *
	(1 - tspjik) * (1 - tspijl);
	if (VA == static_cast<flt_t>(0.0)) continue;

	flt_t dt1dik = (rik2i) - (dctik * sink2i * cos321);
	flt_t dt1djk = (-dctjk * sink2i * cos321);
	flt_t dt1djl = (rjl2i) - (dctjl * sinl2i * cos234);
	flt_t dt1dil = (-dctil * sinl2i * cos234);
	flt_t dt1dij = (2 / (r23mag * r23mag)) -
	(dctij * sink2i * cos321) -
	(dctji * sinl2i * cos234);

	flt_t dt2dikx = (-r23z * cross234y) + (r23y * cross234z);
	flt_t dt2diky = (-r23x * cross234z) + (r23z * cross234x);
	flt_t dt2dikz = (-r23y * cross234x) + (r23x * cross234y);

	flt_t dt2djlx = (-r23y * cross321z) + (r23z * cross321y);
	flt_t dt2djly = (-r23z * cross321x) + (r23x * cross321z);
	flt_t dt2djlz = (-r23x * cross321y) + (r23y * cross321x);

	flt_t dt2dijx = (r21z * cross234y) - (r34z * cross321y) -
	flt_t (r21y * cross234z) + (r34y * cross321z);
	flt_t dt2dijy = (r21x * cross234z) - (r34x * cross321z) -
	flt_t (r21z * cross234x) + (r34z * cross321x);
	flt_t dt2dijz = (r21y * cross234x) - (r34y * cross321x) -
	flt_t (r21x * cross234y) + (r34x * cross321y);

	flt_t aa = (prefactor * 2 * cw / cwnom) * w21 * w34 *
	(1 - tspjik) * (1 - tspijl);
	flt_t aaa1 = -prefactor * (1 - (om1234 * om1234)) *
	(1 - tspjik) * (1 - tspijl);
	flt_t aaa2 = -prefactor * (1 - (om1234 * om1234)) * w21 * w34;
	flt_t at2 = aa * cwnum;

	flt_t fcijpc = (-dt1dij * at2) +
	(aaa2 * dtsjik * dctij * (1 - tspijl)) +
	(aaa2 * dtsijl * dctji * (1 - tspjik));
	flt_t fcikpc = (-dt1dik * at2) +
	(aaa2 * dtsjik * dctik * (1 - tspijl));
	flt_t fcjlpc = (-dt1djl * at2) +
	(aaa2 * dtsijl * dctjl * (1 - tspjik));
	flt_t fcjkpc = (-dt1djk * at2) +
	(aaa2 * dtsjik * dctjk * (1 - tspijl));
	flt_t fcilpc = (-dt1dil * at2) +
	(aaa2 * dtsijl * dctil * (1 - tspjik));

	flt_t F23x = (fcijpc * r23x) + (aa * dt2dijx);
	flt_t F23y = (fcijpc * r23y) + (aa * dt2dijy);
	flt_t F23z = (fcijpc * r23z) + (aa * dt2dijz);

	flt_t F12x = (fcikpc * r21x) + (aa * dt2dikx);
	flt_t F12y = (fcikpc * r21y) + (aa * dt2diky);
	flt_t F12z = (fcikpc * r21z) + (aa * dt2dikz);

	flt_t F34x = (fcjlpc * r34x) + (aa * dt2djlx);
	flt_t F34y = (fcjlpc * r34y) + (aa * dt2djly);
	flt_t F34z = (fcjlpc * r34z) + (aa * dt2djlz);

	flt_t F31x = (fcjkpc * r31x);
	flt_t F31y = (fcjkpc * r31y);
	flt_t F31z = (fcjkpc * r31z);

	flt_t F24x = (fcilpc * r24x);
	flt_t F24y = (fcilpc * r24y);
	flt_t F24z = (fcilpc * r24z);

	flt_t f1x = -F12x - F31x;
	flt_t f1y = -F12y - F31y;
	flt_t f1z = -F12z - F31z;
	flt_t f2x = F12x + F31x;
	flt_t f2y = F12y + F31y;
	flt_t f2z = F12z + F31z;
	flt_t f3x = F34x + F24x;
	flt_t f3y = F34y + F24y;
	flt_t f3z = F34z + F24z;
	flt_t f4x = -F34x - F24x;
	flt_t f4y = -F34y - F24y;
	flt_t f4z = -F34z - F24z;

	fij[0] += F23x + F24x - F31x;
	fij[1] += F23y + F24y - F31y;
	fij[2] += F23z + F24z - F31z;

	// coordination forces

	flt_t tmp20 = VA * ((1 - (om1234 * om1234))) *
	(1 - tspjik) * (1 - tspijl) * dw21 * w34 / r21mag;
	f2x -= tmp20 * r21x;
	f2y -= tmp20 * r21y;
	f2z -= tmp20 * r21z;
	f1x += tmp20 * r21x;
	f1y += tmp20 * r21y;
	f1z += tmp20 * r21z;

	flt_t tmp21 = VA * ((1 - (om1234 * om1234))) *
	(1 - tspjik) * (1 - tspijl) * w21 * dw34 / r34mag;
	f3x -= tmp21 * r34x;
	f3y -= tmp21 * r34y;
	f3z -= tmp21 * r34z;
	f4x += tmp21 * r34x;
	f4y += tmp21 * r34y;
	f4z += tmp21 * r34z;

	result_f[a1].x += f1x;
	result_f[a1].y += f1y;
	result_f[a1].z += f1z;
	result_f[a2].x += f2x;
	result_f[a2].y += f2y;
	result_f[a2].z += f2z;
	result_f[a3].x += f3x;
	result_f[a3].y += f3y;
	result_f[a3].z += f3z;
	result_f[a4].x += f4x;
	result_f[a4].y += f4y;
	result_f[a4].z += f4z;
	}
	}
	return sum_omega;
	}

	/*
	* Implements a scalar implementation the force update due to splines.
	* It is used for both pi^rc_ij and T_ij.
	* Occurs four times in each bondorder and bondorderLJ.
	*/
	template<typename flt_t, typename acc_t>
	inline void frebo_N_spline_force(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i,
	int j, flt_t VA, flt_t dN, flt_t dNconj, flt_t Nconj) {
	int * map = ka->map;
	AtomAIREBOT<flt_t> * x = ka->x;
	ResultForceT<acc_t> * result_f = ka->result_f;
	flt_t * nC = ka->nC;
	flt_t * nH = ka->nH;
	flt_t Nmin = ka->params.Nmin;
	flt_t Nmax = ka->params.Nmax;
	int itype = map[x[i].w];
	int * neighs = ka->neigh_rebo.entries + ka->neigh_rebo.offset[i];
	int knum = ka->neigh_rebo.num[i];
	int kk;
	for (kk = 0; kk < knum; kk++) {
	int k = neighs[kk];
	if (k == j) continue;
	flt_t rikx = x[i].x - x[k].x;
	flt_t riky = x[i].y - x[k].y;
	flt_t rikz = x[i].z - x[k].z;
	flt_t rikmag = overloaded::sqrt(rikx * rikx + riky * riky + rikz * rikz);
	int ktype = map[x[k].w];
	flt_t rcminik = ka->params.rcmin[itype][ktype];
	flt_t rcmaxik = ka->params.rcmax[itype][ktype];
	flt_t dwik;
	flt_t wik = Sp(rikmag, rcminik, rcmaxik, &dwik);
	flt_t Nki = nC[k] + nH[k] - wik;
	flt_t dNki;
	flt_t SpN = Sp(Nki, Nmin, Nmax, &dNki);
	flt_t fdN = VA * dN * dwik / rikmag;
	flt_t fdNconj = VA * dNconj * 2 * Nconj * dwik * SpN / rikmag;
	flt_t ffactor = fdN;
	if (ktype == 0) ffactor += fdNconj;
	flt_t fkx = ffactor * rikx;
	flt_t fky = ffactor * riky;
	flt_t fkz = ffactor * rikz;
	result_f[i].x -= fkx;
	result_f[i].y -= fky;
	result_f[i].z -= fkz;
	result_f[k].x += fkx;
	result_f[k].y += fky;
	result_f[k].z += fkz;
	if (ktype != 0 \|\| fabs(dNki) <= TOL) continue;
	int * neighs_k = ka->neigh_rebo.entries + ka->neigh_rebo.offset[k];
	int nnum = ka->neigh_rebo.num[k];
	int nn;
	for (nn = 0; nn < nnum; nn++) {
	int n = neighs_k[nn];
	if (n == i) continue;
	flt_t rknx = x[k].x - x[n].x;
	flt_t rkny = x[k].y - x[n].y;
	flt_t rknz = x[k].z - x[n].z;
	flt_t rknmag = overloaded::sqrt(rknx * rknx + rkny * rkny + rknz * rknz);
	int ntype = map[x[n].w];
	flt_t rcminkn = ka->params.rcmin[ktype][ntype];
	flt_t rcmaxkn = ka->params.rcmax[ktype][ntype];
	flt_t dwkn;
	Sp(rknmag, rcminkn, rcmaxkn, &dwkn);
	flt_t ffactor = VA * dNconj * 2 * Nconj * wik * dNki * dwkn / rknmag;
	result_f[k].x -= ffactor * rknx;
	result_f[k].y -= ffactor * rkny;
	result_f[k].z -= ffactor * rknz;
	result_f[n].x += ffactor * rknx;
	result_f[n].y += ffactor * rkny;
	result_f[n].z += ffactor * rknz;
	}
	}
	}

	/*
	* This data-structure contains the result of a search through neighbor-lists.
	* It is used to calculate C_ij and the corresponding force updates.
	*/
	template<typename flt_t>
	struct LennardJonesPathAIREBOT {
	AtomAIREBOT<flt_t> del[3];
	int num;
	flt_t w[3];
	flt_t dw[3];
	flt_t r[3];
	int idx[4];
	};

	/*
	* Checks a candidate path stored in idxs whether it is better than *path
	* and updates *path accordingly.
	*/
	template<typename flt_t, typename acc_t>
	inline flt_t ref_lennard_jones_test_path_single(
	KernelArgsAIREBOT<flt_t,acc_t> * ka, flt_t best, int num, int * idxs,
	LennardJonesPathAIREBOT<flt_t> * path) {
	LennardJonesPathAIREBOT<flt_t> result;
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	result.num = num;
	flt_t combined = 1;
	for (int i = num - 2; i >= 0; i--) {
	int a0 = idxs[i+0];
	int a1 = idxs[i+1];
	flt_t delx = x[a1].x - x[a0].x;
	flt_t dely = x[a1].y - x[a0].y;
	flt_t delz = x[a1].z - x[a0].z;
	flt_t rsq = delx * delx + dely * dely + delz * delz;
	int type0 = map[x[a0].w];
	int type1 = map[x[a1].w];
	if (rsq >= ka->params.rcmaxsq[type0][type1]) return best;
	flt_t r = overloaded::sqrt(rsq);
	flt_t dw, w = Sp<flt_t>(r, ka->params.rcmin[type0][type1],
	ka->params.rcmax[type0][type1], &dw);
	if (w == 0) return best;
	combined *= w;
	if (combined <= best) return best;
	result.idx[i] = a0;
	result.del[i].x = delx;
	result.del[i].y = dely;
	result.del[i].z = delz;
	result.r[i] = r;
	result.w[i] = w;
	result.dw[i] = dw;
	}
	result.idx[num - 1] = idxs[num - 1];
	*path = result;
	return combined;
	}

	/*
	* Test through all paths surrounding i and j to find the corresponding
	* best path. Uses the same iteration ordering as FLJ() does.
	* Note that an optimization would use the j neighlist instead in the inner
	* loop.
	*/
	template<typename flt_t, typename acc_t>
	inline flt_t ref_lennard_jones_test_path(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	int i, int j, flt_t rij, flt_t rcmax,
	LennardJonesPathAIREBOT<flt_t> * path) {
	int idxs[4];
	idxs[0] = i;
	idxs[1] = j;
	flt_t best = 0;
	if (rij <= rcmax) {
	best = ref_lennard_jones_test_path_single(ka, best, 2, idxs, path);
	if (best == static_cast<flt_t>(1.0)) return 0;
	}
	for (int kk = 0; kk < ka->neigh_rebo.num[i]; kk++) {
	int k = ka->neigh_rebo.entries[ka->neigh_rebo.offset[i] + kk];
	if (k == j) continue;
	idxs[1] = k;
	idxs[2] = j;
	best = ref_lennard_jones_test_path_single(ka, best, 3, idxs, path);
	if (best == static_cast<flt_t>(1.0)) return 0;
	for (int mm = 0; mm < ka->neigh_rebo.num[k]; mm++) {
	int m = ka->neigh_rebo.entries[ka->neigh_rebo.offset[k] + mm];
	if (m == i \|\| m == j) continue;
	idxs[2] = m;
	idxs[3] = j;
	best = ref_lennard_jones_test_path_single(ka, best, 4, idxs, path);
	if (best == static_cast<flt_t>(1.0)) return 0;
	}
	}
	return 1 - best;
	}

	/*
	* Conducts the force update due to C_ij, given the active path.
	*/
	template<typename flt_t, typename acc_t>
	inline void ref_lennard_jones_force_path(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	flt_t dC, LennardJonesPathAIREBOT<flt_t> * path) {
	AtomAIREBOT<flt_t> * x = ka->x;
	ResultForceT<acc_t> * result_f = ka->result_f;
	for (int i = 0; i < path->num - 1; i++) {
	flt_t fpair = dC * path->dw[i] / path->r[i];
	for (int j = 0; j < path->num - 1; j++) {
	if (i != j) fpair *= path->w[j];
	}
	result_f[path->idx[i+0]].x -= fpair * path->del[i].x;
	result_f[path->idx[i+0]].y -= fpair * path->del[i].y;
	result_f[path->idx[i+0]].z -= fpair * path->del[i].z;
	result_f[path->idx[i+1]].x += fpair * path->del[i].x;
	result_f[path->idx[i+1]].y += fpair * path->del[i].y;
	result_f[path->idx[i+1]].z += fpair * path->del[i].z;
	}
	}

	/*
	* Calculate the bondorderLJ term.
	*/
	template<typename flt_t, typename acc_t>
	inline flt_t ref_lennard_jones_bondorder(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	int i, int j, flt_t VA, acc_t fij[3]) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	ResultForceT<acc_t> * result_f = ka->result_f;

	int itype = map[x[i].w];
	int jtype = map[x[j].w];

	flt_t delx = x[i].x - x[j].x;
	flt_t dely = x[i].y - x[j].y;
	flt_t delz = x[i].z - x[j].z;
	flt_t rsq = delx * delx + dely * dely + delz * delz;
	flt_t rij = overloaded::sqrt(rsq);

	flt_t rcminij = ka->params.rcmin[itype][jtype];
	flt_t rcmaxij = ka->params.rcmax[itype][jtype];
	flt_t dwij;
	flt_t wij = Sp(rij, rcminij, rcmaxij, &dwij);

	flt_t the_r = ka->params.rcmin[itype][jtype];
	flt_t scale = the_r / rij;
	flt_t Nij = ka->nH[i] + ka->nC[i] - wij;
	flt_t Nji = ka->nH[j] + ka->nC[j] - wij;
	flt_t NconjtmpI;
	acc_t fijc[3] = {0}, fjic[3] = {0};
	flt_t pij = frebo_pij<flt_t,acc_t>(ka, i, j, delx * scale, dely * scale,
	delz * scale, the_r, wij, 0.0, &NconjtmpI, fijc);
	flt_t NconjtmpJ;
	flt_t pji = frebo_pij<flt_t,acc_t>(ka, j, i, -delx * scale, -dely * scale,
	-delz * scale, the_r, wij, 0.0, &NconjtmpJ, fjic);
	flt_t Nijconj = 1.0 + (NconjtmpI * NconjtmpI) + (NconjtmpJ * NconjtmpJ);
	flt_t dN3_pi_rc[3];
	flt_t pi_rc = frebo_pi_rc<flt_t,acc_t>(ka, itype, jtype, Nij, Nji, Nijconj,
	dN3_pi_rc);
	flt_t dN3_Tij[3];
	flt_t Tij = frebo_Tij<flt_t,acc_t>(ka, itype, jtype, Nij, Nji, Nijconj,
	dN3_Tij);
	flt_t sum_omega = 0;
	if (fabs(Tij) > TOL) {
	sum_omega = frebo_sum_omega<flt_t,acc_t>(ka, i, j, delx * scale, dely *
	scale, delz * scale, the_r, 0.0,
	fijc);
	}
	flt_t pi_dh = Tij * sum_omega;
	flt_t bij = 0.5 * (pij + pji) + pi_rc + pi_dh;
	flt_t dStb;
	flt_t Stb = Sp2<flt_t>(bij, ka->params.bLJmin[itype][jtype],
	ka->params.bLJmax[itype][jtype], &dStb);
	if (dStb != 0) {
	flt_t pij_reverse = frebo_pij<flt_t,acc_t>(ka, i, j, delx * scale,
	dely * scale, delz * scale, the_r, wij, VA * dStb, &NconjtmpI, fijc);
	flt_t pji_reverse = frebo_pij<flt_t,acc_t>(ka, j, i, -delx * scale,
	-dely * scale, -delz * scale, the_r, wij, VA * dStb, &NconjtmpJ, fjic);
	fijc[0] -= fjic[0];
	fijc[1] -= fjic[1];
	fijc[2] -= fjic[2];
	frebo_N_spline_force<flt_t,acc_t>(ka, i, j, VA * dStb, dN3_pi_rc[0],
	dN3_pi_rc[2], NconjtmpI);
	frebo_N_spline_force<flt_t,acc_t>(ka, j, i, VA * dStb, dN3_pi_rc[1],
	dN3_pi_rc[2], NconjtmpJ);
	if (fabs(Tij) > TOL) {
	flt_t sum_omega_reverse = frebo_sum_omega<flt_t,acc_t>(ka, i, j,
	delx * scale, dely * scale, delz * scale, the_r, VA * dStb * Tij, fijc);
	frebo_N_spline_force(ka, i, j, VA * dStb * sum_omega, dN3_Tij[0],
	dN3_Tij[2], NconjtmpI);
	frebo_N_spline_force(ka, j, i, VA * dStb * sum_omega, dN3_Tij[1],
	dN3_Tij[2], NconjtmpJ);
	}
	assert(fij[0] == 0);
	assert(fij[1] == 0);
	assert(fij[2] == 0);
	fij[0] = scale * (fijc[0] - (delx * delx * fijc[0] + dely * delx *
	fijc[1] + delz * delx * fijc[2]) / rsq);
	fij[1] = scale * (fijc[1] - (delx * dely * fijc[0] + dely * dely *
	fijc[1] + delz * dely * fijc[2]) / rsq);
	fij[2] = scale * (fijc[2] - (delx * delz * fijc[0] + dely * delz *
	fijc[1] + delz * delz * fijc[2]) / rsq);
	}
	return Stb;
	}

	/*
	* Scalar reference implementation of neighbor routine.
	*/
	template<typename flt_t, typename acc_t>
	void ref_rebo_neigh(KernelArgsAIREBOT<flt_t,acc_t> * ka) {
	int offset = ka->neigh_from_atom * ka->num_neighs_per_atom;
	for (int i = ka->neigh_from_atom; i < ka->neigh_to_atom; i++) {
	ka->neigh_rebo.offset[i] = offset;
	int itype = ka->map[ka->x[i].w];
	int n = 0;
	ka->nC[i] = 0;
	ka->nH[i] = 0;
	for (int j = 0; j < ka->neigh_lmp.num[i]; j++) {
	int ji = ka->neigh_lmp.entries[ka->neigh_lmp.offset[i] + j];
	flt_t delx = ka->x[i].x - ka->x[ji].x;
	flt_t dely = ka->x[i].y - ka->x[ji].y;
	flt_t delz = ka->x[i].z - ka->x[ji].z;
	flt_t rsq = delx * delx + dely * dely + delz * delz;
	int jtype = ka->map[ka->x[ji].w];
	if (rsq < ka->params.rcmaxsq[itype][jtype]) {
	ka->neigh_rebo.entries[offset + n++] = ji;
	flt_t rcmin = ka->params.rcmin[itype][jtype];
	flt_t rcmax = ka->params.rcmax[itype][jtype];
	if (jtype == CARBON)
	ka->nC[i] += Sp<flt_t>(overloaded::sqrt(rsq), rcmin, rcmax, NULL);
	else
	ka->nH[i] += Sp<flt_t>(overloaded::sqrt(rsq), rcmin, rcmax, NULL);
	}
	}
	ka->neigh_rebo.num[i] = n;
	offset += n;
	}
	}

	template<typename flt_t, typename acc_t>
	void ref_torsion_single_interaction(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i,
	int j) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	ResultForceT<acc_t> * f = ka->result_f;
	flt_t (*rcmin)[2] = ka->params.rcmin;
	flt_t (*rcmax)[2] = ka->params.rcmax;
	flt_t (*epsilonT)[2] = ka->params.epsilonT;
	flt_t thmin = ka->params.thmin;
	flt_t thmax = ka->params.thmax;
	int itype = map[x[i].w];
	flt_t xtmp = x[i].x;
	flt_t ytmp = x[i].y;
	flt_t ztmp = x[i].z;
	int * REBO_neighs_i = &ka->neigh_rebo.entries[ka->neigh_rebo.offset[i]];
	int jnum = ka->neigh_rebo.num[i];
	int jtype = map[x[j].w];

	flt_t del32x = x[j].x-x[i].x;
	flt_t del32y = x[j].y-x[i].y;
	flt_t del32z = x[j].z-x[i].z;
	flt_t rsq = del32xdel32x + del32ydel32y + del32z*del32z;
	flt_t r32 = overloaded::sqrt(rsq);
	flt_t del23x = -del32x;
	flt_t del23y = -del32y;
	flt_t del23z = -del32z;
	flt_t r23 = r32;
	flt_t dw23, w23 = Sp<flt_t>(r23,rcmin[itype][jtype],rcmax[itype][jtype],
	&dw23);

	assert(itype == 0);
	assert(jtype == 0);

	for (int kk = 0; kk < jnum; kk++) {
	int k = REBO_neighs_i[kk];
	int ktype = map[x[k].w];
	if (k == j) continue;
	flt_t del21x = x[i].x-x[k].x;
	flt_t del21y = x[i].y-x[k].y;
	flt_t del21z = x[i].z-x[k].z;
	flt_t rsq = del21xdel21x + del21ydel21y + del21z*del21z;
	flt_t r21 = overloaded::sqrt(rsq);
	flt_t cos321 = - ((del21xdel32x) + (del21ydel32y) +
	(del21zdel32z)) / (r21r32);
	cos321 = fmin(cos321,1);
	cos321 = fmax(cos321,-1);
	flt_t sin321 = overloaded::sqrt(1 - cos321*cos321);
	if (sin321 < TOL) continue;

	flt_t deljkx = del21x-del23x;
	flt_t deljky = del21y-del23y;
	flt_t deljkz = del21z-del23z;
	flt_t rjk2 = deljkxdeljkx + deljkydeljky + deljkz*deljkz;
	flt_t rjk = overloaded::sqrt(rjk2);
	flt_t rik2 = r21*r21;
	flt_t dw21, w21 = Sp<flt_t>(r21,rcmin[itype][ktype],rcmax[itype][ktype],
	&dw21);

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

	int * REBO_neighs_j = &ka->neigh_rebo.entries[ka->neigh_rebo.offset[j]];
	int lnum = ka->neigh_rebo.num[j];
	for (int ll = 0; ll < lnum; ll++) {
	int l = REBO_neighs_j[ll];
	int ltype = map[x[l].w];
	if (l == i \|\| l == k) continue;
	flt_t del34x = x[j].x-x[l].x;
	flt_t del34y = x[j].y-x[l].y;
	flt_t del34z = x[j].z-x[l].z;
	flt_t rsq = del34xdel34x + del34ydel34y + del34z*del34z;
	flt_t r34 = overloaded::sqrt(rsq);
	flt_t cos234 = (del32xdel34x + del32ydel34y +
	del32zdel34z) / (r32r34);
	cos234 = fmin(cos234,1);
	cos234 = fmax(cos234,-1);
	flt_t sin234 = overloaded::sqrt(1 - cos234*cos234);
	if (sin234 < TOL) continue;
	flt_t dw34, w34 = Sp<flt_t>(r34,rcmin[jtype][ltype],rcmax[jtype][ltype],
	&dw34);
	flt_t delilx = del23x + del34x;
	flt_t delily = del23y + del34y;
	flt_t delilz = del23z + del34z;
	flt_t ril2 = delilxdelilx + delilydelily + delilz*delilz;
	flt_t ril = overloaded::sqrt(ril2);
	flt_t rjl2 = r34*r34;

	flt_t rjl = r34;
	flt_t costmp = static_cast<flt_t>(0.5)*(rij2+rjl2-ril2)/rij/rjl;
	flt_t dtsijl, tspijl = Sp2<flt_t>(costmp,thmin,thmax,&dtsijl);
	dtsijl = -dtsijl; //need minus sign
	flt_t cross321x = (del32ydel21z)-(del32zdel21y);
	flt_t cross321y = (del32zdel21x)-(del32xdel21z);
	flt_t cross321z = (del32xdel21y)-(del32ydel21x);
	flt_t cross321mag = overloaded::sqrt(cross321x*cross321x+
	cross321ycross321y + cross321zcross321z);
	flt_t cross234x = (del23ydel34z)-(del23zdel34y);
	flt_t cross234y = (del23zdel34x)-(del23xdel34z);
	flt_t cross234z = (del23xdel34y)-(del23ydel34x);
	flt_t cross234mag = overloaded::sqrt(cross234x*cross234x+
	cross234ycross234y + cross234zcross234z);
	flt_t cwnum = (cross321x*cross234x) +
	(cross321ycross234y)+(cross321zcross234z);
	flt_t cwnom = r21r34r32r32sin321*sin234;
	flt_t cw = cwnum/cwnom;

	flt_t cw2 = (static_cast<flt_t>(.5)*(1-cw));
	flt_t ekijl = epsilonT[ktype][ltype];
	flt_t Ec = 256*ekijl/405;
	flt_t Vtors = (Ec*(overloaded::pow(cw2,5)))-(ekijl/10);

	ka->result_eng += Vtorsw21w23w34(1-tspjik)*(1-tspijl);

	flt_t dndijx = (cross234ydel21z)-(cross234zdel21y);
	flt_t dndijy = (cross234zdel21x)-(cross234xdel21z);
	flt_t dndijz = (cross234xdel21y)-(cross234ydel21x);

	flt_t tmpvecx = (del34ycross321z)-(del34zcross321y);
	flt_t tmpvecy = (del34zcross321x)-(del34xcross321z);
	flt_t tmpvecz = (del34xcross321y)-(del34ycross321x);

	dndijx = dndijx+tmpvecx;
	dndijy = dndijy+tmpvecy;
	dndijz = dndijz+tmpvecz;

	flt_t dndikx = (del23ycross234z)-(del23zcross234y);
	flt_t dndiky = (del23zcross234x)-(del23xcross234z);
	flt_t dndikz = (del23xcross234y)-(del23ycross234x);

	flt_t dndjlx = (cross321ydel23z)-(cross321zdel23y);
	flt_t dndjly = (cross321zdel23x)-(cross321xdel23z);
	flt_t dndjlz = (cross321xdel23y)-(cross321ydel23x);

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

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

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

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

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

	flt_t ddndij = (dxidijcross234mag)+(cross321magdxjdji);
	flt_t ddndik = dxidik*cross234mag;
	flt_t ddndjk = dxidjk*cross234mag;
	flt_t ddndjl = cross321mag*dxjdjl;
	flt_t ddndil = cross321mag*dxjdil;
	flt_t dcwddn = -cwnum/(cwnom*cwnom);
	flt_t dcwdn = 1/cwnom;
	flt_t dvpdcw = (-1)Ecstatic_cast<flt_t>(-0.5)5overloaded::pow(cw2,4)*
	w23w21w34(1-tspjik)(1-tspijl);

	flt_t Ftmpx = dvpdcw((dcwdndndijx)+(dcwddnddndijdel23x/r23));
	flt_t Ftmpy = dvpdcw((dcwdndndijy)+(dcwddnddndijdel23y/r23));
	flt_t Ftmpz = dvpdcw((dcwdndndijz)+(dcwddnddndijdel23z/r23));
	flt_t fix = Ftmpx;
	flt_t fiy = Ftmpy;
	flt_t fiz = Ftmpz;
	flt_t fjx = -Ftmpx;
	flt_t fjy = -Ftmpy;
	flt_t fjz = -Ftmpz;

	Ftmpx = dvpdcw((dcwdndndikx)+(dcwddnddndikdel21x/r21));
	Ftmpy = dvpdcw((dcwdndndiky)+(dcwddnddndikdel21y/r21));
	Ftmpz = dvpdcw((dcwdndndikz)+(dcwddnddndikdel21z/r21));
	fix += Ftmpx;
	fiy += Ftmpy;
	fiz += Ftmpz;
	flt_t fkx = -Ftmpx;
	flt_t fky = -Ftmpy;
	flt_t fkz = -Ftmpz;

	Ftmpx = (dvpdcwdcwddnddndjk*deljkx)/rjk;
	Ftmpy = (dvpdcwdcwddnddndjk*deljky)/rjk;
	Ftmpz = (dvpdcwdcwddnddndjk*deljkz)/rjk;
	fjx += Ftmpx;
	fjy += Ftmpy;
	fjz += Ftmpz;
	fkx -= Ftmpx;
	fky -= Ftmpy;
	fkz -= Ftmpz;

	Ftmpx = dvpdcw((dcwdndndjlx)+(dcwddnddndjldel34x/r34));
	Ftmpy = dvpdcw((dcwdndndjly)+(dcwddnddndjldel34y/r34));
	Ftmpz = dvpdcw((dcwdndndjlz)+(dcwddnddndjldel34z/r34));
	fjx += Ftmpx;
	fjy += Ftmpy;
	fjz += Ftmpz;
	flt_t flx = -Ftmpx;
	flt_t fly = -Ftmpy;
	flt_t flz = -Ftmpz;

	Ftmpx = (dvpdcwdcwddnddndil*delilx)/ril;
	Ftmpy = (dvpdcwdcwddnddndil*delily)/ril;
	Ftmpz = (dvpdcwdcwddnddndil*delilz)/ril;
	fix += Ftmpx;
	fiy += Ftmpy;
	fiz += Ftmpz;
	flx -= Ftmpx;
	fly -= Ftmpy;
	flz -= Ftmpz;

	// coordination forces

	flt_t fpair = Vtorsdw21w23w34(1-tspjik)*(1-tspijl) / r21;
	fix -= del21x*fpair;
	fiy -= del21y*fpair;
	fiz -= del21z*fpair;
	fkx += del21x*fpair;
	fky += del21y*fpair;
	fkz += del21z*fpair;

	fpair = Vtorsw21dw23w34(1-tspjik)*(1-tspijl) / r23;
	fix -= del23x*fpair;
	fiy -= del23y*fpair;
	fiz -= del23z*fpair;
	fjx += del23x*fpair;
	fjy += del23y*fpair;
	fjz += del23z*fpair;

	fpair = Vtorsw21w23dw34(1-tspjik)*(1-tspijl) / r34;
	fjx -= del34x*fpair;
	fjy -= del34y*fpair;
	fjz -= del34z*fpair;
	flx += del34x*fpair;
	fly += del34y*fpair;
	flz += del34z*fpair;

	// additional cut off function forces

	flt_t fcpc = -Vtorsw21w23w34dtsjik*(1-tspijl);
	fpair = fcpc*dcidij/rij;
	fix += fpair*del23x;
	fiy += fpair*del23y;
	fiz += fpair*del23z;
	fjx -= fpair*del23x;
	fjy -= fpair*del23y;
	fjz -= fpair*del23z;

	fpair = fcpc*dcidik/rik;
	fix += fpair*del21x;
	fiy += fpair*del21y;
	fiz += fpair*del21z;
	fkx -= fpair*del21x;
	fky -= fpair*del21y;
	fkz -= fpair*del21z;

	fpair = fcpc*dcidjk/rjk;
	fjx += fpair*deljkx;
	fjy += fpair*deljky;
	fjz += fpair*deljkz;
	fkx -= fpair*deljkx;
	fky -= fpair*deljky;
	fkz -= fpair*deljkz;

	fcpc = -Vtorsw21w23w34(1-tspjik)*dtsijl;
	fpair = fcpc*dcjdji/rij;
	fix += fpair*del23x;
	fiy += fpair*del23y;
	fiz += fpair*del23z;
	fjx -= fpair*del23x;
	fjy -= fpair*del23y;
	fjz -= fpair*del23z;

	fpair = fcpc*dcjdjl/rjl;
	fjx += fpair*del34x;
	fjy += fpair*del34y;
	fjz += fpair*del34z;
	flx -= fpair*del34x;
	fly -= fpair*del34y;
	flz -= fpair*del34z;

	fpair = fcpc*dcjdil/ril;
	fix += fpair*delilx;
	fiy += fpair*delily;
	fiz += fpair*delilz;
	flx -= fpair*delilx;
	fly -= fpair*delily;
	flz -= fpair*delilz;

	// sum per-atom forces into atom force array

	f[i].x += fix; f[i].y += fiy; f[i].z += fiz;
	f[j].x += fjx; f[j].y += fjy; f[j].z += fjz;
	f[k].x += fkx; f[k].y += fky; f[k].z += fkz;
	f[l].x += flx; f[l].y += fly; f[l].z += flz;
	}
	}
	}

	template<typename flt_t, typename acc_t>
	void ref_torsion(KernelArgsAIREBOT<flt_t,acc_t> * ka) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	int * tag = ka->tag;
	for (int ii = ka->frebo_from_atom; ii < ka->frebo_to_atom; ii++) {
	int i = ii;
	int itag = tag[i];
	int itype = map[x[i].w];
	if (itype != 0) continue;
	flt_t xtmp = x[i].x;
	flt_t ytmp = x[i].y;
	flt_t ztmp = x[i].z;
	int * REBO_neighs_i = &ka->neigh_rebo.entries[ka->neigh_rebo.offset[i]];
	int jnum = ka->neigh_rebo.num[i];
	for (int jj = 0; jj < jnum; jj++) {
	int j = REBO_neighs_i[jj];
	int jtag = tag[j];

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

	int jtype = map[x[j].w];
	if (jtype != 0) continue;
	ref_torsion_single_interaction(ka, i, j);
	}
	}
	}

	/*
	* Calculate single REBO interaction.
	* Corresponds to FREBO method. Note that the bondorder() function is
	* inlined.
	*/
	template<typename flt_t, typename acc_t>
	void ref_frebo_single_interaction(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i,
	int j) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	ResultForceT<acc_t> * result_f = ka->result_f;
	int jj;
	int itype = map[x[i].w];
	flt_t x_i = x[i].x;
	flt_t y_i = x[i].y;
	flt_t z_i = x[i].z;
	int jtype = map[x[j].w];
	flt_t delx = x[i].x - x[j].x;
	flt_t dely = x[i].y - x[j].y;
	flt_t delz = x[i].z - x[j].z;
	flt_t rsq = delx * delx + dely * dely + delz * delz;
	flt_t rij = overloaded::sqrt(rsq);
	flt_t rcminij = ka->params.rcmin[itype][jtype];
	flt_t rcmaxij = ka->params.rcmax[itype][jtype];
	flt_t dwij;
	flt_t wij = Sp(rij, rcminij, rcmaxij, &dwij);
	if (wij <= TOL) return;

	flt_t Qij = ka->params.Q[itype][jtype];
	flt_t Aij = ka->params.A[itype][jtype];
	flt_t alphaij = ka->params.alpha[itype][jtype];

	flt_t exp_alphar = exp(-alphaij * rij);
	flt_t VR_by_wij = (1.0 + (Qij / rij)) * Aij * exp_alphar;
	flt_t VR = wij * VR_by_wij;
	flt_t pre = wij * Aij * exp_alphar;
	flt_t dVRdi = pre * ((-alphaij) - (Qij / rsq) - (Qij * alphaij / rij));
	dVRdi += VR_by_wij * dwij;

	flt_t VA_by_wij = 0, dVA = 0;
	for (int k = 0; k < 3; k++) {
	flt_t BIJc = ka->params.BIJc[itype][jtype][k];
	flt_t Betaij = ka->params.Beta[itype][jtype][k];
	flt_t term = -BIJc * overloaded::exp(-Betaij * rij);
	VA_by_wij += term;
	dVA += -Betaij * wij * term;
	}
	dVA += VA_by_wij * dwij;
	flt_t VA = VA_by_wij * wij;

	acc_t fij[3] = {0};
	flt_t Nij = ka->nH[i] + ka->nC[i] - wij;
	flt_t Nji = ka->nH[j] + ka->nC[j] - wij;
	flt_t NconjtmpI;
	flt_t pij = frebo_pij(ka, i, j, delx, dely, delz, rij, wij, VA, &NconjtmpI,
	fij);
	flt_t NconjtmpJ;
	acc_t fji[3] = {0};
	flt_t pji = frebo_pij(ka, j, i, -delx, -dely, -delz, rij, wij, VA,
	&NconjtmpJ, fji);
	fij[0] -= fji[0]; fij[1] -= fji[1]; fij[2] -= fji[2];
	flt_t Nijconj = 1.0 + (NconjtmpI * NconjtmpI) + (NconjtmpJ * NconjtmpJ);
	flt_t dN3[3];
	flt_t pi_rc = frebo_pi_rc(ka, itype, jtype, Nij, Nji, Nijconj, dN3);
	frebo_N_spline_force(ka, i, j, VA, dN3[0], dN3[2], NconjtmpI);
	frebo_N_spline_force(ka, j, i, VA, dN3[1], dN3[2], NconjtmpJ);
	flt_t Tij = frebo_Tij(ka, itype, jtype, Nij, Nji, Nijconj, dN3);
	flt_t sum_omega = 0.0;
	if (fabs(Tij) > TOL) {
	sum_omega = frebo_sum_omega(ka, i, j, delx, dely, delz, rij, VA * Tij, fij);
	frebo_N_spline_force(ka, i, j, VA * sum_omega, dN3[0], dN3[2], NconjtmpI);
	frebo_N_spline_force(ka, j, i, VA * sum_omega, dN3[1], dN3[2], NconjtmpJ);
	}
	flt_t pi_dh = Tij * sum_omega;
	flt_t bij = static_cast<flt_t>(0.5) * (pij + pji) + pi_rc + pi_dh;
	flt_t dVAdi = bij * dVA;
	flt_t fpair = -(dVRdi + dVAdi) / rij;

	result_f[i].x += fpair * delx + fij[0];
	result_f[i].y += fpair * dely + fij[1];
	result_f[i].z += fpair * delz + fij[2];
	result_f[j].x -= fpair * delx + fij[0];
	result_f[j].y -= fpair * dely + fij[1];
	result_f[j].z -= fpair * delz + fij[2];

	flt_t evdwl = VR + bij * VA;
	ka->result_eng += evdwl;
	result_f[i].w += 0.5 * evdwl;
	result_f[j].w += 0.5 * evdwl;
	}


	template<typename flt_t, typename acc_t>
	inline void ref_frebo_single_atom(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * tag = ka->tag;
	int jj;
	int itag = tag[i];
	flt_t x_i = x[i].x;
	flt_t y_i = x[i].y;
	flt_t z_i = x[i].z;
	int * neighs = ka->neigh_rebo.entries + ka->neigh_rebo.offset[i];
	int jnum = ka->neigh_rebo.num[i];
	for (jj = 0; jj < jnum; jj++) {
	int j = neighs[jj];
	int jtag = tag[j];
	if (itag > jtag) {
	if (((itag + jtag) & 1) == 0)
	continue;
	} else if (itag < jtag) {
	if (((itag + jtag) & 1) == 1)
	continue;
	} else {
	if (x[j].z < z_i)
	continue;
	if (x[j].z == z_i && x[j].y < y_i)
	continue;
	if (x[j].z == z_i && x[j].y == y_i && x[j].x < x_i)
	continue;
	}
	ref_frebo_single_interaction(ka, i, j);
	}
	}


	template<typename flt_t, typename acc_t>
	void ref_frebo(KernelArgsAIREBOT<flt_t,acc_t> * ka, int torflag) {
	for (int i = ka->frebo_from_atom; i < ka->frebo_to_atom; i++) {
	ref_frebo_single_atom(ka, i);
	}
	if (torflag) ref_torsion(ka);
	}

	template<typename flt_t, typename acc_t>
	void ref_lennard_jones_single_interaction(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	int i, int j, int morseflag) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	ResultForceT<acc_t> * result_f = ka->result_f;

	int itype = map[x[i].w];
	int jtype = map[x[j].w];

	flt_t delx = x[i].x - x[j].x;
	flt_t dely = x[i].y - x[j].y;
	flt_t delz = x[i].z - x[j].z;
	flt_t rsq = delx * delx + dely * dely + delz * delz;

	if (rsq >= ka->params.cutljsq[itype][jtype]) { return; }
	flt_t rij = overloaded::sqrt(rsq);

	LennardJonesPathAIREBOT<flt_t> testpath;
	flt_t cij = 1.0;
	if (rij < ka->params.cut3rebo) {
	#pragma noinline
	cij = ref_lennard_jones_test_path<flt_t,acc_t>(ka, i, j, rij,
	ka->params.rcmax[itype][jtype], &testpath);
	}
	if (cij == 0) {
	return;
	}

	flt_t sigcut = ka->params.sigcut;
	flt_t sigmin = ka->params.sigmin;
	flt_t sigma = ka->params.sigma[itype][jtype];
	flt_t rljmax = sigcut * sigma;
	flt_t rljmin = sigmin * sigma;

	flt_t dslw, slw = Sp2(rij, rljmin, rljmax, &dslw);

	flt_t vdw, dvdw;
	if (morseflag) {
	const flt_t exr = exp(-rij * ka->params.lj4[itype][jtype]);
	vdw = ka->params.lj1[itype][jtype] * exr *
	(ka->params.lj2[itype][jtype]*exr - 2);
	dvdw = ka->params.lj3[itype][jtype] * exr *
	(1 - ka->params.lj2[itype][jtype]*exr);
	} else {
	flt_t r2inv = 1 / rsq;
	flt_t r6inv = r2inv * r2inv * r2inv;

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

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

	flt_t dStr, Str = Sp2<flt_t>(rij, ka->params.rcLJmin[itype][jtype],
	ka->params.rcLJmax[itype][jtype], &dStr);
	flt_t VA = Str * cij * VLJ;
	flt_t Stb = 0;
	acc_t fij[3] = {0};
	if (Str > 0) {
	#pragma noinline
	Stb = ref_lennard_jones_bondorder(ka, i, j, VA, fij);
	}
	flt_t fpair = -(dStr * (Stb * cij * VLJ - cij * VLJ) +
	dVLJ * (Str * Stb * cij + cij - Str * cij)) / rij;
	flt_t evdwl = VA * Stb + (1 - Str) * cij * VLJ;
	result_f[i].x += fpair * delx + fij[0];
	result_f[i].y += fpair * dely + fij[1];
	result_f[i].z += fpair * delz + fij[2];
	result_f[j].x -= fpair * delx + fij[0];
	result_f[j].y -= fpair * dely + fij[1];
	result_f[j].z -= fpair * delz + fij[2];
	ka->result_eng += evdwl;

	if (cij < 1) {
	#pragma noinline
	ref_lennard_jones_force_path(ka, Str * Stb * VLJ + (1 - Str) * VLJ,
	&testpath);
	}
	}

	template<typename flt_t, typename acc_t>
	void ref_lennard_jones_single_atom(KernelArgsAIREBOT<flt_t,acc_t> * ka, int i,
	int morseflag) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * tag = ka->tag;
	int jj;
	int itag = tag[i];
	int * neighs = ka->neigh_lmp.entries + ka->neigh_lmp.offset[i];
	int jnum = ka->neigh_lmp.num_half[i];
	for (jj = 0; jj < jnum; jj++) {
	int j = neighs[jj];
	ref_lennard_jones_single_interaction(ka, i, j, morseflag);
	}
	}

	template<typename flt_t, typename acc_t>
	void ref_lennard_jones(KernelArgsAIREBOT<flt_t,acc_t> * ka, int morseflag) {
	for (int i = ka->frebo_from_atom; i < ka->frebo_to_atom; i++) {
	#pragma noinline
	ref_lennard_jones_single_atom(ka, i, morseflag);
	}
	}

	/* ----------------------------------------------------------------------
	Vectorized AIREBO implementation, standalone, using caching to reduce
	memory access.
	---------------------------------------------------------------------- */

	template<typename flt_t, typename acc_t>
	struct aut_wrap {

	typedef typename intr_types<flt_t, acc_t>::fvec fvec;
	typedef typename intr_types<flt_t, acc_t>::avec avec;
	typedef typename intr_types<flt_t, acc_t>::ivec ivec;
	typedef typename intr_types<flt_t, acc_t>::bvec bvec;

	VEC_INLINE inline
	static void aut_loadatoms_vec(
	AtomAIREBOT<flt_t> * atoms, ivec j_vec,
	fvec x, fvec y, fvec * z, bvec * type_mask, int * map, ivec map_i,
	ivec c_1
	) {
	const ivec c_4 = ivec::set1(4);
	ivec j_vec_4 = ivec::mullo(c_4, j_vec);
	fvec w;
	fvec::gather_4_adjacent(j_vec_4, &atoms[0].x, sizeof(flt_t), x, y, z, &w);
	ivec jtype = fvec::unpackloepi32(w);
	jtype = ivec::srlv(map_i, jtype); //_mm512_castpd_si512(w));
	jtype = ivec::the_and(c_1, jtype);
	bvec jtype_mask = ivec::cmpneq(jtype, ivec::setzero());
	*type_mask = jtype_mask;
	}

	VEC_INLINE inline
	static void aut_loadatoms_vec_notype(
	AtomAIREBOT<flt_t> * atoms, ivec j_vec,
	fvec x, fvec y, fvec * z
	) {
	const ivec c_4 = ivec::set1(4);
	ivec j_vec_4 = ivec::mullo(c_4, j_vec);
	fvec::gather_3_adjacent(j_vec_4, &atoms[0].x, sizeof(flt_t), x, y, z);
	}

	static fvec aut_Sp2_deriv(fvec r, fvec lo, fvec hi, fvec * d) {
	fvec c_1 = fvec::set1(1);
	fvec c_2 = fvec::set1(2);
	fvec c_3 = fvec::set1(3);
	fvec c_6 = fvec::set1(6);
	bvec m_lo = fvec::cmple(r, lo);
	bvec m_hi = fvec::cmpnlt(r, hi); // nlt == ge
	bvec m_tr = bvec::kandn(m_lo, ~ m_hi);
	fvec ret = c_1;
	ret = fvec::mask_blend(m_hi, ret, fvec::setzero());
	fvec der = fvec::setzero();
	if (bvec::test_any_set(m_tr)) {
	fvec diff = hi - lo;
	fvec rcp = fvec::recip(diff);
	fvec t = (r - lo) * rcp;
	fvec v = c_1 - t * t * ( c_3 - c_2 * t);
	ret = fvec::mask_blend(m_tr, ret, v);
	fvec dv = c_6 * rcp * ( t * t - t);
	der = fvec::mask_blend(m_tr, der, dv);
	}
	*d = der;
	return ret;
	}

	static fvec aut_Sp_deriv(fvec r, fvec lo, fvec hi, fvec * d) {
	fvec c_1 = fvec::set1(1);
	fvec c_0_5 = fvec::set1(0.5);
	fvec c_m0_5 = fvec::set1(-0.5);
	fvec c_PI = fvec::set1(M_PI);
	bvec m_lo = fvec::cmple(r, lo);
	bvec m_hi = fvec::cmpnlt(r, hi); // nlt == ge
	bvec m_tr = bvec::kandn(m_lo, ~ m_hi);
	fvec ret = c_1;
	ret = fvec::mask_blend(m_hi, ret, fvec::setzero());
	fvec der = fvec::setzero();
	if (bvec::test_any_set(m_tr)) {
	fvec diff = hi - lo;
	fvec rcp = fvec::mask_recip(c_1, m_tr, diff);
	fvec t = (r - lo) / diff;
	fvec sinval, cosval;
	sinval = fvec::mask_sincos(&cosval, fvec::setzero(), c_1, m_tr, c_PI * t);
	fvec v = c_0_5 * ( c_1 + cosval);
	ret = fvec::mask_blend(m_tr, ret, v);
	fvec dv = c_PI * c_m0_5 * rcp * sinval;
	der = fvec::mask_blend(m_tr, der, dv);
	}
	*d = der;
	return ret;
	}

	static fvec aut_mask_Sp(bvec mask, fvec r, fvec lo, fvec hi) {
	fvec c_1 = fvec::set1(1);
	fvec c_0_5 = fvec::set1(0.5);
	fvec c_PI = fvec::set1(M_PI);
	bvec m_lo = fvec::mask_cmple(mask, r, lo);
	bvec m_hi = fvec::mask_cmpnlt(mask, r, hi); // nlt == ge
	bvec m_tr = bvec::kandn(m_lo, bvec::kandn(m_hi, mask));
	fvec ret = c_1;
	ret = fvec::mask_blend(m_hi, ret, fvec::setzero());
	if (bvec::test_any_set(m_tr)) {
	fvec rcp = fvec::mask_recip(c_1, m_tr, hi - lo);
	fvec t = (r - lo) * rcp;
	fvec v = c_0_5 * ( c_1 + fvec::mask_cos(c_1, m_tr, c_PI * t));
	ret = fvec::mask_blend(m_tr, ret, v);
	}
	return ret;
	}

	static void aut_rebo_neigh(KernelArgsAIREBOT<flt_t,acc_t> * ka) {
	int offset = ka->neigh_from_atom * ka->num_neighs_per_atom;
	ivec c_CARBON = ivec::setzero();
	int map_i = 0;
	int i;
	for (i = 1; i < ka->num_types; i++) {
	if (ka->map[i])
	map_i \|= (1 << i);
	}
	ivec c_i1 = ivec::set1(1);
	ivec c_im = ivec::set1(map_i);
	AtomAIREBOT<flt_t> * _noalias x = ka->x;

	for (i = ka->neigh_from_atom; i < ka->neigh_to_atom; i++) {

	fvec x_i = fvec::set1(x[i].x);
	fvec y_i = fvec::set1(x[i].y);
	fvec z_i = fvec::set1(x[i].z);
	int itype = ka->map[ka->x[i].w];

	fvec rcmaxsq0 = fvec::set1(ka->params.rcmaxsq[itype][0]);
	fvec rcmaxsq1 = fvec::set1(ka->params.rcmaxsq[itype][1]);
	fvec rcmax0 = fvec::set1(ka->params.rcmax[itype][0]);
	fvec rcmax1 = fvec::set1(ka->params.rcmax[itype][1]);
	fvec rcmin0 = fvec::set1(ka->params.rcmin[itype][0]);
	fvec rcmin1 = fvec::set1(ka->params.rcmin[itype][1]);
	fvec rcmaxskinsq0 = fvec::set1(
	(ka->params.rcmax[itype][0] + ka->skin) * (ka->params.rcmax[itype][0] +
	ka->skin));
	fvec rcmaxskinsq1 = fvec::set1(
	(ka->params.rcmax[itype][1] + ka->skin) * (ka->params.rcmax[itype][1] +
	ka->skin));
	fvec nC = fvec::setzero();
	fvec nH = fvec::setzero();

	ka->neigh_rebo.offset[i] = offset;

	int jnum = ka->rebuild_flag ? ka->neigh_lmp.num[i] :
	ka->neigh_rebo.num_half[i];
	int * neighs = ka->rebuild_flag ?
	&ka->neigh_lmp.entries[ka->neigh_lmp.offset[i]] :
	&ka->neigh_rebo.entries[ka->neigh_rebo.offset[i]+jnum];
	int * skin_target = &ka->neigh_rebo.entries[offset+ka->num_neighs_per_atom];
	int n = 0;
	int n_skin = 0;

	int lowest_idx;
	#pragma unroll(4)
	for (lowest_idx = 0; lowest_idx < jnum; lowest_idx += fvec::VL) {
	bvec j_mask = bvec::full();
	if (lowest_idx + fvec::VL > jnum) j_mask = bvec::only(jnum - lowest_idx);

	int * _noalias neighs_l = neighs + lowest_idx;
	fvec x_j, y_j, z_j;
	bvec jtype_mask;
	ivec ji = ivec::maskz_loadu(j_mask, neighs_l);
	aut_loadatoms_vec(x, ji,
	&x_j, &y_j, &z_j, &jtype_mask, ka->map, c_im, c_i1);
	fvec delx = x_i - x_j;
	fvec dely = y_i - y_j;
	fvec delz = z_i - z_j;
	fvec rsq = delx * delx + dely * dely + delz * delz;
	if (ka->rebuild_flag) {
	fvec rcmaxskinsq = fvec::mask_blend(jtype_mask, rcmaxskinsq0,
	rcmaxskinsq1);
	bvec c_mask = fvec::mask_cmplt(j_mask, rsq, rcmaxskinsq);
	ivec::mask_compressstore(c_mask, &skin_target[n_skin], ji);
	n_skin += bvec::popcnt(c_mask);
	}
	fvec rcmaxsq = fvec::mask_blend(jtype_mask, rcmaxsq0, rcmaxsq1);
	bvec c_mask = fvec::mask_cmplt(j_mask, rsq, rcmaxsq);
	if (bvec::test_all_unset(c_mask)) continue;
	ivec::mask_compressstore(c_mask, &ka->neigh_rebo.entries[offset + n], ji);
	n += bvec::popcnt(c_mask);
	fvec rcmax = fvec::mask_blend(jtype_mask, rcmax0, rcmax1);
	fvec rcmin = fvec::mask_blend(jtype_mask, rcmin0, rcmin1);
	fvec sp = aut_mask_Sp(c_mask, fvec::sqrt(rsq), rcmin, rcmax);
	nC = fvec::mask_add(nC, bvec::kandn(jtype_mask, c_mask), nC, sp);
	nH = fvec::mask_add(nH, bvec::kand (jtype_mask, c_mask), nH, sp);
	}
	ka->neigh_rebo.num[i] = n;
	if (ka->rebuild_flag) {
	for (int i = 0; i < n_skin; i++) {
	ka->neigh_rebo.entries[offset+n_skin+i] = skin_target[i];
	}
	}
	if (ka->rebuild_flag) {
	assert(n <= n_skin);
	offset += 2 * n_skin;
	ka->neigh_rebo.num_half[i] = n_skin;
	} else {
	assert(n <= jnum);
	offset += 2 * jnum;
	}
	ka->nC[i] = fvec::reduce_add(nC);
	ka->nH[i] = fvec::reduce_add(nH);
	}
	}


	static fvec aut_eval_poly_lin_pd_2(int n, flt_t * vals, ivec idx, fvec x,
	fvec * deriv) {
	fvec c_1 = fvec::set1(1);
	fvec x_i = c_1;
	fvec x_im1 = fvec::setzero();
	fvec result = fvec::setzero();
	fvec i_v = fvec::setzero();
	*deriv = fvec::setzero();
	int i;
	for (i = 0; i < n; i++) {
	fvec coeff = fvec::gather(idx, vals + i, sizeof(flt_t));
	result = result + coeff * x_i;
	deriv = deriv + coeff * x_im1 * i_v;
	x_im1 = x_i;
	x_i = x_i * x;
	i_v = i_v + c_1;
	}
	return result;
	}

	static fvec aut_mask_gSpline_pd_2(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	bvec active_mask, int itype, fvec cosjik,
	fvec Nij, fvec dgdc, fvec dgdN) {
	int i;
	flt_t * gDom = NULL;
	int nDom = 0;
	ivec offs = ivec::setzero();
	fvec NCmin = fvec::set1(ka->params.NCmin);
	bvec Ngt = fvec::cmpnle(Nij, NCmin); //gt
	if (itype == 0) {
	nDom = 4;
	gDom = &ka->params.gCdom[0];
	offs = ivec::mask_blend(Ngt, offs, ivec::set1(4*6));
	} else {
	nDom = 3;
	gDom = &ka->params.gHdom[0];
	offs = ivec::set1(8 * 6);
	}
	cosjik = fvec::max(fvec::set1(gDom[0]), fvec::min(fvec::set1(gDom[nDom]),
	cosjik));
	ivec index6 = ivec::setzero();
	for (i = 0; i < nDom; i++) {
	bvec cosge = fvec::cmpnlt(cosjik, fvec::set1(gDom[i])); //ge
	bvec cosle = fvec::cmple(cosjik, fvec::set1(gDom[i+1]));
	index6 = ivec::mask_blend(cosge & cosle, index6, ivec::set1(6*i));
	}
	fvec g = aut_eval_poly_lin_pd_2(6, &ka->params.gVal[0], offs + index6,
	cosjik, dgdc);
	*dgdN = fvec::setzero();
	if (itype == 0) {
	fvec NCmax = fvec::set1(ka->params.NCmax);
	bvec Nlt = fvec::cmplt(Nij, NCmax); //gt
	bvec Nmask = Ngt & Nlt;
	if (bvec::test_any_set(Nmask)) {
	fvec dg1;
	fvec g1 = aut_eval_poly_lin_pd_2(6, &ka->params.gVal[0], index6, cosjik,
	&dg1);
	fvec dS;
	fvec cut = aut_Sp_deriv(Nij, NCmin, NCmax, &dS);
	dgdN = fvec::mask_mul(dgdN, Nmask, dS, g1 - g);
	g = fvec::mask_add(g, Nmask, g, cut * ( g1 - g));
	dgdc = fvec::mask_add(dgdc, Nmask, dgdc, cut ( dg1 - *dgdc));
	}
	}
	return g;
	}

	static fvec aut_PijSpline(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, fvec NijC, fvec NijH, fvec *dN2) {
	flt_t ret[fvec::VL] __attribute__((aligned(64)));
	flt_t dN20[fvec::VL] __attribute__((aligned(64)));
	flt_t dN21[fvec::VL] __attribute__((aligned(64)));
	flt_t NijC_[fvec::VL] __attribute__((aligned(64)));
	flt_t NijH_[fvec::VL] __attribute__((aligned(64)));
	flt_t tmp_dN2[2];
	fvec::store(NijC_, NijC);
	fvec::store(NijH_, NijH);
	int i;
	for (i = 0; i < fvec::VL; i++) {
	ret[i] = PijSpline(ka, itype, jtype, NijC_[i], NijH_[i], tmp_dN2);
	dN20[i] = tmp_dN2[0];
	dN21[i] = tmp_dN2[1];
	}
	dN2[0] = fvec::load(dN20);
	dN2[1] = fvec::load(dN21);
	return fvec::load(ret);
	}

	/*
	* aut_frebo_data stores all the short-ranged coordinations
	* and intermediate values that get reused frequently during
	* bondorder calculations.
	* BUF_CAP should rarely exceed 4, so 8 is a very conservative
	* value.
	*/
	static const int BUF_CAP = 8;
	struct aut_frebo_data {
	fvec rikx_buf[BUF_CAP];
	fvec riky_buf[BUF_CAP];
	fvec rikz_buf[BUF_CAP];
	fvec rikmag_buf[BUF_CAP];
	fvec cosjik_buf[BUF_CAP];
	ivec k_buf[BUF_CAP];
	fvec g_buf[BUF_CAP];
	fvec dgdc_buf[BUF_CAP];
	fvec ex_lam_buf[BUF_CAP];
	fvec wik_buf[BUF_CAP];
	fvec dwik_buf[BUF_CAP];
	fvec cutN_buf[BUF_CAP];
	fvec dcutN_buf[BUF_CAP];
	bvec ktype_buf[BUF_CAP];
	bvec mask_buf[BUF_CAP];
	fvec force_k_x_buf[BUF_CAP];
	fvec force_k_y_buf[BUF_CAP];
	fvec force_k_z_buf[BUF_CAP];
	int buf_len;
	fvec x_i;
	fvec y_i;
	fvec z_i;
	fvec x_j;
	fvec y_j;
	fvec z_j;
	fvec nCi;
	fvec nHi;
	fvec force_i_x;
	fvec force_i_y;
	fvec force_i_z;
	fvec force_j_x;
	fvec force_j_y;
	fvec force_j_z;
	};

	/*
	* Initialize values in aut_frebo_data and perform the calculations
	* for p_ij.
	*/
	static fvec aut_frebo_pij_pd_2(
	KernelArgsAIREBOT<flt_t,acc_t> * _noalias ka,
	struct aut_frebo_data * _noalias data,
	int itype, int jtype,
	ivec vi, ivec vj,
	fvec rijx, fvec rijy, fvec rijz, fvec rijmag,
	fvec wij, fvec VA, fvec * sum_N, fvec fij[3]
	) {
	AtomAIREBOT<flt_t> * _noalias x = ka->x;
	int * _noalias map = ka->map;
	flt_t * _noalias nC = ka->nC;
	flt_t * _noalias nH = ka->nH;
	fvec x_i, y_i, z_i;
	fvec x_j, y_j, z_j;
	x_i = data->x_i;
	y_i = data->y_i;
	z_i = data->z_i;
	x_j = data->x_j;
	y_j = data->y_j;
	z_j = data->z_j;
	fvec invrijm = fvec::recip(rijmag);
	fvec invrijm2 = invrijm * invrijm;
	fvec rcminij = fvec::set1(ka->params.rcmin[itype][jtype]);
	fvec rcmaxij = fvec::set1(ka->params.rcmax[itype][jtype]);
	fvec Nmin = fvec::set1(ka->params.Nmin);
	fvec Nmax = fvec::set1(ka->params.Nmax);
	int map_i_scalar = 0;
	{
	int i;
	for (i = 1; i < ka->num_types; i++) {
	if (ka->map[i])
	map_i_scalar \|= (1 << i);
	}
	}
	ivec map_i = ivec::set1(map_i_scalar);
	fvec nCi = data->nCi;
	fvec nHi = data->nHi;
	fvec Nij = nHi + nCi - wij;
	fvec factor_jtype, factor_not_jtype;
	if (jtype) {
	factor_jtype = fvec::set1(1);
	factor_not_jtype = fvec::set1(0);
	} else {
	factor_jtype = fvec::set1(0);
	factor_not_jtype = fvec::set1(1);
	}
	fvec NijC = nCi - wij * factor_not_jtype;
	fvec NijH = nHi - wij * factor_jtype;
	fvec sum_pij = fvec::setzero();
	fvec sum_dpij_dN = fvec::setzero();
	fvec dN2[2];
	ivec offseti = ivec::mask_gather(ivec::setzero(), bvec::full(), vi,
	ka->neigh_rebo.offset, sizeof(int));
	int buf_len = 0;
	ivec knum = ivec::mask_gather(ivec::setzero(), bvec::full(), vi,
	ka->neigh_rebo.num, sizeof(int));
	ivec kk = ivec::setzero();
	bvec active_mask = ivec::cmplt(kk, knum);
	ivec c_i1 = ivec::set1(1);
	fvec rho_j = fvec::set1(ka->params.rho[jtype][1]);
	fvec rho_k0 = fvec::set1(ka->params.rho[0][1]);
	fvec rho_k1 = fvec::set1(ka->params.rho[1][1]);
	fvec c_4 = fvec::set1(4);
	fvec c_2_0 = fvec::set1(2.0);
	fvec c_m2_0 = fvec::set1(-2.0);
	fvec c_4_0 = fvec::set1(4.0);
	fvec c_0_5 = fvec::set1(0.5);
	fvec c_m0_5 = fvec::set1(-0.5);
	fvec c_1 = fvec::set1(1);
	fvec c_m1 = fvec::set1(-1);
	fvec factor_itype = itype ? c_1 : fvec::setzero();
	fvec rcmax0 = fvec::set1(ka->params.rcmax[itype][0]);
	fvec rcmax1 = fvec::set1(ka->params.rcmax[itype][1]);
	fvec rcmin0 = fvec::set1(ka->params.rcmin[itype][0]);
	fvec rcmin1 = fvec::set1(ka->params.rcmin[itype][1]);
	fvec result_f_i_x = fvec::setzero();
	fvec result_f_i_y = fvec::setzero();
	fvec result_f_i_z = fvec::setzero();
	fvec result_f_j_x = fvec::setzero();
	fvec result_f_j_y = fvec::setzero();
	fvec result_f_j_z = fvec::setzero();
	*sum_N = fvec::setzero();
	{
	while (bvec::test_any_set(active_mask)) {
	ivec k = ivec::mask_gather(ivec::setzero(), active_mask, kk + offseti,
	ka->neigh_rebo.entries, sizeof(int));
	bvec excluded_mask = ivec::cmpeq(k, vj) & active_mask;
	if (bvec::test_any_set(excluded_mask)) {
	kk = ivec::mask_add(kk, excluded_mask, kk, c_i1);
	active_mask = ivec::cmplt(kk, knum);
	continue;
	}
	fvec x_k, y_k, z_k;
	bvec ktype_mask;
	aut_loadatoms_vec(x, k, &x_k, &y_k, &z_k, &ktype_mask, ka->map, map_i,
	c_i1);
	fvec rikx = x_i - x_k;
	fvec riky = y_i - y_k;
	fvec rikz = z_i - z_k;
	fvec rikmag = fvec::sqrt(rikx * rikx + riky * riky + rikz * rikz);
	fvec rho_k = fvec::mask_blend(ktype_mask, rho_k0, rho_k1);
	fvec lamdajik = c_4 * factor_itype * ( rho_k - rikmag - ( rho_j -
	rijmag));
	fvec ex_lam = fvec::exp(lamdajik);
	fvec rcmax = fvec::mask_blend(ktype_mask, rcmax0, rcmax1);
	fvec rcmin = fvec::mask_blend(ktype_mask, rcmin0, rcmin1);
	fvec dwik;
	fvec wik = aut_Sp_deriv(rikmag, rcmin, rcmax, &dwik);
	fvec Nki = fvec::gather(k, nC, sizeof(flt_t)) +
	fvec::gather(k, nH, sizeof(flt_t)) - wik;
	fvec cosjik = (rijx * rikx + rijy * riky + rijz * rikz) /
	( rijmag * rikmag);
	cosjik = fvec::min(c_1, fvec::max(c_m1, cosjik));
	fvec dgdc, dgdN;
	fvec g = aut_mask_gSpline_pd_2(ka, active_mask, itype, cosjik, Nij,
	&dgdc, &dgdN);
	sum_pij = fvec::mask_add(sum_pij, active_mask, sum_pij, wik * g * ex_lam);
	sum_dpij_dN = fvec::mask_add(sum_dpij_dN, active_mask, sum_dpij_dN,
	wik * ex_lam * dgdN);
	fvec dcutN;
	fvec cutN = aut_Sp_deriv(Nki, Nmin, Nmax, &dcutN);
	sum_N = fvec::mask_add(sum_N, active_mask, *sum_N,
	fvec::mask_blend(ktype_mask, c_1,
	fvec::setzero()) * wik * cutN);
	if (buf_len == BUF_CAP) goto exceed_buffer;
	data->rikx_buf[buf_len] = rikx;
	data->riky_buf[buf_len] = riky;
	data->rikz_buf[buf_len] = rikz;
	data->rikmag_buf[buf_len] = rikmag;
	data->cosjik_buf[buf_len] = cosjik;
	data->ktype_buf[buf_len] = ktype_mask;
	data->k_buf[buf_len] = k;
	data->g_buf[buf_len] = g;
	data->dgdc_buf[buf_len] = dgdc;
	data->ex_lam_buf[buf_len] = ex_lam;
	data->wik_buf[buf_len] = wik;
	data->dwik_buf[buf_len] = dwik;
	data->mask_buf[buf_len] = active_mask;
	data->cutN_buf[buf_len] = cutN;
	data->dcutN_buf[buf_len] = dcutN;
	buf_len += 1;
	kk = ivec::mask_add(kk, active_mask, kk, c_i1);
	active_mask = ivec::cmplt(kk, knum);
	}
	data->buf_len = buf_len;
	fvec PijS = aut_PijSpline(ka, itype, jtype, NijC, NijH, &dN2[0]);
	fvec pij = fvec::invsqrt(c_1 + sum_pij + PijS);
	fvec tmp = c_m0_5 * pij * pij * pij;
	int buf_idx;
	for (buf_idx = 0; buf_idx < buf_len; buf_idx++) {
	fvec rikx = data->rikx_buf[buf_idx];
	fvec riky = data->riky_buf[buf_idx];
	fvec rikz = data->rikz_buf[buf_idx];
	fvec rikmag = data->rikmag_buf[buf_idx];
	fvec cosjik = data->cosjik_buf[buf_idx];
	bvec ktype_mask = data->ktype_buf[buf_idx];
	ivec k = data->k_buf[buf_idx];
	fvec g = data->g_buf[buf_idx];
	fvec dgdc = data->dgdc_buf[buf_idx];
	fvec ex_lam = data->ex_lam_buf[buf_idx];
	fvec wik = data->wik_buf[buf_idx];
	fvec dwik = data->dwik_buf[buf_idx];
	bvec mask = data->mask_buf[buf_idx];
	fvec invrikm = fvec::recip(rikmag);
	fvec rjkx = rikx - rijx;
	fvec rjky = riky - rijy;
	fvec rjkz = rikz - rijz;
	fvec rjkmag = fvec::sqrt(
	rjkx * rjkx + rjky * rjky + rjkz * rjkz);
	fvec rijrik = c_2_0 * rijmag * rikmag;
	fvec rr = rijmag * rijmag - rikmag * rikmag;
	fvec dctdjk = c_m2_0 / rijrik;
	fvec dctdik = (rjkmag * rjkmag - rr) / ( rijrik * rikmag * rikmag);
	fvec dctdij = (rjkmag * rjkmag + rr) / ( rijrik * rijmag * rijmag);
	fvec fi[3], fj[3], fk[3];
	fvec pref = c_0_5 * VA * tmp;
	fvec tmp20 = pref * wik * dgdc * ex_lam;
	fj[0] = fj[1] = fj[2] = fvec::setzero();
	fvec tmpdik = tmp20 * dctdik;
	fi[0] = fvec::setzero() - tmpdik * rikx;
	fi[1] = fvec::setzero() - tmpdik * riky;
	fi[2] = fvec::setzero() - tmpdik * rikz;
	fk[0] = tmpdik * rikx;
	fk[1] = tmpdik * riky;
	fk[2] = tmpdik * rikz;

	fvec tmpdij = tmp20 * dctdij;
	fij[0] = fvec::mask_sub(fij[0], mask, fij[0], tmpdij * rijx);
	fij[1] = fvec::mask_sub(fij[1], mask, fij[1], tmpdij * rijy);
	fij[2] = fvec::mask_sub(fij[2], mask, fij[2], tmpdij * rijz);

	fvec tmpdjk = tmp20 * dctdjk;
	fi[0] = fi[0] - tmpdjk * rjkx;
	fi[1] = fi[1] - tmpdjk * rjky;
	fi[2] = fi[2] - tmpdjk * rjkz;
	fk[0] = fk[0] + tmpdjk * rjkx;
	fk[1] = fk[1] + tmpdjk * rjky;
	fk[2] = fk[2] + tmpdjk * rjkz;
	fij[0] = fvec::mask_add(fij[0], mask, fij[0], tmpdjk * rjkx);
	fij[1] = fvec::mask_add(fij[1], mask, fij[1], tmpdjk * rjky);
	fij[2] = fvec::mask_add(fij[2], mask, fij[2], tmpdjk * rjkz);

	if (itype) {
	fvec tmp21 = pref * wik * g * ex_lam * c_4_0;
	fvec tmp21ij = tmp21 * invrijm;
	fij[0] = fvec::mask_sub(fij[0], mask, fij[0], tmp21ij * rijx);
	fij[1] = fvec::mask_sub(fij[1], mask, fij[1], tmp21ij * rijy);
	fij[2] = fvec::mask_sub(fij[2], mask, fij[2], tmp21ij * rijz);
	fvec tmp21ik = tmp21 * invrikm;
	fi[0] = fi[0] + tmp21ik * rikx;
	fi[1] = fi[1] + tmp21ik * riky;
	fi[2] = fi[2] + tmp21ik * rikz;
	fk[0] = fk[0] - tmp21ik * rikx;
	fk[1] = fk[1] - tmp21ik * riky;
	fk[2] = fk[2] - tmp21ik * rikz;
	}

	// coordination forces

	// dwik forces
	fvec tmp22 = pref * dwik * g * ex_lam * invrikm;
	fi[0] = fi[0] - tmp22 * rikx;
	fi[1] = fi[1] - tmp22 * riky;
	fi[2] = fi[2] - tmp22 * rikz;
	fk[0] = fk[0] + tmp22 * rikx;
	fk[1] = fk[1] + tmp22 * riky;
	fk[2] = fk[2] + tmp22 * rikz;

	// PIJ forces
	fvec dN2ktype = fvec::mask_blend(ktype_mask, dN2[0], dN2[1]);
	fvec tmp23 = pref * dN2ktype * dwik * invrikm;
	fi[0] = fi[0] - tmp23 * rikx;
	fi[1] = fi[1] - tmp23 * riky;
	fi[2] = fi[2] - tmp23 * rikz;
	fk[0] = fk[0] + tmp23 * rikx;
	fk[1] = fk[1] + tmp23 * riky;
	fk[2] = fk[2] + tmp23 * rikz;

	// dgdN forces
	fvec tmp24 = pref * sum_dpij_dN * dwik * invrikm;
	fi[0] = fi[0] - tmp24 * rikx;
	fi[1] = fi[1] - tmp24 * riky;
	fi[2] = fi[2] - tmp24 * rikz;
	fk[0] = fk[0] + tmp24 * rikx;
	fk[1] = fk[1] + tmp24 * riky;
	fk[2] = fk[2] + tmp24 * rikz;

	result_f_i_x = fvec::mask_add(result_f_i_x, mask, result_f_i_x, fi[0]);
	result_f_i_y = fvec::mask_add(result_f_i_y, mask, result_f_i_y, fi[1]);
	result_f_i_z = fvec::mask_add(result_f_i_z, mask, result_f_i_z, fi[2]);
	result_f_j_x = fvec::mask_add(result_f_j_x, mask, result_f_j_x, fj[0]);
	result_f_j_y = fvec::mask_add(result_f_j_y, mask, result_f_j_y, fj[1]);
	result_f_j_z = fvec::mask_add(result_f_j_z, mask, result_f_j_z, fj[2]);

	data->force_k_x_buf[buf_idx] = fk[0];
	data->force_k_y_buf[buf_idx] = fk[1];
	data->force_k_z_buf[buf_idx] = fk[2];
	}
	data->force_i_x = result_f_i_x;
	data->force_i_y = result_f_i_y;
	data->force_i_z = result_f_i_z;
	data->force_j_x = result_f_j_x;
	data->force_j_y = result_f_j_y;
	data->force_j_z = result_f_j_z;
	return pij;
	}
	exceed_buffer:
	data->buf_len = -1;
	return fvec::setzero();
	}

	/*
	* Apply the force values stored iin aut_frebo_data to
	* the respective neighbors.
	*/
	static void aut_frebo_data_writeback(
	KernelArgsAIREBOT<flt_t,acc_t> * _noalias ka,
	struct aut_frebo_data * _noalias data) {
	ResultForceT<acc_t> * _noalias result_f = ka->result_f;
	flt_t fk_x_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fk_y_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fk_z_buf[fvec::VL] __attribute__((aligned(64)));
	int fk_k_buf[ivec::VL] __attribute__((aligned(64)));
	int buf_idx;
	for (buf_idx = 0; buf_idx < data->buf_len; buf_idx++) {
	ivec k = data->k_buf[buf_idx];
	bvec active_mask = data->mask_buf[buf_idx];

	fvec::store(fk_x_buf, data->force_k_x_buf[buf_idx]);
	fvec::store(fk_y_buf, data->force_k_y_buf[buf_idx]);
	fvec::store(fk_z_buf, data->force_k_z_buf[buf_idx]);
	ivec::store(fk_k_buf, k);

	int lane;
	for (lane = 0; lane < fvec::VL; lane++) {
	if (bvec::test_at(active_mask, lane)) {} else continue;
	int kk = fk_k_buf[lane];
	result_f[kk].x += fk_x_buf[lane];
	result_f[kk].y += fk_y_buf[lane];
	result_f[kk].z += fk_z_buf[lane];
	}
	}
	}

	static void aut_frebo_N_spline_force(
	KernelArgsAIREBOT<flt_t,acc_t> * _noalias ka,
	struct aut_frebo_data * _noalias data, int itype, int jtype, ivec vi,
	ivec vj, fvec VA, fvec dN, fvec dNconj, fvec Nconj) {
	ivec c_i1 = ivec::set1(1);
	fvec c_2 = fvec::set1(2);
	fvec c_TOL = fvec::set1(TOL);
	ResultForceT<acc_t> * _noalias result_f = ka->result_f;
	AtomAIREBOT<flt_t> * _noalias x = ka->x;
	int * _noalias map = ka->map;
	flt_t * _noalias nC = ka->nC;
	flt_t * _noalias nH = ka->nH;
	fvec x_i, y_i, z_i;
	x_i = data->x_i;
	y_i = data->y_i;
	z_i = data->z_i;
	fvec Nmin = fvec::set1(ka->params.Nmin);
	fvec Nmax = fvec::set1(ka->params.Nmax);
	int map_i_scalar = 0;
	{
	int i;
	for (i = 1; i < ka->num_types; i++) {
	if (ka->map[i])
	map_i_scalar \|= (1 << i);
	}
	}
	ivec map_i = ivec::set1(map_i_scalar);
	fvec dN2[2];
	ivec kk = ivec::setzero();
	fvec rcmax0 = fvec::set1(ka->params.rcmax[itype][0]);
	fvec rcmax1 = fvec::set1(ka->params.rcmax[itype][1]);
	fvec rcmin0 = fvec::set1(ka->params.rcmin[itype][0]);
	fvec rcmin1 = fvec::set1(ka->params.rcmin[itype][1]);
	fvec result_f_i_x = fvec::setzero();
	fvec result_f_i_y = fvec::setzero();
	fvec result_f_i_z = fvec::setzero();
	int buf_idx;
	for (buf_idx = 0; buf_idx < data->buf_len; buf_idx++) {
	ivec k = data->k_buf[buf_idx];
	bvec active_mask = data->mask_buf[buf_idx];
	fvec rikx = data->rikx_buf[buf_idx];
	fvec riky = data->riky_buf[buf_idx];
	fvec rikz = data->rikz_buf[buf_idx];
	fvec rikmag = data->rikmag_buf[buf_idx];
	bvec ktype_mask = data->ktype_buf[buf_idx];

	fvec dwik = data->dwik_buf[buf_idx];
	fvec wik = data->wik_buf[buf_idx];

	fvec dNki = data->dcutN_buf[buf_idx];
	fvec SpN = data->cutN_buf[buf_idx];

	fvec invrikmag = fvec::recip(rikmag);
	fvec pref = VA * dwik * invrikmag;
	fvec fdN = dN * pref;
	fvec fdNconj = pref * SpN * c_2 * dNconj * Nconj;
	fvec ffactor = fdN;
	bvec ktype_is_C = ~ ktype_mask;
	ffactor = fvec::mask_add(ffactor, ktype_is_C, ffactor, fdNconj);

	fvec fkx = ffactor * rikx;
	fvec fky = ffactor * riky;
	fvec fkz = ffactor * rikz;

	data->force_k_x_buf[buf_idx] = data->force_k_x_buf[buf_idx] + fkx;
	data->force_k_y_buf[buf_idx] = data->force_k_y_buf[buf_idx] + fky;
	data->force_k_z_buf[buf_idx] = data->force_k_z_buf[buf_idx] + fkz;

	result_f_i_x = fvec::mask_sub(result_f_i_x, active_mask, result_f_i_x, fkx);
	result_f_i_y = fvec::mask_sub(result_f_i_y, active_mask, result_f_i_y, fky);
	result_f_i_z = fvec::mask_sub(result_f_i_z, active_mask, result_f_i_z, fkz);

	bvec need_k_neighs = fvec::mask_cmpnle(active_mask, fvec::abs(dNki), c_TOL)
	& ktype_is_C;
	if (bvec::test_any_set(need_k_neighs)) {
	int lane;
	for (lane = 0; lane < fvec::VL; lane++) {
	if (! bvec::test_at(need_k_neighs, lane)) continue;
	int kk = ivec::at(k, lane);
	int k = kk;
	int ktype = map[x[k].w];
	int i = ivec::at(vi, lane);
	fvec oldVA = VA;
	double VA = fvec::at(oldVA, lane);
	fvec oldwik = wik;
	double wik = fvec::at(oldwik, lane);
	fvec olddNconj = dNconj;
	double dNconj = fvec::at(olddNconj, lane);
	fvec oldNconj = Nconj;
	double Nconj = fvec::at(oldNconj, lane);
	fvec olddNki = dNki;
	double dNki = fvec::at(olddNki, lane);
	int * neighs_k = ka->neigh_rebo.entries + ka->neigh_rebo.offset[k];
	int nnum = ka->neigh_rebo.num[k];
	int nn;
	for (nn = 0; nn < nnum; nn++) {
	int n = neighs_k[nn];
	if (n == i) continue;
	double rknx = x[k].x - x[n].x;
	double rkny = x[k].y - x[n].y;
	double rknz = x[k].z - x[n].z;
	double rknmag = sqrt(rknx * rknx + rkny * rkny + rknz * rknz);
	int ntype = map[x[n].w];
	double rcminkn = ka->params.rcmin[ktype][ntype];
	double rcmaxkn = ka->params.rcmax[ktype][ntype];
	double dwkn;
	Sp(rknmag, rcminkn, rcmaxkn, &dwkn);
	double ffactor = VA * dNconj * 2 * Nconj * wik * dNki * dwkn / rknmag;
	result_f[k].x -= ffactor * rknx;
	result_f[k].y -= ffactor * rkny;
	result_f[k].z -= ffactor * rknz;
	result_f[n].x += ffactor * rknx;
	result_f[n].y += ffactor * rkny;
	result_f[n].z += ffactor * rknz;
	}
	}
	}
	}
	data->force_i_x = data->force_i_x + result_f_i_x;
	data->force_i_y = data->force_i_y + result_f_i_y;
	data->force_i_z = data->force_i_z + result_f_i_z;
	}

	static fvec aut_frebo_pi_rc_pd(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, fvec Nij, fvec Nji, fvec Nijconj,
	fvec * dN3) {
	flt_t ret[fvec::VL] __attribute__((aligned(64)));
	flt_t dN3ret[3][fvec::VL] __attribute__((aligned(64)));
	int i;
	for (i = 0; i < fvec::VL; i++) {
	flt_t dN3tmp[3];
	ret[i] = frebo_pi_rc(ka, itype, jtype, fvec::at(Nij, i), fvec::at(Nji, i),
	fvec::at(Nijconj, i), &dN3tmp[0]);
	dN3ret[0][i] = dN3tmp[0];
	dN3ret[1][i] = dN3tmp[1];
	dN3ret[2][i] = dN3tmp[2];
	}
	dN3[0] = fvec::load(&dN3ret[0][0]);
	dN3[1] = fvec::load(&dN3ret[1][0]);
	dN3[2] = fvec::load(&dN3ret[2][0]);
	return fvec::load(&ret[0]);
	}

	static fvec aut_frebo_Tij(KernelArgsAIREBOT<flt_t,acc_t> * ka, int itype,
	int jtype, fvec Nij, fvec Nji, fvec Nijconj,
	fvec * dN3) {
	flt_t ret[fvec::VL] __attribute__((aligned(64)));
	flt_t dN3ret[3][fvec::VL] __attribute__((aligned(64)));
	int i;
	for (i = 0; i < fvec::VL; i++) {
	flt_t dN3tmp[3];
	ret[i] = frebo_Tij(ka, itype, jtype, fvec::at(Nij, i), fvec::at(Nji, i),
	fvec::at(Nijconj, i), &dN3tmp[0]);
	dN3ret[0][i] = dN3tmp[0];
	dN3ret[1][i] = dN3tmp[1];
	dN3ret[2][i] = dN3tmp[2];
	}
	dN3[0] = fvec::load(&dN3ret[0][0]);
	dN3[1] = fvec::load(&dN3ret[1][0]);
	dN3[2] = fvec::load(&dN3ret[2][0]);
	return fvec::load(&ret[0]);
	}

	static fvec aut_frebo_sum_omega(
	KernelArgsAIREBOT<flt_t,acc_t> * _noalias ka,
	struct aut_frebo_data * _noalias i_data,
	struct aut_frebo_data * _noalias j_data,
	int itype, int jtype,
	ivec vi, ivec vj,
	fvec r23x, fvec r23y, fvec r23z, fvec r23mag,
	fvec VA, fvec fij[3]
	) {
	fvec c_1 = fvec::set1(1);
	fvec c_m1 = fvec::set1(-1);
	fvec c_2 = fvec::set1(2);
	fvec c_m2 = fvec::set1(-2);
	fvec sum_omega = fvec::setzero();
	fvec thmin = fvec::set1(ka->params.thmin);
	fvec thmax = fvec::set1(ka->params.thmax);
	// 2 == i, 3 == j
	fvec r32x = fvec::setzero() - r23x;
	fvec r32y = fvec::setzero() - r23y;
	fvec r32z = fvec::setzero() - r23z;
	int buf_idx_i, buf_idx_j;
	for (buf_idx_i = 0; buf_idx_i < i_data->buf_len; buf_idx_i++) {
	// a1 == k == buf_idx_i
	bvec mask_start = i_data->mask_buf[buf_idx_i];
	fvec r21x = i_data->rikx_buf[buf_idx_i]; // a2 - a1 -> i - k
	fvec r21y = i_data->riky_buf[buf_idx_i]; // a2 - a1 -> i - k
	fvec r21z = i_data->rikz_buf[buf_idx_i]; // a2 - a1 -> i - k
	fvec r21mag = i_data->rikmag_buf[buf_idx_i];
	// TODO use buffered cosjik
	fvec cos321 = (
	r23x * r21x + r23y * r21y + r23z * r21z) / ( r23mag * r21mag);
	cos321 = fvec::min(c_1, fvec::max(c_m1, cos321));
	fvec sin321 = fvec::sqrt(c_1 - cos321 * cos321);
	bvec mask_outer = fvec::cmpneq(fvec::setzero(), sin321) & mask_start;
	// add "continue"
	fvec sink2i = fvec::mask_recip(fvec::undefined(), mask_outer,
	sin321 * sin321);
	fvec rik2i = fvec::mask_recip(fvec::undefined(), mask_outer,
	r21mag * r21mag);
	fvec rr = r23mag * r23mag - r21mag * r21mag;
	fvec r31x = r21x - r23x;
	fvec r31y = r21y - r23y;
	fvec r31z = r21z - r23z;
	fvec r31mag2 = r31x * r31x + r31y * r31y + r31z * r31z;
	fvec rijrik = c_2 * r23mag * r21mag;
	fvec r21mag2 = r21mag * r21mag;
	fvec dctik = fvec::mask_div(fvec::undefined(), mask_outer, r31mag2 - rr,
	rijrik * r21mag2);
	fvec dctij = fvec::mask_div(fvec::undefined(), mask_outer, r31mag2 + rr,
	rijrik * r23mag * r23mag);
	fvec dctjk = fvec::mask_div(fvec::undefined(), mask_outer, c_m2, rijrik);
	fvec dw21 = i_data->dwik_buf[buf_idx_i];
	fvec w21 = i_data->wik_buf[buf_idx_i];
	fvec dtsjik;
	fvec tspjik = aut_Sp2_deriv(cos321, thmin, thmax, &dtsjik);
	dtsjik = fvec::setzero() - dtsjik; // todo replace by appropriate xor.
	ivec k = i_data->k_buf[buf_idx_i];
	for (buf_idx_j = 0; buf_idx_j < j_data->buf_len; buf_idx_j++) {
	// check l == k in second loop.
	// l == a4 == buf_idx_j
	ivec l = j_data->k_buf[buf_idx_j];
	bvec mask_inner_0 = ivec::mask_cmpneq(mask_outer, k, l) &
	j_data->mask_buf[buf_idx_j];
	// add "continue"
	fvec r34x = j_data->rikx_buf[buf_idx_j];
	fvec r34y = j_data->riky_buf[buf_idx_j];
	fvec r34z = j_data->rikz_buf[buf_idx_j];
	fvec r34mag = j_data->rikmag_buf[buf_idx_j];
	fvec cos234 = fvec::mask_div(fvec::undefined(), mask_inner_0,
	r32x * r34x + r32y * r34y + r32z * r34z,
	r23mag * r34mag);
	cos234 = fvec::min(c_1, fvec::max(c_m1, cos234));
	fvec sin234 = fvec::mask_sqrt(fvec::undefined(), mask_inner_0,
	c_1 - cos234 * cos234);
	bvec mask_inner_1 = fvec::mask_cmpneq(mask_inner_0, sin234,
	fvec::setzero());
	// add "continue"
	fvec sinl2i = fvec::mask_recip(fvec::undefined(), mask_inner_1,
	sin234 * sin234);
	fvec rjl2i = fvec::mask_recip(fvec::undefined(), mask_inner_1,
	r34mag * r34mag);
	fvec dw34 = j_data->dwik_buf[buf_idx_j];
	fvec w34 = j_data->wik_buf[buf_idx_j];
	fvec rr = r23mag * r23mag - r34mag * r34mag;
	fvec r24x = r23x + r34x;
	fvec r24y = r23y + r34y;
	fvec r24z = r23z + r34z;
	fvec r242 = r24x * r24x + r24y * r24y + r24z * r24z;
	fvec rijrjl = c_2 * r23mag * r34mag;
	fvec rjl2 = r34mag * r34mag;
	fvec dctjl = fvec::mask_div(fvec::undefined(), mask_inner_1, r242 - rr,
	rijrjl * rjl2);
	fvec dctji = fvec::mask_div(fvec::undefined(), mask_inner_1, r242 + rr,
	rijrjl * r23mag * r23mag);
	fvec dctil = fvec::mask_div(fvec::undefined(), mask_inner_1, c_m2,
	rijrjl);
	fvec dtsijl;
	fvec tspijl = aut_Sp2_deriv(cos234, thmin, thmax, &dtsijl);
	dtsijl = fvec::setzero() - dtsijl;
	fvec prefactor = VA;

	fvec cross321x = r32y * r21z - r32z * r21y;
	fvec cross321y = r32z * r21x - r32x * r21z;
	fvec cross321z = r32x * r21y - r32y * r21x;
	fvec cross234x = r23y * r34z - r23z * r34y;
	fvec cross234y = r23z * r34x - r23x * r34z;
	fvec cross234z = r23x * r34y - r23y * r34x;

	fvec cwnum = cross321x * cross234x + cross321y * cross234y + cross321z *
	cross234z;
	fvec cwnom = r21mag * r34mag * r23mag * r23mag * sin321 * sin234;
	fvec om1234 = fvec::mask_div(fvec::undefined(), mask_inner_1, cwnum,
	cwnom);
	fvec cw = om1234;
	fvec sum_omega_contrib = (c_1 - om1234 * om1234) * w21 * w34 *
	(c_1 - tspjik) * ( c_1 - tspijl);
	sum_omega = fvec::mask_add(sum_omega, mask_inner_1, sum_omega,
	sum_omega_contrib);
	fvec dt1dik = rik2i - dctik * sink2i * cos321;
	fvec dt1djk = fvec::setzero() - dctjk * sink2i * cos321;
	fvec dt1djl = rjl2i - dctjl * sinl2i * cos234;
	fvec dt1dil = fvec::setzero() - dctil * sinl2i * cos234;
	fvec dt1dij = fvec::mask_div(fvec::undefined(), mask_inner_1, c_2,
	r23mag * r23mag) -
	dctij * sink2i * cos321 - dctji * sinl2i * cos234;

	fvec dt2dikx = r23y * cross234z - r23z * cross234y;
	fvec dt2diky = r23z * cross234x - r23x * cross234z;
	fvec dt2dikz = r23x * cross234y - r23y * cross234x;

	fvec dt2djlx = r23z * cross321y - r23y * cross321z;
	fvec dt2djly = r23x * cross321z - r23z * cross321x;
	fvec dt2djlz = r23y * cross321x - r23x * cross321y;

	fvec dt2dijx = r21z * cross234y + r34y * cross321z -
	( r34z * cross321y + r21y * cross234z);
	fvec dt2dijy = r21x * cross234z + r34z * cross321x -
	( r34x * cross321z + r21z * cross234x);
	fvec dt2dijz = r21y * cross234x + r34x * cross321y -
	( r34y * cross321x + r21x * cross234y);

	fvec aa = prefactor * c_2 * fvec::mask_div(fvec::undefined(),
	mask_inner_1, cw, cwnom) *
	w21 * w34 * (c_1 - tspjik) * ( c_1 - tspijl);
	fvec aaa1 = (fvec::setzero() - prefactor) * (c_1 - om1234 * om1234) *
	(c_1 - tspjik) * (c_1 - tspijl);
	fvec aaa2 = (fvec::setzero() - prefactor) * (c_1 - om1234 * om1234) *
	w21 * w34;
	fvec at2 = aa * cwnum;

	fvec fcijpc = aaa2 * dtsjik * dctij * (c_1 - tspijl) + aaa2 * dtsijl *
	dctji * (c_1 - tspjik) - dt1dij * at2;
	fvec fcikpc = aaa2 * dtsjik * dctik * (c_1 - tspijl) - dt1dik * at2;
	fvec fcjlpc = aaa2 * dtsijl * dctjl * (c_1 - tspjik) - dt1djl * at2;
	fvec fcjkpc = aaa2 * dtsjik * dctjk * (c_1 - tspijl) - dt1djk * at2;
	fvec fcilpc = aaa2 * dtsijl * dctil * (c_1 - tspjik) - dt1dil * at2;

	fvec F23x = fcijpc * r23x + aa * dt2dijx;
	fvec F23y = fcijpc * r23y + aa * dt2dijy;
	fvec F23z = fcijpc * r23z + aa * dt2dijz;

	fvec F12x = fcikpc * r21x + aa * dt2dikx;
	fvec F12y = fcikpc * r21y + aa * dt2diky;
	fvec F12z = fcikpc * r21z + aa * dt2dikz;

	fvec F34x = fcjlpc * r34x + aa * dt2djlx;
	fvec F34y = fcjlpc * r34y + aa * dt2djly;
	fvec F34z = fcjlpc * r34z + aa * dt2djlz;

	fvec F31x = fcjkpc * r31x;
	fvec F31y = fcjkpc * r31y;
	fvec F31z = fcjkpc * r31z;

	fvec F24x = fcilpc * r24x;
	fvec F24y = fcilpc * r24y;
	fvec F24z = fcilpc * r24z;

	fvec f1x = fvec::setzero() - ( F12x + F31x);
	fvec f1y = fvec::setzero() - ( F12y + F31y);
	fvec f1z = fvec::setzero() - ( F12z + F31z);
	fvec f2x = F12x + F31x;
	fvec f2y = F12y + F31y;
	fvec f2z = F12z + F31z;
	fvec f3x = F34x + F24x;
	fvec f3y = F34y + F24y;
	fvec f3z = F34z + F24z;
	fvec f4x = fvec::setzero() - ( F34x + F24x);
	fvec f4y = fvec::setzero() - ( F34y + F24y);
	fvec f4z = fvec::setzero() - ( F34z + F24z);

	fij[0] = fvec::mask_add(fij[0], mask_inner_1, fij[0],
	F23x + F24x - F31x);
	fij[1] = fvec::mask_add(fij[1], mask_inner_1, fij[1],
	F23y + F24y - F31y);
	fij[2] = fvec::mask_add(fij[2], mask_inner_1, fij[2],
	F23z + F24z - F31z);

	fvec tmp20 = VA * (c_1 - om1234 * om1234) * (c_1 - tspjik) *
	(c_1 - tspijl) * dw21 * w34 * fvec::mask_recip(fvec::undefined(),
	mask_inner_1, r21mag);
	f2x = f2x - tmp20 * r21x;
	f2y = f2y - tmp20 * r21y;
	f2z = f2z - tmp20 * r21z;
	f1x = f1x + tmp20 * r21x;
	f1y = f1y + tmp20 * r21y;
	f1z = f1z + tmp20 * r21z;

	fvec tmp21 = VA * (c_1 - om1234 * om1234) * (c_1 - tspjik) *
	(c_1 - tspijl) * w21 * dw34 * fvec::mask_recip(fvec::undefined(),
	mask_inner_1, r34mag);
	f3x = f3x - tmp21 * r34x;
	f3y = f3y - tmp21 * r34y;
	f3z = f3z - tmp21 * r34z;
	f4x = f4x + tmp21 * r34x;
	f4y = f4y + tmp21 * r34y;
	f4z = f4z + tmp21 * r34z;

	// 1 == buf_idx_i, 2 == i, 3 == j, 4 == buf_idx_j
	i_data->force_k_x_buf[buf_idx_i] =
	fvec::mask_add(i_data->force_k_x_buf[buf_idx_i],
	mask_inner_1, i_data->force_k_x_buf[buf_idx_i], f1x);
	i_data->force_k_y_buf[buf_idx_i] =
	fvec::mask_add(i_data->force_k_y_buf[buf_idx_i], mask_inner_1,
	i_data->force_k_y_buf[buf_idx_i], f1y);
	i_data->force_k_z_buf[buf_idx_i] =
	fvec::mask_add(i_data->force_k_z_buf[buf_idx_i], mask_inner_1,
	i_data->force_k_z_buf[buf_idx_i], f1z);
	i_data->force_i_x =
	fvec::mask_add(i_data->force_i_x, mask_inner_1, i_data->force_i_x, f2x);
	i_data->force_i_y =
	fvec::mask_add(i_data->force_i_y, mask_inner_1, i_data->force_i_y, f2y);
	i_data->force_i_z =
	fvec::mask_add(i_data->force_i_z, mask_inner_1, i_data->force_i_z, f2z);
	j_data->force_i_x =
	fvec::mask_add(j_data->force_i_x, mask_inner_1, j_data->force_i_x, f3x);
	j_data->force_i_y =
	fvec::mask_add(j_data->force_i_y, mask_inner_1, j_data->force_i_y, f3y);
	j_data->force_i_z =
	fvec::mask_add(j_data->force_i_z, mask_inner_1, j_data->force_i_z, f3z);
	j_data->force_k_x_buf[buf_idx_j] =
	fvec::mask_add(j_data->force_k_x_buf[buf_idx_j], mask_inner_1,
	j_data->force_k_x_buf[buf_idx_j], f4x);
	j_data->force_k_y_buf[buf_idx_j] =
	fvec::mask_add(j_data->force_k_y_buf[buf_idx_j], mask_inner_1,
	j_data->force_k_y_buf[buf_idx_j], f4y);
	j_data->force_k_z_buf[buf_idx_j] =
	fvec::mask_add(j_data->force_k_z_buf[buf_idx_j], mask_inner_1,
	j_data->force_k_z_buf[buf_idx_j], f4z);
	}
	}
	return sum_omega;
	}

	static fvec aut_frebo_pi_dh(
	KernelArgsAIREBOT<flt_t,acc_t> * _noalias ka,
	struct aut_frebo_data * _noalias i_data,
	struct aut_frebo_data * _noalias j_data,
	int itype, int jtype, ivec vi, ivec vj,
	fvec r23x, fvec r23y, fvec r23z, fvec r23mag,
	fvec VA,
	fvec Nij, fvec Nji, fvec Nijconj, fvec NconjtmpI, fvec NconjtmpJ,
	fvec fij[3]
	) {
	fvec c_TOL = fvec::set1(TOL);
	fvec dN3[3];
	fvec Tij = aut_frebo_Tij(ka, itype, jtype, Nij, Nji, Nijconj, &dN3[0]);
	bvec TijgtTOLmask = fvec::cmpnle(fvec::abs(Tij), c_TOL);
	fvec sum_omega = fvec::setzero();
	if (bvec::test_any_set(TijgtTOLmask)) {
	sum_omega = aut_frebo_sum_omega(
	ka, i_data, j_data, itype, jtype, vi, vj,
	r23x, r23y, r23z, r23mag, VA * Tij, fij);
	sum_omega = fvec::mask_blend(TijgtTOLmask, fvec::setzero(), sum_omega);
	aut_frebo_N_spline_force(ka, i_data, itype, jtype, vi, vj, VA * sum_omega,
	dN3[0], dN3[2], NconjtmpI);
	aut_frebo_N_spline_force(ka, j_data, jtype, itype, vj, vi, VA * sum_omega,
	dN3[1], dN3[2], NconjtmpJ);
	}
	return Tij * sum_omega;
	}

	/*
	We can reuse the aut_frebo_data buffers here to do this calculation very
	cheaply.
	*/
	static void aut_torsion_vec(
	KernelArgsAIREBOT<flt_t,acc_t> * ka,
	struct aut_frebo_data * i_data,
	struct aut_frebo_data * j_data,
	ivec i, ivec j, fvec wij, fvec dwij
	) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	flt_t (*epsilonT)[2] = ka->params.epsilonT;
	fvec epsilonT00 = fvec::set1(epsilonT[0][0]);
	fvec epsilonT01 = fvec::set1(epsilonT[0][1]);
	fvec epsilonT10 = fvec::set1(epsilonT[1][0]);
	fvec epsilonT11 = fvec::set1(epsilonT[1][1]);
	fvec thmin = fvec::set1(ka->params.thmin);
	fvec thmax = fvec::set1(ka->params.thmax);

	const fvec c_1_0 = fvec::set1(1.0);
	const fvec c_0_5 = fvec::set1(0.5);
	const fvec c_0_1 = fvec::set1(0.1);
	const fvec c_2_0 = fvec::set1(2.0);
	const fvec c_2_5 = fvec::set1(2.5);
	const fvec c_256_405 = fvec::set1(256.0/405.0);

	fvec del32x = j_data->x_i - i_data->x_i;
	fvec del32y = j_data->y_i - i_data->y_i;
	fvec del32z = j_data->z_i - i_data->z_i;
	fvec rsq = del32x * del32x + del32y * del32y + del32z * del32z;
	fvec r32 = fvec::sqrt(rsq);
	fvec del23x = fvec::setzero() - del32x;
	fvec del23y = fvec::setzero() - del32y;
	fvec del23z = fvec::setzero() - del32z;
	fvec r23 = r32;
	fvec w23 = wij;
	fvec dw23 = dwij;

	for (int buf_idx_i = 0; buf_idx_i < i_data->buf_len; buf_idx_i++) {
	bvec mask_start = i_data->mask_buf[buf_idx_i];
	fvec del21x = i_data->rikx_buf[buf_idx_i]; // a2 - a1 -> i - k
	fvec del21y = i_data->riky_buf[buf_idx_i]; // a2 - a1 -> i - k
	fvec del21z = i_data->rikz_buf[buf_idx_i]; // a2 - a1 -> i - k
	fvec r21 = i_data->rikmag_buf[buf_idx_i];
	fvec cos321 = i_data->cosjik_buf[buf_idx_i];
	fvec sin321 = fvec::sqrt(c_1_0 - cos321 * cos321);
	// strictly equivalent to sin321 < TOL
	mask_start = fvec::mask_cmpneq(mask_start, fvec::setzero(), sin321);
	if (! bvec::test_any_set(mask_start)) continue;

	fvec deljkx = del21x - del23x;
	fvec deljky = del21y - del23y;
	fvec deljkz = del21z - del23z;
	fvec rjk2 = deljkx * deljkx + deljky * deljky + deljkz * deljkz;
	fvec rjk = fvec::sqrt(rjk2);
	fvec rik2 = r21 * r21;
	fvec w21 = i_data->wik_buf[buf_idx_i];
	fvec dw21 = i_data->dwik_buf[buf_idx_i];

	fvec rij = r32;
	fvec rik = r21;
	fvec rij2 = r32 * r32;
	fvec dtsjik;
	fvec tspjik = aut_Sp2_deriv(cos321, thmin, thmax, &dtsjik);
	dtsjik = fvec::setzero() - dtsjik;

	bvec ktype_mask = i_data->ktype_buf[buf_idx_i];
	fvec epsilonT0 = fvec::mask_blend(ktype_mask, epsilonT00, epsilonT10);
	fvec epsilonT1 = fvec::mask_blend(ktype_mask, epsilonT01, epsilonT11);

	ivec k = i_data->k_buf[buf_idx_i];
	for (int buf_idx_j = 0; buf_idx_j < j_data->buf_len; buf_idx_j++) {
	ivec l = j_data->k_buf[buf_idx_j];
	bvec mask_inner_0 = ivec::mask_cmpneq(mask_start, k, l) &
	j_data->mask_buf[buf_idx_j];
	if (! bvec::test_any_set(mask_inner_0)) continue;
	fvec del34x = j_data->rikx_buf[buf_idx_j];
	fvec del34y = j_data->riky_buf[buf_idx_j];
	fvec del34z = j_data->rikz_buf[buf_idx_j];
	fvec r34 = j_data->rikmag_buf[buf_idx_j];
	bvec ltype_mask = j_data->ktype_buf[buf_idx_j];
	fvec cos234 = j_data->cosjik_buf[buf_idx_j];
	fvec sin234 = fvec::sqrt(c_1_0 - cos234 * cos234);
	// strictly equivalent to sin234 < TOL
	mask_inner_0 = fvec::mask_cmpneq(mask_inner_0, sin234, fvec::setzero());
	if (! bvec::test_any_set(mask_inner_0)) continue;
	fvec dw34 = j_data->dwik_buf[buf_idx_j];
	fvec w34 = j_data->wik_buf[buf_idx_j];
	fvec delilx = del23x + del34x;
	fvec delily = del23y + del34y;
	fvec delilz = del23z + del34z;
	fvec ril2 = delilx * delilx + delily * delily + delilz * delilz;
	fvec ril = fvec::sqrt(ril2);
	fvec rjl2 = r34 * r34;

	fvec rjl = r34;
	fvec dtsijl;
	fvec tspijl = aut_Sp2_deriv(cos234, thmin, thmax, &dtsijl);
	dtsijl = fvec::setzero() - dtsijl;
	fvec cross321x = del32y * del21z - del32z * del21y;
	fvec cross321y = del32z * del21x - del32x * del21z;
	fvec cross321z = del32x * del21y - del32y * del21x;
	fvec cross321mag = fvec::sqrt(cross321x * cross321x +
	cross321y * cross321y +
	cross321z * cross321z);
	fvec cross234x = del23y * del34z - del23z * del34y;
	fvec cross234y = del23z * del34x - del23x * del34z;
	fvec cross234z = del23x * del34y - del23y * del34x;
	fvec cross234mag = fvec::sqrt(cross234x * cross234x +
	cross234y * cross234y +
	cross234z * cross234z);
	fvec cwnum = cross321x * cross234x + cross321y * cross234y +
	cross321z * cross234z;
	fvec cwnom = r21 * r34 * r32 * r32 * sin321 * sin234;
	fvec cw = cwnum / cwnom;

	fvec cw2 = c_0_5 * ( c_1_0 - cw);
	fvec ekijl = fvec::mask_blend(ltype_mask, epsilonT0, epsilonT1);
	fvec Ec = c_256_405 * ekijl;
	fvec cw2_5 = cw2 * cw2 * cw2 * cw2 * cw2;
	fvec Vtors = Ec * cw2_5 - ekijl * c_0_1;

	fvec evdwl = Vtors * w21 * w23 * w34 * (c_1_0-tspjik) * (c_1_0-tspijl);
	ka->result_eng += fvec::mask_reduce_add(mask_inner_0, evdwl);

	fvec dndijx = cross234y * del21z - cross234z * del21y;
	fvec dndijy = cross234z * del21x - cross234x * del21z;
	fvec dndijz = cross234x * del21y - cross234y * del21x;

	fvec tmpvecx = del34y * cross321z - del34z * cross321y;
	fvec tmpvecy = del34z * cross321x - del34x * cross321z;
	fvec tmpvecz = del34x * cross321y - del34y * cross321x;

	dndijx = dndijx + tmpvecx;
	dndijy = dndijy + tmpvecy;
	dndijz = dndijz + tmpvecz;

	fvec dndikx = del23y * cross234z - del23z * cross234y;
	fvec dndiky = del23z * cross234x - del23x * cross234z;
	fvec dndikz = del23x * cross234y - del23y * cross234x;

	fvec dndjlx = cross321y * del23z - cross321z * del23y;
	fvec dndjly = cross321z * del23x - cross321x * del23z;
	fvec dndjlz = cross321x * del23y - cross321y * del23x;

	fvec r23sq = r23 * r23;
	fvec r21sq = r21 * r21;
	fvec r34sq = r34 * r34;
	fvec rjksq = rjk * rjk;
	fvec rilsq = ril * ril;
	fvec dcidij = (r23sq - r21sq + rjksq) / ( c_2_0 * r23sq * r21);
	fvec dcidik = (r21sq - r23sq + rjksq) / ( c_2_0 * r21sq * r23);
	fvec dcidjk = fvec::setzero() - rjk / ( r23 * r21);
	fvec dcjdji = (r23sq - r34sq + rilsq) / ( c_2_0 * r23sq * r34);
	fvec dcjdjl = (r34sq - r23sq + rilsq) / ( c_2_0 * r34sq * r23);
	fvec dcjdil = fvec::setzero() - ril / ( r23 * r34);

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

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

	fvec dxidij = r21 * sin321 + r23 * r21 * dsidij;
	fvec dxidik = r23 * sin321 + r23 * r21 * dsidik;
	fvec dxidjk = r23 * r21 * dsidjk;

	fvec dxjdji = r34 * sin234 + r23 * r34 * dsjdji;
	fvec dxjdjl = r23 * sin234 + r23 * r34 * dsjdjl;
	fvec dxjdil = r23 * r34 * dsjdil;

	fvec ddndij = dxidij * cross234mag + cross321mag * dxjdji;
	fvec ddndik = dxidik * cross234mag;
	fvec ddndjk = dxidjk * cross234mag;
	fvec ddndjl = cross321mag * dxjdjl;
	fvec ddndil = cross321mag * dxjdil;
	fvec dcwddn = fvec::setzero() - cwnum / ( cwnom * cwnom);
	fvec dcwdn = fvec::recip(cwnom);
	fvec cw2_4 = cw2 * cw2 * cw2 * cw2;
	fvec dvpdcw = c_2_5 * Ec * cw2_4 * w23 * w21 * w34 * (c_1_0 - tspjik) *
	(c_1_0 - tspijl);

	fvec Ftmpx = dvpdcw * (dcwdn * dndijx + dcwddn * ddndij * del23x / r23);
	fvec Ftmpy = dvpdcw * (dcwdn * dndijy + dcwddn * ddndij * del23y / r23);
	fvec Ftmpz = dvpdcw * (dcwdn * dndijz + dcwddn * ddndij * del23z / r23);
	fvec fix = Ftmpx;
	fvec fiy = Ftmpy;
	fvec fiz = Ftmpz;
	fvec fjx = fvec::setzero() - Ftmpx;
	fvec fjy = fvec::setzero() - Ftmpy;
	fvec fjz = fvec::setzero() - Ftmpz;

	Ftmpx = dvpdcw * (dcwdn * dndikx + dcwddn * ddndik * del21x / r21);
	Ftmpy = dvpdcw * (dcwdn * dndiky + dcwddn * ddndik * del21y / r21);
	Ftmpz = dvpdcw * (dcwdn * dndikz + dcwddn * ddndik * del21z / r21);
	fix = fix + Ftmpx;
	fiy = fiy + Ftmpy;
	fiz = fiz + Ftmpz;
	fvec fkx = fvec::setzero() - Ftmpx;
	fvec fky = fvec::setzero() - Ftmpy;
	fvec fkz = fvec::setzero() - Ftmpz;

	Ftmpx = dvpdcw * dcwddn * ddndjk * deljkx / rjk;
	Ftmpy = dvpdcw * dcwddn * ddndjk * deljky / rjk;
	Ftmpz = dvpdcw * dcwddn * ddndjk * deljkz / rjk;
	fjx = fjx + Ftmpx;
	fjy = fjy + Ftmpy;
	fjz = fjz + Ftmpz;
	fkx = fkx - Ftmpx;
	fky = fky - Ftmpy;
	fkz = fkz - Ftmpz;

	Ftmpx = dvpdcw * (dcwdn * dndjlx + dcwddn * ddndjl * del34x / r34);
	Ftmpy = dvpdcw * (dcwdn * dndjly + dcwddn * ddndjl * del34y / r34);
	Ftmpz = dvpdcw * (dcwdn * dndjlz + dcwddn * ddndjl * del34z / r34);
	fjx = fjx + Ftmpx;
	fjy = fjy + Ftmpy;
	fjz = fjz + Ftmpz;
	fvec flx = fvec::setzero() - Ftmpx;
	fvec fly = fvec::setzero() - Ftmpy;
	fvec flz = fvec::setzero() - Ftmpz;

	Ftmpx = dvpdcw * dcwddn * ddndil * delilx / ril;
	Ftmpy = dvpdcw * dcwddn * ddndil * delily / ril;
	Ftmpz = dvpdcw * dcwddn * ddndil * delilz / ril;
	fix = fix + Ftmpx;
	fiy = fiy + Ftmpy;
	fiz = fiz + Ftmpz;
	flx = flx - Ftmpx;
	fly = fly - Ftmpy;
	flz = flz - Ftmpz;

	// coordination forces

	fvec fpair = Vtors * dw21 * w23 * w34 * (c_1_0 - tspjik) *
	(c_1_0 - tspijl) / r21;
	fix = fix - del21x * fpair;
	fiy = fiy - del21y * fpair;
	fiz = fiz - del21z * fpair;
	fkx = fkx + del21x * fpair;
	fky = fky + del21y * fpair;
	fkz = fkz + del21z * fpair;

	fpair = Vtors * w21 * dw23 * w34 * (c_1_0 - tspjik) * (c_1_0 - tspijl) /
	r23;
	fix = fix - del23x * fpair;
	fiy = fiy - del23y * fpair;
	fiz = fiz - del23z * fpair;
	fjx = fjx + del23x * fpair;
	fjy = fjy + del23y * fpair;
	fjz = fjz + del23z * fpair;

	fpair = Vtors * w21 * w23 * dw34 * (c_1_0 - tspjik) * (c_1_0 - tspijl) /
	r34;
	fjx = fjx - del34x * fpair;
	fjy = fjy - del34y * fpair;
	fjz = fjz - del34z * fpair;
	flx = flx + del34x * fpair;
	fly = fly + del34y * fpair;
	flz = flz + del34z * fpair;

	// additional cut off function forces

	fvec fcpc = fvec::setzero() - Vtors * w21 * w23 * w34 * dtsjik * (c_1_0 -
	tspijl);
	fpair = fcpc * dcidij / rij;
	fix = fix + fpair * del23x;
	fiy = fiy + fpair * del23y;
	fiz = fiz + fpair * del23z;
	fjx = fjx - fpair * del23x;
	fjy = fjy - fpair * del23y;
	fjz = fjz - fpair * del23z;

	fpair = fcpc * dcidik / rik;
	fix = fix + fpair * del21x;
	fiy = fiy + fpair * del21y;
	fiz = fiz + fpair * del21z;
	fkx = fkx - fpair * del21x;
	fky = fky - fpair * del21y;
	fkz = fkz - fpair * del21z;

	fpair = fcpc * dcidjk / rjk;
	fjx = fjx + fpair * deljkx;
	fjy = fjy + fpair * deljky;
	fjz = fjz + fpair * deljkz;
	fkx = fkx - fpair * deljkx;
	fky = fky - fpair * deljky;
	fkz = fkz - fpair * deljkz;

	fcpc = fvec::setzero() - Vtors * w21 * w23 * w34 * (c_1_0 - tspjik) *
	dtsijl;
	fpair = fcpc * dcjdji / rij;
	fix = fix + fpair * del23x;
	fiy = fiy + fpair * del23y;
	fiz = fiz + fpair * del23z;
	fjx = fjx - fpair * del23x;
	fjy = fjy - fpair * del23y;
	fjz = fjz - fpair * del23z;

	fpair = fcpc * dcjdjl / rjl;
	fjx = fjx + fpair * del34x;
	fjy = fjy + fpair * del34y;
	fjz = fjz + fpair * del34z;
	flx = flx - fpair * del34x;
	fly = fly - fpair * del34y;
	flz = flz - fpair * del34z;

	fpair = fcpc * dcjdil / ril;
	fix = fix + fpair * delilx;
	fiy = fiy + fpair * delily;
	fiz = fiz + fpair * delilz;
	flx = flx - fpair * delilx;
	fly = fly - fpair * delily;
	flz = flz - fpair * delilz;

	// sum per-atom forces into atom force array

	i_data->force_i_x = fvec::mask_add(i_data->force_i_x, mask_inner_0,
	i_data->force_i_x, fix);
	i_data->force_i_y = fvec::mask_add(i_data->force_i_y, mask_inner_0,
	i_data->force_i_y, fiy);
	i_data->force_i_z = fvec::mask_add(i_data->force_i_z, mask_inner_0,
	i_data->force_i_z, fiz);
	i_data->force_j_x = fvec::mask_add(i_data->force_j_x, mask_inner_0,
	i_data->force_j_x, fjx);
	i_data->force_j_y = fvec::mask_add(i_data->force_j_y, mask_inner_0,
	i_data->force_j_y, fjy);
	i_data->force_j_z = fvec::mask_add(i_data->force_j_z, mask_inner_0,
	i_data->force_j_z, fjz);
	i_data->force_k_x_buf[buf_idx_i] =
	fvec::mask_add(i_data->force_k_x_buf[buf_idx_i], mask_inner_0,
	i_data->force_k_x_buf[buf_idx_i], fkx);
	i_data->force_k_y_buf[buf_idx_i] =
	fvec::mask_add(i_data->force_k_y_buf[buf_idx_i], mask_inner_0,
	i_data->force_k_y_buf[buf_idx_i], fky);
	i_data->force_k_z_buf[buf_idx_i] =
	fvec::mask_add(i_data->force_k_z_buf[buf_idx_i], mask_inner_0,
	i_data->force_k_z_buf[buf_idx_i], fkz);
	j_data->force_k_x_buf[buf_idx_j] =
	fvec::mask_add(j_data->force_k_x_buf[buf_idx_j], mask_inner_0,
	j_data->force_k_x_buf[buf_idx_j], flx);
	j_data->force_k_y_buf[buf_idx_j] =
	fvec::mask_add(j_data->force_k_y_buf[buf_idx_j], mask_inner_0,
	j_data->force_k_y_buf[buf_idx_j], fly);
	j_data->force_k_z_buf[buf_idx_j] =
	fvec::mask_add(j_data->force_k_z_buf[buf_idx_j], mask_inner_0,
	j_data->force_k_z_buf[buf_idx_j], flz);
	}
	}
	}

	/*
	* Processes VL elements of the same type itype/jtype for REBO and TORSION
	* interactions. This allows us to reuse the aut_frebo_data buffes in the
	* torsion calculaltion.
	*/
	static void aut_frebo_batch_of_kind(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	int torflag, int itype, int jtype,
	int * i_buf, int * j_buf) {
	{ // jump-scope for exceed_limits
	AtomAIREBOT<flt_t> * x = ka->x;
	int * tag = ka->tag;
	int * map = ka->map;
	ResultForceT<acc_t> * result_f = ka->result_f;
	flt_t rcminij = ka->params.rcmin[itype][jtype];
	flt_t rcmaxij = ka->params.rcmax[itype][jtype];
	flt_t Qij = ka->params.Q[itype][jtype];
	flt_t Aij = ka->params.A[itype][jtype];
	flt_t alphaij = ka->params.alpha[itype][jtype];
	fvec vrcminij = fvec::set1(ka->params.rcmin[itype][jtype]);
	fvec vrcmaxij = fvec::set1(ka->params.rcmax[itype][jtype]);
	fvec vQij = fvec::set1(ka->params.Q[itype][jtype]);
	fvec vAij = fvec::set1(ka->params.A[itype][jtype]);
	fvec malphaij = fvec::set1(-ka->params.alpha[itype][jtype]);
	fvec c_1_0 = fvec::set1(1);
	fvec c_0_5 = fvec::set1(0.5);
	fvec c_TOL = fvec::set1(1e-9);
	struct aut_frebo_data i_data, j_data;

	fvec evdwl_vacc = fvec::setzero();
	ivec vi = ivec::maskz_loadu(bvec::full(), i_buf);
	int tmp;
	ivec vj = ivec::maskz_loadu(bvec::full(), j_buf);
	fvec x_i, y_i, z_i;
	fvec x_j, y_j, z_j;
	aut_loadatoms_vec_notype(x, vi, &x_i, &y_i, &z_i);
	aut_loadatoms_vec_notype(x, vj, &x_j, &y_j, &z_j);
	i_data.x_i = x_i;
	i_data.y_i = y_i;
	i_data.z_i = z_i;
	i_data.x_j = x_j;
	i_data.y_j = y_j;
	i_data.z_j = z_j;
	j_data.x_i = x_j;
	j_data.y_i = y_j;
	j_data.z_i = z_j;
	j_data.x_j = x_i;
	j_data.y_j = y_i;
	j_data.z_j = z_i;
	fvec delx = x_i - x_j;
	fvec dely = y_i - y_j;
	fvec delz = z_i - z_j;
	fvec rsq = delx * delx + dely * dely + delz * delz;
	fvec rij = fvec::sqrt(rsq);
	fvec dwij;
	fvec wij = aut_Sp_deriv(rij, vrcminij, vrcmaxij, &dwij);

	fvec exp_alphar = fvec::exp(malphaij * rij);
	fvec Qij_over_rij = vQij / rij;
	fvec Qij_over_rsq = vQij / rsq;
	fvec VR_by_wij = ( c_1_0 + Qij_over_rij) * vAij * exp_alphar;
	fvec VR = wij * VR_by_wij;
	fvec pre = wij * vAij * exp_alphar;
	fvec dVRdi = pre * ( malphaij + malphaij * Qij_over_rij - Qij_over_rsq);
	dVRdi = dVRdi + VR_by_wij * dwij;

	fvec VA_by_wij = fvec::setzero();
	fvec dVA = fvec::setzero();

	int k;
	for (k = 0; k < 3; k++) {
	fvec mBIJc = fvec::set1(-ka->params.BIJc[itype][jtype][k]);
	fvec mBetaij = fvec::set1(-ka->params.Beta[itype][jtype][k]);
	fvec term = mBIJc * fvec::exp(mBetaij * rij);
	VA_by_wij = VA_by_wij + term;
	dVA = dVA + mBetaij * wij * term;
	}

	dVA = dVA + dwij * VA_by_wij;
	fvec VA = wij * VA_by_wij;

	bvec tol_check = fvec::cmplt(wij, c_TOL);
	VA = fvec::mask_blend(tol_check, VA, fvec::setzero());
	dVA = fvec::mask_blend(tol_check, dVA, fvec::setzero());
	VR = fvec::mask_blend(tol_check, VR, fvec::setzero());
	dVRdi = fvec::mask_blend(tol_check, dVRdi, fvec::setzero());

	fvec nHi = fvec::gather(vi, ka->nH, sizeof(flt_t));
	fvec nCi = fvec::gather(vi, ka->nC, sizeof(flt_t));
	fvec nHj = fvec::gather(vj, ka->nH, sizeof(flt_t));
	fvec nCj = fvec::gather(vj, ka->nC, sizeof(flt_t));
	fvec Nij = (nHi + nCi) - wij;
	fvec Nji = (nHj + nCj) - wij;
	i_data.nHi = nHi;
	i_data.nCi = nCi;
	j_data.nHi = nHj;
	j_data.nCi = nCj;
	fvec fij[3], fji[3];
	fij[0] = fvec::setzero(); fij[1] = fvec::setzero();
	fij[2] = fvec::setzero();
	fji[0] = fvec::setzero(); fji[1] = fvec::setzero();
	fji[2] = fvec::setzero();

	fvec NconjtmpI;
	fvec pij = aut_frebo_pij_pd_2(
	ka, &i_data, itype, jtype, vi, vj,
	delx, dely, delz, rij, wij, VA, &NconjtmpI, fij);

	if (i_data.buf_len < 0) goto exceed_limits;

	fvec NconjtmpJ;
	fvec rjix = fvec::setzero() - delx;
	fvec rjiy = fvec::setzero() - dely;
	fvec rjiz = fvec::setzero() - delz;
	fvec pji = aut_frebo_pij_pd_2(
	ka, &j_data, jtype, itype, vj, vi,
	rjix, rjiy, rjiz, rij, wij, VA, &NconjtmpJ, fji);
	fij[0] = fij[0] - fji[0];
	fij[1] = fij[1] - fji[1];
	fij[2] = fij[2] - fji[2];

	if (j_data.buf_len < 0) goto exceed_limits;

	if (torflag && itype == 0 && jtype == 0)
	aut_torsion_vec(ka, &i_data, &j_data, vi, vj, wij, dwij);

	fvec Nijconj = c_1_0 + NconjtmpI * NconjtmpI + NconjtmpJ * NconjtmpJ;
	fvec dN3[3];
	fvec pi_rc = aut_frebo_pi_rc_pd(ka, itype, jtype, Nij, Nji, Nijconj, dN3);
	aut_frebo_N_spline_force(ka, &i_data, itype, jtype, vi, vj, VA, dN3[0],
	dN3[2], NconjtmpI);
	aut_frebo_N_spline_force(ka, &j_data, jtype, itype, vj, vi, VA, dN3[1],
	dN3[2], NconjtmpJ);
	fvec pi_dh = aut_frebo_pi_dh(ka, &i_data, &j_data, itype, jtype, vi, vj,
	delx, dely, delz, rij, VA, Nij, Nji, Nijconj,
	NconjtmpI, NconjtmpJ, fij);

	fvec bij = c_0_5 * ( pij + pji) + pi_rc + pi_dh;
	fvec dVAdi = bij * dVA;
	fvec fpair = (dVAdi + dVRdi) * fvec::recip(rij);
	fvec result_f_j_x = fpair * delx - fij[0];
	fvec result_f_j_y = fpair * dely - fij[1];
	fvec result_f_j_z = fpair * delz - fij[2];
	fvec result_f_i_x = fvec::setzero() - result_f_j_x;
	fvec result_f_i_y = fvec::setzero() - result_f_j_y;
	fvec result_f_i_z = fvec::setzero() - result_f_j_z;
	fvec evdwl = VR + bij * VA;
	evdwl_vacc = evdwl_vacc + evdwl;

	aut_frebo_data_writeback(ka, &i_data);
	aut_frebo_data_writeback(ka, &j_data);

	flt_t fi_x_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fi_y_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fi_z_buf[fvec::VL] __attribute__((aligned(64)));
	int fi_i_buf[ivec::VL] __attribute__((aligned(64)));
	flt_t fj_x_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fj_y_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fj_z_buf[fvec::VL] __attribute__((aligned(64)));
	int fj_j_buf[ivec::VL] __attribute__((aligned(64)));
	flt_t evdwl_buf[fvec::VL] __attribute__((aligned(64)));

	result_f_i_x = i_data.force_i_x + result_f_i_x;
	result_f_i_y = i_data.force_i_y + result_f_i_y;
	result_f_i_z = i_data.force_i_z + result_f_i_z;
	result_f_j_x = i_data.force_j_x + result_f_j_x;
	result_f_j_y = i_data.force_j_y + result_f_j_y;
	result_f_j_z = i_data.force_j_z + result_f_j_z;

	result_f_i_x = j_data.force_j_x + result_f_i_x;
	result_f_i_y = j_data.force_j_y + result_f_i_y;
	result_f_i_z = j_data.force_j_z + result_f_i_z;
	result_f_j_x = j_data.force_i_x + result_f_j_x;
	result_f_j_y = j_data.force_i_y + result_f_j_y;
	result_f_j_z = j_data.force_i_z + result_f_j_z;

	fvec::store(fi_x_buf, result_f_i_x);
	fvec::store(fi_y_buf, result_f_i_y);
	fvec::store(fi_z_buf, result_f_i_z);
	ivec::store(fi_i_buf, vi);
	fvec::store(fj_x_buf, result_f_j_x);
	fvec::store(fj_y_buf, result_f_j_y);
	fvec::store(fj_z_buf, result_f_j_z);
	ivec::store(fj_j_buf, vj);
	fvec::store(evdwl_buf, evdwl);

	int lane;
	for (lane = 0; lane < fvec::VL; lane++) {
	int ii = fi_i_buf[lane];
	result_f[ii].x += fi_x_buf[lane];
	result_f[ii].y += fi_y_buf[lane];
	result_f[ii].z += fi_z_buf[lane];
	result_f[ii].w += 0.5 * evdwl_buf[lane];
	int jj = fj_j_buf[lane];
	result_f[jj].x += fj_x_buf[lane];
	result_f[jj].y += fj_y_buf[lane];
	result_f[jj].z += fj_z_buf[lane];
	result_f[jj].w += 0.5 * evdwl_buf[lane];
	}
	ka->result_eng += fvec::reduce_add(evdwl_vacc);
	return;
	}
	exceed_limits:
	for (int l = 0; l < fvec::VL; l++) {
	int i = i_buf[l];
	int j = j_buf[l];
	ref_frebo_single_interaction(ka, i, j);
	if (torflag && itype == 0 && jtype == 0)
	ref_torsion_single_interaction(ka, i, j);
	}
	}

	/*
	Orders the interactions by itype and jtype and passes chunks to the above
	method.
	*/
	static void aut_frebo(KernelArgsAIREBOT<flt_t,acc_t> * ka, int torflag) {
	AtomAIREBOT<flt_t> * _noalias x = ka->x;
	int * _noalias tag = ka->tag;
	int * _noalias map = ka->map;
	int i_buf[2][2][fvec::VL];
	int j_buf[2][2][fvec::VL];
	int n_buf[2][2] = {0};
	for (int i = ka->frebo_from_atom; i < ka->frebo_to_atom; i++) {
	int itag = tag[i];
	int itype = map[x[i].w];
	flt_t x_i = x[i].x;
	flt_t y_i = x[i].y;
	flt_t z_i = x[i].z;
	int * neighs = ka->neigh_rebo.entries + ka->neigh_rebo.offset[i];
	int jnum = ka->neigh_rebo.num[i];
	for (int jj = 0; jj < jnum; jj++) {
	int j = neighs[jj];
	int jtag = tag[j];
	if (itag > jtag) {
	if (((itag + jtag) & 1) == 0)
	continue;
	} else if (itag < jtag) {
	if (((itag + jtag) & 1) == 1)
	continue;
	} else {
	if (x[j].z < z_i)
	continue;
	if (x[j].z == z_i && x[j].y < y_i)
	continue;
	if (x[j].z == z_i && x[j].y == y_i && x[j].x < x_i)
	continue;
	}
	int jtype = map[x[j].w];
	int ins = n_buf[itype][jtype];
	i_buf[itype][jtype][ins] = i;
	j_buf[itype][jtype][ins] = j;
	n_buf[itype][jtype] += 1;
	if (n_buf[itype][jtype] == fvec::VL) {
	aut_frebo_batch_of_kind(ka, torflag, itype, jtype,
	i_buf[itype][jtype], j_buf[itype][jtype]);
	n_buf[itype][jtype] = 0;
	}
	}
	}
	for (int itype = 0; itype < 2; itype++) {
	for (int jtype = 0; jtype < 2; jtype++) {
	for (int l = 0; l < n_buf[itype][jtype]; l++) {
	int i = i_buf[itype][jtype][l];
	int j = j_buf[itype][jtype][l];
	ref_frebo_single_interaction(ka, i, j);
	if (torflag && itype == 0 && jtype == 0)
	ref_torsion_single_interaction(ka, i, j);
	}
	}
	}
	}

	/*
	* Apply paths in scalar fashion, not crucial for performance.
	*/
	static void aut_airebo_lj_force_path(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	bvec mask, fvec dC, LennardJonesPathAIREBOT<flt_t> path[fvec::VL]) {
	for (int i = 0; i < fvec::VL; i++) {
	if (bvec::test_at(mask, i)) {
	ref_lennard_jones_force_path(ka, fvec::at(dC, i), &path[i]);
	}
	}
	}

	/*
	* Hash-Map for efficient calculation of C_ij.
	* Can have up to ITEMS entries with associated paths, as well as
	* 1024 entries. Open addressing, invalidation by using a different i.
	* Only needs to be reset once per timestep.
	*/
	static const int OPT_TEST_PATH_SIZE = 1024;
	static const int OPT_TEST_PATH_ITEMS = 128;
	struct aut_airebo_lj_test_path_result_data {
	LennardJonesPathAIREBOT<flt_t> testpath[OPT_TEST_PATH_ITEMS];
	int i[OPT_TEST_PATH_SIZE];
	int j[OPT_TEST_PATH_SIZE];
	flt_t cij[OPT_TEST_PATH_SIZE];
	int testpath_idx[OPT_TEST_PATH_SIZE];
	};
	static const unsigned int OPT_TEST_PATH_HASH = 2654435761;

	static int aut_lj_tap_hash_fn(int j, int attempt) {
	uint32_t result = j;
	result *= (uint32_t) OPT_TEST_PATH_HASH;
	result += (uint32_t) attempt;
	result %= (uint32_t) OPT_TEST_PATH_SIZE;
	return result;
	}

	static ivec aut_airebo_lj_tap_hash_fn_vec(ivec val, ivec attempt) {
	const ivec golden = ivec::set1(OPT_TEST_PATH_HASH);
	const ivec mask = ivec::set1(OPT_TEST_PATH_SIZE - 1);
	ivec a = ivec::mullo(golden, val);
	ivec b = a + attempt;
	ivec c = ivec::the_and(b, mask);
	return c;
	}

	/*
	* Enter all those (potential) neighbors of i (including 2nd and 3rd degree)
	* into the hash-map. There is no good way to vectorize this, and it does not
	* seem time-critical.
	*/
	static bool aut_airebo_lj_test_all_paths(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	int i, struct aut_airebo_lj_test_path_result_data * result) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * map = ka->map;
	flt_t (*rcmin)[2] = &ka->params.rcmin[0];
	flt_t (*rcmax)[2] = &ka->params.rcmax[0];
	flt_t rcminsq[2][2];
	rcminsq[0][0] = rcmin[0][0] * rcmin[0][0];
	rcminsq[0][1] = rcmin[0][1] * rcmin[0][1];
	rcminsq[1][0] = rcmin[1][0] * rcmin[1][0];
	rcminsq[1][1] = rcmin[1][1] * rcmin[1][1];
	int * neighs_i = &ka->neigh_rebo.entries[ka->neigh_rebo.offset[i]];
	int itype = map[x[i].w];
	int path_insert_pos = 0;
	for (int jj = 0; jj < ka->neigh_rebo.num[i]; jj++) {
	int j = neighs_i[jj];
	int jtype = map[x[j].w];
	flt_t dijx = x[j].x - x[i].x;
	flt_t dijy = x[j].y - x[i].y;
	flt_t dijz = x[j].z - x[i].z;
	flt_t rijsq = dijx * dijx + dijy * dijy + dijz * dijz;
	flt_t wj = 1, dwj = 0;
	flt_t rij = 0;
	if (rijsq >= rcminsq[itype][jtype]) {
	rij = overloaded::sqrt(rijsq);
	wj = Sp(rij, rcmin[itype][jtype], rcmax[itype][jtype], &dwj);
	}
	int attempt = 0;
	int start_hash_slot = aut_lj_tap_hash_fn(j, attempt);
	int hash_slot = start_hash_slot;
	while (result->i[hash_slot] == i && result->j[hash_slot] != j &&
	attempt < OPT_TEST_PATH_SIZE) {
	hash_slot = aut_lj_tap_hash_fn(j, ++attempt);
	}
	if (attempt >= OPT_TEST_PATH_SIZE) goto exceed_limits;
	bool init_slot = result->i[hash_slot] != i;
	if (init_slot \|\| (1 - wj < result->cij[hash_slot])) {
	result->i[hash_slot] = i;
	result->j[hash_slot] = j;
	result->cij[hash_slot] = 1 - wj;
	if (wj != 1.0) {
	if (path_insert_pos >= OPT_TEST_PATH_ITEMS) goto exceed_limits;
	result->testpath_idx[hash_slot] = path_insert_pos;
	LennardJonesPathAIREBOT<flt_t> *path =
	&result->testpath[path_insert_pos++];
	path->num = 2;
	path->del[0].x = dijx;
	path->del[0].y = dijy;
	path->del[0].z = dijz;
	if (rij == 0) rij = sqrt(rijsq);
	path->r[0] = rij;
	path->w[0] = wj;
	path->dw[0] = dwj;
	path->idx[0] = i;
	path->idx[1] = j;
	}
	}
	int * neighs_j = &ka->neigh_rebo.entries[ka->neigh_rebo.offset[j]];
	for (int kk = 0; kk < ka->neigh_rebo.num[j]; kk++) {
	int k = neighs_j[kk];
	if (k == i) continue;
	int ktype = map[x[k].w];
	flt_t djkx = x[k].x - x[j].x;
	flt_t djky = x[k].y - x[j].y;
	flt_t djkz = x[k].z - x[j].z;
	flt_t rjksq = djkx * djkx + djky * djky + djkz * djkz;
	flt_t wk = 1, dwk = 0;
	flt_t rjk = 0;
	if (rjksq >= rcminsq[jtype][ktype]) {
	rjk = overloaded::sqrt(rjksq);
	wk = Sp(rjk, rcmin[jtype][ktype], rcmax[jtype][ktype], &dwk);
	}
	int attempt = 0;
	int start_hash_slot = aut_lj_tap_hash_fn(k, attempt);
	int hash_slot = start_hash_slot;
	while (result->i[hash_slot] == i && result->j[hash_slot] != k &&
	attempt < OPT_TEST_PATH_SIZE) {
	hash_slot = aut_lj_tap_hash_fn(k, ++attempt);
	}
	if (attempt >= OPT_TEST_PATH_SIZE) goto exceed_limits;
	bool init_slot = result->i[hash_slot] != i;
	if (init_slot \|\| (1 - wj * wk < result->cij[hash_slot])) {
	result->i[hash_slot] = i;
	result->j[hash_slot] = k;
	result->cij[hash_slot] = 1 - wj * wk;
	if (wj * wk != 1.0) {
	if (path_insert_pos >= OPT_TEST_PATH_ITEMS) goto exceed_limits;
	result->testpath_idx[hash_slot] = path_insert_pos;
	LennardJonesPathAIREBOT<flt_t> *path =
	&result->testpath[path_insert_pos++];
	path->num = 3;
	path->del[0].x = dijx;
	path->del[0].y = dijy;
	path->del[0].z = dijz;
	if (rij == 0) rij = sqrt(rijsq);
	path->r[0] = rij;
	path->del[1].x = djkx;
	path->del[1].y = djky;
	path->del[1].z = djkz;
	if (rjk == 0) rjk = sqrt(rjksq);
	path->r[1] = rjk;
	path->w[0] = wj;
	path->dw[0] = dwj;
	path->w[1] = wk;
	path->dw[1] = dwk;
	path->idx[0] = i;
	path->idx[1] = j;
	path->idx[2] = k;
	}
	}
	int * neighs_k = &ka->neigh_rebo.entries[ka->neigh_rebo.offset[k]];
	for (int ll = 0; ll < ka->neigh_rebo.num[k]; ll++) {
	int l = neighs_k[ll];
	if ((l == i) \|\| (l == j)) continue;
	int ltype = map[x[l].w];
	flt_t dklx = x[l].x - x[k].x;
	flt_t dkly = x[l].y - x[k].y;
	flt_t dklz = x[l].z - x[k].z;
	flt_t rklsq = dklx * dklx + dkly * dkly + dklz * dklz;
	flt_t wl = 1, dwl = 0;
	flt_t rkl = 0;
	if (rklsq >= rcminsq[ktype][ltype]) {
	rkl = overloaded::sqrt(rklsq);
	wl = Sp(rkl, rcmin[ktype][ltype], rcmax[ktype][ltype], &dwl);
	}
	int attempt = 0;
	int start_hash_slot = aut_lj_tap_hash_fn(l, attempt);
	int hash_slot = start_hash_slot;
	while (result->i[hash_slot] == i && result->j[hash_slot] != l &&
	attempt < OPT_TEST_PATH_SIZE) {
	hash_slot = aut_lj_tap_hash_fn(l, ++attempt);
	}
	if (attempt >= OPT_TEST_PATH_SIZE) goto exceed_limits;
	bool init_slot = result->i[hash_slot] != i;
	if (init_slot \|\| (1 - wj * wk * wl < result->cij[hash_slot])) {
	result->i[hash_slot] = i;
	result->j[hash_slot] = l;
	result->cij[hash_slot] = 1 - wj * wk * wl;
	if (wj * wk * wl != 1.0) {
	if (path_insert_pos >= OPT_TEST_PATH_ITEMS) goto exceed_limits;
	result->testpath_idx[hash_slot] = path_insert_pos;
	LennardJonesPathAIREBOT<flt_t> *path =
	&result->testpath[path_insert_pos++];
	path->num = 4;
	path->del[0].x = dijx;
	path->del[0].y = dijy;
	path->del[0].z = dijz;
	if (rij == 0) rij = sqrt(rijsq);
	path->r[0] = rij;
	path->del[1].x = djkx;
	path->del[1].y = djky;
	path->del[1].z = djkz;
	if (rjk == 0) rjk = sqrt(rjksq);
	path->r[1] = rjk;
	path->del[2].x = dklx;
	path->del[2].y = dkly;
	path->del[2].z = dklz;
	if (rkl == 0) rkl = sqrt(rklsq);
	path->r[2] = rkl;
	path->w[0] = wj;
	path->dw[0] = dwj;
	path->w[1] = wk;
	path->dw[1] = dwk;
	path->w[2] = wl;
	path->dw[2] = dwl;
	path->idx[0] = i;
	path->idx[1] = j;
	path->idx[2] = k;
	path->idx[3] = l;
	}
	}
	}
	}
	}
	return true;
	exceed_limits:
	return false;
	}

	/*
	* Attempt to look up an element in the hash-map.
	*/
	static fvec aut_airebo_lj_tap_test_path(KernelArgsAIREBOT<flt_t,acc_t> * ka,
	struct aut_airebo_lj_test_path_result_data * test_path_result,
	bvec need_search, ivec i_bc, ivec j,
	LennardJonesPathAIREBOT<flt_t> path[fvec::VL]
	) {
	const ivec c_i1 = ivec::set1(1);
	fvec cij = fvec::set1(1.0);
	// first round: hash all j
	// lookup i/j in hash list.
	// if i matches and j matches: congrats
	// if i matches and j does not: look up attempts
	// if attempts > current_attempts:
	// do another round of hashing
	// for all those found:

	// fill in the path
	// -----------------------------------------------
	// find all the correct hash slots, and a mask of where found.
	ivec attempt = ivec::setzero();
	ivec hash_slot = aut_airebo_lj_tap_hash_fn_vec(j, attempt);
	ivec lookup_i = ivec::mask_gather(ivec::undefined(), need_search, hash_slot,
	&test_path_result->i[0], sizeof(int));
	bvec correct_i = ivec::mask_cmpeq(need_search, lookup_i, i_bc);
	ivec lookup_j = ivec::mask_gather(ivec::undefined(), correct_i, hash_slot,
	&test_path_result->j[0], sizeof(int));
	bvec found_items = ivec::mask_cmpeq(correct_i, lookup_j, j);
	bvec another_attempt = correct_i & ~ found_items;
	while (bvec::test_any_set(another_attempt)) {
	attempt = ivec::mask_add(attempt, another_attempt, attempt, c_i1);
	hash_slot = aut_airebo_lj_tap_hash_fn_vec(j, attempt);
	ivec lookup_i_2 = ivec::mask_gather(lookup_i, another_attempt, hash_slot,
	&test_path_result->i[0], sizeof(int));
	lookup_i = lookup_i_2;
	correct_i = ivec::mask_cmpeq(need_search, lookup_i, i_bc);
	lookup_j = ivec::mask_gather(lookup_j, another_attempt, hash_slot,
	&test_path_result->j[0], sizeof(int));
	found_items = ivec::mask_cmpeq(correct_i, lookup_j, j);
	another_attempt = correct_i & ~ found_items;
	}
	cij = fvec::mask_gather(cij, found_items, hash_slot,
	&test_path_result->cij[0], sizeof(flt_t));
	bvec need_testpath = fvec::mask_cmplt(found_items, fvec::setzero(), cij);
	if (bvec::test_any_set(need_testpath)) {
	for (int i = 0; i < fvec::VL; i++) {
	if (bvec::test_at(need_testpath, i)) {
	int testpath_idx =
	test_path_result->testpath_idx[ivec::at(hash_slot, i)];
	path[i] = test_path_result->testpath[testpath_idx];
	}
	}
	}
	return cij;
	}

	/*
	* This function calculates the Lennard-Jones interaciton for those
	* elements that require a bond-order calculation.
	* It is similarly structured as the aut_frebo_batch_of_kind function.
	* The forces due to bondorders are calculated speculatively and later
	* updated with the correct outer derivative.
	*/
	template<int MORSEFLAG>
	static void aut_lj_with_bo(
	KernelArgsAIREBOT<flt_t,acc_t> * ka,
	int itype, int jtype,
	ivec i, ivec j,
	fvec cij, LennardJonesPathAIREBOT<flt_t> testpath[fvec::VL]
	) {
	{ // jump-scope for exceed_limits
	AtomAIREBOT<flt_t> * _noalias x = ka->x;
	ResultForceT<acc_t> * result_f = ka->result_f;

	ivec c_i4 = ivec::set1(4);
	fvec c_1_0 = fvec::set1(1.0);
	fvec c_2_0 = fvec::set1(2.0);
	fvec c_0_5 = fvec::set1(0.5);

	fvec x_i, y_i, z_i;
	aut_loadatoms_vec_notype(x, i, &x_i, &y_i, &z_i);
	fvec x_j, y_j, z_j;
	aut_loadatoms_vec_notype(x, j, &x_j, &y_j, &z_j);
	fvec delx = x_i - x_j;
	fvec dely = y_i - y_j;
	fvec delz = z_i - z_j;
	fvec rsq = delx * delx + dely * dely + delz * delz;

	fvec rij = fvec::sqrt(rsq);
	bvec need_path_force = fvec::cmplt(cij, c_1_0);
	flt_t sigcut = ka->params.sigcut;
	flt_t sigmin = ka->params.sigmin;
	flt_t sigma = ka->params.sigma[itype][jtype];
	flt_t rljmax = sigcut * sigma;
	flt_t rljmin = sigmin * sigma;
	fvec p_rljmin = fvec::set1(rljmin);
	fvec p_rljmax = fvec::set1(rljmax);

	fvec dslw, slw = aut_Sp2_deriv(rij, p_rljmin, p_rljmax, &dslw);

	fvec p_lj1 = fvec::set1(ka->params.lj1[itype][jtype]);
	fvec p_lj2 = fvec::set1(ka->params.lj2[itype][jtype]);
	fvec p_lj3 = fvec::set1(ka->params.lj3[itype][jtype]);
	fvec p_lj4 = fvec::set1(ka->params.lj4[itype][jtype]);

	fvec r2inv = fvec::recip(rsq);

	fvec vdw, dvdw;
	if (MORSEFLAG) {
	fvec exr = fvec::exp(fvec::setzero() - rij * p_lj4);
	vdw = p_lj1 * exr * (p_lj2 * exr - c_2_0);
	dvdw = p_lj3 * exr * (c_1_0 - p_lj2 * exr);
	} else {
	fvec r6inv = r2inv * r2inv * r2inv;

	vdw = r6inv * ( p_lj3 * r6inv - p_lj4);
	fvec r7inv = r6inv * rij * r2inv;
	dvdw = r7inv * ( p_lj2 - p_lj1 * r6inv);
	}

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

	fvec p_rcLJmin = fvec::set1(ka->params.rcLJmin[itype][jtype]);
	fvec p_rcLJmax = fvec::set1(ka->params.rcLJmax[itype][jtype]);
	fvec dStr, Str = aut_Sp2_deriv(rij, p_rcLJmin, p_rcLJmax, &dStr);
	fvec VA = cij * VLJ * Str;

	fvec fij[3], fji[3];
	fij[0] = fvec::setzero(); fij[1] = fvec::setzero();
	fij[2] = fvec::setzero();
	fji[0] = fvec::setzero(); fji[1] = fvec::setzero();
	fji[2] = fvec::setzero();

	ivec vi = i;
	ivec vj = j;

	struct aut_frebo_data i_data, j_data;
	i_data.x_i = x_i;
	i_data.y_i = y_i;
	i_data.z_i = z_i;
	i_data.x_j = x_j;
	i_data.y_j = y_j;
	i_data.z_j = z_j;
	j_data.x_i = x_j;
	j_data.y_i = y_j;
	j_data.z_i = z_j;
	j_data.x_j = x_i;
	j_data.y_j = y_i;
	j_data.z_j = z_i;

	fvec p_rcmin = fvec::set1(ka->params.rcmin[itype][jtype]);
	fvec p_rcmax = fvec::set1(ka->params.rcmax[itype][jtype]);
	fvec dwij;
	fvec wij = aut_Sp_deriv(rij, p_rcmin, p_rcmax, &dwij);

	fvec nHi = fvec::gather(vi, ka->nH, sizeof(flt_t));
	fvec nCi = fvec::gather(vi, ka->nC, sizeof(flt_t));
	fvec nHj = fvec::gather(vj, ka->nH, sizeof(flt_t));
	fvec nCj = fvec::gather(vj, ka->nC, sizeof(flt_t));
	fvec Nij = nHi + nCi - wij;
	fvec Nji = nHj + nCj - wij;
	i_data.nHi = nHi;
	i_data.nCi = nCi;
	j_data.nHi = nHj;
	j_data.nCi = nCj;

	fvec the_r = fvec::set1(ka->params.rcmin[itype][jtype]);
	fvec scale = the_r / rij;

	fvec NconjtmpI;
	fvec pij = aut_frebo_pij_pd_2(ka, &i_data, itype, jtype, vi, vj,
	delx * scale, dely * scale, delz * scale,
	the_r, wij, VA, &NconjtmpI, fij);

	if (i_data.buf_len < 0) goto exceed_limits;

	fvec NconjtmpJ;
	fvec rjix = fvec::setzero() - delx;
	fvec rjiy = fvec::setzero() - dely;
	fvec rjiz = fvec::setzero() - delz;
	fvec pji = aut_frebo_pij_pd_2(ka, &j_data, jtype, itype, vj, vi,
	rjix * scale, rjiy * scale, rjiz * scale,
	the_r, wij, VA, &NconjtmpJ, fji);
	fij[0] = fij[0] - fji[0];
	fij[1] = fij[1] - fji[1];
	fij[2] = fij[2] - fji[2];

	if (j_data.buf_len < 0) goto exceed_limits;

	fvec Nijconj = c_1_0 + NconjtmpI * NconjtmpI + NconjtmpJ * NconjtmpJ;
	fvec dN3[3];
	fvec pi_rc = aut_frebo_pi_rc_pd(ka, itype, jtype, Nij, Nji, Nijconj, dN3);

	fvec c_TOL = fvec::set1(TOL);
	fvec dN3_dh[3];
	fvec Tij = aut_frebo_Tij(ka, itype, jtype, Nij, Nji, Nijconj, &dN3_dh[0]);
	bvec TijgtTOLmask = fvec::cmpnle(fvec::abs(Tij), c_TOL);
	fvec sum_omega = fvec::setzero();
	if (bvec::test_any_set(TijgtTOLmask)) {
	sum_omega = aut_frebo_sum_omega(
	ka, &i_data, &j_data, itype, jtype, vi, vj,
	delx * scale, dely * scale, delz * scale, the_r, VA * Tij, fij);
	sum_omega = fvec::mask_blend(TijgtTOLmask, fvec::setzero(), sum_omega);
	}
	fvec pi_dh = Tij * sum_omega;

	fvec bij = c_0_5 * ( pij + pji) + pi_rc + pi_dh;

	fvec p_bLJmin = fvec::set1(ka->params.bLJmin[itype][jtype]);
	fvec p_bLJmax = fvec::set1(ka->params.bLJmax[itype][jtype]);
	fvec dStb, Stb = aut_Sp2_deriv(bij, p_bLJmin, p_bLJmax, &dStb);

	bvec need_bo_deriv = fvec::cmpneq(dStb, fvec::setzero());
	// fix up j_data, i_data, fij:
	// multiply each by dStb
	if (bvec::test_any_set(need_bo_deriv)) {
	i_data.force_i_x = dStb * i_data.force_i_x;
	i_data.force_i_y = dStb * i_data.force_i_y;
	i_data.force_i_z = dStb * i_data.force_i_z;
	i_data.force_j_x = dStb * i_data.force_j_x;
	i_data.force_j_y = dStb * i_data.force_j_y;
	i_data.force_j_z = dStb * i_data.force_j_z;
	j_data.force_i_x = dStb * j_data.force_i_x;
	j_data.force_i_y = dStb * j_data.force_i_y;
	j_data.force_i_z = dStb * j_data.force_i_z;
	j_data.force_j_x = dStb * j_data.force_j_x;
	j_data.force_j_y = dStb * j_data.force_j_y;
	j_data.force_j_z = dStb * j_data.force_j_z;
	for (int k = 0; k < i_data.buf_len; k++) {
	i_data.force_k_x_buf[k] = dStb * i_data.force_k_x_buf[k];
	i_data.force_k_y_buf[k] = dStb * i_data.force_k_y_buf[k];
	i_data.force_k_z_buf[k] = dStb * i_data.force_k_z_buf[k];
	}
	for (int k = 0; k < j_data.buf_len; k++) {
	j_data.force_k_x_buf[k] = dStb * j_data.force_k_x_buf[k];
	j_data.force_k_y_buf[k] = dStb * j_data.force_k_y_buf[k];
	j_data.force_k_z_buf[k] = dStb * j_data.force_k_z_buf[k];
	}
	fvec fijc[3];
	fijc[0] = dStb * fij[0];
	fijc[1] = dStb * fij[1];
	fijc[2] = dStb * fij[2];
	fij[0] = scale * (fijc[0] - (delx * delx * fijc[0] + dely * delx *
	fijc[1] + delz * delx * fijc[2]) / rsq);
	fij[1] = scale * (fijc[1] - (delx * dely * fijc[0] + dely * dely *
	fijc[1] + delz * dely * fijc[2]) / rsq);
	fij[2] = scale * (fijc[2] - (delx * delz * fijc[0] + dely * delz *
	fijc[1] + delz * delz * fijc[2]) / rsq);

	aut_frebo_N_spline_force(ka, &i_data, itype, jtype, vi, vj, dStb * VA,
	dN3[0], dN3[2], NconjtmpI);
	aut_frebo_N_spline_force(ka, &j_data, jtype, itype, vj, vi, dStb * VA,
	dN3[1], dN3[2], NconjtmpJ);
	if (bvec::test_any_set(TijgtTOLmask)) {
	aut_frebo_N_spline_force(ka, &i_data, itype, jtype, vi, vj,
	dStb * VA * sum_omega, dN3_dh[0], dN3_dh[2],
	NconjtmpI);
	aut_frebo_N_spline_force(ka, &j_data, jtype, itype, vj, vi,
	dStb * VA * sum_omega, dN3_dh[1], dN3_dh[2],
	NconjtmpJ);
	}

	aut_frebo_data_writeback(ka, &i_data);
	aut_frebo_data_writeback(ka, &j_data);
	} else {
	fij[0] = fvec::setzero();
	fij[1] = fvec::setzero();
	fij[2] = fvec::setzero();
	}

	fvec fpdVLJ = cij * dVLJ * ( c_1_0 + Str * ( Stb - c_1_0));
	fvec fpdStr = dStr * cij * ( Stb * VLJ - VLJ);
	fvec fpair = r2inv * rij * ( fvec::setzero() - ( fpdVLJ + fpdStr));
	fvec evdwl = VA * Stb + cij * VLJ * ( c_1_0 - Str);

	fvec result_f_i_x = fpair * delx + fij[0];
	fvec result_f_i_y = fpair * dely + fij[1];
	fvec result_f_i_z = fpair * delz + fij[2];
	fvec result_f_j_x = fvec::setzero() - result_f_i_x;
	fvec result_f_j_y = fvec::setzero() - result_f_i_y;
	fvec result_f_j_z = fvec::setzero() - result_f_i_z;

	flt_t fi_x_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fi_y_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fi_z_buf[fvec::VL] __attribute__((aligned(64)));
	int fi_i_buf[ivec::VL] __attribute__((aligned(64)));
	flt_t fj_x_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fj_y_buf[fvec::VL] __attribute__((aligned(64)));
	flt_t fj_z_buf[fvec::VL] __attribute__((aligned(64)));
	int fj_j_buf[ivec::VL] __attribute__((aligned(64)));
	flt_t evdwl_buf[fvec::VL] __attribute__((aligned(64)));

	if (bvec::test_any_set(need_bo_deriv)) {
	result_f_i_x = i_data.force_i_x + result_f_i_x;
	result_f_i_y = i_data.force_i_y + result_f_i_y;
	result_f_i_z = i_data.force_i_z + result_f_i_z;
	result_f_j_x = i_data.force_j_x + result_f_j_x;
	result_f_j_y = i_data.force_j_y + result_f_j_y;
	result_f_j_z = i_data.force_j_z + result_f_j_z;

	result_f_i_x = j_data.force_j_x + result_f_i_x;
	result_f_i_y = j_data.force_j_y + result_f_i_y;
	result_f_i_z = j_data.force_j_z + result_f_i_z;
	result_f_j_x = j_data.force_i_x + result_f_j_x;
	result_f_j_y = j_data.force_i_y + result_f_j_y;
	result_f_j_z = j_data.force_i_z + result_f_j_z;
	}

	fvec::store(fi_x_buf, result_f_i_x);
	fvec::store(fi_y_buf, result_f_i_y);
	fvec::store(fi_z_buf, result_f_i_z);
	ivec::store(fi_i_buf, vi);
	fvec::store(fj_x_buf, result_f_j_x);
	fvec::store(fj_y_buf, result_f_j_y);
	fvec::store(fj_z_buf, result_f_j_z);
	ivec::store(fj_j_buf, vj);
	fvec::store(evdwl_buf, evdwl);

	int lane;
	for (lane = 0; lane < fvec::VL; lane++) {
	int ii = fi_i_buf[lane];
	result_f[ii].x += fi_x_buf[lane];
	result_f[ii].y += fi_y_buf[lane];
	result_f[ii].z += fi_z_buf[lane];
	result_f[ii].w += 0.5 * evdwl_buf[lane];
	int jj = fj_j_buf[lane];
	result_f[jj].x += fj_x_buf[lane];
	result_f[jj].y += fj_y_buf[lane];
	result_f[jj].z += fj_z_buf[lane];
	result_f[jj].w += 0.5 * evdwl_buf[lane];
	}
	ka->result_eng += fvec::reduce_add(evdwl);

	if (bvec::test_any_set(need_path_force)) {
	fvec dC = VLJ * ( Str * Stb + c_1_0 - Str);
	aut_airebo_lj_force_path(ka, need_path_force, dC, testpath);
	}
	return;
	}
	exceed_limits:
	for (int l = 0; l < fvec::VL; l++) {
	ref_lennard_jones_single_interaction(ka, ivec::at(i, l), ivec::at(j, l),
	MORSEFLAG);
	}
	return;
	}

	/*
	* Calculate the lennard-jones interaction.
	* Uses the above hash-map, and outlines the calculation if the bondorder is
	* needed.
	* Agressively compresses to get the most values calculated.
	*/
	template<int MORSEFLAG>
	static void aut_lennard_jones(KernelArgsAIREBOT<flt_t,acc_t> * ka) {
	AtomAIREBOT<flt_t> * x = ka->x;
	int * tag = ka->tag;
	int * map = ka->map;
	ResultForceT<acc_t> * result_f = ka->result_f;
	ivec c_i1 = ivec::set1(1);
	ivec c_i4 = ivec::set1(4);
	fvec c_1_0 = fvec::set1(1.0);
	fvec c_2_0 = fvec::set1(2.0);
	fvec c_0_0 = fvec::set1(0.0);
	int map_i_scalar = 0;
	{
	int i;
	for (i = 1; i < ka->num_types; i++) {
	if (ka->map[i])
	map_i_scalar \|= (1 << i);
	}
	}
	ivec map_i = ivec::set1(map_i_scalar);
	fvec result_eng = fvec::setzero();

	struct aut_airebo_lj_test_path_result_data test_path_result;
	for (int i = 0; i < OPT_TEST_PATH_SIZE; i++) {
	test_path_result.i[i] = -1;
	}

	ivec i_bo[2][2];
	ivec j_bo[2][2];
	fvec cij_bo[2][2];
	LennardJonesPathAIREBOT<flt_t> testpath_bo[2][2][fvec::VL];
	int num_bo[2][2] = {0};

	for (int i = ka->frebo_from_atom; i < ka->frebo_to_atom; i++) {
	ivec itag_bc = ivec::set1(tag[i]);
	int itype = map[x[i].w];
	fvec x_i = fvec::set1(x[i].x);
	fvec y_i = fvec::set1(x[i].y);
	fvec z_i = fvec::set1(x[i].z);
	ivec i_bc = ivec::set1(i);

	fvec cutljsq0 = fvec::set1(ka->params.cutljsq[itype][0]);
	fvec cutljsq1 = fvec::set1(ka->params.cutljsq[itype][1]);
	fvec p_rcmax0 = fvec::set1(ka->params.rcmax[itype][0]);
	fvec p_rcmax1 = fvec::set1(ka->params.rcmax[itype][1]);
	flt_t sigcut = ka->params.sigcut;
	flt_t sigmin = ka->params.sigmin;
	flt_t sigma0 = ka->params.sigma[itype][0];
	flt_t rljmax0 = sigcut * sigma0;
	flt_t rljmin0 = sigmin * sigma0;
	flt_t sigma1 = ka->params.sigma[itype][1];
	flt_t rljmax1 = sigcut * sigma1;
	flt_t rljmin1 = sigmin * sigma1;
	fvec p_rljmax0 = fvec::set1(rljmax0);
	fvec p_rljmax1 = fvec::set1(rljmax1);
	fvec p_rljmin0 = fvec::set1(rljmin0);
	fvec p_rljmin1 = fvec::set1(rljmin1);
	fvec p_rcLJmax0 = fvec::set1(ka->params.rcLJmax[itype][0]);
	fvec p_rcLJmax1 = fvec::set1(ka->params.rcLJmax[itype][1]);
	fvec p_rcLJmin0 = fvec::set1(ka->params.rcLJmin[itype][0]);
	fvec p_rcLJmin1 = fvec::set1(ka->params.rcLJmin[itype][1]);
	fvec p_lj10 = fvec::set1(ka->params.lj1[itype][0]);
	fvec p_lj11 = fvec::set1(ka->params.lj1[itype][1]);
	fvec p_lj20 = fvec::set1(ka->params.lj2[itype][0]);
	fvec p_lj21 = fvec::set1(ka->params.lj2[itype][1]);
	fvec p_lj30 = fvec::set1(ka->params.lj3[itype][0]);
	fvec p_lj31 = fvec::set1(ka->params.lj3[itype][1]);
	fvec p_lj40 = fvec::set1(ka->params.lj4[itype][0]);
	fvec p_lj41 = fvec::set1(ka->params.lj4[itype][1]);

	int * neighs = ka->neigh_lmp.entries + ka->neigh_lmp.offset[i];
	int jnum = ka->neigh_lmp.num_half[i];

	bool tap_success = aut_airebo_lj_test_all_paths(ka, i, &test_path_result);
	if (! tap_success) {
	for (int jj = 0; jj < jnum; jj++) {
	ref_lennard_jones_single_interaction(ka, i, neighs[jj], MORSEFLAG);
	}
	continue;
	}

	ivec j_2;
	fvec delx_2, dely_2, delz_2, rsq_2;
	bvec jtype_mask_2;
	int num_2 = 0;

	fvec result_f_i_x = fvec::setzero();
	fvec result_f_i_y = fvec::setzero();
	fvec result_f_i_z = fvec::setzero();

	int jj = 0;
	bool rest_j = jj < jnum;
	bool rest_2 = fvec::fast_compress();
	#pragma forceinline recursive
	while (rest_j \|\| rest_2) {
	fvec delx, dely, delz, rsq;
	bvec jtype_mask, within_cutoff;
	ivec j;
	if (rest_j) {
	bvec mask_0 = bvec::full();
	//0xFF >> (8 - (jnum - jj));
	if (jj + (fvec::VL - 1) >= jnum) mask_0 = bvec::only(jnum - jj);
	j = ivec::maskz_loadu(mask_0, &neighs[jj]);
	fvec x_j, y_j, z_j;
	aut_loadatoms_vec(x, j, &x_j, &y_j, &z_j, &jtype_mask, map, map_i,
	c_i1);
	fvec::gather_prefetch0(ivec::mullo(c_i4,
	ivec::maskz_loadu(bvec::full(), &neighs[jj + fvec::VL])), x);
	_mm_prefetch((const char)&neighs[jj + 2 fvec::VL], _MM_HINT_T0);
	delx = x_i - x_j;
	dely = y_i - y_j;
	delz = z_i - z_j;
	rsq = delx * delx + dely * dely + delz * delz;
	fvec cutoff_sq = fvec::mask_blend(jtype_mask, cutljsq0, cutljsq1);
	within_cutoff = fvec::mask_cmplt(mask_0, rsq, cutoff_sq);

	if (fvec::fast_compress()) {
	j = ivec::masku_compress(within_cutoff, j);
	delx = fvec::masku_compress(within_cutoff, delx);
	dely = fvec::masku_compress(within_cutoff, dely);
	delz = fvec::masku_compress(within_cutoff, delz);
	rsq = fvec::masku_compress(within_cutoff, rsq);
	jtype_mask = bvec::masku_compress(within_cutoff, jtype_mask);
	//within_cutoff = 0xFF >> (8 - _cc_popcnt(within_cutoff));

	bvec mask_2 = bvec::after(num_2);//0xFF << num_2;
	j_2 = ivec::mask_expand(j_2, mask_2, j);
	delx_2 = fvec::mask_expand(delx_2, mask_2, delx);
	dely_2 = fvec::mask_expand(dely_2, mask_2, dely);
	delz_2 = fvec::mask_expand(delz_2, mask_2, delz);
	rsq_2 = fvec::mask_expand(rsq_2, mask_2, rsq);
	jtype_mask_2 = bvec::mask_expand(jtype_mask_2, mask_2, jtype_mask);
	num_2 = num_2 + bvec::popcnt(within_cutoff);
	if (num_2 < fvec::VL) {
	jj += fvec::VL;
	rest_j = jj < jnum;
	continue;
	}

	num_2 -= fvec::VL;
	//(0xFF >> (8 - num_2)) << (_cc_popcnt(within_cutoff) - num_2);
	mask_2 = bvec::onlyafter(num_2, bvec::popcnt(within_cutoff) - num_2);
	{
	ivec tmp_j = j_2;
	j_2 = ivec::masku_compress(mask_2, j);
	j = tmp_j;
	fvec tmp_delx = delx_2;
	delx_2 = fvec::masku_compress(mask_2, delx);
	delx = tmp_delx;
	fvec tmp_dely = dely_2;
	dely_2 = fvec::masku_compress(mask_2, dely);
	dely = tmp_dely;
	fvec tmp_delz = delz_2;
	delz_2 = fvec::masku_compress(mask_2, delz);
	delz = tmp_delz;
	fvec tmp_rsq = rsq_2;
	rsq_2 = fvec::masku_compress(mask_2, rsq);
	rsq = tmp_rsq;
	bvec tmp_jtype_mask = jtype_mask_2;
	jtype_mask_2 = bvec::masku_compress(mask_2, jtype_mask);
	jtype_mask = tmp_jtype_mask;
	within_cutoff = bvec::full();
	}
	}
	} else if (rest_2) {
	rest_2 = false;
	j = j_2;
	delx = delx_2;
	dely = dely_2;
	delz = delz_2;
	rsq = rsq_2;
	jtype_mask = jtype_mask_2;
	within_cutoff = bvec::only(num_2);
	num_2 = 0;
	}

	bvec current_mask = within_cutoff;
	if (bvec::test_all_unset(current_mask)) {
	jj += fvec::VL;
	rest_j = jj < jnum;
	continue;
	}

	fvec rij = fvec::sqrt(rsq);
	LennardJonesPathAIREBOT<flt_t> testpath[fvec::VL];
	fvec cij = c_1_0;
	fvec p_cut3rebo = fvec::set1(ka->params.cut3rebo);
	bvec need_search = fvec::mask_cmplt(current_mask, rij, p_cut3rebo);
	if (bvec::test_any_set(need_search)) {
	fvec p_rcmax = fvec::mask_blend(jtype_mask, p_rcmax0, p_rcmax1);
	#pragma noinline
	cij = aut_airebo_lj_tap_test_path(ka, &test_path_result, need_search,
	i_bc, j, testpath);
	}
	current_mask = fvec::mask_cmplt(current_mask, c_0_0, cij);
	if (bvec::test_all_unset(current_mask)) {
	jj += fvec::VL;
	rest_j = jj < jnum;
	continue;
	}
	bvec need_path_force = fvec::mask_cmplt(current_mask, cij, c_1_0);

	fvec p_rljmax = fvec::mask_blend(jtype_mask, p_rljmax0, p_rljmax1);
	fvec p_rljmin = fvec::mask_blend(jtype_mask, p_rljmin0, p_rljmin1);

	fvec dslw, slw = aut_Sp2_deriv(rij, p_rljmin, p_rljmax, &dslw);

	fvec p_lj1 = fvec::mask_blend(jtype_mask, p_lj10, p_lj11);
	fvec p_lj2 = fvec::mask_blend(jtype_mask, p_lj20, p_lj21);
	fvec p_lj3 = fvec::mask_blend(jtype_mask, p_lj30, p_lj31);
	fvec p_lj4 = fvec::mask_blend(jtype_mask, p_lj40, p_lj41);

	fvec vdw, dvdw;

	fvec r2inv = fvec::recip(rsq);

	if (MORSEFLAG) {
	fvec exr = fvec::exp(fvec::setzero() - rij * p_lj4);
	vdw = p_lj1 * exr * (p_lj2 * exr - c_2_0);
	dvdw = p_lj3 * exr * (c_1_0 - p_lj2 * exr);
	} else {
	fvec r6inv = r2inv * r2inv * r2inv;

	vdw = r6inv * ( p_lj3 * r6inv - p_lj4);
	fvec r7inv = r6inv * rij * r2inv;
	dvdw = r7inv * ( p_lj2 - p_lj1 * r6inv);
	}

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

	fvec p_rcLJmin = fvec::mask_blend(jtype_mask, p_rcLJmin0, p_rcLJmin1);
	fvec p_rcLJmax = fvec::mask_blend(jtype_mask, p_rcLJmax0, p_rcLJmax1);
	fvec dStr, Str = aut_Sp2_deriv(rij, p_rcLJmin, p_rcLJmax, &dStr);
	fvec VA = cij * VLJ * Str;
	bvec need_bondorder = fvec::mask_cmplt(current_mask, c_0_0, Str);
	fvec Stb = fvec::setzero();
	fvec fij[3];
	fij[0] = fvec::setzero();
	fij[1] = fvec::setzero();
	fij[2] = fvec::setzero();
	if (bvec::test_any_set(need_bondorder)) {
	for (int jtype = 0; jtype < 2; jtype++) {
	bvec need_bo_with_jtype = need_bondorder;
	if (jtype) need_bo_with_jtype = need_bo_with_jtype & jtype_mask;
	else need_bo_with_jtype = need_bo_with_jtype & ~ jtype_mask;
	ivec jtmp = ivec::masku_compress(need_bo_with_jtype, j);
	ivec itmp = ivec::masku_compress(need_bo_with_jtype, ivec::set1(i));
	fvec cijtmp = fvec::masku_compress(need_bo_with_jtype, cij);
	bvec insert_mask = bvec::after(num_bo[itype][jtype]);
	i_bo[itype][jtype] = ivec::mask_expand(i_bo[itype][jtype],
	insert_mask, itmp);
	j_bo[itype][jtype] = ivec::mask_expand(j_bo[itype][jtype],
	insert_mask, jtmp);
	cij_bo[itype][jtype] = fvec::mask_expand(cij_bo[itype][jtype],
	insert_mask, cijtmp);
	bvec need_path_force_with_jtype = need_bo_with_jtype &
	need_path_force;
	int testpath_end = fvec::VL;
	if (bvec::test_any_set(need_path_force_with_jtype)) {
	int pos = num_bo[itype][jtype];
	for (int l = 0; l < fvec::VL; l++) {
	if (pos >= fvec::VL) {
	testpath_end = l;
	break;
	}
	if (bvec::test_at(need_path_force_with_jtype, l)) {
	testpath_bo[itype][jtype][pos] = testpath[l];
	}
	if (bvec::test_at(need_bo_with_jtype, l)) {
	pos += 1;
	}
	}
	}
	num_bo[itype][jtype] = num_bo[itype][jtype] +
	bvec::popcnt(need_bo_with_jtype);
	if (num_bo[itype][jtype] >= fvec::VL) {
	#pragma noinline
	aut_lj_with_bo<MORSEFLAG>(ka, itype, jtype, i_bo[itype][jtype],
	j_bo[itype][jtype], cij_bo[itype][jtype],
	testpath_bo[itype][jtype]);
	num_bo[itype][jtype] -= fvec::VL;
	insert_mask = bvec::onlyafter(num_bo[itype][jtype],
	bvec::popcnt(need_bo_with_jtype) -
	num_bo[itype][jtype]);
	i_bo[itype][jtype] = ivec::masku_compress(insert_mask, itmp);
	j_bo[itype][jtype] = ivec::masku_compress(insert_mask, jtmp);
	cij_bo[itype][jtype] = fvec::masku_compress(insert_mask, cijtmp);
	if (bvec::test_any_set(need_path_force_with_jtype)) {
	int pos = 0;
	for (int l = testpath_end; l < fvec::VL; l++) {
	if (bvec::test_at(need_path_force_with_jtype, l)) {
	testpath_bo[itype][jtype][pos] = testpath[l];
	}
	if (bvec::test_at(need_bo_with_jtype, l)) {
	pos += 1;
	}
	}
	}
	}
	}
	current_mask = current_mask & ~ need_bondorder;
	need_path_force = need_path_force & ~ need_bondorder;
	}

	fvec fpdVLJ = cij * dVLJ * ( c_1_0 + Str * ( Stb - c_1_0));
	fvec fpdStr = dStr * cij * ( Stb * VLJ - VLJ);
	fvec fpair = r2inv * rij * ( fvec::setzero() - ( fpdVLJ + fpdStr));
	fvec evdwl = VA * Stb + cij * VLJ * ( c_1_0 - Str);

	fvec fix = fpair * delx + fij[0];
	fvec fiy = fpair * dely + fij[1];
	fvec fiz = fpair * delz + fij[2];
	result_f_i_x = fvec::mask_add(result_f_i_x, current_mask, result_f_i_x,
	fix);
	result_f_i_y = fvec::mask_add(result_f_i_y, current_mask, result_f_i_y,
	fiy);
	result_f_i_z = fvec::mask_add(result_f_i_z, current_mask, result_f_i_z,
	fiz);
	result_eng = fvec::mask_add(result_eng, current_mask, result_eng, evdwl);

	ivec j_dbl_idx = ivec::mullo(j, c_i4);
	avec fjx = avec::mask_gather(avec::undefined(), current_mask, j_dbl_idx,
	&ka->result_f[0].x, sizeof(acc_t));
	avec fjy = avec::mask_gather(avec::undefined(), current_mask, j_dbl_idx,
	&ka->result_f[0].y, sizeof(acc_t));
	avec fjz = avec::mask_gather(avec::undefined(), current_mask, j_dbl_idx,
	&ka->result_f[0].z, sizeof(acc_t));

	fjx = fjx - fix;
	fjy = fjy - fiy;
	fjz = fjz - fiz;
	avec::mask_i32loscatter(&ka->result_f[0].x, current_mask, j_dbl_idx, fjx,
	sizeof(acc_t));
	avec::mask_i32loscatter(&ka->result_f[0].y, current_mask, j_dbl_idx, fjy,
	sizeof(acc_t));
	avec::mask_i32loscatter(&ka->result_f[0].z, current_mask, j_dbl_idx, fjz,
	sizeof(acc_t));

	if (bvec::test_any_set(need_path_force)) {
	fvec dC = VLJ * ( Str * Stb + c_1_0 - Str);
	#pragma noinline
	aut_airebo_lj_force_path(ka, need_path_force, dC, testpath);
	}
	jj += fvec::VL;
	rest_j = jj < jnum;
	}
	ka->result_f[i].x += fvec::reduce_add(result_f_i_x);
	ka->result_f[i].y += fvec::reduce_add(result_f_i_y);
	ka->result_f[i].z += fvec::reduce_add(result_f_i_z);
	}
	for (int itype = 0; itype < 2; itype++) {
	for (int jtype = 0; jtype < 2; jtype++) {
	for (int l = 0; l < num_bo[itype][jtype]; l++) {
	ref_lennard_jones_single_interaction(ka,ivec::at(i_bo[itype][jtype],l),
	ivec::at(j_bo[itype][jtype], l),
	MORSEFLAG);
	}
	}
	}
	ka->result_eng += fvec::reduce_add(result_eng);
	}

	};

	template<typename flt_t, typename acc_t>
	void aut_lennard_jones(KernelArgsAIREBOT<flt_t,acc_t> * ka, int morseflag) {
	#ifdef LMP_INTEL_AIREBO_REF
	ref_lennard_jones(ka, morseflag);
	#else
	if (morseflag) {
	aut_wrap<flt_t,acc_t>::template aut_lennard_jones<1>(ka);
	} else {
	aut_wrap<flt_t,acc_t>::template aut_lennard_jones<0>(ka);
	}
	#endif
	}

	template<typename flt_t, typename acc_t>
	void aut_rebo_neigh(KernelArgsAIREBOT<flt_t,acc_t> * ka) {
	#ifdef LMP_INTEL_AIREBO_REF
	ref_rebo_neigh(ka);
	#else
	aut_wrap<flt_t,acc_t>::aut_rebo_neigh(ka);
	#endif
	}

	template<typename flt_t, typename acc_t>
	void aut_frebo(KernelArgsAIREBOT<flt_t,acc_t> * ka, int torsion_flag) {
	#ifdef LMP_INTEL_AIREBO_REF
	ref_frebo(ka, torsion_flag);
	#else
	aut_wrap<flt_t,acc_t>::aut_frebo(ka, torsion_flag);
	#endif
	}

	#ifdef __INTEL_OFFLOAD
	#pragma offload_attribute(pop)
	#endif

	}

pair_airebo_intel.cppNo OneTemporaryActions

File Metadata

pair_airebo_intel.cppView Options

Event Timeline

pair_airebo_intel.cpp
No OneTemporary
Actions

pair_airebo_intel.cpp
View Options