// **************************************************************************
//                               tersoff_zbl.cu
//                             -------------------
//                              Trung Dac Nguyen
//
//  Device code for acceleration of the tersoff/zbl pair style
//
// __________________________________________________________________________
//    This file is part of the LAMMPS Accelerator Library (LAMMPS_AL)
// __________________________________________________________________________
//
//       begin                :
//       email                : ndactrung@gmail.com
// ***************************************************************************/

#ifdef NV_KERNEL
#include "lal_tersoff_zbl_extra.h"

#ifndef _DOUBLE_DOUBLE
texture<float4> pos_tex;
texture<float4> ts1_tex;
texture<float4> ts2_tex;
texture<float4> ts3_tex;
texture<float4> ts4_tex;
texture<float4> ts5_tex;
texture<float4> ts6_tex;
#else
texture<int4,1> pos_tex;
texture<int4> ts1_tex;
texture<int4> ts2_tex;
texture<int4> ts3_tex;
texture<int4> ts4_tex;
texture<int4> ts5_tex;
texture<int4> ts6_tex;
#endif

#else
#define pos_tex x_
#define ts1_tex ts1
#define ts2_tex ts2
#define ts3_tex ts3
#define ts4_tex ts4
#define ts5_tex ts5
#define ts6_tex ts6
#endif

//#define THREE_CONCURRENT

#define TWOTHIRD (numtyp)0.66666666666666666667

#define zeta_idx(nbor_mem, packed_mem, nbor_pitch, n_stride, t_per_atom,    \
                 i, nbor_j, offset_j, idx)                                  \
  if (nbor_mem==packed_mem) {                                               \
    int jj = (nbor_j-offset_j-2*nbor_pitch)/n_stride;                       \
    idx = jj*n_stride + i*t_per_atom + offset_j;                            \
  } else {                                                                  \
    idx = nbor_j;                                                           \
  }

#if (ARCH < 300)

#define store_answers_p(f, energy, virial, ii, inum, tid, t_per_atom,       \
                        offset, eflag, vflag, ans, engv)                    \
  if (t_per_atom>1) {                                                       \
    __local acctyp red_acc[6][BLOCK_PAIR];                                  \
    red_acc[0][tid]=f.x;                                                    \
    red_acc[1][tid]=f.y;                                                    \
    red_acc[2][tid]=f.z;                                                    \
    red_acc[3][tid]=energy;                                                 \
    for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
      if (offset < s) {                                                     \
        for (int r=0; r<4; r++)                                             \
          red_acc[r][tid] += red_acc[r][tid+s];                             \
      }                                                                     \
    }                                                                       \
    f.x=red_acc[0][tid];                                                    \
    f.y=red_acc[1][tid];                                                    \
    f.z=red_acc[2][tid];                                                    \
    energy=red_acc[3][tid];                                                 \
    if (vflag>0) {                                                          \
      for (int r=0; r<6; r++)                                               \
        red_acc[r][tid]=virial[r];                                          \
      for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                       \
        if (offset < s) {                                                   \
          for (int r=0; r<6; r++)                                           \
            red_acc[r][tid] += red_acc[r][tid+s];                           \
        }                                                                   \
      }                                                                     \
      for (int r=0; r<6; r++)                                               \
        virial[r]=red_acc[r][tid];                                          \
    }                                                                       \
  }                                                                         \
  if (offset==0) {                                                          \
    int ei=ii;                                                              \
    if (eflag>0) {                                                          \
      engv[ei]+=energy*(acctyp)0.5;                                         \
      ei+=inum;                                                             \
    }                                                                       \
    if (vflag>0) {                                                          \
      for (int i=0; i<6; i++) {                                             \
        engv[ei]+=virial[i]*(acctyp)0.5;                                    \
        ei+=inum;                                                           \
      }                                                                     \
    }                                                                       \
    acctyp4 old=ans[ii];                                                    \
    old.x+=f.x;                                                             \
    old.y+=f.y;                                                             \
    old.z+=f.z;                                                             \
    ans[ii]=old;                                                            \
  }

#define acc_zeta(z, tid, t_per_atom, offset)                                \
  if (t_per_atom>1) {                                                       \
    __local acctyp red_acc[BLOCK_PAIR];                                     \
    red_acc[tid]=z;                                                         \
    for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
      if (offset < s) {                                                     \
        red_acc[tid] += red_acc[tid+s];                                     \
      }                                                                     \
    }                                                                       \
    z=red_acc[tid];                                                         \
  }

#else

#define store_answers_p(f, energy, virial, ii, inum, tid, t_per_atom,       \
                        offset, eflag, vflag, ans, engv)                    \
  if (t_per_atom>1) {                                                       \
    for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
      f.x += shfl_xor(f.x, s, t_per_atom);                                  \
      f.y += shfl_xor(f.y, s, t_per_atom);                                  \
      f.z += shfl_xor(f.z, s, t_per_atom);                                  \
      energy += shfl_xor(energy, s, t_per_atom);                            \
    }                                                                       \
    if (vflag>0) {                                                          \
      for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                       \
          for (int r=0; r<6; r++)                                           \
            virial[r] += shfl_xor(virial[r], s, t_per_atom);                \
      }                                                                     \
    }                                                                       \
  }                                                                         \
  if (offset==0) {                                                          \
    int ei=ii;                                                              \
    if (eflag>0) {                                                          \
      engv[ei]+=energy*(acctyp)0.5;                                         \
      ei+=inum;                                                             \
    }                                                                       \
    if (vflag>0) {                                                          \
      for (int i=0; i<6; i++) {                                             \
        engv[ei]+=virial[i]*(acctyp)0.5;                                    \
        ei+=inum;                                                           \
      }                                                                     \
    }                                                                       \
    acctyp4 old=ans[ii];                                                    \
    old.x+=f.x;                                                             \
    old.y+=f.y;                                                             \
    old.z+=f.z;                                                             \
    ans[ii]=old;                                                            \
  }

#define acc_zeta(z, tid, t_per_atom, offset)                                \
  if (t_per_atom>1) {                                                       \
    for (unsigned int s=t_per_atom/2; s>0; s>>=1) {                         \
      z += shfl_xor(z, s, t_per_atom);                                      \
    }                                                                       \
  }

