NAMD
Public Member Functions | List of all members
LjPmeCompute Class Reference

#include <LjPmeCompute.h>

Public Member Functions

 LjPmeCompute ()
 
 ~LjPmeCompute ()
 
void initialize (const double *pmeCoord, const double *parameter, const double *parameter14, const int *type, double *forceNonbond, double *forceSlow, const int &ljTableDim, const int &nAtoms)
 
void computeLJpotential (const double &alphaLJ, const double &cutoff, const Lattice &lattice, double virialNonbonded[][3], double virialSlow[][3], double &energyNonbond, double &energySlow, bool doEnergy, bool doSlow)
 
template<bool doSlow, bool doEnergy>
void computeNonbonded (const double &alphaLJ, const double &cutoff, const Lattice &lattice, double virialNonbonded[][3], double virialSlow[][3], double &energyNonbond, double &energySlow)
 

Detailed Description

Definition at line 20 of file LjPmeCompute.h.

Constructor & Destructor Documentation

◆ LjPmeCompute()

LjPmeCompute::LjPmeCompute ( )
inline

Definition at line 22 of file LjPmeCompute.h.

22  :
23  ljPmeCoord(NULL), ljParameter(NULL),
24  ljParameter14(NULL), atomType(NULL),
25  ljForceNonbond(NULL), ljForceSlow(NULL)
26  {
27  initialized = false;
28  //atomChanged = true;
29  numAtoms = 0;
30  }

◆ ~LjPmeCompute()

LjPmeCompute::~LjPmeCompute ( )
inline

Definition at line 32 of file LjPmeCompute.h.

32  {
33  ljPmeCoord = NULL; ljParameter = NULL;
34  ljParameter14 = NULL; atomType = NULL;
35  ljForceNonbond = NULL; ljForceSlow = NULL;
36  numAtoms = 0;
37  initialized = false;
38  }

Member Function Documentation

◆ computeLJpotential()

void LjPmeCompute::computeLJpotential ( const double &  alphaLJ,
const double &  cutoff,
const Lattice &  lattice,
double  virialNonbonded[][3],
double  virialSlow[][3],
double &  energyNonbond,
double &  energySlow,
bool  doEnergy,
bool  doSlow 
)
Parameters
alphaLJalpha/beta value in LJ-PME
cutoffrcut distance
latticeCell lattice
virialNonbondedNonbonded virial 3x3
virialSlowSlow virial 3x3

Definition at line 34 of file LjPmeCompute.C.

References LjPmeMgr::computeLongRange(), and NAMD_die().

Referenced by ComputeLjPmeSerialMgr::recvCoord().

