NAMD
ComputeImpropers.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeImpropers.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, Improper** structarray)
26 {
27 #ifdef MEM_OPT_VERSION
28  NAMD_die("Should not be called in ImproperElem::getMoleculePointers in memory optimized version!");
29 #else
30  *count = mol->numImpropers;
31  *byatom = mol->impropersByAtom;
32  *structarray = mol->impropers;
33 #endif
34 }
35 
37  *v = p->improper_array;
38 }
39 
40 void ImproperElem::computeForce(ImproperElem *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 ImproperElem &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 ImproperValue * const(&value)(tup.value);
65 
66  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
67  << localIndex[1] << " " << localIndex[2] << " " <<
68  localIndex[3] << std::endl);
69 
70  // Vector r12, r23, r34; // vector between atoms
71  Vector A,B,C; // cross products
72  BigReal rA, rB, rC; // length of vectors A, B, and C
73  BigReal energy=0; // energy from the angle
74  BigReal phi; // angle between the plans
75  double cos_phi; // cos(phi)
76  double sin_phi; // sin(phi)
77  Vector dcosdA; // Derivative d(cos(phi))/dA
78  Vector dcosdB; // Derivative d(cos(phi))/dB
79  Vector dsindC; // Derivative d(sin(phi))/dC
80  Vector dsindB; // Derivative d(sin(phi))/dB
81  BigReal K,K1; // energy constants
82  BigReal diff; // for periodicity
83  Force f1(0,0,0),f2(0,0,0),f3(0,0,0); // force components
84 
85  //DebugM(3, "::computeForce() -- starting with improper type " << improperType << std::endl);
86 
87  // get the improper information
88  int multiplicity = value->multiplicity;
89 
90  // Calculate the vectors between atoms
91  const Position & pos0 = p[0]->x[localIndex[0]].position;
92  const Position & pos1 = p[1]->x[localIndex[1]].position;
93  const Position & pos2 = p[2]->x[localIndex[2]].position;
94  const Position & pos3 = p[3]->x[localIndex[3]].position;
95  const Vector r12 = lattice.delta(pos0,pos1);
96  const Vector r23 = lattice.delta(pos1,pos2);
97  const Vector r34 = lattice.delta(pos2,pos3);
98 
99  // Calculate the cross products
100  A = cross(r12,r23);
101  B = cross(r23,r34);
102  C = cross(r23,A);
103 
104  // Calculate the distances
105  rA = A.length();
106  rB = B.length();
107  rC = C.length();
108 
109  // Calculate the sin and cos
110  cos_phi = A*B/(rA*rB);
111  sin_phi = C*B/(rC*rB);
112 
113  // Normalize B
114  rB = 1.0/rB;
115  B *= rB;
116 
117  phi= -atan2(sin_phi,cos_phi);
118 
119  if (fabs(sin_phi) > 0.1)
120  {
121  // Normalize A
122  rA = 1.0/rA;
123  A *= rA;
124  dcosdA = rA*(cos_phi*A-B);
125  dcosdB = rB*(cos_phi*B-A);
126  }
127  else
128  {
129  // Normalize C
130  rC = 1.0/rC;
131  C *= rC;
132  dsindC = rC*(sin_phi*C-B);
133  dsindB = rB*(sin_phi*B-C);
134  }
135 
136  // Loop through the multiple parameter sets for this
137  // bond. We will only loop more than once if this
138  // has multiple parameter sets from Charmm22
139  for (int mult_num=0; mult_num<multiplicity; mult_num++)
140  {
141  /* get angle information */
142  Real k = value->values[mult_num].k * scale;
143  Real delta = value->values[mult_num].delta;
144  int n = value->values[mult_num].n;
145 
146  // Calculate the energy
147  if (n)
148  {
149  // Periodicity is greater than 0, so use cos form
150  K = k*(1+cos(n*phi - delta));
151  K1 = -n*k*sin(n*phi - delta);
152  }
153  else
154  {
155  // Periodicity is 0, so just use the harmonic form
156  diff = phi-delta;
157  if (diff < -PI) diff += TWOPI;
158  else if (diff > PI) diff -= TWOPI;
159 
160  K = k*diff*diff;
161  K1 = 2.0*k*diff;
162  }
163 
164  // Add the energy from this improper to the total energy
165  energy += K;
166 
167  // Next, we want to calculate the forces. In order
168  // to do that, we first need to figure out whether the
169  // sin or cos form will be more stable. For this,
170  // just look at the value of phi
171  if (fabs(sin_phi) > 0.1)
172  {
173  // use the sin version to avoid 1/cos terms
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  K1 = -K1/cos_phi;
196 
197  f1.x += K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
198  - r23.x*r23.y*dsindC.y
199  - r23.x*r23.z*dsindC.z);
200  f1.y += K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
201  - r23.y*r23.z*dsindC.z
202  - r23.y*r23.x*dsindC.x);
203  f1.z += K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
204  - r23.z*r23.x*dsindC.x
205  - r23.z*r23.y*dsindC.y);
206 
207  f3 += cross(K1,dsindB,r23);
208 
209  f2.x += K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
210  +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
211  +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
212  +dsindB.z*r34.y - dsindB.y*r34.z);
213  f2.y += K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
214  +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
215  +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
216  +dsindB.x*r34.z - dsindB.z*r34.x);
217  f2.z += K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
218  +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
219  +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
220  +dsindB.y*r34.x - dsindB.x*r34.y);
221  }
222  } /* for multiplicity */
223 
224  //fepb - BKR scaling of alchemical bonded terms
225  // NB: TI derivative is the _unscaled_ energy.
226  if ( simParams->alchOn && !simParams->singleTopology) {
227  switch ( mol->get_fep_bonded_type(atomID, 4) ) {
228  case 1:
229  reduction[improperEnergyIndex_ti_1] += energy;
230  reduction[improperEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1) *
231  energy;
232  energy *= bond_lambda_1;
233  f1 *= bond_lambda_1;
234  f2 *= bond_lambda_1;
235  f3 *= bond_lambda_1;
236  break;
237  case 2:
238  reduction[improperEnergyIndex_ti_2] += energy;
239  reduction[improperEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2) *
240  energy;
241  energy *= bond_lambda_2;
242  f1 *= bond_lambda_2;
243  f2 *= bond_lambda_2;
244  f3 *= bond_lambda_2;
245  break;
246  }
247  }
248  //fepe
249 
250  /* store the forces */
251  p[0]->f[localIndex[0]] += f1;
252  p[1]->f[localIndex[1]] += f2 - f1;
253  p[2]->f[localIndex[2]] += f3 - f2;
254  p[3]->f[localIndex[3]] += -f3;
255 
256  DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
257  reduction[improperEnergyIndex] += energy;
258  reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
259  reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
260  reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
261  reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
262  reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
263  reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
264  reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
265  reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
266  reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
267 
268  if (pressureProfileData) {
269  BigReal z1 = p[0]->x[localIndex[0]].position.z;
270  BigReal z2 = p[1]->x[localIndex[1]].position.z;
271  BigReal z3 = p[2]->x[localIndex[2]].position.z;
272  BigReal z4 = p[3]->x[localIndex[3]].position.z;
273  int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
274  int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
275  int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
276  int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
281  int p1 = p[0]->x[localIndex[0]].partition;
282  int p2 = p[1]->x[localIndex[1]].partition;
283  int p3 = p[2]->x[localIndex[2]].partition;
284  int p4 = p[3]->x[localIndex[3]].partition;
285  int pn = pressureProfileAtomTypes;
287  p1, p2, pn,
288  f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
289  pressureProfileData);
291  p2, p3, pn,
292  f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
293  pressureProfileData);
295  p3, p4, pn,
296  f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
297  pressureProfileData);
298  }
299 
300  }
301 }
302 
304 {
309  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
310 }
311 
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
static int pressureProfileAtomTypes
static void submitReductionData(BigReal *, SubmitReduction *)
static void getMoleculePointers(Molecule *, int *, int32 ***, Improper **)
TuplePatchElem * p[size]
const ImproperValue * value
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)
AtomID atomID[size]
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
int get_fep_bonded_type(const int *atomID, unsigned int order) const
Definition: Molecule.h:1385
BigReal length(void) const
Definition: Vector.h:169
static void getParameterPointers(Parameters *, const ImproperValue **)
static BigReal pressureProfileMin
Flags flags
Definition: Patch.h:127
#define PI
Definition: common.h:81
void pp_clamp(int &n, int nslabs)
static int pressureProfileSlabs
BigReal getBondLambda(const BigReal)
BigReal getCurrentLambda2(const int)
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
BigReal x
Definition: Vector.h:66
ImproperValue * improper_array
Definition: Parameters.h:244
void NAMD_die(const char *err_msg)
Definition: common.C:83
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:114
#define TWOPI
Definition: common.h:85
int numImpropers
Definition: Molecule.h:566
#define simParams
Definition: Output.C:127
static void computeForce(ImproperElem *, int, BigReal *, BigReal *)
int localIndex[size]
BigReal y
Definition: Vector.h:66
const BigReal B
BigReal getCurrentLambda(const int)
static BigReal pressureProfileThickness
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:112
int step
Definition: PatchTypes.h:16