ComputeLjPmeSerial.C

Go to the documentation of this file.

#include "InfoStream.h"
#include "Node.h"
#include "PatchMap.h"
#include "PatchMap.inl"
#include "AtomMap.h"
#include "ComputeLjPmeSerial.h"
#include "ComputeLjPmeSerialMgr.decl.h"
#include "PatchMgr.h"
#include "Molecule.h"
#include "ReductionMgr.h"
#include "ComputeMgr.h"
#include "ComputeMgr.decl.h"
// #define DEBUGM
#define MIN_DEBUG_LEVEL 3
#include "Debug.h"
#include "SimParameters.h"
#include "Parameters.h"
#include "LJTable.h"
#include "WorkDistrib.h"
#include "varsizemsg.h"
#include "LjPmeCompute.h"
#include <stdlib.h>
#include <stdio.h>
#include <errno.h>


struct LjPmeSerialAtom {
  Position position;
  short vdwType;
  int id;
};

typedef Force LjPmeSerialForce;

class LjPmeSerialCoordMsg : public CMessage_LjPmeSerialCoordMsg {
public:
  int sourceNode;
  int numAtoms;
  Lattice lattice;
  LjPmeSerialAtom *coord;
};

class LjPmeSerialForceMsg : public CMessage_LjPmeSerialForceMsg {
public:
  BigReal energyNonbond;
  BigReal energySlow;
  BigReal virialNonbond[3][3], virialSlow[3][3];
  LjPmeSerialForce *forceNonbond;
  LjPmeSerialForce *forceSlow;
};

class ComputeLjPmeSerialMgr : public CBase_ComputeLjPmeSerialMgr {
public:
  ComputeLjPmeSerialMgr();
  ~ComputeLjPmeSerialMgr();

  void setCompute(ComputeLjPmeSerial *c) { ljPmeCompute = c; }

  void recvCoord(LjPmeSerialCoordMsg *);
  void recvForce(LjPmeSerialForceMsg *);
  void getLJparameters(const int &i, const int &j,
        Real &A, Real &B, Real &A14, Real &B14);
        void get_vdw_parameters(Real &sigma, Real &epsilon, Real &sigma14, 
                            Real &epsilon14, Index index);
  void setLJparameters();

private:
  CProxy_ComputeLjPmeSerialMgr ljPmeProxy;
  ComputeLjPmeSerial *ljPmeCompute;

  int numSources;
  int numArrived;
  LjPmeSerialCoordMsg **coordMsgs;
  int numAtoms;
  int table_dim; // dimension of LJ table
  LjPmeSerialAtom *coord;
  LjPmeSerialForce *forceNonbond, *forceSlow;
  LjPmeSerialForceMsg *oldmsg;

  LjPmeCompute LjPme;
  double *ljPmeCoord; // stores atomic coordinate and equivalent of q in LJ-PME
  double *ljForceSlow; // slow LJ-PME force
  double *ljForceNonbond; // fast LJ-PME force
  double *ljParameter; // eps & sigma for nonbonded LJ paramters
  double *ljParameter14; // eps & sigma for 1-4 LJ paramters
  double oldNonbondedScaling; // check if nonbonded scaling has changed
  double oldSoluteScaling; // check if solute VDW scaling has changed
  int *atomType;
  Bool initialized;
};

ComputeLjPmeSerialMgr::ComputeLjPmeSerialMgr() :
  ljPmeProxy(thisgroup), ljPmeCompute(0), numSources(0), numArrived(0),
  coordMsgs(0), coord(0), forceNonbond(0), forceSlow(0), oldmsg(0),
  numAtoms(0), ljPmeCoord(0), ljForceSlow(0), ljForceNonbond(0),
  ljParameter(0), ljParameter14(0), atomType(0)
{
  initialized = false;
  CkpvAccess(BOCclass_group).computeLjPmeSerialMgr = thisgroup;
}

ComputeLjPmeSerialMgr::~ComputeLjPmeSerialMgr()
{
  initialized = false;
  for (int i=0;  i < numSources;  i++)  delete coordMsgs[i];
  delete [] coordMsgs;
  delete [] coord;
  delete [] forceNonbond;
  delete [] forceSlow;
  delete oldmsg;
  // extra
  if (ljPmeCoord) delete[] ljPmeCoord;
  if (ljForceSlow) delete[] ljForceSlow;
  if (ljForceNonbond) delete[] ljForceNonbond;
  if (ljParameter) delete[] ljParameter;
  if (ljParameter14) delete[] ljParameter14;
  if (atomType) delete[] atomType;
}

