NAMD
ComputeDihedrals.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeDihedrals.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, Dihedral** structarray)
26 {
27 #ifdef MEM_OPT_VERSION
28  NAMD_die("Should not be called in DihedralElem::getMoleculePointers in memory optimized version!");
29 #else
30  *count = mol->numDihedrals;
31  *byatom = mol->dihedralsByAtom;
32  *structarray = mol->dihedrals;
33 #endif
34 }
35 
37  *v = p->dihedral_array;
38 }
39 
40 void DihedralElem::computeForce(DihedralElem *tuples, int ntuple, BigReal *reduction,
41  BigReal *pressureProfileData)
42 {
43  const Lattice & lattice = tuples[0].p[0]->p->lattice;
44 
45  //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 DihedralElem &tup = tuples[ituple];
59  enum { size = 4 };
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 DihedralValue * const(&value)(tup.value);
65 
66  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
67  << localIndex[1] << " " << localIndex[2] << std::endl);
68 
69  // Calculate the vectors between atoms
70  const Position & pos0 = p[0]->x[localIndex[0]].position;
71  const Position & pos1 = p[1]->x[localIndex[1]].position;
72  const Vector r12 = lattice.delta(pos0,pos1);
73  const Position & pos2 = p[2]->x[localIndex[2]].position;
74  const Vector r23 = lattice.delta(pos1,pos2);
75  const Position & pos3 = p[3]->x[localIndex[3]].position;
76  const Vector r34 = lattice.delta(pos2,pos3);
77 
78  // Calculate the cross products and distances
79  Vector A = cross(r12,r23);
80  register BigReal rAinv = A.rlength();
81  Vector B = cross(r23,r34);
82  register BigReal rBinv = B.rlength();
83  Vector C = cross(r23,A);
84  register BigReal rCinv = C.rlength();
85 
86  // Calculate the sin and cos
87  BigReal cos_phi = (A*B)*(rAinv*rBinv);
88  BigReal sin_phi = (C*B)*(rCinv*rBinv);
89 
90  BigReal phi= -atan2(sin_phi,cos_phi);
91 
92  BigReal K=0; // energy
93  BigReal K1=0; // force
94 
95  // get the dihedral information
96  int multiplicity = value->multiplicity;
97 
98  // Loop through the multiple parameter sets for this
99  // bond. We will only loop more than once if this
100  // has multiple parameter sets from Charmm22
101  for (int mult_num=0; mult_num<multiplicity; mult_num++)
102  {
103  /* get angle information */
104  Real k = value->values[mult_num].k * scale;
105  Real delta = value->values[mult_num].delta;
106  int n = value->values[mult_num].n;
107 
108  // Calculate the energy
109  if (n)
110  {
111  // Periodicity is greater than 0, so use cos form
112  K += k*(1+cos(n*phi - delta));
113  K1 += -n*k*sin(n*phi - delta);
114  }
115  else
116  {
117  // Periodicity is 0, so just use the harmonic form
118  BigReal diff = phi-delta;
119  if (diff < -PI) diff += TWOPI;
120  else if (diff > PI) diff -= TWOPI;
121 
122  K += k*diff*diff;
123  K1 += 2.0*k*diff;
124  }
125  } /* for multiplicity */
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, 4) ) {
131  case 1:
132  reduction[dihedralEnergyIndex_ti_1] += K;
133  reduction[dihedralEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*K;
134  K *= bond_lambda_1;
135  K1 *= bond_lambda_1;
136  break;
137  case 2:
138  reduction[dihedralEnergyIndex_ti_2] += K;
139  reduction[dihedralEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*K;
140  K *= bond_lambda_2;
141  K1 *= bond_lambda_2;
142  break;
143  }
144  }
145  //fepe
146 
147  Force f1,f2,f3;
148 
149  // Normalize B
150  //rB = 1.0/rB;
151  B *= rBinv;
152 
153  // Next, we want to calculate the forces. In order
154  // to do that, we first need to figure out whether the
155  // sin or cos form will be more stable. For this,
156  // just look at the value of phi
157  if (fabs(sin_phi) > 0.1)
158  {
159  // use the sin version to avoid 1/cos terms
160 
161  // Normalize A
162  A *= rAinv;
163  Vector dcosdA;
164  Vector dcosdB;
165 
166  dcosdA.x = rAinv*(cos_phi*A.x-B.x);
167  dcosdA.y = rAinv*(cos_phi*A.y-B.y);
168  dcosdA.z = rAinv*(cos_phi*A.z-B.z);
169 
170  dcosdB.x = rBinv*(cos_phi*B.x-A.x);
171  dcosdB.y = rBinv*(cos_phi*B.y-A.y);
172  dcosdB.z = rBinv*(cos_phi*B.z-A.z);
173 
174  K1 = K1/sin_phi;
175 
176  f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
177  f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
178  f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
179 
180  f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
181  f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
182  f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
183 
184  f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z
185  + r34.y*dcosdB.z - r34.z*dcosdB.y);
186  f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x
187  + r34.z*dcosdB.x - r34.x*dcosdB.z);
188  f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y
189  + r34.x*dcosdB.y - r34.y*dcosdB.x);
190  }
191  else
192  {
193  // This angle is closer to 0 or 180 than it is to
194  // 90, so use the cos version to avoid 1/sin terms
195 
196  // Normalize C
197  // rC = 1.0/rC;
198  C *= rCinv;
199 
200  Vector dsindC;
201  Vector dsindB;
202 
203  dsindC.x = rCinv*(sin_phi*C.x-B.x);
204  dsindC.y = rCinv*(sin_phi*C.y-B.y);
205  dsindC.z = rCinv*(sin_phi*C.z-B.z);
206 
207  dsindB.x = rBinv*(sin_phi*B.x-C.x);
208  dsindB.y = rBinv*(sin_phi*B.y-C.y);
209  dsindB.z = rBinv*(sin_phi*B.z-C.z);
210 
211  K1 = -K1/cos_phi;
212 
213  f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
214  - r23.x*r23.y*dsindC.y
215  - r23.x*r23.z*dsindC.z);
216  f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
217  - r23.y*r23.z*dsindC.z
218  - r23.y*r23.x*dsindC.x);
219  f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
220  - r23.z*r23.x*dsindC.x
221  - r23.z*r23.y*dsindC.y);
222 
223  f3 = cross(K1,dsindB,r23);
224 
225  f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
226  +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
227  +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
228  +dsindB.z*r34.y - dsindB.y*r34.z);
229  f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
230  +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
231  +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
232  +dsindB.x*r34.z - dsindB.z*r34.x);
233  f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
234  +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
235  +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
236  +dsindB.y*r34.x - dsindB.x*r34.y);
237  }
238 
239  /* store the forces */
240  // p[0]->f[localIndex[0]] += f1;
241  // p[1]->f[localIndex[1]] += f2 - f1;
242  // p[2]->f[localIndex[2]] += f3 - f2;
243  // p[3]->f[localIndex[3]] += -f3;
244 
245  p[0]->f[localIndex[0]].x += f1.x;
246  p[0]->f[localIndex[0]].y += f1.y;
247  p[0]->f[localIndex[0]].z += f1.z;
248 
249  p[1]->f[localIndex[1]].x += f2.x - f1.x;
250  p[1]->f[localIndex[1]].y += f2.y - f1.y;
251  p[1]->f[localIndex[1]].z += f2.z - f1.z;
252 
253  p[2]->f[localIndex[2]].x += f3.x - f2.x;
254  p[2]->f[localIndex[2]].y += f3.y - f2.y;
255  p[2]->f[localIndex[2]].z += f3.z - f2.z;
256 
257  p[3]->f[localIndex[3]].x += -f3.x;
258  p[3]->f[localIndex[3]].y += -f3.y;
259  p[3]->f[localIndex[3]].z += -f3.z;
260 
261  /* store the force for dihedral-only accelMD */
262  if ( p[0]->af ) {
263  p[0]->af[localIndex[0]].x += f1.x;
264  p[0]->af[localIndex[0]].y += f1.y;
265  p[0]->af[localIndex[0]].z += f1.z;
266 
267  p[1]->af[localIndex[1]].x += f2.x - f1.x;
268  p[1]->af[localIndex[1]].y += f2.y - f1.y;
269  p[1]->af[localIndex[1]].z += f2.z - f1.z;
270 
271  p[2]->af[localIndex[2]].x += f3.x - f2.x;
272  p[2]->af[localIndex[2]].y += f3.y - f2.y;
273  p[2]->af[localIndex[2]].z += f3.z - f2.z;
274 
275  p[3]->af[localIndex[3]].x += -f3.x;
276  p[3]->af[localIndex[3]].y += -f3.y;
277  p[3]->af[localIndex[3]].z += -f3.z;
278  }
279 
280  DebugM(3, "::computeForce() -- ending with delta energy " << K << std::endl);
281  reduction[dihedralEnergyIndex] += K;
282  reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
283  reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
284  reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
285  reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
286  reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
287  reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
288  reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
289  reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
290  reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
291 
292  if (pressureProfileData) {
293  BigReal z1 = p[0]->x[localIndex[0]].position.z;
294  BigReal z2 = p[1]->x[localIndex[1]].position.z;
295  BigReal z3 = p[2]->x[localIndex[2]].position.z;
296  BigReal z4 = p[3]->x[localIndex[3]].position.z;
297  int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
298  int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
299  int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
300  int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
305  int p1 = p[0]->x[localIndex[0]].partition;
306  int p2 = p[1]->x[localIndex[1]].partition;
307  int p3 = p[2]->x[localIndex[2]].partition;
308  int p4 = p[3]->x[localIndex[3]].partition;
309  int pn = pressureProfileAtomTypes;
311  p1, p2, pn,
312  f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
313  pressureProfileData);
315  p2, p3, pn,
316  f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
317  pressureProfileData);
319  p3, p4, pn,
320  f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
321  pressureProfileData);
322  }
323 
324  }
325 }
326 
327 
329 {
334  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
335  ADD_TENSOR(reduction,REDUCTION_VIRIAL_AMD_DIHE,data,virialIndex);
336 }
337 
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)
const BigReal A
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
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
DihedralValue * dihedral_array
Definition: Parameters.h:243
static void computeForce(DihedralElem *, int, BigReal *, BigReal *)
int get_fep_bonded_type(const int *atomID, unsigned int order) const
Definition: Molecule.h:1385
static BigReal pressureProfileMin
int localIndex[size]
Flags flags
Definition: Patch.h:127
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:108
static int pressureProfileAtomTypes
#define PI
Definition: common.h:81
void pp_clamp(int &n, int nslabs)
static BigReal pressureProfileThickness
BigReal rlength(void)
Definition: Vector.h:177
const DihedralValue * value
BigReal getBondLambda(const BigReal)
AtomID atomID[size]
BigReal getCurrentLambda2(const int)
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
int numDihedrals
Definition: Molecule.h:561
#define TWOPI
Definition: common.h:85
#define simParams
Definition: Output.C:127
BigReal y
Definition: Vector.h:66
const BigReal B
static void submitReductionData(BigReal *, SubmitReduction *)
BigReal getCurrentLambda(const int)
static int pressureProfileSlabs
TuplePatchElem * p[size]
Molecule * molecule
Definition: Node.h:176
static void getParameterPointers(Parameters *, const DihedralValue **)
static void getMoleculePointers(Molecule *, int *, int32 ***, Dihedral **)
double BigReal
Definition: common.h:112
int step
Definition: PatchTypes.h:16