#endif

__kernel void k_tersoff_zbl_short_nbor(const __global numtyp4 *restrict x_,
                                   const __global numtyp *restrict cutsq,
                                   const __global int *restrict map,
                                   const __global int *restrict elem2param,
                                   const int nelements, const int nparams,
                                   const __global int * dev_nbor,
                                   const __global int * dev_packed,
                                   __global int * dev_short_nbor,
                                   const int inum, const int nbor_pitch,
                                   const int t_per_atom) {
  __local int n_stride;
  int tid, ii, offset;
  atom_info(t_per_atom,ii,tid,offset);

  if (ii<inum) {
    int nbor, nbor_end;
    int i, numj;
    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
              n_stride,nbor_end,nbor);

    numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
    int itype=ix.w;
    itype=map[itype];

    int ncount = 0;
    int m = nbor;
    dev_short_nbor[m] = 0;
    int nbor_short = nbor+n_stride;

    for ( ; nbor<nbor_end; nbor+=n_stride) {

      int j=dev_packed[nbor];
      int nj = j;
      j &= NEIGHMASK;

      numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
      int jtype=jx.w;
      jtype=map[jtype];
      int ijparam=elem2param[itype*nelements*nelements+jtype*nelements+jtype];

      // Compute r12
      numtyp delx = ix.x-jx.x;
      numtyp dely = ix.y-jx.y;
      numtyp delz = ix.z-jx.z;
      numtyp rsq = delx*delx+dely*dely+delz*delz;

      if (rsq<cutsq[ijparam]) {
        dev_short_nbor[nbor_short] = nj;
        nbor_short += n_stride;
        ncount++;
      }
    } // for nbor

    // store the number of neighbors for each thread
    dev_short_nbor[m] = ncount;

  } // if ii
}

// Tersoff is currently used for 3 elements at most: 3*3*3 = 27 entries
// while the block size should never be less than 32.
// SHARED_SIZE = 32 for now to reduce the pressure on the shared memory per block
// must be increased if there will be more than 3 elements in the future.

#define SHARED_SIZE 32

__kernel void k_tersoff_zbl_zeta(const __global numtyp4 *restrict x_,
                             const __global numtyp4 *restrict ts1_in,
                             const __global numtyp4 *restrict ts2_in,
                             const __global numtyp4 *restrict ts3_in,
                             const __global numtyp4 *restrict ts4_in,
                             const __global numtyp4 *restrict ts5_in,
                             const __global numtyp4 *restrict ts6_in,
                             const __global numtyp *restrict cutsq,
                             const __global int *restrict map,
                             const __global int *restrict elem2param,
                             const int nelements, const int nparams,
                             __global acctyp4 * zetaij,
                             const __global int * dev_nbor,
                             const __global int * dev_packed,
                             const __global int * dev_short_nbor,
                             const int eflag, const int inum,
                             const int nbor_pitch, const int t_per_atom) {
  __local int tpa_sq,n_stride;
  tpa_sq = fast_mul(t_per_atom,t_per_atom);

  int tid, ii, offset;
  atom_info(tpa_sq,ii,tid,offset);

  // must be increased if there will be more than 3 elements in the future.
  __local numtyp4 ts1[SHARED_SIZE];
  __local numtyp4 ts2[SHARED_SIZE];
  __local numtyp4 ts3[SHARED_SIZE];
  __local numtyp4 ts4[SHARED_SIZE];
  __local numtyp4 ts5[SHARED_SIZE];
  __local numtyp4 ts6[SHARED_SIZE];
  if (tid<nparams) {
    ts1[tid]=ts1_in[tid];
    ts2[tid]=ts2_in[tid];
    ts3[tid]=ts3_in[tid];
    ts4[tid]=ts4_in[tid];
    ts5[tid]=ts5_in[tid];
    ts6[tid]=ts6_in[tid];
  }

  acctyp z = (acctyp)0;

  __syncthreads();

  if (ii<inum) {
    int nbor_j, nbor_end, i, numj;
    const __global int* nbor_mem=dev_packed;
    int offset_j=offset/t_per_atom;
    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
              n_stride,nbor_end,nbor_j);
    int offset_k=tid & (t_per_atom-1);

    numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
    int itype=ix.w;
    itype=map[itype];

    // recalculate numj and nbor_end for use of the short nbor list
    if (dev_packed==dev_nbor) {
      numj = dev_short_nbor[nbor_j];
      nbor_j += n_stride;
      nbor_end = nbor_j+fast_mul(numj,n_stride);
      nbor_mem = dev_short_nbor;
    }
    int nborj_start = nbor_j;

    for ( ; nbor_j<nbor_end; nbor_j+=n_stride) {

      int j=nbor_mem[nbor_j];
      j &= NEIGHMASK;

      numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
      int jtype=jx.w;
      jtype=map[jtype];
      int ijparam=elem2param[itype*nelements*nelements+jtype*nelements+jtype];

      // Compute rij
      numtyp4 delr1, delr2;
      delr1.x = jx.x-ix.x;
      delr1.y = jx.y-ix.y;
      delr1.z = jx.z-ix.z;
      numtyp rsq1 = delr1.x*delr1.x+delr1.y*delr1.y+delr1.z*delr1.z;

      // compute zeta_ij
      z = (acctyp)0;

      int nbor_k = nborj_start-offset_j+offset_k;
      int k_end = nbor_end;
      if (dev_packed==dev_nbor) {
        int numk = dev_short_nbor[nbor_k-n_stride];
        k_end = nbor_k+fast_mul(numk,n_stride);
      }

      for ( ; nbor_k < k_end; nbor_k+=n_stride) {
        int k=nbor_mem[nbor_k];
        k &= NEIGHMASK;

        if (k == j) continue;

        numtyp4 kx; fetch4(kx,k,pos_tex); //x_[k];
        int ktype=kx.w;
        ktype=map[ktype];
        int ijkparam=elem2param[itype*nelements*nelements+jtype*nelements+ktype];

        // Compute rik
        delr2.x = kx.x-ix.x;
        delr2.y = kx.y-ix.y;
        delr2.z = kx.z-ix.z;
        numtyp rsq2 = delr2.x*delr2.x+delr2.y*delr2.y+delr2.z*delr2.z;

        if (rsq2 > cutsq[ijkparam]) continue;

        numtyp4 ts1_ijkparam = ts1[ijkparam]; //fetch4(ts1_ijkparam,ijkparam,ts1_tex);
        numtyp ijkparam_lam3 = ts1_ijkparam.z;
        numtyp ijkparam_powermint = ts1_ijkparam.w;
        numtyp4 ts2_ijkparam = ts2[ijkparam]; //fetch4(ts2_ijkparam,ijkparam,ts2_tex);
        numtyp ijkparam_bigr = ts2_ijkparam.z;
        numtyp ijkparam_bigd = ts2_ijkparam.w;
        numtyp4 ts4_ijkparam = ts4[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex);
        numtyp ijkparam_c = ts4_ijkparam.x;
        numtyp ijkparam_d = ts4_ijkparam.y;
        numtyp ijkparam_h = ts4_ijkparam.z;
        numtyp ijkparam_gamma = ts4_ijkparam.w;
        z += zeta(ijkparam_powermint, ijkparam_lam3, ijkparam_bigr, ijkparam_bigd,
                  ijkparam_c, ijkparam_d, ijkparam_h, ijkparam_gamma,
                  rsq1, rsq2, delr1, delr2);
      }

      // idx to zetaij is shifted by n_stride relative to nbor_j in dev_short_nbor
      int idx = nbor_j;
      if (dev_packed==dev_nbor) idx -= n_stride;
      acc_zeta(z, tid, t_per_atom, offset_k);

      numtyp4 ts1_ijparam = ts1[ijparam]; //fetch4(ts1_ijparam,ijparam,ts1_tex);
      numtyp ijparam_lam2 = ts1_ijparam.y;
      numtyp4 ts2_ijparam = ts2[ijparam]; //fetch4(ts2_ijparam,ijparam,ts2_tex);
      numtyp ijparam_bigb = ts2_ijparam.y;
      numtyp ijparam_bigr = ts2_ijparam.z;
      numtyp ijparam_bigd = ts2_ijparam.w;
      numtyp4 ts3_ijparam = ts3[ijparam]; //fetch4(ts3_ijparam,ijparam,ts3_tex);
      numtyp ijparam_c1 = ts3_ijparam.x;
      numtyp ijparam_c2 = ts3_ijparam.y;
      numtyp ijparam_c3 = ts3_ijparam.z;
      numtyp ijparam_c4 = ts3_ijparam.w;
      numtyp4 ts5_ijparam = ts5[ijparam]; //fetch4(ts5_ijparam,ijparam,ts5_tex);
      numtyp ijparam_beta = ts5_ijparam.x;
      numtyp ijparam_powern = ts5_ijparam.y;
      numtyp4 ts6_ijparam = ts6[ijparam]; //fetch4(ts6_ijparam,ijparam,ts6_tex);
      numtyp ijparam_ZBLcut = ts6_ijparam.z;
      numtyp ijparam_ZBLexpscale = ts6_ijparam.w;

      if (offset_k == 0) {
        numtyp fpfeng[4];
        force_zeta(ijparam_bigb, ijparam_bigr, ijparam_bigd, ijparam_lam2,
                   ijparam_beta, ijparam_powern, ijparam_c1, ijparam_c2, ijparam_c3,
                   ijparam_c4, ijparam_ZBLcut, ijparam_ZBLexpscale, rsq1, z, eflag, fpfeng);
        acctyp4 zij;
        zij.x = fpfeng[0];
        zij.y = fpfeng[1];
        zij.z = fpfeng[2];
        zij.w = z;
        zetaij[idx] = zij;
      }

    } // for nbor
  } // if ii
}