ComputeLjPmeSerial::~ComputeLjPmeSerial() {}

ComputeLjPmeSerial::ComputeLjPmeSerial(ComputeID c) :
  ComputeHomePatches(c)
{
  CProxy_ComputeLjPmeSerialMgr::ckLocalBranch(
        CkpvAccess(BOCclass_group).computeLjPmeSerialMgr)->setCompute(this);
  reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
}

void ComputeLjPmeSerialMgr::getLJparameters(const int &i, const int &j,
  Real &A, Real &B, Real &A14, Real &B14) {

  Parameters *params = Node::Object()->parameters;
  SimParameters *simParams = Node::Object()->simParameters;
  int useGeom = simParams->vdwGeometricSigma;
  Bool tabulatedEnergies = simParams->tabulatedEnergies;
  Bool soluteScalingOn = simParams->soluteScalingOn;
  BigReal soluteScalingFactor = simParams->soluteScalingFactorVdw;
  unsigned int table_dim_org = params->get_num_vdw_params();
  int ss_dim = Node::Object()->molecule->ss_num_vdw_params;
  int K = -1;
  int *ss_vdw_type = Node::Object()->molecule->ss_vdw_type;
  // BigReal sigma_max;
  //  We need the A and B parameters for the Van der Waals.  These can
  //  be explicitly be specified for this pair or calculated from the
  //  sigma and epsilon values for the two atom types
  //  printf("Looking at interaction of  %i with %i\n", i, j);
  if ( tabulatedEnergies && params->get_table_pair_params(i,j,&K)) {

    if ( K < 0 ) { NAMD_bug(
        "LJPME::compute_vdw_params: energy table index is negative");
    }
    A = -1 - K;
    B = 0;
    A14 = -1 - K;
    B14 = 0;
  } else if (params->get_vdw_pair_params(i,j, &A, &B, &A14, &B14)) {
    if ( tabulatedEnergies && ( A < 0 || A14 < 0 ) ) {
      NAMD_die("LJ A is negative with tabulatedEnergies enabled");
    }
  } else {
    //  We didn't find explicit parameters for this pair. So instead,
    //  get the parameters for each atom type separately and use them
    //  to calculate the values we need
    Real sigma_i, sigma_i14, epsilon_i, epsilon_i14;
    Real sigma_j, sigma_j14, epsilon_j, epsilon_j14;

    if (!soluteScalingOn) {
      params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
                &epsilon_i14,i);
      params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14, 
                &epsilon_j14,j);
    } else {
      int i_type = (i >= table_dim_org)? ss_vdw_type[i-table_dim_org]:i;
      int j_type = (j >= table_dim_org)? ss_vdw_type[j-table_dim_org]:j;
      params->get_vdw_params(&sigma_i, &epsilon_i, &sigma_i14,
                                          &epsilon_i14,i_type);
      params->get_vdw_params(&sigma_j, &epsilon_j, &sigma_j14,
                                          &epsilon_j14,j_type);
    }   
  
    BigReal sigma_ij =
        useGeom ? sqrt(sigma_i*sigma_j) : 0.5*(sigma_i+sigma_j);
    BigReal sigma_ij14 =
        useGeom ? sqrt(sigma_i14*sigma_j14) : 0.5 * (sigma_i14+sigma_j14);
    BigReal epsilon_ij = sqrt(epsilon_i*epsilon_j);
    BigReal epsilon_ij14 = sqrt(epsilon_i14*epsilon_j14);

    //  Calculate sigma^6
    sigma_ij *= sigma_ij*sigma_ij;
    sigma_ij *= sigma_ij;
    sigma_ij14 *= sigma_ij14*sigma_ij14;
    sigma_ij14 *= sigma_ij14;
      
    //  Calculate LJ constants A & B
    B = 4.0 * sigma_ij * epsilon_ij;
    A = B * sigma_ij;
    B14 = 4.0 * sigma_ij14 * epsilon_ij14;
    A14 = B14 * sigma_ij14;

    if (soluteScalingOn) {
      if (i >= table_dim_org && i < (table_dim_org+ss_dim) && j < table_dim_org) {
        A *= sqrt(soluteScalingFactor);
        B *= sqrt(soluteScalingFactor);
        A14 *= sqrt(soluteScalingFactor);
        B14 *= sqrt(soluteScalingFactor);
      }
      if (i < table_dim_org && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
        A *= sqrt(soluteScalingFactor);
        B *= sqrt(soluteScalingFactor);
        A14 *= sqrt(soluteScalingFactor);
        B14 *= sqrt(soluteScalingFactor);
      }
      if (i >=table_dim_org && i < (table_dim_org+ss_dim) && j >= table_dim_org && j < (table_dim_org+ss_dim)) {
        A *= soluteScalingFactor;
        B *= soluteScalingFactor;
        A14 *= soluteScalingFactor;
        B14 *= soluteScalingFactor;
      }
    }

    if ( tabulatedEnergies && ( A < 0 || A14 < 0 ) ) {
      NAMD_die("LJPME: LJ A is negative with tabulatedEnergies enabled");
    }
  }
}

