NAMD
ComputeGromacsPair.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeGromacsPair.h"
9 #include "Molecule.h"
10 #include "Parameters.h"
11 #include "LJTable.h"
12 #include "Node.h"
13 #include "ReductionMgr.h"
14 #include "Lattice.h"
15 #include "PressureProfile.h"
16 #include "Debug.h"
17 #include <iostream>
18 
19 // static initialization
24 
25 
27  (Molecule* mol, int* count, int32*** byatom, GromacsPair** structarray)
28 {
29 #ifdef MEM_OPT_VERSION
30  NAMD_die("Should not be called in GromacsPairElem::getMoleculePointers in memory optimized version!");
31 #else
32  *count = mol->numLJPair;
33  *byatom = mol->gromacsPairByAtom;
34  *structarray = mol->gromacsPair;
35 #endif
36 }
37 
38 //void GromacsPairElem::getParameterPointers(Parameters *p, const GromacsPairValue **v) {
40  // JLai
41  *v=p->gromacsPair_array;
42 }
43 
44 void GromacsPairElem::computeForce(GromacsPairElem *tuples, int ntuple, BigReal *reduction,
45  BigReal *pressureProfileData)
46 {
47  const Lattice & lattice = tuples[0].p[0]->p->lattice;
48 
49  for ( int ituple=0; ituple<ntuple; ++ituple ) {
50  const GromacsPairElem &tup = tuples[ituple];
51  enum { size = 2 };
52  const AtomID (&atomID)[size](tup.atomID);
53  const int (&localIndex)[size](tup.localIndex);
54  TuplePatchElem * const(&p)[size](tup.p);
55  const Real (&scale)(tup.scale);
56  const GromacsPairValue * const(&value)(tup.value);
57 
58  DebugM(1, "::computeforce() localIndex = " << localIndex [0] << " "
59  << localIndex[1] << std::endl);
60 
61  // compute vectors between atoms and their distances
62  const Vector r12 = lattice.delta(p[0]->x[localIndex[0]].position,
63  p[1]->x[localIndex[1]].position);
64 
65  if ( p[0]->patchID == p[1]->patchID && localIndex[0] == localIndex[1] ) {
66  continue;
67  }
69  BigReal cutoff = simParams->cutoff;
70  BigReal r12_len = r12.length2();
71  if ( r12_len == 0 ) continue;
72  //if ( r12_len > cutoff) {
73  //continue;
74  //}
75  BigReal ri2 = 1.0/r12_len;
76  BigReal ri = sqrt(ri2);
77  BigReal ri6 = ri2*ri2*ri2;
78  BigReal ri7 = ri*ri6;
79  BigReal ri8 = ri2*ri2*ri2*ri2;
80  BigReal ri12 = ri6*ri6;
81  BigReal ri13 = ri12*ri;
82  BigReal ri14 = ri12*ri2;
83 
84  BigReal energy = 0;
85  BigReal diff = 0;
86  BigReal pairC12 = value->pairC12;
87  BigReal pairC6 = value->pairC6;
88 
89  // Add the energy for the 12-6 LJ interaction to the total energy
90  energy = (pairC12*ri12) - (pairC6*ri6);
91  // This is a dirty hack; currently the code is looping over N^2 instead of N(N-1)/2
92  // This is happening because the LJ list is 2x long as it has to be
93  //energy *= 0.5;
94 
95  // Determine the magnitude of the force
96  diff = ((12*pairC12*ri14) - (6*pairC6*ri8));
97  // This is a dirty hack; currently the code is looping over N^2 instead of N(N-1)/2
98  // This is happening because the LJ list is 2x long as it has to be
99  //diff *= 0.5;
100  //std::cout << "Force: " << diff << " " << pairC12 << " " << pairC6 << " " << r12.length() << "\n";
101 
102  //Scale the force vector accordingly
103  const Force f12 = diff * r12;
104  //std::cout << "Atoms: " << localIndex[0] << " " << localIndex[1] << " " << f12.length() << " " << f12 << "\n";
105  //std::cout << "Force2: " << f12 << "\n";
106 
107  // Now add the forces to each force vector
108  p[0]->f[localIndex[0]] += f12;
109  p[1]->f[localIndex[1]] -= f12;
110 
111  DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
112  reduction[gromacsPairEnergyIndex] += energy;
113  reduction[virialIndex_XX] += f12.x * r12.x;
114  reduction[virialIndex_XY] += f12.x * r12.y;
115  reduction[virialIndex_XZ] += f12.x * r12.z;
116  reduction[virialIndex_YX] += f12.y * r12.x;
117  reduction[virialIndex_YY] += f12.y * r12.y;
118  reduction[virialIndex_YZ] += f12.y * r12.z;
119  reduction[virialIndex_ZX] += f12.z * r12.x;
120  reduction[virialIndex_ZY] += f12.z * r12.y;
121  reduction[virialIndex_ZZ] += f12.z * r12.z;
122 
123  if (pressureProfileData) {
124  BigReal z1 = p[0]->x[localIndex[0]].position.z;
125  BigReal z2 = p[1]->x[localIndex[1]].position.z;
126  int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
127  int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
130  int p1 = p[0]->x[localIndex[0]].partition;
131  int p2 = p[1]->x[localIndex[1]].partition;
132  int pn = pressureProfileAtomTypes;
134  n1, n2,
135  p1, p2, pn,
136  f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
137  pressureProfileData);
138  }
139 
140  }
141 }
142 
144 {
145  // JLai
147  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
148  return;
149 }
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
GromacsPairValue * gromacsPair_array
Definition: Parameters.h:247
static void getParameterPointers(Parameters *, const GromacsPairValue **)
TuplePatchElem * p[size]
static int pressureProfileSlabs
const GromacsPairValue * value
short int32
Definition: dumpdcd.c:24
int AtomID
Definition: NamdTypes.h:29
static void computeForce(GromacsPairElem *, int, BigReal *, BigReal *)
Lattice & lattice
Definition: Patch.h:126
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:32
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
float Real
Definition: common.h:107
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
static void submitReductionData(BigReal *, SubmitReduction *)
void pp_clamp(int &n, int nslabs)
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:83
AtomID atomID[size]
BigReal length2(void) const
Definition: Vector.h:173
#define simParams
Definition: Output.C:127
BigReal y
Definition: Vector.h:66
static int pressureProfileAtomTypes
static void getMoleculePointers(Molecule *, int *, int32 ***, GromacsPair **)
static BigReal pressureProfileThickness
gridSize x
static BigReal pressureProfileMin
double BigReal
Definition: common.h:112
int numLJPair
Definition: Molecule.h:661