__kernel void k_tersoff_zbl_repulsive(const __global numtyp4 *restrict x_,
                                  const __global numtyp4 *restrict ts1_in,
                                  const __global numtyp4 *restrict ts2_in,
                                  const __global numtyp4 *restrict ts6_in,
                                  const numtyp global_e, const numtyp global_a_0,
                                  const numtyp global_epsilon_0,
                                  const __global numtyp *restrict cutsq,
                                  const __global int *restrict map,
                                  const __global int *restrict elem2param,
                                  const int nelements, const int nparams,
                                  const __global int * dev_nbor,
                                  const __global int * dev_packed,
                                  const __global int * dev_short_nbor,
                                  __global acctyp4 *restrict ans,
                                  __global acctyp *restrict engv,
                                  const int eflag, const int vflag,
                                  const int inum, const int nbor_pitch,
                                  const int t_per_atom) {
  __local int n_stride;
  int tid, ii, offset;
  atom_info(t_per_atom,ii,tid,offset);

  __local numtyp4 ts1[SHARED_SIZE];
  __local numtyp4 ts2[SHARED_SIZE];
  __local numtyp4 ts6[SHARED_SIZE];
  if (tid<nparams) {
    ts1[tid]=ts1_in[tid];
    ts2[tid]=ts2_in[tid];
    ts6[tid]=ts6_in[tid];
  }

  acctyp energy=(acctyp)0;
  acctyp4 f;
  f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
  acctyp virial[6];
  for (int i=0; i<6; i++)
    virial[i]=(acctyp)0;

  __syncthreads();

  if (ii<inum) {
    int nbor, nbor_end, i, numj;
    const __global int* nbor_mem=dev_packed;
    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset,i,numj,
              n_stride,nbor_end,nbor);

    numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
    int itype=ix.w;
    itype=map[itype];

    // recalculate numj and nbor_end for use of the short nbor list
    if (dev_packed==dev_nbor) {
      numj = dev_short_nbor[nbor];
      nbor += n_stride;
      nbor_end = nbor+fast_mul(numj,n_stride);
      nbor_mem = dev_short_nbor;
    }

    for ( ; nbor<nbor_end; nbor+=n_stride) {

      int j=nbor_mem[nbor];
      j &= NEIGHMASK;

      numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
      int jtype=jx.w;
      jtype=map[jtype];
      int ijparam=elem2param[itype*nelements*nelements+jtype*nelements+jtype];

      // Compute r12

      numtyp delx = ix.x-jx.x;
      numtyp dely = ix.y-jx.y;
      numtyp delz = ix.z-jx.z;
      numtyp rsq = delx*delx+dely*dely+delz*delz;

      // rsq<cutsq[ijparam]
      numtyp feng[2];
      numtyp ijparam_lam1 = ts1[ijparam].x;
      numtyp4 ts2_ijparam = ts2[ijparam];
      numtyp ijparam_biga = ts2_ijparam.x;
      numtyp ijparam_bigr = ts2_ijparam.z;
      numtyp ijparam_bigd = ts2_ijparam.w;
      numtyp4 ts6_ijparam = ts6[ijparam];
      numtyp ijparam_Z_i = ts6_ijparam.x;
      numtyp ijparam_Z_j = ts6_ijparam.y;
      numtyp ijparam_ZBLcut = ts6_ijparam.z;
      numtyp ijparam_ZBLexpscale = ts6_ijparam.w;

      repulsive(ijparam_bigr, ijparam_bigd, ijparam_lam1, ijparam_biga,
                ijparam_Z_i, ijparam_Z_j, ijparam_ZBLcut, ijparam_ZBLexpscale,
                global_e, global_a_0, global_epsilon_0, rsq, eflag, feng);

      numtyp force = feng[0];
      f.x+=delx*force;
      f.y+=dely*force;
      f.z+=delz*force;

      if (eflag>0)
        energy+=feng[1];
      if (vflag>0) {
        virial[0] += delx*delx*force;
        virial[1] += dely*dely*force;
        virial[2] += delz*delz*force;
        virial[3] += delx*dely*force;
        virial[4] += delx*delz*force;
        virial[5] += dely*delz*force;
      }
    } // for nbor

    store_answers(f,energy,virial,ii,inum,tid,t_per_atom,offset,eflag,vflag,
                  ans,engv);
  } // if ii

}