void ComputeLjPmeSerialMgr::get_vdw_parameters(Real &sigma, Real &epsilon, Real &sigma14, 
        Real &epsilon14, Index index) {
  // Get the index index interaction
  Real A, B, A14, B14;
  this->getLJparameters(index, index, A, B, A14, B14);
  // Copied from get_vdw_params in Paramters.h
  sigma = sqrt(cbrt(A)) / sqrt(cbrt(B));
  epsilon = B*B / (4*A);

  sigma14 = sqrt(cbrt(A14)) / sqrt(cbrt(B14));
  epsilon14 = B14*B14 / (4*A14);
}

void ComputeLjPmeSerialMgr::setLJparameters() {
  SimParameters *simParams = Node::Object()->simParameters;
  Real scaling = simParams->nonbondedScaling;
  // These won't change
  for (int i = 0; i < table_dim; i++) {
    for (int j = i; j < table_dim; j++) {
      int index_ij = 2*(i*table_dim+j);
      int index_ji = 2*(j*table_dim+i);
      Real A, B, A14, B14;
      this->getLJparameters(i, j, A, B, A14, B14);
      ljParameter[index_ij] = scaling * A; // repulsion
      ljParameter[index_ij+1] = scaling * B; // attraction
      ljParameter14[index_ij] = scaling * A14; // 1-4 repulsion
      ljParameter14[index_ij+1] = scaling * B14; // 1-4 attraction
      // Copy to transpose entry
      ljParameter[index_ji] = ljParameter[index_ij];
      ljParameter[index_ji+1] = ljParameter[index_ij+1];
      ljParameter14[index_ji] = ljParameter14[index_ij];
      ljParameter14[index_ji+1] = ljParameter14[index_ij+1];
  
      //printf("Table Dim: %d, i: %d, j: %d \n", table_dim, i, j);
      //printf("A: %3.6f, B: %3.6f, A14: %3.6f, B14: %3.6f \n", A, B, A14, B14);
    }
  }

  // Store the dispersion part of LJ using Geometric sigma combining rules
  for (int i = 0;  i < numAtoms;  i++) {
    Real sigma, sigma_14, epsilon, epsilon_14;
    atomType[i] = coord[i].vdwType;
    this->get_vdw_parameters(sigma, epsilon, sigma_14,
        epsilon_14, coord[i].vdwType);
    // printf("VDWtype: %d, eps: %3.6f, sigma: %3.6f, eps14: %3.6f, sigma14: %3.6f \n",
    //         atomType[i], epsilon, sigma, epsilon_14, sigma_14);

    // store sqrt(4*eps*sig^6) = 2*sqrt(eps)*sigma^3
    ljPmeCoord[4*i+3] = 2.0*sigma*sigma*sigma*sqrt(scaling * epsilon);
  }
}

