NAMD
Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | List of all members
TholeElem Class Reference

#include <ComputeThole.h>

Public Types

enum  { size = 4 }
 
enum  {
  tholeEnergyIndex, tholeEnergyIndex_f, tholeEnergyIndex_ti_1, tholeEnergyIndex_ti_2,
  TENSOR =(virialIndex), reductionDataSize
}
 
enum  { reductionChecksumLabel = REDUCTION_THOLE_CHECKSUM }
 

Public Member Functions

int hash () const
 
 TholeElem ()
 
 TholeElem (AtomID atom0, const TupleSignature *sig, const TholeValue *v)
 
 TholeElem (const Thole *a, const TholeValue *v)
 
 TholeElem (AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3)
 
 ~TholeElem ()
 
int operator== (const TholeElem &a) const
 
int operator< (const TholeElem &a) const
 

Static Public Member Functions

static void computeForce (TholeElem *, int, BigReal *, BigReal *)
 
static void getMoleculePointers (Molecule *, int *, int32 ***, Thole **)
 
static void getParameterPointers (Parameters *, const TholeValue **)
 
static void getTupleInfo (AtomSignature *sig, int *count, TupleSignature **t)
 
static void submitReductionData (BigReal *, SubmitReduction *)
 

Public Attributes

AtomID atomID [size]
 
int localIndex [size]
 
TuplePatchElemp [size]
 
Real scale
 
const TholeValuevalue
 

Static Public Attributes

static int pressureProfileSlabs = 0
 
static int pressureProfileAtomTypes = 1
 
static BigReal pressureProfileThickness = 0
 
static BigReal pressureProfileMin = 0
 

Detailed Description

Definition at line 18 of file ComputeThole.h.

Member Enumeration Documentation

anonymous enum
Enumerator
size 

Definition at line 21 of file ComputeThole.h.

21 { size = 4 };
anonymous enum
Enumerator
tholeEnergyIndex 
tholeEnergyIndex_f 
tholeEnergyIndex_ti_1 
tholeEnergyIndex_ti_2 
TENSOR 
reductionDataSize 

Definition at line 49 of file ComputeThole.h.

anonymous enum
Enumerator
reductionChecksumLabel 

Definition at line 51 of file ComputeThole.h.

Constructor & Destructor Documentation

TholeElem::TholeElem ( )
inline

Definition at line 12 of file ComputeThole.inl.

12 { ; }
TholeElem::TholeElem ( AtomID  atom0,
const TupleSignature sig,
const TholeValue v 
)
inline

Definition at line 14 of file ComputeThole.inl.

References NAMD_die().

14  {
15  NAMD_die("Can't use Thole with memory optimized version of NAMD.");
16  // atomID[0] = atom0;
17  // atomID[1] = atom0 + sig->offset[0];
18  // atomID[2] = atom0 + sig->offset[1];
19  // atomID[3] = atom0 + sig->offset[2];
20  // value = &v[sig->tupleParamType];
21 }
void NAMD_die(const char *err_msg)
Definition: common.C:85
TholeElem::TholeElem ( const Thole a,
const TholeValue v 
)
inline

Definition at line 23 of file ComputeThole.inl.

References thole::atom1, thole::atom2, thole::atom3, thole::atom4, atomID, and value.

24  {
25  atomID[0] = a->atom1;
26  atomID[1] = a->atom2;
27  atomID[2] = a->atom3;
28  atomID[3] = a->atom4;
29  value = a; // expect v to be NULL
30  }
int32 atom3
Definition: structures.h:145
int32 atom4
Definition: structures.h:146
int32 atom2
Definition: structures.h:144
int32 atom1
Definition: structures.h:143
const TholeValue * value
Definition: ComputeThole.h:43
AtomID atomID[size]
Definition: ComputeThole.h:22
TholeElem::TholeElem ( AtomID  atom0,
AtomID  atom1,
AtomID  atom2,
AtomID  atom3 
)
inline

Definition at line 32 of file ComputeThole.inl.

References atomID.

34  {
35  // atoms arranged: HEAVY DRUDE HEAVY DRUDE
36  if (atom0 > atom2) { // Swap heavy atoms so lowest is first!
37  AtomID tmp = atom2; atom2 = atom0; atom0 = tmp;
38  tmp = atom1; atom1 = atom3; atom3 = tmp;
39  }
40  atomID[0] = atom0;
41  atomID[1] = atom1;
42  atomID[2] = atom2;
43  atomID[3] = atom3;
44  }
int AtomID
Definition: NamdTypes.h:29
AtomID atomID[size]
Definition: ComputeThole.h:22
TholeElem::~TholeElem ( )
inline

Definition at line 58 of file ComputeThole.h.

58 {};

Member Function Documentation

void TholeElem::computeForce ( TholeElem tuples,
int  ntuple,
BigReal reduction,
BigReal pressureProfileData 
)
static

Definition at line 41 of file ComputeThole.C.