__kernel void k_tersoff_zbl_three_center(const __global numtyp4 *restrict x_,
                                     const __global numtyp4 *restrict ts1_in,
                                     const __global numtyp4 *restrict ts2_in,
                                     const __global numtyp4 *restrict ts4_in,
                                     const __global numtyp *restrict cutsq,
                                     const __global int *restrict map,
                                     const __global int *restrict elem2param,
                                     const int nelements, const int nparams,
                                     const __global acctyp4 *restrict zetaij,
                                     const __global int * dev_nbor,
                                     const __global int * dev_packed,
                                     const __global int * dev_short_nbor,
                                     __global acctyp4 *restrict ans,
                                     __global acctyp *restrict engv,
                                     const int eflag, const int vflag,
                                     const int inum,  const int nbor_pitch,
                                     const int t_per_atom, const int evatom) {
  __local int tpa_sq, n_stride;
  tpa_sq=fast_mul(t_per_atom,t_per_atom);
  numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;

  int tid, ii, offset;
  atom_info(tpa_sq,ii,tid,offset); // offset ranges from 0 to tpa_sq-1

  __local numtyp4 ts1[SHARED_SIZE];
  __local numtyp4 ts2[SHARED_SIZE];
  __local numtyp4 ts4[SHARED_SIZE];
  if (tid<nparams) {
    ts1[tid]=ts1_in[tid];
    ts2[tid]=ts2_in[tid];
    ts4[tid]=ts4_in[tid];
  }

  acctyp energy=(acctyp)0;
  acctyp4 f;
  f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
  acctyp virial[6];
  for (int i=0; i<6; i++)
    virial[i]=(acctyp)0;
  numtyp tpainv = ucl_recip((numtyp)t_per_atom);

  __syncthreads();

  if (ii<inum) {
    int i, numj, nbor_j, nbor_end;
    const __global int* nbor_mem=dev_packed;
    int offset_j=offset/t_per_atom;
    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
              n_stride,nbor_end,nbor_j);
    int offset_k=tid & (t_per_atom-1);

    numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
    int itype=ix.w;
    itype=map[itype];

    // recalculate numj and nbor_end for use of the short nbor list
    if (dev_packed==dev_nbor) {
      numj = dev_short_nbor[nbor_j];
      nbor_j += n_stride;
      nbor_end = nbor_j+fast_mul(numj,n_stride);
      nbor_mem = dev_short_nbor;
    }
    int nborj_start = nbor_j;

    for ( ; nbor_j<nbor_end; nbor_j+=n_stride) {

      int j=nbor_mem[nbor_j];
      j &= NEIGHMASK;

      numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
      int jtype=jx.w;
      jtype=map[jtype];
      int ijparam=elem2param[itype*nelements*nelements+jtype*nelements+jtype];

      // Compute r12
      numtyp delr1[3];
      delr1[0] = jx.x-ix.x;
      delr1[1] = jx.y-ix.y;
      delr1[2] = jx.z-ix.z;
      numtyp rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];

      numtyp r1 = ucl_sqrt(rsq1);
      numtyp r1inv = ucl_rsqrt(rsq1);

      // look up for zeta_ij
      // idx to zetaij is shifted by n_stride relative to nbor_j in dev_short_nbor
      int idx = nbor_j;
      if (dev_packed==dev_nbor) idx -= n_stride;
      acctyp4 zeta_ij = zetaij[idx]; // fetch(zeta_ij,idx,zeta_tex);
      numtyp force = zeta_ij.x*tpainv;
      numtyp prefactor = zeta_ij.y;
      f.x += delr1[0]*force;
      f.y += delr1[1]*force;
      f.z += delr1[2]*force;

      if (eflag>0) {
        energy+=zeta_ij.z*tpainv;
      }
      if (vflag>0) {
        numtyp mforce = -force;
        virial[0] += delr1[0]*delr1[0]*mforce;
        virial[1] += delr1[1]*delr1[1]*mforce;
        virial[2] += delr1[2]*delr1[2]*mforce;
        virial[3] += delr1[0]*delr1[1]*mforce;
        virial[4] += delr1[0]*delr1[2]*mforce;
        virial[5] += delr1[1]*delr1[2]*mforce;
      }

      int nbor_k = nborj_start-offset_j+offset_k;
      int k_end = nbor_end;
      if (dev_packed==dev_nbor) {
        int numk = dev_short_nbor[nbor_k-n_stride];
        k_end = nbor_k+fast_mul(numk,n_stride);
      }

      for ( ; nbor_k<k_end; nbor_k+=n_stride) {
        int k=nbor_mem[nbor_k];
        k &= NEIGHMASK;

        if (j == k) continue;

        numtyp4 kx; fetch4(kx,k,pos_tex);
        int ktype=kx.w;
        ktype=map[ktype];
        int ijkparam=elem2param[itype*nelements*nelements+jtype*nelements+ktype];

        numtyp delr2[3];
        delr2[0] = kx.x-ix.x;
        delr2[1] = kx.y-ix.y;
        delr2[2] = kx.z-ix.z;
        numtyp rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];

        if (rsq2 > cutsq[ijkparam]) continue;
        numtyp r2 = ucl_sqrt(rsq2);
        numtyp r2inv = ucl_rsqrt(rsq2);

        numtyp fi[3], fj[3], fk[3];
        numtyp4 ts1_ijkparam = ts1[ijkparam]; //fetch4(ts1_ijkparam,ijkparam,ts1_tex);
        lam3 = ts1_ijkparam.z;
        powermint = ts1_ijkparam.w;
        numtyp4 ts2_ijkparam = ts2[ijkparam]; //fetch4(ts2_ijkparam,ijkparam,ts2_tex);
        bigr = ts2_ijkparam.z;
        bigd = ts2_ijkparam.w;
        numtyp4 ts4_ijkparam = ts4[ijkparam]; //fetch4(ts4_ijkparam,ijkparam,ts4_tex);
        c = ts4_ijkparam.x;
        d = ts4_ijkparam.y;
        h = ts4_ijkparam.z;
        gamma = ts4_ijkparam.w;
        if (vflag>0)
          attractive(bigr, bigd, powermint, lam3, c, d, h, gamma,
                     prefactor, r1, r1inv, r2, r2inv, delr1, delr2, fi, fj, fk);
        else
          attractive_fi(bigr, bigd, powermint, lam3, c, d, h, gamma,
                        prefactor, r1, r1inv, r2, r2inv, delr1, delr2, fi);
        f.x += fi[0];
        f.y += fi[1];
        f.z += fi[2];

        if (vflag>0) {
          acctyp v[6];
          numtyp pre = (numtyp)2.0;
          if (evatom==1) pre = TWOTHIRD;
          v[0] = pre*(delr1[0]*fj[0] + delr2[0]*fk[0]);
          v[1] = pre*(delr1[1]*fj[1] + delr2[1]*fk[1]);
          v[2] = pre*(delr1[2]*fj[2] + delr2[2]*fk[2]);
          v[3] = pre*(delr1[0]*fj[1] + delr2[0]*fk[1]);
          v[4] = pre*(delr1[0]*fj[2] + delr2[0]*fk[2]);
          v[5] = pre*(delr1[1]*fj[2] + delr2[1]*fk[2]);

          virial[0] += v[0]; virial[1] += v[1]; virial[2] += v[2];
          virial[3] += v[3]; virial[4] += v[4]; virial[5] += v[5];
        }
      } // nbor_k
    } // for nbor_j

    store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,
                     offset,eflag,vflag,ans,engv);
  } // if ii
}