38  {
39  if(!initialized) {
40  NAMD_die("LjPmeCompute has been called without initialization!");
41  }
42  // Initialize the energy and virial values to zero
43  energySlow = energyNonbond = 0.0;
44  for(int i = 0; i < 3; ++i) {
45  for(int j = 0; j < 3; ++j) {
46  virialNonbonded[i][j] = 0.0;
47  virialSlow[i][j] = 0.0;
48  }
49  }
50 
51  // Calculate the short range and long range Nonbonded LJ potential
52  if(doSlow) {
53  if (doEnergy) {
54  // Calculate the LJ nonBonded term
55  this->computeNonbonded<true, true>(alphaLJ, cutoff,
56  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
57  } else {
58  energySlow = 0.0;
59  this->computeNonbonded<true, false>(alphaLJ, cutoff,
60  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
61  }
62  // Compute reciprocal force, virial, energy, and self energy
63  ljPmeMgr.computeLongRange(ljPmeCoord, lattice, alphaLJ,
64  ljForceSlow, energySlow, virialSlow, doEnergy);
65  } else {
66  if (doEnergy) {
67  this->computeNonbonded<false, true>(alphaLJ, cutoff,
68  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
69  } else {
70  this->computeNonbonded<false, false>(alphaLJ, cutoff,
71  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
72  }
73  }
74 }
LjPmeMgr::computeLongRange
void computeLongRange(const double *ljPmeCoord, const Lattice &lattice, const double &alphaLJ, double *force, double &energy, double virial[][3], bool doEnergy)
Definition: LjPmeMgr.C:158
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:147

◆ computeNonbonded()

template<bool doSlow, bool doEnergy>
void LjPmeCompute::computeNonbonded ( const double &  alphaLJ,
const double &  cutoff,
const Lattice &  lattice,
double  virialNonbonded[][3],
double  virialSlow[][3],
double &  energyNonbond,
double &  energySlow 
)

< distance square
< Flag to check if the pair is excluded (1), scaled1-4(2), or not (0)

Parameters
alphaLJalpha/beta value in LJ-PME
cutoffrcut distance
latticeCell lattice
virialNonbondedNonbonded virial 3x3
virialSlowSlow virial 3x3

Definition at line 77 of file LjPmeCompute.C.

References Molecule::checkexcl(), Lattice::delta(), Vector::length2(), Node::molecule, Node::Object(), Vector::x, Vector::y, and Vector::z.

80  {
81  Molecule *molecule = Node::Object()->molecule;
82  double dispersion = 0.0; double repulsion = 0.0;
83  double nonBondEnergy = 0.0;
84  double slowEnergy = 0.0;
85  double r2 = 0.0;
86  double rcut2 = cutoff * cutoff;
87  double rcut2inv = 1.0 / rcut2;
88  double rcut6inv = rcut2inv * rcut2inv * rcut2inv;
89  double rcut12inv = rcut6inv * rcut6inv;
90  double aRc = alphaLJ * cutoff;
91  double aRc2 = aRc * aRc;
92  int exclusion = 0;
93  // Screening function at cutoff distance (1 - g(alpha*cutoff))
94  double screen_rc = (1.0 - (1 + aRc2 + 0.5*aRc2*aRc2)*exp(-aRc2));
95  Vector iCoord, jCoord, distance;
96  for(int i = 0; i < numAtoms; ++i) {
97  for(int j = i+1; j < numAtoms; ++j) {
98  iCoord = Vector(ljPmeCoord[4*i], ljPmeCoord[4*i+1], ljPmeCoord[4*i+2]);
99  jCoord = Vector(ljPmeCoord[4*j], ljPmeCoord[4*j+1], ljPmeCoord[4*j+2]);
100  // Shortest vector from jCoords to iCoords (iCoord - jCoords)
101  distance = lattice.delta(iCoord, jCoord);
102  // Distance square
103  r2 = distance.length2();
104 
105  if (r2 < rcut2) {
106 
107  // Check exclusion flag; 0: no exclusion. 1: full exclusion. 2: 1-4 scaling
108  #ifdef MEM_OPT_VERSION
109  // For now we dont support MEM_OPT
110  exclusion = molecule->checkExclByIdx(i, i, j); // Not sure if this is correct!
111  #else
112  exclusion = molecule->checkexcl(i, j);
113  #endif
114 
115  // Check if the pair of atoms have full exclusion (1)
116  if(exclusion == 1) {
117  // Remove the reciprocal contribution for excluded pairs
118  if(doSlow ) {
119  double r2inv = 1.0/r2;
120  double r6inv = r2inv * r2inv * r2inv;
121  double r12inv = r6inv * r6inv;
122  double aR2 = r2 * alphaLJ * alphaLJ;
123  double aR4 = aR2 * aR2;
124  // Get the C6 term (4*eps*sig^6) using geometric combining rule
125  double c6 = ljPmeCoord[4*i+3] * ljPmeCoord[4*j+3];
126  // Derivative of Screening function (1 - g'(alpha*r))
127  double screen_dr = (1.0 - (1.0 + aR2 + 0.5*aR4 + aR4*aR2/6.0)*exp(-aR2));
128  // Subtract LJ dispersion force that explicitly included in reciprocal sum
129  // Check! I think this should be positive! OpenMM uses negative sign!
130  double slowForce = 6.0 * c6 * screen_dr * r6inv * r2inv;
131 
132  if(doEnergy) {
133  // Screening function (1 - g(alpha*r))
134  double screen_r = (1.0 - (1.0 + aR2 + 0.5*aR4)*exp(-aR2));
135  // Subtract LJ dispersion energy that explicitly included in reciprocal sum
136  // Since dispersion is negative value, subtraction turns to addition
137  slowEnergy += c6 * screen_r * r6inv;
138  }
139 
140  // Store the force
141  Vector force_r = slowForce * distance;
142  ljForceSlow[3*i] += force_r.x;
143  ljForceSlow[3*i+1] += force_r.y;
144  ljForceSlow[3*i+2] += force_r.z;
145  ljForceSlow[3*j] -= force_r.x;
146  ljForceSlow[3*j+1] -= force_r.y;
147  ljForceSlow[3*j+2] -= force_r.z;
148  // Store virial
149  virialSlow[0][0] += force_r.x * distance.x;
150  virialSlow[0][1] += force_r.x * distance.y;
151  virialSlow[0][2] += force_r.x * distance.z;
152  virialSlow[1][0] += force_r.y * distance.x;
153  virialSlow[1][1] += force_r.y * distance.y;
154  virialSlow[1][2] += force_r.y * distance.z;
155  virialSlow[2][0] += force_r.z * distance.x;
156  virialSlow[2][1] += force_r.z * distance.y;
157  virialSlow[2][2] += force_r.z * distance.z;
158  }
159 
160  } else /* if(r2 < rcut2) */ { // Make sure the distance is within the cutoff if pairs are not excluded (0, 2)
161  int iType = atomType[i];
162  int jType = atomType[j];
163  int typeIndex = 2*(iType * tableDim + jType);
164  double r2inv = 1.0/r2;
165  double r6inv = r2inv * r2inv * r2inv;
166  double r12inv = r6inv * r6inv;
167  if (exclusion == 2) { // Use 1-4 parameters
168  repulsion = ljParameter14[typeIndex];
169  dispersion = ljParameter14[typeIndex + 1];
170  } else { // No full and 1-4 exclusion
171  repulsion = ljParameter[typeIndex];
172  dispersion = ljParameter[typeIndex + 1];
173  }
174  // LJ energy
175  if(doEnergy) {
176  // Add standard LJ energy term
177  nonBondEnergy += (repulsion*r12inv - dispersion*r6inv);
178  }
179  // LJ force
180  double nonBondForce = (12.0*repulsion*r12inv - 6.0*dispersion*r6inv)*r2inv;
181 
182  // LJ dispersion approximation contribution to force and energy
183  if(doSlow) {
184  double aR2 = r2 * alphaLJ * alphaLJ;
185  double aR4 = aR2 * aR2;
186  // Get the C6 term (4*eps*sig^6) using geometric combining rule
187  double c6 = ljPmeCoord[4*i+3] * ljPmeCoord[4*j+3];
188  // Derivative of Screening function (1 - g'(alpha*r))
189  double screen_dr = (1.0 - (1.0 + aR2 + 0.5*aR4 + aR4*aR2/6.0)*exp(-aR2));
190  // Add LJ dispersion approximation contribution to force
191  nonBondForce += 6.0 * c6 * screen_dr * r6inv * r2inv;
192 
193  if(doEnergy) {
194  // Screening function (1 - g(alpha*r))
195  double screen_r = (1.0 - (1.0 + aR2 + 0.5*aR4)*exp(-aR2));
196  // Add LJ dispersion approximation contribution to energy
197  nonBondEnergy += c6 * screen_r * r6inv;
198  /* Calculate the shifted energy */
199  // LJ energy at cutoff distance
200  double potentialShift = (repulsion*rcut12inv - dispersion*rcut6inv);
201  // Add the LJ dispersion approximation energy at cutoff
202  potentialShift += c6 * screen_rc * rcut6inv;
203  // Subtract the potentialShift to have zero energy at cutoff
204  nonBondEnergy -= potentialShift;
205  }
206  }
207 
208  // Store the force
209  Vector force_r = nonBondForce * distance;
210  ljForceNonbond[3*i] += force_r.x;
211  ljForceNonbond[3*i+1] += force_r.y;
212  ljForceNonbond[3*i+2] += force_r.z;
213  ljForceNonbond[3*j] -= force_r.x;
214  ljForceNonbond[3*j+1] -= force_r.y;
215  ljForceNonbond[3*j+2] -= force_r.z;
216 
217  // Store virial
218  virialNonbonded[0][0] += force_r.x * distance.x;
219  virialNonbonded[0][1] += force_r.x * distance.y;
220  virialNonbonded[0][2] += force_r.x * distance.z;
221  virialNonbonded[1][0] += force_r.y * distance.x;
222  virialNonbonded[1][1] += force_r.y * distance.y;
223  virialNonbonded[1][2] += force_r.y * distance.z;
224  virialNonbonded[2][0] += force_r.z * distance.x;
225  virialNonbonded[2][1] += force_r.z * distance.y;
226  virialNonbonded[2][2] += force_r.z * distance.z;
227  }
228 
229  }
230  }
231  }
232  energyNonbond += nonBondEnergy;
233  energySlow += slowEnergy;
234 }
Node::Object
static Node * Object()
Definition: Node.h:86
Molecule::checkexcl
int checkexcl(int atom1, int atom2) const
Vector
Definition: Vector.h:72
Vector::z
BigReal z
Definition: Vector.h:74
Molecule
Molecule stores the structural information for the system.
Definition: Molecule.h:175
Vector::length2
NAMD_HOST_DEVICE BigReal length2(void) const
Definition: Vector.h:206
Vector::x
BigReal x
Definition: Vector.h:74
Vector::y
BigReal y
Definition: Vector.h:74
Node::molecule
Molecule * molecule
Definition: Node.h:179
Lattice::delta
NAMD_HOST_DEVICE Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:149

◆ initialize()

void LjPmeCompute::initialize ( const double *  pmeCoord,
const double *  parameter,
const double *  parameter14,
const int *  type,
double *  forceNonbond,
double *  forceSlow,
const int &  ljTableDim,
const int &  nAtoms 
)

Setting up the data needed for LJ calculation.

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 13 of file LjPmeCompute.C.

References LjPmeMgr::initialize(), NAMD_die(), Node::Object(), and Node::simParameters.

Referenced by ComputeLjPmeSerialMgr::recvCoord().

16  {
17  if (! initialized) {
18  ljPmeCoord = pmeCoord;
19  ljParameter = parameter;
20  ljParameter14 = parameter14;
21  atomType = type;
22  ljForceNonbond = forceNonbond;
23  ljForceSlow = forceSlow;
24  tableDim = ljTableDim;
25  numAtoms = nAtoms;
26  SimParameters *simParm = Node::Object()->simParameters;
27  ljPmeMgr.initialize(simParm, numAtoms);
28  initialized = true;
29  } else {
30  NAMD_die("LjPmeCompute has already been initialized!");
31  }
32 }
Node::Object
static Node * Object()
Definition: Node.h:86
Node::simParameters
SimParameters * simParameters
Definition: Node.h:181
LjPmeMgr::initialize
void initialize(const SimParameters *simParams, const int nAtoms)
Definition: LjPmeMgr.C:14
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:147
The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.14