void ComputeLjPmeSerial::doWork()
{
  ResizeArrayIter<PatchElem> ap(patchList);

  // For now we do the LJ-PME calculation every steps

  // Skip computations if nothing to do.
  // if ( ! patchList[0].p->flags.doFullElectrostatics )
  // {
  //   for (ap = ap.begin(); ap != ap.end(); ap++) {
  //     CompAtom *x = (*ap).positionBox->open();
  //     Results *r = (*ap).forceBox->open();
  //     (*ap).positionBox->close(&x);
  //     (*ap).forceBox->close(&r);
  //   }
  //   reduction->submit();
  //   return;
  // }

  // allocate message
  int numLocalAtoms = 0;
  for (ap = ap.begin(); ap != ap.end(); ap++) {
    numLocalAtoms += (*ap).p->getNumAtoms();
  }

  LjPmeSerialCoordMsg *msg = new (numLocalAtoms, 0) LjPmeSerialCoordMsg;
  msg->sourceNode = CkMyPe();
  msg->numAtoms = numLocalAtoms;
  msg->lattice = patchList[0].p->flags.lattice;
  LjPmeSerialAtom *data_ptr = msg->coord;

  // get positions
  for (ap = ap.begin();  ap != ap.end();  ap++) {
    CompAtom *x = (*ap).positionBox->open();
    CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
    if ( patchList[0].p->flags.doMolly ) {
      (*ap).positionBox->close(&x);
      x = (*ap).avgPositionBox->open();
    }
    int numAtoms = (*ap).p->getNumAtoms();

    for(int i=0;  i < numAtoms;  i++)
    {
      data_ptr->position = x[i].position;
      data_ptr->vdwType = x[i].vdwType;
      data_ptr->id = xExt[i].id;
      ++data_ptr;
    }

    if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
    else { (*ap).positionBox->close(&x); }
  }

  CProxy_ComputeLjPmeSerialMgr ljPmeProxy(
      CkpvAccess(BOCclass_group).computeLjPmeSerialMgr);
  ljPmeProxy[0].recvCoord(msg);
}


