NAMD
ComputeThole.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeThole.h"
9 #include "Molecule.h"
10 #include "Parameters.h"
11 #include "Node.h"
12 #include "ReductionMgr.h"
13 #include "Lattice.h"
14 #include "PressureProfile.h"
15 #include "Debug.h"
16 
17 #define CALCULATE_THOLE_CORRECTION
18 
19 // static initialization
24 
26  (Molecule* mol, int* count, int32*** byatom, Thole** structarray)
27 {
28 #ifdef MEM_OPT_VERSION
29  NAMD_die("Should not be called in TholeElem::getMoleculePointers in memory optimized version!");
30 #else
31  *count = mol->numTholes;
32  *byatom = mol->tholesByAtom;
33  *structarray = mol->tholes;
34 #endif
35 }
36 
38  *v = NULL; // parameters are stored in the structure
39 }
40 
41 void TholeElem::computeForce(TholeElem *tuples, int ntuple, BigReal *reduction,
42  BigReal *pressureProfileData)
43 {
44  const Lattice & lattice = tuples[0].p[0]->p->lattice;
45 
46  //fepb BKR
48  const int step = tuples[0].p[0]->p->flags.step;
49  const BigReal alchLambda = simParams->getCurrentLambda(step);
50  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
51  const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
52  const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
53  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda2);
54  const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda2);
55  Molecule *const mol = Node::Object()->molecule;
56  //fepe
57 
58 
59  for ( int ituple=0; ituple<ntuple; ++ituple ) {
60  const TholeElem &tup = tuples[ituple];
61  enum { size = 4 };
62  const AtomID (&atomID)[size](tup.atomID);
63  const int (&localIndex)[size](tup.localIndex);
64  TuplePatchElem * const(&p)[size](tup.p);
65  const Real (&scale)(tup.scale);
66  const TholeValue * const(&value)(tup.value);
67 
68  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
69  << localIndex[1] << " " << localIndex[2] << " "
70  << localIndex[3] << std::endl);
71 
72 #ifdef CALCULATE_THOLE_CORRECTION
73  const BigReal aa = value->aa;
74  const BigReal qq = value->qq;
75 
76  // Calculate the vectors between atoms
77  const Position & rai = p[0]->x[localIndex[0]].position; // atom i
78  const Position & rdi = p[1]->x[localIndex[1]].position; // atom i's Drude
79  const Position & raj = p[2]->x[localIndex[2]].position; // atom j
80  const Position & rdj = p[3]->x[localIndex[3]].position; // atom j's Drude
81 
82  // r_ij = r_i - r_j
83  Vector raa = lattice.delta(rai,raj); // shortest vector image: rai - raj
84  Vector rad = lattice.delta(rai,rdj); // shortest vector image: rai - rdj
85  Vector rda = lattice.delta(rdi,raj); // shortest vector image: rdi - raj
86  Vector rdd = lattice.delta(rdi,rdj); // shortest vector image: rdi - rdj
87 
88  // 1/r, r = |r_ij|
89  BigReal raa_invlen = raa.rlength(); // reciprocal of length
90  BigReal rad_invlen = rad.rlength();
91  BigReal rda_invlen = rda.rlength();
92  BigReal rdd_invlen = rdd.rlength();
93 
94  // ar
95  BigReal auaa = aa / raa_invlen;
96  BigReal auad = aa / rad_invlen;
97  BigReal auda = aa / rda_invlen;
98  BigReal audd = aa / rdd_invlen;
99 
100  // exp(-ar)
101  BigReal expauaa = exp(-auaa);
102  BigReal expauad = exp(-auad);
103  BigReal expauda = exp(-auda);
104  BigReal expaudd = exp(-audd);
105 
106  // (1 + ar/2)
107  BigReal polyauaa = 1 + 0.5*auaa;
108  BigReal polyauad = 1 + 0.5*auad;
109  BigReal polyauda = 1 + 0.5*auda;
110  BigReal polyaudd = 1 + 0.5*audd;
111 
112  // U(r) = qq/r (1 - (1 + ar/2) exp(-ar))
113  BigReal ethole = 0;
114  ethole += qq * raa_invlen * (1 - polyauaa * expauaa);
115  ethole += -qq * rad_invlen * (1 - polyauad * expauad);
116  ethole += -qq * rda_invlen * (1 - polyauda * expauda);
117  ethole += qq * rdd_invlen * (1 - polyaudd * expaudd);
118 
119  polyauaa = 1 + auaa*polyauaa;
120  polyauad = 1 + auad*polyauad;
121  polyauda = 1 + auda*polyauda;
122  polyaudd = 1 + audd*polyaudd;
123 
124  BigReal raa_invlen3 = raa_invlen * raa_invlen * raa_invlen;
125  BigReal rad_invlen3 = rad_invlen * rad_invlen * rad_invlen;
126  BigReal rda_invlen3 = rda_invlen * rda_invlen * rda_invlen;
127  BigReal rdd_invlen3 = rdd_invlen * rdd_invlen * rdd_invlen;
128 
129  // df = (1/r) (dU/dr)
130  BigReal dfaa = qq * raa_invlen3 * (polyauaa*expauaa - 1);
131  BigReal dfad = -qq * rad_invlen3 * (polyauad*expauad - 1);
132  BigReal dfda = -qq * rda_invlen3 * (polyauda*expauda - 1);
133  BigReal dfdd = qq * rdd_invlen3 * (polyaudd*expaudd - 1);
134 
135  Vector faa = -dfaa * raa;
136  Vector fad = -dfad * rad;
137  Vector fda = -dfda * rda;
138  Vector fdd = -dfdd * rdd;
139 
140  //fepb - BKR scaling of alchemical drude terms
141  // NB: TI derivative is the _unscaled_ energy.
142  if ( simParams->alchOn ) {
143  // THIS ASSUMES THAT AN ATOM AND ITS DRUDE PARTICLE ARE ALWAYS IN THE SAME
144  // PARTITION!
145  int typeSum = 0;
146  typeSum += (mol->get_fep_type(atomID[0]) == 2 ? -1 :\
147  mol->get_fep_type(atomID[0]));
148  typeSum += (mol->get_fep_type(atomID[2]) == 2 ? -1 :\
149  mol->get_fep_type(atomID[2]));
150  int order = (!simParams->alchDecouple ? 3 : 2);
151 
152  if ( 0 < typeSum && typeSum < order ) {
153  reduction[tholeEnergyIndex_ti_1] += ethole;
154  reduction[tholeEnergyIndex_f] += elec_lambda_12*ethole;
155  ethole *= elec_lambda_1;
156  faa *= elec_lambda_1;
157  fad *= elec_lambda_1;
158  fda *= elec_lambda_1;
159  fdd *= elec_lambda_1;
160  } else if ( 0 > typeSum && typeSum > -order ) {
161  reduction[tholeEnergyIndex_ti_2] += ethole;
162  reduction[tholeEnergyIndex_f] += elec_lambda_22*ethole;
163  ethole *= elec_lambda_2;
164  faa *= elec_lambda_2;
165  fad *= elec_lambda_2;
166  fda *= elec_lambda_2;
167  fdd *= elec_lambda_2;
168  }
169  }
170  //fepe
171 
172  p[0]->f[localIndex[0]] += faa + fad;
173  p[1]->f[localIndex[1]] += fda + fdd;
174  p[2]->f[localIndex[2]] -= faa + fda;
175  p[3]->f[localIndex[3]] -= fad + fdd;
176 
177  DebugM(3, "::computeForce() -- ending with delta energy " << ethole
178  << std::endl);
179  reduction[tholeEnergyIndex] += ethole;
180 
181  reduction[virialIndex_XX] += faa.x * raa.x + fad.x * rad.x
182  + fda.x * rda.x + fdd.x * rdd.x;
183  reduction[virialIndex_XY] += faa.x * raa.y + fad.x * rad.y
184  + fda.x * rda.y + fdd.x * rdd.y;
185  reduction[virialIndex_XZ] += faa.x * raa.z + fad.x * rad.z
186  + fda.x * rda.z + fdd.x * rdd.z;
187  reduction[virialIndex_YX] += faa.y * raa.x + fad.y * rad.x
188  + fda.y * rda.x + fdd.y * rdd.x;
189  reduction[virialIndex_YY] += faa.y * raa.y + fad.y * rad.y
190  + fda.y * rda.y + fdd.y * rdd.y;
191  reduction[virialIndex_YZ] += faa.y * raa.z + fad.y * rad.z
192  + fda.y * rda.z + fdd.y * rdd.z;
193  reduction[virialIndex_ZX] += faa.z * raa.x + fad.z * rad.x
194  + fda.z * rda.x + fdd.z * rdd.x;
195  reduction[virialIndex_ZY] += faa.z * raa.y + fad.z * rad.y
196  + fda.z * rda.y + fdd.z * rdd.y;
197  reduction[virialIndex_ZZ] += faa.z * raa.z + fad.z * rad.z
198  + fda.z * rda.z + fdd.z * rdd.z;
199 
200  if (pressureProfileData) {
201  BigReal zai = p[0]->x[localIndex[0]].position.z;
202  BigReal zdi = p[1]->x[localIndex[1]].position.z;
203  BigReal zaj = p[2]->x[localIndex[2]].position.z;
204  BigReal zdj = p[3]->x[localIndex[3]].position.z;
205  int nai = (int)floor((zai-pressureProfileMin)/pressureProfileThickness);
206  int ndi = (int)floor((zdi-pressureProfileMin)/pressureProfileThickness);
207  int naj = (int)floor((zaj-pressureProfileMin)/pressureProfileThickness);
208  int ndj = (int)floor((zdj-pressureProfileMin)/pressureProfileThickness);
213  int pai = p[0]->x[localIndex[0]].partition;
214  int pdi = p[1]->x[localIndex[1]].partition;
215  int paj = p[2]->x[localIndex[2]].partition;
216  int pdj = p[3]->x[localIndex[3]].partition;
217  int pn = pressureProfileAtomTypes;
219  pai, paj, pn, faa.x * raa.x, faa.y * raa.y, faa.z * raa.z,
220  pressureProfileData);
222  pai, pdj, pn, fad.x * rad.x, fad.y * rad.y, fad.z * rad.z,
223  pressureProfileData);
225  pdi, paj, pn, fda.x * rda.x, fda.y * rda.y, fda.z * rda.z,
226  pressureProfileData);
228  pdi, pdj, pn, fdd.x * rdd.x, fdd.y * rdd.y, fdd.z * rdd.z,
229  pressureProfileData);
230  }
231 #endif
232 
233  }
234 }
235 
236 
237 // The energy from the screened Coulomb correction of Thole is
238 // accumulated into the electrostatic potential energy.
240 {
241  reduction->item(REDUCTION_ELECT_ENERGY) += data[tholeEnergyIndex];
242  reduction->item(REDUCTION_ELECT_ENERGY_F) += data[tholeEnergyIndex_f];
245  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
246 }
247 
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
short int32
Definition: dumpdcd.c:24
int AtomID
Definition: NamdTypes.h:29
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 int pressureProfileSlabs
Definition: ComputeThole.h:37
Flags flags
Definition: Patch.h:127
static void getMoleculePointers(Molecule *, int *, int32 ***, Thole **)
Definition: ComputeThole.C:26
static void computeForce(TholeElem *, int, BigReal *, BigReal *)
Definition: ComputeThole.C:41
void pp_clamp(int &n, int nslabs)
#define order
Definition: PmeRealSpace.C:235
BigReal rlength(void)
Definition: Vector.h:177
static void getParameterPointers(Parameters *, const TholeValue **)
Definition: ComputeThole.C:37
BigReal getCurrentLambda2(const int)
static void submitReductionData(BigReal *, SubmitReduction *)
Definition: ComputeThole.C:239
int numTholes
Number of Thole terms.
Definition: Molecule.h:582
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
static BigReal pressureProfileMin
Definition: ComputeThole.h:40
void NAMD_die(const char *err_msg)
Definition: common.C:83
Real qq
Definition: structures.h:148
const TholeValue * value
Definition: ComputeThole.h:43
unsigned char get_fep_type(int anum) const
Definition: Molecule.h:1345
Real aa
Definition: structures.h:147
#define simParams
Definition: Output.C:127
Real scale
Definition: ComputeThole.h:25
BigReal getCurrentLambda(const int)
static int pressureProfileAtomTypes
Definition: ComputeThole.h:38
int localIndex[size]
Definition: ComputeThole.h:23
TuplePatchElem * p[size]
Definition: ComputeThole.h:24
static BigReal pressureProfileThickness
Definition: ComputeThole.h:39
Molecule * molecule
Definition: Node.h:176
AtomID atomID[size]
Definition: ComputeThole.h:22
BigReal getElecLambda(const BigReal)
double BigReal
Definition: common.h:112
int step
Definition: PatchTypes.h:16