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(), SimParameters::LJPMESerial, NAMD_die(), Node::Object(), and Node::simParameters.

Referenced by ComputeLjPmeSerialMgr::recvCoord().

38  {
39  if(!initialized) {
40  NAMD_die("LjPmeCompute has been called without initialization!");
41  }
42  SimParameters *simParm = Node::Object()->simParameters;
43  // Initialize the energy and virial values to zero
44  energySlow = energyNonbond = 0.0;
45  for(int i = 0; i < 3; ++i) {
46  for(int j = 0; j < 3; ++j) {
47  virialNonbonded[i][j] = 0.0;
48  virialSlow[i][j] = 0.0;
49  }
50  }
51 
52  // Calculate the short range and long range Nonbonded LJ potential
53  if(doSlow) {
54  if (doEnergy) {
55  if (simParm->LJPMESerial) {
56  // Calculate the LJ nonBonded term
57  this->computeNonbonded<true, true>(alphaLJ, cutoff,
58  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
59  }
60  } else {
61  energySlow = 0.0;
62  if (simParm->LJPMESerial) {
63  this->computeNonbonded<true, false>(alphaLJ, cutoff,
64  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
65  }
66  }
67  // Compute reciprocal force, virial, energy, and self energy
68  ljPmeMgr.computeLongRange(ljPmeCoord, lattice, alphaLJ,
69  ljForceSlow, energySlow, virialSlow, doEnergy);
70  } else {
71  if (doEnergy) {
72  this->computeNonbonded<false, true>(alphaLJ, cutoff,
73  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
74  } else {
75  this->computeNonbonded<false, false>(alphaLJ, cutoff,
76  lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
77  }
78  }
79 }
Node::Object
static Node * Object()
Definition: Node.h:86
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
Node::simParameters
SimParameters * simParameters
Definition: Node.h:181
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 82 of file LjPmeCompute.C.

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

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