__kernel void k_tersoff_zbl_three_end(const __global numtyp4 *restrict x_,
                                  const __global numtyp4 *restrict ts1_in,
                                  const __global numtyp4 *restrict ts2_in,
                                  const __global numtyp4 *restrict ts4_in,
                                  const __global numtyp *restrict cutsq,
                                  const __global int *restrict map,
                                  const __global int *restrict elem2param,
                                  const int nelements, const int nparams,
                                  const __global acctyp4 *restrict zetaij,
                                  const __global int * dev_nbor,
                                  const __global int * dev_packed,
                                  const __global int * dev_ilist,
                                  const __global int * dev_short_nbor,
                                  __global acctyp4 *restrict ans,
                                  __global acctyp *restrict engv,
                                  const int eflag, const int vflag,
                                  const int inum,  const int nbor_pitch,
                                  const int t_per_atom, const int gpu_nbor) {
  __local int tpa_sq, n_stride;
  tpa_sq=fast_mul(t_per_atom,t_per_atom);
  numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;

  int tid, ii, offset;
  atom_info(tpa_sq,ii,tid,offset);

  __local numtyp4 ts1[SHARED_SIZE];
  __local numtyp4 ts2[SHARED_SIZE];
  __local numtyp4 ts4[SHARED_SIZE];
  if (tid<nparams) {
    ts1[tid]=ts1_in[tid];
    ts2[tid]=ts2_in[tid];
    ts4[tid]=ts4_in[tid];
  }

  acctyp energy=(acctyp)0;
  acctyp4 f;
  f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
  acctyp virial[6];
  for (int i=0; i<6; i++)
    virial[i]=(acctyp)0;

  __local int red_acc[2*BLOCK_PAIR];

  __syncthreads();

  if (ii<inum) {
    int i, numj, nbor_j, nbor_end, k_end;
    const __global int* nbor_mem=dev_packed;
    int offset_j=offset/t_per_atom;
    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
              n_stride,nbor_end,nbor_j);
    int offset_k=tid & (t_per_atom-1);

    numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
    int itype=ix.w;
    itype=map[itype];

    numtyp tpainv = ucl_recip((numtyp)t_per_atom);

    // recalculate numj and nbor_end for use of the short nbor list
    if (dev_packed==dev_nbor) {
      numj = dev_short_nbor[nbor_j];
      nbor_j += n_stride;
      nbor_end = nbor_j+fast_mul(numj,n_stride);
      nbor_mem = dev_short_nbor;
    }

    for ( ; nbor_j<nbor_end; nbor_j+=n_stride) {

      int j=nbor_mem[nbor_j];
      j &= NEIGHMASK;

      numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
      int jtype=jx.w;
      jtype=map[jtype];
      int ijparam=elem2param[itype*nelements*nelements+jtype*nelements+jtype];

      // Compute r12
      numtyp delr1[3];
      delr1[0] = jx.x-ix.x;
      delr1[1] = jx.y-ix.y;
      delr1[2] = jx.z-ix.z;
      numtyp rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];

      numtyp mdelr1[3];
      mdelr1[0] = -delr1[0];
      mdelr1[1] = -delr1[1];
      mdelr1[2] = -delr1[2];

      int nbor_k,numk;
      if (dev_nbor==dev_packed) {
        if (gpu_nbor) nbor_k=j+nbor_pitch;
        else nbor_k=dev_ilist[j]+nbor_pitch;
        numk=dev_nbor[nbor_k];
        nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
        k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
        nbor_k+=offset_k;
      } else {
        nbor_k=dev_ilist[j]+nbor_pitch;
        numk=dev_nbor[nbor_k];
        nbor_k+=nbor_pitch;
        nbor_k=dev_nbor[nbor_k];
        k_end=nbor_k+numk;
        nbor_k+=offset_k;
      }

      // recalculate numk and k_end for the use of short neighbor list
      if (dev_packed==dev_nbor) {
        numk = dev_short_nbor[nbor_k];
        nbor_k += n_stride;
        k_end = nbor_k+fast_mul(numk,n_stride);
      }
      int nbork_start = nbor_k;

      // look up for zeta_ji: find i in the j's neighbor list
      int m = tid / t_per_atom;
      int ijnum = -1;
      for ( ; nbor_k<k_end; nbor_k+=n_stride) {
        int k=nbor_mem[nbor_k];
        k &= NEIGHMASK;
        if (k == i) {
          ijnum = nbor_k;
          red_acc[2*m+0] = ijnum;
          red_acc[2*m+1] = offset_k;
          break;
        }
      }

      numtyp r1 = ucl_sqrt(rsq1);
      numtyp r1inv = ucl_rsqrt(rsq1);
      int offset_kf;
      if (ijnum >= 0) {
        offset_kf = offset_k;
      } else {
        ijnum = red_acc[2*m+0];
        offset_kf = red_acc[2*m+1];
      }

      // idx to zetaij is shifted by n_stride relative to ijnum in dev_short_nbor
      int idx = ijnum;
      if (dev_packed==dev_nbor) idx -= n_stride;
      acctyp4 zeta_ji = zetaij[idx]; // fetch(zeta_ji,idx,zeta_tex);
      numtyp force = zeta_ji.x*tpainv;
      numtyp prefactor_ji = zeta_ji.y;
      f.x += delr1[0]*force;
      f.y += delr1[1]*force;
      f.z += delr1[2]*force;

      if (eflag>0) {
        energy+=zeta_ji.z*tpainv;
      }
      if (vflag>0) {
        numtyp mforce = -force;
        virial[0] += mdelr1[0]*mdelr1[0]*mforce;
        virial[1] += mdelr1[1]*mdelr1[1]*mforce;
        virial[2] += mdelr1[2]*mdelr1[2]*mforce;
        virial[3] += mdelr1[0]*mdelr1[1]*mforce;
        virial[4] += mdelr1[0]*mdelr1[2]*mforce;
        virial[5] += mdelr1[1]*mdelr1[2]*mforce;
      }

      // attractive forces
      for (nbor_k = nbork_start ; nbor_k<k_end; nbor_k+=n_stride) {
        int k=nbor_mem[nbor_k];
        k &= NEIGHMASK;

        if (k == i) continue;

        numtyp4 kx; fetch4(kx,k,pos_tex);
        int ktype=kx.w;
        ktype=map[ktype];
        int jikparam=elem2param[jtype*nelements*nelements+itype*nelements+ktype];

        numtyp delr2[3];
        delr2[0] = kx.x-jx.x;
        delr2[1] = kx.y-jx.y;
        delr2[2] = kx.z-jx.z;
        numtyp rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];

        if (rsq2 > cutsq[jikparam]) continue;
        numtyp r2 = ucl_sqrt(rsq2);
        numtyp r2inv = ucl_rsqrt(rsq2);
        numtyp4 ts1_param, ts2_param, ts4_param;
        numtyp fi[3];

        ts1_param = ts1[jikparam]; //fetch4(ts1_jikparam,jikparam,ts1_tex);
        lam3 = ts1_param.z;
        powermint = ts1_param.w;
        ts2_param = ts2[jikparam]; //fetch4(ts2_jikparam,jikparam,ts2_tex);
        bigr = ts2_param.z;
        bigd = ts2_param.w;
        ts4_param = ts4[jikparam]; //fetch4(ts4_jikparam,jikparam,ts4_tex);
        c = ts4_param.x;
        d = ts4_param.y;
        h = ts4_param.z;
        gamma = ts4_param.w;
        attractive_fj(bigr, bigd, powermint, lam3, c, d, h, gamma,
                      prefactor_ji, r1, r1inv, r2, r2inv, mdelr1, delr2, fi);
        f.x += fi[0];
        f.y += fi[1];
        f.z += fi[2];

        // idx to zetaij is shifted by n_stride relative to nbor_k in dev_short_nbor
        int idx = nbor_k;
        if (dev_packed==dev_nbor) idx -= n_stride;
        acctyp4 zeta_jk = zetaij[idx]; // fetch(zeta_jk,idx,zeta_tex);
        numtyp prefactor_jk = zeta_jk.y;
        int jkiparam=elem2param[jtype*nelements*nelements+ktype*nelements+itype];
        ts1_param = ts1[jkiparam]; //fetch4(ts1_jkiparam,jkiparam,ts1_tex);
        lam3 = ts1_param.z;
        powermint = ts1_param.w;
        ts2_param = ts2[jkiparam]; //fetch4(ts2_jkiparam,jkiparam,ts2_tex);
        bigr = ts2_param.z;
        bigd = ts2_param.w;
        ts4_param = ts4[jkiparam]; //fetch4(ts4_jkiparam,jkiparam,ts4_tex);
        c = ts4_param.x;
        d = ts4_param.y;
        h = ts4_param.z;
        gamma = ts4_param.w;
        attractive_fk(bigr, bigd, powermint, lam3, c, d, h, gamma,
                      prefactor_jk, r2, r2inv, r1, r1inv, delr2, mdelr1, fi);
        f.x += fi[0];
        f.y += fi[1];
        f.z += fi[2];
      } // for nbor_k
    } // for nbor_j

    #ifdef THREE_CONCURRENT
    store_answers(f,energy,virial,ii,inum,tid,tpa_sq,offset,
                  eflag,vflag,ans,engv);
    #else
    store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,offset,
                    eflag,vflag,ans,engv);
    #endif
  } // if ii
}

__kernel void k_tersoff_zbl_three_end_vatom(const __global numtyp4 *restrict x_,
                                        const __global numtyp4 *restrict ts1_in,
                                        const __global numtyp4 *restrict ts2_in,
                                        const __global numtyp4 *restrict ts4_in,
                                        const __global numtyp *restrict cutsq,
                                        const __global int *restrict map,
                                        const __global int *restrict elem2param,
                                        const int nelements, const int nparams,
                                        const __global acctyp4 *restrict zetaij,
                                        const __global int * dev_nbor,
                                        const __global int * dev_packed,
                                        const __global int * dev_ilist,
                                        const __global int * dev_short_nbor,
                                        __global acctyp4 *restrict ans,
                                        __global acctyp *restrict engv,
                                        const int eflag, const int vflag,
                                        const int inum,  const int nbor_pitch,
                                        const int t_per_atom, const int gpu_nbor) {
  __local int tpa_sq, n_stride;
  tpa_sq=fast_mul(t_per_atom,t_per_atom);
  numtyp lam3, powermint, bigr, bigd, c, d, h, gamma;

  int tid, ii, offset;
  atom_info(tpa_sq,ii,tid,offset);

  __local numtyp4 ts1[SHARED_SIZE];
  __local numtyp4 ts2[SHARED_SIZE];
  __local numtyp4 ts4[SHARED_SIZE];
  if (tid<nparams) {
    ts1[tid]=ts1_in[tid];
    ts2[tid]=ts2_in[tid];
    ts4[tid]=ts4_in[tid];
  }

  acctyp energy=(acctyp)0;
  acctyp4 f;
  f.x=(acctyp)0; f.y=(acctyp)0; f.z=(acctyp)0;
  acctyp virial[6];
  for (int i=0; i<6; i++)
    virial[i]=(acctyp)0;

  __local int red_acc[2*BLOCK_PAIR];

  __syncthreads();

  if (ii<inum) {
    int i, numj, nbor_j, nbor_end, k_end;
    const __global int* nbor_mem = dev_packed;
    int offset_j=offset/t_per_atom;
    nbor_info(dev_nbor,dev_packed,nbor_pitch,t_per_atom,ii,offset_j,i,numj,
              n_stride,nbor_end,nbor_j);
    int offset_k=tid & (t_per_atom-1);

    numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
    int itype=ix.w;
    itype=map[itype];

    numtyp tpainv = ucl_recip((numtyp)t_per_atom);

    // recalculate numj and nbor_end for use of the short nbor list
    if (dev_packed==dev_nbor) {
      numj = dev_short_nbor[nbor_j];
      nbor_j += n_stride;
      nbor_end = nbor_j+fast_mul(numj,n_stride);
      nbor_mem = dev_short_nbor;
    }

    for ( ; nbor_j<nbor_end; nbor_j+=n_stride) {

      int j=nbor_mem[nbor_j];
      j &= NEIGHMASK;

      numtyp4 jx; fetch4(jx,j,pos_tex); //x_[j];
      int jtype=jx.w;
      jtype=map[jtype];
      int ijparam=elem2param[itype*nelements*nelements+jtype*nelements+jtype];

      // Compute r12
      numtyp delr1[3];
      delr1[0] = jx.x-ix.x;
      delr1[1] = jx.y-ix.y;
      delr1[2] = jx.z-ix.z;
      numtyp rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];

      numtyp mdelr1[3];
      mdelr1[0] = -delr1[0];
      mdelr1[1] = -delr1[1];
      mdelr1[2] = -delr1[2];

      int nbor_k,numk;
      if (dev_nbor==dev_packed) {
        if (gpu_nbor) nbor_k=j+nbor_pitch;
        else nbor_k=dev_ilist[j]+nbor_pitch;
        numk=dev_nbor[nbor_k];
        nbor_k+=nbor_pitch+fast_mul(j,t_per_atom-1);
        k_end=nbor_k+fast_mul(numk/t_per_atom,n_stride)+(numk & (t_per_atom-1));
        nbor_k+=offset_k;
      } else {
        nbor_k=dev_ilist[j]+nbor_pitch;
        numk=dev_nbor[nbor_k];
        nbor_k+=nbor_pitch;
        nbor_k=dev_nbor[nbor_k];
        k_end=nbor_k+numk;
        nbor_k+=offset_k;
      }

      // recalculate numk and k_end for the use of short neighbor list
      if (dev_packed==dev_nbor) {
        numk = dev_short_nbor[nbor_k];
        nbor_k += n_stride;
        k_end = nbor_k+fast_mul(numk,n_stride);
      }
      int nbork_start = nbor_k;

      // look up for zeta_ji
      int m = tid / t_per_atom;
      int ijnum = -1;
      for ( ; nbor_k<k_end; nbor_k+=n_stride) {
        int k=nbor_mem[nbor_k];
        k &= NEIGHMASK;
        if (k == i) {
          ijnum = nbor_k;
          red_acc[2*m+0] = ijnum;
          red_acc[2*m+1] = offset_k;
          break;
        }
      }

      numtyp r1 = ucl_sqrt(rsq1);
      numtyp r1inv = ucl_rsqrt(rsq1);
      int offset_kf;
      if (ijnum >= 0) {
        offset_kf = offset_k;
      } else {
        ijnum = red_acc[2*m+0];
        offset_kf = red_acc[2*m+1];
      }

      // idx to zetaij is shifted by n_stride relative to ijnum in dev_short_nbor
      int idx = ijnum;
      if (dev_packed==dev_nbor) idx -= n_stride;
      acctyp4 zeta_ji = zetaij[idx]; //  fetch(zeta_ji,idx,zeta_tex);
      numtyp force = zeta_ji.x*tpainv;
      numtyp prefactor_ji = zeta_ji.y;
      f.x += delr1[0]*force;
      f.y += delr1[1]*force;
      f.z += delr1[2]*force;

      if (eflag>0) {
        energy+=zeta_ji.z*tpainv;
      }
      if (vflag>0) {
        numtyp mforce = -force;
        virial[0] += mdelr1[0]*mdelr1[0]*mforce;
        virial[1] += mdelr1[1]*mdelr1[1]*mforce;
        virial[2] += mdelr1[2]*mdelr1[2]*mforce;
        virial[3] += mdelr1[0]*mdelr1[1]*mforce;
        virial[4] += mdelr1[0]*mdelr1[2]*mforce;
        virial[5] += mdelr1[1]*mdelr1[2]*mforce;
      }

      // attractive forces
      for (nbor_k = nbork_start; nbor_k<k_end; nbor_k+=n_stride) {
        int k=nbor_mem[nbor_k];
        k &= NEIGHMASK;

        if (k == i) continue;

        numtyp4 kx; fetch4(kx,k,pos_tex);
        int ktype=kx.w;
        ktype=map[ktype];
        int jikparam=elem2param[jtype*nelements*nelements+itype*nelements+ktype];

        numtyp delr2[3];
        delr2[0] = kx.x-jx.x;
        delr2[1] = kx.y-jx.y;
        delr2[2] = kx.z-jx.z;
        numtyp rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];

        if (rsq2 > cutsq[jikparam]) continue;
        numtyp r2 = ucl_sqrt(rsq2);
        numtyp r2inv = ucl_rsqrt(rsq2);

        numtyp fi[3], fj[3], fk[3];
        numtyp4 ts1_param, ts2_param, ts4_param;
        ts1_param = ts1[jikparam]; //fetch4(ts1_jikparam,jikparam,ts1_tex);
        lam3 = ts1_param.z;
        powermint = ts1_param.w;
        ts2_param = ts2[jikparam]; //fetch4(ts2_jikparam,jikparam,ts2_tex);
        bigr = ts2_param.z;
        bigd = ts2_param.w;
        ts4_param = ts4[jikparam]; //fetch4(ts4_jikparam,jikparam,ts4_tex);
        c = ts4_param.x;
        d = ts4_param.y;
        h = ts4_param.z;
        gamma = ts4_param.w;
        attractive(bigr, bigd, powermint, lam3, c, d, h, gamma,
                   prefactor_ji, r1, r1inv, r2, r2inv, mdelr1, delr2, fi, fj, fk);
        f.x += fj[0];
        f.y += fj[1];
        f.z += fj[2];

        virial[0] += TWOTHIRD*(mdelr1[0]*fj[0] + delr2[0]*fk[0]);
        virial[1] += TWOTHIRD*(mdelr1[1]*fj[1] + delr2[1]*fk[1]);
        virial[2] += TWOTHIRD*(mdelr1[2]*fj[2] + delr2[2]*fk[2]);
        virial[3] += TWOTHIRD*(mdelr1[0]*fj[1] + delr2[0]*fk[1]);
        virial[4] += TWOTHIRD*(mdelr1[0]*fj[2] + delr2[0]*fk[2]);
        virial[5] += TWOTHIRD*(mdelr1[1]*fj[2] + delr2[1]*fk[2]);

        // idx to zetaij is shifted by n_stride relative to nbor_k in dev_short_nbor
        int idx = nbor_k;
        if (dev_packed==dev_nbor) idx -= n_stride;
        acctyp4 zeta_jk = zetaij[idx]; // fetch(zeta_jk,idx,zeta_tex);
        numtyp prefactor_jk = zeta_jk.y;

        int jkiparam=elem2param[jtype*nelements*nelements+ktype*nelements+itype];
        ts1_param = ts1[jkiparam]; //fetch4(ts1_jkiparam,jkiparam,ts1_tex);
        lam3 = ts1_param.z;
        powermint = ts1_param.w;
        ts2_param = ts2[jkiparam]; //fetch4(ts2_jkiparam,jkiparam,ts2_tex);
        bigr = ts2_param.z;
        bigd = ts2_param.w;
        ts4_param = ts4[jkiparam]; //fetch4(ts4_jkiparam,jkiparam,ts4_tex);
        c = ts4_param.x;
        d = ts4_param.y;
        h = ts4_param.z;
        gamma = ts4_param.w;
        attractive(bigr, bigd, powermint, lam3, c, d, h, gamma,
                   prefactor_jk, r2, r2inv, r1, r1inv, delr2, mdelr1, fi, fj, fk);
        f.x += fk[0];
        f.y += fk[1];
        f.z += fk[2];

        virial[0] += TWOTHIRD*(delr2[0]*fj[0] + mdelr1[0]*fk[0]);
        virial[1] += TWOTHIRD*(delr2[1]*fj[1] + mdelr1[1]*fk[1]);
        virial[2] += TWOTHIRD*(delr2[2]*fj[2] + mdelr1[2]*fk[2]);
        virial[3] += TWOTHIRD*(delr2[0]*fj[1] + mdelr1[0]*fk[1]);
        virial[4] += TWOTHIRD*(delr2[0]*fj[2] + mdelr1[0]*fk[2]);
        virial[5] += TWOTHIRD*(delr2[1]*fj[2] + mdelr1[1]*fk[2]);
      }
    } // for nbor

    #ifdef THREE_CONCURRENT
    store_answers(f,energy,virial,ii,inum,tid,tpa_sq,offset,
                  eflag,vflag,ans,engv);
    #else
    store_answers_p(f,energy,virial,ii,inum,tid,tpa_sq,offset,
                    eflag,vflag,ans,engv);
    #endif
  } // if ii
}