void ComputeLjPmeSerialMgr::recvCoord(LjPmeSerialCoordMsg *msg) {
  SimParameters *simParams = Node::Object()->simParameters;
  Parameters *params = Node::Object()->parameters;
  if ( ! numSources ) {
    numSources = (PatchMap::Object())->numNodesWithPatches();
    coordMsgs = new LjPmeSerialCoordMsg*[numSources];
    for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
    numArrived = 0;
    numAtoms = Node::Object()->molecule->numAtoms;
    // Allocate memory to receive coords and send force
    coord = new LjPmeSerialAtom[numAtoms];
    forceNonbond = new LjPmeSerialForce[numAtoms];
    forceSlow = new LjPmeSerialForce[numAtoms];
  }

  // Receive the message at store using global atom index
  for (int i=0; i < msg->numAtoms; ++i) {
    coord[msg->coord[i].id] = msg->coord[i];
  }
  coordMsgs[numArrived] = msg;
  ++numArrived;

  if ( numArrived < numSources ) return;
  numArrived = 0;

  //###########################################
  // ALL DATA ARRIVED --- CALCULATE FORCES
  Lattice lattice = msg->lattice;

  if(!initialized) {
    initialized = true;
    // Allocate memory for LJ-PME
    table_dim = params->get_num_vdw_params();
    ljParameter = new double[2*table_dim*table_dim];
    ljParameter14 = new double[2*table_dim*table_dim];
    ljPmeCoord = new double[4*numAtoms];
    ljForceSlow = new double[3*numAtoms];
    ljForceNonbond = new double[3*numAtoms];
    atomType = new int[numAtoms];

    if (ljPmeCoord==0 || ljForceSlow==0 || ljForceNonbond==0 || ljParameter==0 
        || ljParameter==0 || atomType == 0) {
      NAMD_die("can't allocate LJ-PME atom buffers");
    }

    // Set the data for LJPME compute. To save memory, all arrays in LjPmeCompute
    // are pointing to arrayes in ComputeLjPmeSerialMgr class
    LjPme.initialize(ljPmeCoord, ljParameter, ljParameter14, atomType,
                    ljForceNonbond, ljForceSlow, table_dim, numAtoms);
    // Set LJ parameters
    this->setLJparameters();
    oldNonbondedScaling = simParams->nonbondedScaling;
    oldSoluteScaling = simParams->soluteScalingFactorVdw;
  }

  // Nonbonded scaling or solute VDW scaling has changed, so update the LJ parameters
  if((simParams->nonbondedScaling != oldNonbondedScaling) ||
    (simParams->soluteScalingOn && (simParams->soluteScalingFactorVdw != oldSoluteScaling))) {
    this->setLJparameters();
    oldNonbondedScaling = simParams->nonbondedScaling;
    oldSoluteScaling = simParams->soluteScalingFactorVdw;
  }

  // Update the coordinates and initialize the forces to zero
  for (int i = 0;  i < numAtoms;  i++) {
    ljPmeCoord[4*i  ] = coord[i].position.x;
    ljPmeCoord[4*i+1] = coord[i].position.y;
    ljPmeCoord[4*i+2] = coord[i].position.z;
    // Initialize the force to zero
    ljForceSlow[3*i  ] = 0.0;
    ljForceSlow[3*i+1] = 0.0;
    ljForceSlow[3*i+2] = 0.0;
    ljForceNonbond[3*i  ] = 0.0;
    ljForceNonbond[3*i+1] = 0.0;
    ljForceNonbond[3*i+2] = 0.0;
  }

  double energySlow, energyNonbond;
  double virialSlow[3][3], virialNonbond[3][3];
  LjPme.computeLJpotential(simParams->LJPMEEwaldCoefficient,
        simParams->cutoff,
        lattice,
        virialNonbond,
        virialSlow,
        energyNonbond,
        energySlow,
        true,
        true);
  
  if (0) {
    printf("LJ-PME: E_nb: %6.12f, E_rec: %6.12f, E_sum: %6.12f! \n", 
    energyNonbond, energySlow, energyNonbond+energySlow);
    printf("LJ-PME: virial_nb: [%3.6f, %3.6f, %3.6f, %3.6f, %3.6f, %3.6f] \n",
        virialNonbond[0][0], virialNonbond[0][1], virialNonbond[0][2],
        virialNonbond[1][1], virialNonbond[1][2], virialNonbond[2][2]);
    printf("LJ-PME: virial_slow: [%3.6f, %3.6f, %3.6f, %3.6f, %3.6f, %3.6f] \n",
        virialSlow[0][0], virialSlow[0][1], virialSlow[0][2],
        virialSlow[1][1], virialSlow[1][2], virialSlow[2][2]);
  }

  for (int i = 0;  i < numAtoms;  i++) {
    forceSlow[i].x = ljForceSlow[3*i  ];
    forceSlow[i].y = ljForceSlow[3*i+1];
    forceSlow[i].z = ljForceSlow[3*i+2];

    forceNonbond[i].x = ljForceNonbond[3*i  ];
    forceNonbond[i].y = ljForceNonbond[3*i+1];
    forceNonbond[i].z = ljForceNonbond[3*i+2];

    if (0) {
      forceNonbond[i].x = 0.0;
      forceNonbond[i].y = 0.0;
      forceNonbond[i].z = 0.0;
      double fx = forceNonbond[i].x + forceSlow[i].x;
      double fy = forceNonbond[i].y + forceSlow[i].y;
      double fz = forceNonbond[i].z + forceSlow[i].z;
      printf("LJ-PME send: atom %5d<-> nb: (%3.6f,%3.6f,%3.6f), slow:(%3.6f,%3.6f,%3.6f), sum:(%3.6f,%3.6f,%3.6f)! \n", i,
      forceNonbond[i].x, forceNonbond[i].y, forceNonbond[i].z,
      forceSlow[i].x, forceSlow[i].y, forceSlow[i].z,
      fx, fy, fz);
    }
  }

  // distribute forces
  for (int j=0;  j < numSources;  j++) {
    LjPmeSerialCoordMsg *cmsg = coordMsgs[j];
    coordMsgs[j] = 0;
    LjPmeSerialForceMsg *fmsg = new (cmsg->numAtoms, cmsg->numAtoms, 0) LjPmeSerialForceMsg;

    for (int i=0;  i < cmsg->numAtoms;  i++) {
      fmsg->forceSlow[i] = forceSlow[cmsg->coord[i].id];
      fmsg->forceNonbond[i] = forceNonbond[cmsg->coord[i].id];
    }

    if ( ! j ) {  // set virial and energy only for first message
      fmsg->energyNonbond = energyNonbond;
      fmsg->energySlow = energySlow;
      for (int k=0;  k < 3;  k++) {
        for (int l=0;  l < 3;  l++) {
          fmsg->virialSlow[k][l] = virialSlow[k][l];
          fmsg->virialNonbond[k][l] = virialNonbond[k][l];
        }
      }
    }
    else {  // set other messages to zero, add into reduction only once
      fmsg->energyNonbond = 0;
      fmsg->energySlow = 0;
      for (int k=0;  k < 3;  k++) {
        for (int l=0;  l < 3;  l++) {
          fmsg->virialSlow[k][l] = 0;
          fmsg->virialNonbond[k][l] = 0;
        }
      }
    }

    ljPmeProxy[cmsg->sourceNode].recvForce(fmsg);
    delete cmsg;
  }
}

void ComputeLjPmeSerialMgr::recvForce(LjPmeSerialForceMsg *msg) {
  ljPmeCompute->saveResults(msg);
  delete oldmsg;
  oldmsg = msg;
}

void ComputeLjPmeSerial::saveResults(LjPmeSerialForceMsg *msg)
{
  ResizeArrayIter<PatchElem> ap(patchList);

  LjPmeSerialForce *results_Nonbond_ptr = msg->forceNonbond;
  LjPmeSerialForce *results_Slow_ptr = msg->forceSlow;

  // add in forces
  for (ap = ap.begin(); ap != ap.end(); ap++) {
    Results *r = (*ap).forceBox->open();
    Force *fSlow = r->f[Results::slow];
    Force *fNonbond = r->f[Results::nbond];
    int numAtoms = (*ap).p->getNumAtoms();

    for(int i=0; i<numAtoms; ++i) {
      fSlow[i].x += results_Slow_ptr->x;
      fSlow[i].y += results_Slow_ptr->y;
      fSlow[i].z += results_Slow_ptr->z;
      fNonbond[i].x += results_Nonbond_ptr->x;
      fNonbond[i].y += results_Nonbond_ptr->y;
      fNonbond[i].z += results_Nonbond_ptr->z;
      ++results_Slow_ptr;
      ++results_Nonbond_ptr;
    }
  
    (*ap).forceBox->close(&r);
  }

  // reduction->item(REDUCTION_LJ_ENERGY) += msg->energyNonbond;
  // reduction->item(REDUCTION_ELECT_ENERGY_SLOW) +=  msg->energySlow;

  // reduction->item(REDUCTION_VIRIAL_SLOW_XX) += msg->virialSlow[0][0];
  // reduction->item(REDUCTION_VIRIAL_SLOW_XY) += msg->virialSlow[0][1];
  // reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += msg->virialSlow[0][2];
  // reduction->item(REDUCTION_VIRIAL_SLOW_YX) += msg->virialSlow[1][0];
  // reduction->item(REDUCTION_VIRIAL_SLOW_YY) += msg->virialSlow[1][1];
  // reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += msg->virialSlow[1][2];
  // reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += msg->virialSlow[2][0];
  // reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += msg->virialSlow[2][1];
  // reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->virialSlow[2][2];

  // reduction->item(REDUCTION_VIRIAL_NBOND_XX) += msg->virialNonbond[0][0];
  // reduction->item(REDUCTION_VIRIAL_NBOND_XY) += msg->virialNonbond[0][1];
  // reduction->item(REDUCTION_VIRIAL_NBOND_XZ) += msg->virialNonbond[0][2];
  // reduction->item(REDUCTION_VIRIAL_NBOND_YX) += msg->virialNonbond[1][0];
  // reduction->item(REDUCTION_VIRIAL_NBOND_YY) += msg->virialNonbond[1][1];
  // reduction->item(REDUCTION_VIRIAL_NBOND_YZ) += msg->virialNonbond[1][2];
  // reduction->item(REDUCTION_VIRIAL_NBOND_ZX) += msg->virialNonbond[2][0];
  // reduction->item(REDUCTION_VIRIAL_NBOND_ZY) += msg->virialNonbond[2][1];
  // reduction->item(REDUCTION_VIRIAL_NBOND_ZZ) += msg->virialNonbond[2][2];

  reduction->item(REDUCTION_LJ_ENERGY) += msg->energyNonbond + msg->energySlow;

  reduction->item(REDUCTION_VIRIAL_NBOND_XX) += msg->virialNonbond[0][0] + msg->virialSlow[0][0];
  reduction->item(REDUCTION_VIRIAL_NBOND_XY) += msg->virialNonbond[0][1] + msg->virialSlow[0][1];
  reduction->item(REDUCTION_VIRIAL_NBOND_XZ) += msg->virialNonbond[0][2] + msg->virialSlow[0][2];
  reduction->item(REDUCTION_VIRIAL_NBOND_YX) += msg->virialNonbond[1][0] + msg->virialSlow[1][0];
  reduction->item(REDUCTION_VIRIAL_NBOND_YY) += msg->virialNonbond[1][1] + msg->virialSlow[1][1];
  reduction->item(REDUCTION_VIRIAL_NBOND_YZ) += msg->virialNonbond[1][2] + msg->virialSlow[1][2];
  reduction->item(REDUCTION_VIRIAL_NBOND_ZX) += msg->virialNonbond[2][0] + msg->virialSlow[2][0];
  reduction->item(REDUCTION_VIRIAL_NBOND_ZY) += msg->virialNonbond[2][1] + msg->virialSlow[2][1];
  reduction->item(REDUCTION_VIRIAL_NBOND_ZZ) += msg->virialNonbond[2][2] + msg->virialSlow[2][2];
  reduction->submit();
}

#include "ComputeLjPmeSerialMgr.def.h"