References thole::aa, SimParameters::alchDecouple, SimParameters::alchOn, atomID, DebugM, Lattice::delta(), TuplePatchElem::f, Patch::flags, Molecule::get_fep_type(), SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), SimParameters::getElecLambda(), Patch::lattice, localIndex, Node::molecule, Node::Object(), order, p, TuplePatchElem::p, CompAtom::partition, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, thole::qq, Vector::rlength(), scale, Node::simParameters, simParams, size, Flags::step, tholeEnergyIndex, tholeEnergyIndex_f, tholeEnergyIndex_ti_1, tholeEnergyIndex_ti_2, value, TuplePatchElem::x, and Vector::z.

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 }
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
int AtomID
Definition: NamdTypes.h:29
Lattice & lattice
Definition: Patch.h:126
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:109
#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
void pp_clamp(int &n, int nslabs)
#define order
Definition: PmeRealSpace.C:235
BigReal rlength(void)
Definition: Vector.h:177
BigReal getCurrentLambda2(const int)
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
static BigReal pressureProfileMin
Definition: ComputeThole.h:40
Real qq
Definition: structures.h:148
const TholeValue * value
Definition: ComputeThole.h:43
unsigned char get_fep_type(int anum) const
Definition: Molecule.h:1349
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:114
int step
Definition: PatchTypes.h:16
void TholeElem::getMoleculePointers ( Molecule mol,
int *  count,
int32 ***  byatom,
Thole **  structarray 
)
static

Definition at line 26 of file ComputeThole.C.

References NAMD_die(), and Molecule::numTholes.

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 }
int numTholes
Number of Thole terms.
Definition: Molecule.h:583
void NAMD_die(const char *err_msg)
Definition: common.C:85
void TholeElem::getParameterPointers ( Parameters p,
const TholeValue **  v 
)
static

Definition at line 37 of file ComputeThole.C.

37  {
38  *v = NULL; // parameters are stored in the structure
39 }
static void TholeElem::getTupleInfo ( AtomSignature sig,
int *  count,
TupleSignature **  t 
)
inlinestatic

Definition at line 30 of file ComputeThole.h.

References NAMD_die().

30  {
31  NAMD_die("Can't use Thole with memory optimized version of NAMD.");
32  // *count = sig->tholeCnt;
33  // *t = sig->tholeSigs;
34  }
void NAMD_die(const char *err_msg)
Definition: common.C:85
int TholeElem::hash ( void  ) const
inline

Definition at line 45 of file ComputeThole.h.

References atomID.

45  {
46  return 0x7FFFFFFF &((atomID[0]<<24) + (atomID[1]<<16) + (atomID[2]<<8) + atomID[3]);
47  }
AtomID atomID[size]
Definition: ComputeThole.h:22
int TholeElem::operator< ( const TholeElem a) const
inline

Definition at line 52 of file ComputeThole.inl.

References atomID.

53  {
54  return (atomID[0] < a.atomID[0] ||
55  (atomID[0] == a.atomID[0] &&
56  (atomID[1] < a.atomID[1] ||
57  (atomID[1] == a.atomID[1] &&
58  (atomID[2] < a.atomID[2] ||
59  (atomID[2] == a.atomID[2] &&
60  atomID[3] < a.atomID[3]
61  ))))));
62  }
AtomID atomID[size]
Definition: ComputeThole.h:22
int TholeElem::operator== ( const TholeElem a) const
inline

Definition at line 46 of file ComputeThole.inl.

References atomID.

47  {
48  return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
49  a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3]);
50  }
AtomID atomID[size]
Definition: ComputeThole.h:22
void TholeElem::submitReductionData ( BigReal data,
SubmitReduction reduction 
)
static

Member Data Documentation

AtomID TholeElem::atomID[size]

Definition at line 22 of file ComputeThole.h.

Referenced by computeForce(), hash(), operator<(), operator==(), and TholeElem().

int TholeElem::localIndex[size]

Definition at line 23 of file ComputeThole.h.

Referenced by computeForce().

TuplePatchElem* TholeElem::p[size]

Definition at line 24 of file ComputeThole.h.

Referenced by computeForce().

int TholeElem::pressureProfileAtomTypes = 1
static

Definition at line 38 of file ComputeThole.h.

Referenced by computeForce().

BigReal TholeElem::pressureProfileMin = 0
static

Definition at line 40 of file ComputeThole.h.

Referenced by computeForce().

int TholeElem::pressureProfileSlabs = 0
static

Definition at line 37 of file ComputeThole.h.

Referenced by computeForce().

BigReal TholeElem::pressureProfileThickness = 0
static

Definition at line 39 of file ComputeThole.h.

Referenced by computeForce().

Real TholeElem::scale

Definition at line 25 of file ComputeThole.h.

Referenced by computeForce().

const TholeValue* TholeElem::value

Definition at line 43 of file ComputeThole.h.

Referenced by computeForce(), and TholeElem().


The documentation for this class was generated from the following files: