NAMD
ComputeBonds.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeBonds.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 
18 // static initialization
23 
25  (Molecule* mol, int* count, int32*** byatom, Bond** structarray)
26 {
27 #ifdef MEM_OPT_VERSION
28  NAMD_die("Should not be called in BondElem::getMoleculePointers in memory optimized version!");
29 #else
30  *count = mol->numBonds;
31  *byatom = mol->bondsByAtom;
32  *structarray = mol->bonds;
33 #endif
34 }
35 
37  *v = p->bond_array;
38 }
39 
40 void BondElem::computeForce(BondElem *tuples, int ntuple, BigReal *reduction,
41  BigReal *pressureProfileData)
42 {
44  const Lattice & lattice = tuples[0].p[0]->p->lattice;
45 
46  //fepb BKR
47  const int step = tuples[0].p[0]->p->flags.step;
48  const BigReal alchLambda = simParams->getCurrentLambda(step);
49  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
50  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
51  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
52  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
53  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
54  Molecule *const mol = Node::Object()->molecule;
55  //fepe
56 
57  for ( int ituple=0; ituple<ntuple; ++ituple ) {
58  const BondElem &tup = tuples[ituple];
59  enum { size = 2 };
60  const AtomID (&atomID)[size](tup.atomID);
61  const int (&localIndex)[size](tup.localIndex);
62  TuplePatchElem * const(&p)[size](tup.p);
63  const Real (&scale)(tup.scale);
64  const BondValue * const(&value)(tup.value);
65 
66  DebugM(1, "::computeForce() localIndex = " << localIndex[0] << " "
67  << localIndex[1] << std::endl);
68 
69  // skip Lonepair bonds (other k=0. bonds have been filtered out)
70  if (0. == value->k) continue;
71 
72  BigReal scal; // force scaling
73  BigReal energy; // energy from the bond
74 
75  //DebugM(3, "::computeForce() -- starting with bond type " << bondType << std::endl);
76 
77  // get the bond information
78  Real k = value->k * scale;
79  Real x0 = value->x0;
80  Real x1 = value->x1;
81 
82  // compute vectors between atoms and their distances
83  const Vector r12 = lattice.delta(p[0]->x[localIndex[0]].position,
84  p[1]->x[localIndex[1]].position);
85 
86  if (0. == x0) {
87  BigReal r2 = r12.length2();
88  scal = -2.0*k; // bond interaction for equilibrium length 0
89  energy = k*r2;
90  // TODO: Instead flag by mass for Drude particles, otherwise mixing and
91  // matching Drude and alchemical "twin" particles will likely not work as
92  // expected.
93  if( simParams->drudeOn && !simParams->drudeHardWallOn ) { // for Drude bonds
94  BigReal drudeBondLen = simParams->drudeBondLen;
95  BigReal drudeBondConst = simParams->drudeBondConst;
96  if ( drudeBondConst > 0 && r2 > drudeBondLen*drudeBondLen ) {
97  // add a quartic restraining potential to keep Drude bond short
98  BigReal r = sqrt(r2);
99  BigReal diff = r - drudeBondLen;
100  BigReal diff2 = diff*diff;
101 
102  scal += -4*drudeBondConst * diff2 * diff / r;
103  energy += drudeBondConst * diff2 * diff2;
104  }
105  }
106  }
107  else {
108  BigReal r = r12.length(); // Distance between atoms
109  BigReal diff = r - x0; // Compare it to the rest bond
110 
111  if (x1) {
112  // in this case, the bond represents a harmonic wall potential
113  // where x0 is the lower wall and x1 is the upper
114  diff = (r > x1 ? r - x1 : (r > x0 ? 0 : diff));
115  }
116 
117  // Add the energy from this bond to the total energy
118  energy = k*diff*diff;
119 
120  // Determine the magnitude of the force
121  diff *= -2.0*k;
122 
123  // Scale the force vector accordingly
124  scal = (diff/r);
125  }
126 
127  //fepb - BKR scaling of alchemical bonded terms
128  // NB: TI derivative is the _unscaled_ energy.
129  if ( simParams->alchOn && !simParams->singleTopology) {
130  switch ( mol->get_fep_bonded_type(atomID, 2) ) {
131  case 1:
132  reduction[bondEnergyIndex_ti_1] += energy;
133  reduction[bondEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
134  energy *= bond_lambda_1;
135  scal *= bond_lambda_1;
136  break;
137  case 2:
138  reduction[bondEnergyIndex_ti_2] += energy;
139  reduction[bondEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
140  energy *= bond_lambda_2;
141  scal *= bond_lambda_2;
142  break;
143  }
144  }
145  //fepe
146  const Force f12 = scal * r12;
147 
148 
149  // Now add the forces to each force vector
150  p[0]->f[localIndex[0]] += f12;
151  p[1]->f[localIndex[1]] -= f12;
152 
153  DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
154  reduction[bondEnergyIndex] += energy;
155  reduction[virialIndex_XX] += f12.x * r12.x;
156  reduction[virialIndex_XY] += f12.x * r12.y;
157  reduction[virialIndex_XZ] += f12.x * r12.z;
158  reduction[virialIndex_YX] += f12.y * r12.x;
159  reduction[virialIndex_YY] += f12.y * r12.y;
160  reduction[virialIndex_YZ] += f12.y * r12.z;
161  reduction[virialIndex_ZX] += f12.z * r12.x;
162  reduction[virialIndex_ZY] += f12.z * r12.y;
163  reduction[virialIndex_ZZ] += f12.z * r12.z;
164 
165  if (pressureProfileData) {
166  BigReal z1 = p[0]->x[localIndex[0]].position.z;
167  BigReal z2 = p[1]->x[localIndex[1]].position.z;
168  int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
169  int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
172  int p1 = p[0]->x[localIndex[0]].partition;
173  int p2 = p[1]->x[localIndex[1]].partition;
174  int pn = pressureProfileAtomTypes;
176  n1, n2,
177  p1, p2, pn,
178  f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
179  pressureProfileData);
180  }
181 
182  }
183 }
184 
186 {
187  reduction->item(REDUCTION_BOND_ENERGY) += data[bondEnergyIndex];
188  reduction->item(REDUCTION_BONDED_ENERGY_F) += data[bondEnergyIndex_f];
191  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
192 }
193 
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
int numBonds
Definition: Molecule.h:559
static void getParameterPointers(Parameters *, const BondValue **)
Definition: ComputeBonds.C:36
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)
static void getMoleculePointers(Molecule *, int *, int32 ***, Bond **)
Definition: ComputeBonds.C:25
Definition: Vector.h:64
static void computeForce(BondElem *, int, BigReal *, BigReal *)
Definition: ComputeBonds.C:40
SimParameters * simParameters
Definition: Node.h:178
float Real
Definition: common.h:109
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 pressureProfileAtomTypes
Definition: ComputeBonds.h:39
int get_fep_bonded_type(const int *atomID, unsigned int order) const
Definition: Molecule.h:1385
BigReal length(void) const
Definition: Vector.h:169
Real x1
Definition: Parameters.h:86
Real scale
Definition: ComputeBonds.h:27
AtomID atomID[size]
Definition: ComputeBonds.h:24
static BigReal pressureProfileMin
Definition: ComputeBonds.h:41
static BigReal pressureProfileThickness
Definition: ComputeBonds.h:40
Flags flags
Definition: Patch.h:127
void pp_clamp(int &n, int nslabs)
int localIndex[size]
Definition: ComputeBonds.h:25
BigReal getBondLambda(const BigReal)
BigReal getCurrentLambda2(const int)
BigReal drudeBondLen
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
static void submitReductionData(BigReal *, SubmitReduction *)
Definition: ComputeBonds.C:185
BigReal drudeBondConst
static int pressureProfileSlabs
Definition: ComputeBonds.h:38
BigReal length2(void) const
Definition: Vector.h:173
#define simParams
Definition: Output.C:127
Real x0
Definition: Parameters.h:85
BondValue * bond_array
Definition: Parameters.h:241
BigReal y
Definition: Vector.h:66
BigReal getCurrentLambda(const int)
const BondValue * value
Definition: ComputeBonds.h:46
TuplePatchElem * p[size]
Definition: ComputeBonds.h:26
gridSize x
Molecule * molecule
Definition: Node.h:176
Real k
Definition: Parameters.h:84
double BigReal
Definition: common.h:114
int step
Definition: PatchTypes.h:16