ComputeOneFourNbTholes.C

Go to the documentation of this file.

#include "ComputeNonbondedUtil.h"
#include "ComputeOneFourNbTholes.h"
#include "Molecule.h"
#include "common.h"
#include "structures.h"
#include "Debug.h"
#include "PressureProfile.h"

// static initialization
int OneFourNbTholeElem::pressureProfileSlabs = 0;
int OneFourNbTholeElem::pressureProfileAtomTypes = 1;
BigReal OneFourNbTholeElem::pressureProfileThickness = 0;
BigReal OneFourNbTholeElem::pressureProfileMin = 0;

void OneFourNbTholeElem::getMoleculePointers (
  Molecule* mol, int* count, int32*** byatom,
  OneFourNbThole** structarray) {
#ifdef MEM_OPT_VERSION
  NAMD_die("Should not be called in OneFourNbTholeElem::getMoleculePointers in memory optimized version!");
#else
  *count = mol->numOneFourNbTholes;
  *byatom = mol->oneFourNbTholesByAtom;
  *structarray = mol->oneFourNbTholes;
#endif
}

void OneFourNbTholeElem::getParameterPointers(
  Parameters *p, const OneFourNbTholeValue **v) {
  *v = NULL;
}

void OneFourNbTholeElem::computeForce(
  OneFourNbTholeElem *tuples, int ntuple,
  BigReal *reduction, BigReal *pressureProfileData) {
  const Lattice & lattice = tuples[0].p[0]->p->lattice;
  const bool doEnergy = tuples[0].p[0]->p->flags.doEnergy;
  const bool doVirial = tuples[0].p[0]->p->flags.doVirial;
  SimParameters *const simParams = Node::Object()->simParameters;
  /*
   * Although CHARMM Drude force field expects the 1-4 scaling factor to be 1.0,
   * I try my best to generalize the case here. The GPU nonbonded kernel differs
   * from the CPU kernel as the former always calculate the full 1-4 interactions
   * without scaling (or just with a scaling factor 1.0), and then the bonded
   * kernel "modifiedExclusionForce" in ComputeBondedCUDAKernel.cu actually
   * computes a complementary part scaled by -(1.0 - ComputeNonbondedUtil::scale14).
   *
   * As a result, the following scaling factor uses the same convention. In the
   * GPU kernel, the NbThole correction of 1-4 interactions is computed by
   * calcForceEnergyNbThole, so the following code should only computes the
   * complementary part.
   *
   * For the CPU kernel, the NbThole correction in ComputeNonbondedBase.h skips
   * modified exclusions, so we just compute the 1-4 NbThole correction with
   * scaling below.
   */
  BigReal s;
#if defined (NAMD_CUDA) || defined(NAMD_HIP)
  s = -(1.0 - ComputeNonbondedUtil::scale14);
#else
  s = ComputeNonbondedUtil::scale14;
#endif
  if (s == 0.0) return;
  const BigReal CC = COULOMB * s * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1;
  if (simParams->alchOn && ntuple > 0) {
    NAMD_die("Alchemical calculation for 1-4 atom pairs with nonbonded "
             "Thole parameters is not supported.");
  }
  for (int ituple=0; ituple<ntuple; ++ituple) {
    const OneFourNbTholeElem &tup = tuples[ituple];
    // const AtomID (&atomID)[OneFourNbTholeElem::size](tup.atomID);
    const int    (&localIndex)[OneFourNbTholeElem::size](tup.localIndex);
    TuplePatchElem * const(&p)[OneFourNbTholeElem::size](tup.p);
    const OneFourNbTholeValue* const(&value)(tup.value);
    DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
               << localIndex[1] << " " << localIndex[2] << " "
               << localIndex[3] << std::endl);
    const BigReal aa = value->aa;
    //  Calculate the vectors between atoms
    const Position & rai = p[0]->x[localIndex[0]].position;  // atom i
    const Position & rdi = p[1]->x[localIndex[1]].position;  // atom i's Drude
    const Position & raj = p[2]->x[localIndex[2]].position;  // atom j
    const Position & rdj = p[3]->x[localIndex[3]].position;  // atom j's Drude

    const BigReal qqaa = CC * p[0]->x[localIndex[0]].charge * p[2]->x[localIndex[2]].charge;
    const BigReal qqad = CC * p[0]->x[localIndex[0]].charge * p[3]->x[localIndex[3]].charge;
    const BigReal qqda = CC * p[1]->x[localIndex[1]].charge * p[2]->x[localIndex[2]].charge;
    const BigReal qqdd = CC * p[1]->x[localIndex[1]].charge * p[3]->x[localIndex[3]].charge;

    // r_ij = r_i - r_j
    const Vector raa = lattice.delta(rai,raj);  // shortest vector image:  rai - raj
    const Vector rad = lattice.delta(rai,rdj);  // shortest vector image:  rai - rdj
    const Vector rda = lattice.delta(rdi,raj);  // shortest vector image:  rdi - raj
    const Vector rdd = lattice.delta(rdi,rdj);  // shortest vector image:  rdi - rdj

    // 1/r, r = |r_ij|
    const BigReal raa_invlen = raa.rlength();  // reciprocal of length
    const BigReal rad_invlen = rad.rlength();
    const BigReal rda_invlen = rda.rlength();
    const BigReal rdd_invlen = rdd.rlength();

    // ar
    const BigReal auaa = aa / raa_invlen;
    const BigReal auad = aa / rad_invlen;
    const BigReal auda = aa / rda_invlen;
    const BigReal audd = aa / rdd_invlen;

    // exp(-ar)
    const BigReal expauaa = exp(-auaa);
    const BigReal expauad = exp(-auad);
    const BigReal expauda = exp(-auda);
    const BigReal expaudd = exp(-audd);

    // (1 + ar/2)
    BigReal polyauaa = 1 + 0.5*auaa;
    BigReal polyauad = 1 + 0.5*auad;
    BigReal polyauda = 1 + 0.5*auda;
    BigReal polyaudd = 1 + 0.5*audd;

    BigReal ethole;
    if (doEnergy) {
      ethole = 0;
      ethole += qqaa * raa_invlen * (-polyauaa * expauaa);
      ethole += qqad * rad_invlen * (-polyauad * expauad);
      ethole += qqda * rda_invlen * (-polyauda * expauda);
      ethole += qqdd * rdd_invlen * (-polyaudd * expaudd);
    }

    polyauaa = 1 + auaa*polyauaa;
    polyauad = 1 + auad*polyauad;
    polyauda = 1 + auda*polyauda;
    polyaudd = 1 + audd*polyaudd;

    const BigReal raa_invlen3 = raa_invlen * raa_invlen * raa_invlen;
    const BigReal rad_invlen3 = rad_invlen * rad_invlen * rad_invlen;
    const BigReal rda_invlen3 = rda_invlen * rda_invlen * rda_invlen;
    const BigReal rdd_invlen3 = rdd_invlen * rdd_invlen * rdd_invlen;

    const BigReal dfaa = qqaa * raa_invlen3 * (polyauaa*expauaa);
    const BigReal dfad = qqad * rad_invlen3 * (polyauad*expauad);
    const BigReal dfda = qqda * rda_invlen3 * (polyauda*expauda);
    const BigReal dfdd = qqdd * rdd_invlen3 * (polyaudd*expaudd);

    const Vector faa = -dfaa * raa;
    const Vector fad = -dfad * rad;
    const Vector fda = -dfda * rda;
    const Vector fdd = -dfdd * rdd;

    p[0]->f[localIndex[0]] += faa + fad;
    p[1]->f[localIndex[1]] += fda + fdd;
    p[2]->f[localIndex[2]] -= faa + fda;
    p[3]->f[localIndex[3]] -= fad + fdd;

    if (doEnergy) {
      reduction[OneFourNbTholeEnergyIndex] += ethole;
    }

    if (doVirial) {
      reduction[virialIndex_XX] += faa.x * raa.x + fad.x * rad.x
        + fda.x * rda.x + fdd.x * rdd.x;
      reduction[virialIndex_XY] += faa.x * raa.y + fad.x * rad.y
        + fda.x * rda.y + fdd.x * rdd.y;
      reduction[virialIndex_XZ] += faa.x * raa.z + fad.x * rad.z
        + fda.x * rda.z + fdd.x * rdd.z;
      reduction[virialIndex_YX] += faa.y * raa.x + fad.y * rad.x
        + fda.y * rda.x + fdd.y * rdd.x;
      reduction[virialIndex_YY] += faa.y * raa.y + fad.y * rad.y
        + fda.y * rda.y + fdd.y * rdd.y;
      reduction[virialIndex_YZ] += faa.y * raa.z + fad.y * rad.z
        + fda.y * rda.z + fdd.y * rdd.z;
      reduction[virialIndex_ZX] += faa.z * raa.x + fad.z * rad.x
        + fda.z * rda.x + fdd.z * rdd.x;
      reduction[virialIndex_ZY] += faa.z * raa.y + fad.z * rad.y
        + fda.z * rda.y + fdd.z * rdd.y;
      reduction[virialIndex_ZZ] += faa.z * raa.z + fad.z * rad.z
        + fda.z * rda.z + fdd.z * rdd.z;
    }

    if (pressureProfileData) {
      BigReal zai = p[0]->x[localIndex[0]].position.z;
      BigReal zdi = p[1]->x[localIndex[1]].position.z;
      BigReal zaj = p[2]->x[localIndex[2]].position.z;
      BigReal zdj = p[3]->x[localIndex[3]].position.z;
      int nai = (int)floor((zai-pressureProfileMin)/pressureProfileThickness);
      int ndi = (int)floor((zdi-pressureProfileMin)/pressureProfileThickness);
      int naj = (int)floor((zaj-pressureProfileMin)/pressureProfileThickness);
      int ndj = (int)floor((zdj-pressureProfileMin)/pressureProfileThickness);
      pp_clamp(nai, pressureProfileSlabs);
      pp_clamp(ndi, pressureProfileSlabs);
      pp_clamp(naj, pressureProfileSlabs);
      pp_clamp(ndj, pressureProfileSlabs);
      int pai = p[0]->x[localIndex[0]].partition;
      int pdi = p[1]->x[localIndex[1]].partition;
      int paj = p[2]->x[localIndex[2]].partition;
      int pdj = p[3]->x[localIndex[3]].partition;
      int pn = pressureProfileAtomTypes;
      pp_reduction(pressureProfileSlabs, nai, naj,
          pai, paj, pn, faa.x * raa.x, faa.y * raa.y, faa.z * raa.z,
          pressureProfileData);
      pp_reduction(pressureProfileSlabs, nai, ndj,
          pai, pdj, pn, fad.x * rad.x, fad.y * rad.y, fad.z * rad.z,
          pressureProfileData);
      pp_reduction(pressureProfileSlabs, ndi, naj,
          pdi, paj, pn, fda.x * rda.x, fda.y * rda.y, fda.z * rda.z,
          pressureProfileData);
      pp_reduction(pressureProfileSlabs, ndi, ndj,
          pdi, pdj, pn, fdd.x * rdd.x, fdd.y * rdd.y, fdd.z * rdd.z,
          pressureProfileData);
    }
  }
}

void OneFourNbTholeElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
{
  reduction->item(REDUCTION_ELECT_ENERGY) += data[OneFourNbTholeEnergyIndex];
  reduction->item(REDUCTION_ELECT_ENERGY_F) += data[OneFourNbTholeEnergyIndex_f];
  reduction->item(REDUCTION_ELECT_ENERGY_TI_1) += data[OneFourNbTholeEnergyIndex_ti_1];
  reduction->item(REDUCTION_ELECT_ENERGY_TI_2) += data[OneFourNbTholeEnergyIndex_ti_2];
  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
}

int OneFourNbTholeElem::hash() const {
  return 0x7FFFFFFF &((atomID[0]<<24) + (atomID[1]<<16) + (atomID[2]<<8) + atomID[3]);
}

OneFourNbTholeElem::OneFourNbTholeElem() {
  scale = 0;
}

OneFourNbTholeElem::OneFourNbTholeElem(AtomID atom0, const TupleSignature *sig, const OneFourNbTholeValue *v) {
  NAMD_die("Can't use OneFourNbThole with memory optimized version of NAMD.");
}

OneFourNbTholeElem::OneFourNbTholeElem(const OneFourNbThole *a, const OneFourNbTholeValue *v) {
  atomID[0] = a->atom1;
  atomID[1] = a->atom2;
  atomID[2] = a->atom3;
  atomID[3] = a->atom4;
  value = a;  // expect v to be NULL
}

OneFourNbTholeElem::OneFourNbTholeElem(
  AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3) {
  // atoms arranged:  HEAVY DRUDE HEAVY DRUDE
  if (atom0 > atom2) {  // Swap heavy atoms so lowest is first!
    std::swap(atom0, atom2);
    std::swap(atom1, atom3);
  }
  atomID[0] = atom0;
  atomID[1] = atom1;
  atomID[2] = atom2;
  atomID[3] = atom3;
}

inline int OneFourNbTholeElem::operator==(const OneFourNbTholeElem &a) const {
  return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
      a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3]);
}

inline int OneFourNbTholeElem::operator<(const OneFourNbTholeElem &a) const {
  return  (atomID[0] < a.atomID[0] ||
          (atomID[0] == a.atomID[0] &&
          (atomID[1] < a.atomID[1] ||
          (atomID[1] == a.atomID[1] &&
          (atomID[2] < a.atomID[2] ||
          (atomID[2] == a.atomID[2] &&
            atomID[3] < a.atomID[3]
          ))))));
}