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

#include <ComputeImpropers.h>

Public Types

enum  { size = 4 }
 
enum  {
  improperEnergyIndex, improperEnergyIndex_f, improperEnergyIndex_ti_1, improperEnergyIndex_ti_2,
  TENSOR =(virialIndex), reductionDataSize
}
 
enum  { reductionChecksumLabel = REDUCTION_IMPROPER_CHECKSUM }
 

Public Member Functions

int hash () const
 
 ImproperElem ()
 
 ImproperElem (AtomID atom0, const TupleSignature *sig, const ImproperValue *v)
 
 ImproperElem (const Improper *a, const ImproperValue *v)
 
 ImproperElem (AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3)
 
 ~ImproperElem ()
 
int operator== (const ImproperElem &a) const
 
int operator< (const ImproperElem &a) const
 

Static Public Member Functions

static void computeForce (ImproperElem *, int, BigReal *, BigReal *)
 
static void getMoleculePointers (Molecule *, int *, int32 ***, Improper **)
 
static void getParameterPointers (Parameters *, const ImproperValue **)
 
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 ImproperValuevalue
 

Static Public Attributes

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

Detailed Description

Definition at line 17 of file ComputeImpropers.h.

Member Enumeration Documentation

anonymous enum
Enumerator
size 

Definition at line 20 of file ComputeImpropers.h.

anonymous enum
anonymous enum

Constructor & Destructor Documentation

ImproperElem::ImproperElem ( )
inline

Definition at line 51 of file ComputeImpropers.h.

51 { ; }
ImproperElem::ImproperElem ( AtomID  atom0,
const TupleSignature sig,
const ImproperValue v 
)
inline

Definition at line 53 of file ComputeImpropers.h.

References atomID, TupleSignature::offset, TupleSignature::tupleParamType, and value.

53  {
54  atomID[0] = atom0;
55  atomID[1] = atom0 + sig->offset[0];
56  atomID[2] = atom0 + sig->offset[1];
57  atomID[3] = atom0 + sig->offset[2];
58  value = &v[sig->tupleParamType];
59  }
const ImproperValue * value
AtomID atomID[size]
Index tupleParamType
Definition: structures.h:202
ImproperElem::ImproperElem ( const Improper a,
const ImproperValue v 
)
inline

Definition at line 61 of file ComputeImpropers.h.

References improper::atom1, improper::atom2, improper::atom3, improper::atom4, atomID, improper::improper_type, and value.

61  {
62  atomID[0] = a->atom1;
63  atomID[1] = a->atom2;
64  atomID[2] = a->atom3;
65  atomID[3] = a->atom4;
66  value = &v[a->improper_type];
67  }
int32 atom4
Definition: structures.h:75
Index improper_type
Definition: structures.h:76
const ImproperValue * value
AtomID atomID[size]
int32 atom1
Definition: structures.h:72
int32 atom2
Definition: structures.h:73
int32 atom3
Definition: structures.h:74
ImproperElem::ImproperElem ( AtomID  atom0,
AtomID  atom1,
AtomID  atom2,
AtomID  atom3 
)
inline

Definition at line 69 of file ComputeImpropers.h.

References atomID.

69  {
70  if (atom0 > atom3) { // Swap end atoms so lowest is first!
71  AtomID tmp = atom3; atom3 = atom0; atom0 = tmp;
72  tmp = atom1; atom1 = atom2; atom2 = tmp;
73  }
74  atomID[0] = atom0;
75  atomID[1] = atom1;
76  atomID[2] = atom2;
77  atomID[3] = atom3;
78  }
int AtomID
Definition: NamdTypes.h:29
AtomID atomID[size]
ImproperElem::~ImproperElem ( )
inline

Definition at line 79 of file ComputeImpropers.h.

79 {};

Member Function Documentation

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

Definition at line 40 of file ComputeImpropers.C.

References A, SimParameters::alchOn, atomID, B, cross(), DebugM, four_body_consts::delta, Lattice::delta(), TuplePatchElem::f, Patch::flags, Molecule::get_fep_bonded_type(), SimParameters::getBondLambda(), SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), improperEnergyIndex, improperEnergyIndex_f, improperEnergyIndex_ti_1, improperEnergyIndex_ti_2, four_body_consts::k, Patch::lattice, Vector::length(), localIndex, Node::molecule, ImproperValue::multiplicity, four_body_consts::n, Node::Object(), p, TuplePatchElem::p, CompAtom::partition, PI, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, scale, Node::simParameters, simParams, SimParameters::singleTopology, size, Flags::step, TWOPI, value, ImproperValue::values, TuplePatchElem::x, Vector::x, Vector::y, and Vector::z.

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 }
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
static int pressureProfileAtomTypes
TuplePatchElem * p[size]
const ImproperValue * value
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)
const BigReal A
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
float Real
Definition: common.h:109
#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 BigReal pressureProfileMin
Flags flags
Definition: Patch.h:127
#define PI
Definition: common.h:83
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
FourBodyConsts values[MAX_MULTIPLICITY]
Definition: Parameters.h:114
#define TWOPI
Definition: common.h:87
#define simParams
Definition: Output.C:127
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:114
int step
Definition: PatchTypes.h:16
void ImproperElem::getMoleculePointers ( Molecule mol,
int *  count,
int32 ***  byatom,
Improper **  structarray 
)
static

Definition at line 25 of file ComputeImpropers.C.

References NAMD_die(), and Molecule::numImpropers.

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 }
void NAMD_die(const char *err_msg)
Definition: common.C:83
int numImpropers
Definition: Molecule.h:566
void ImproperElem::getParameterPointers ( Parameters p,
const ImproperValue **  v 
)
static

Definition at line 36 of file ComputeImpropers.C.

References Parameters::improper_array.

36  {
37  *v = p->improper_array;
38 }
ImproperValue * improper_array
Definition: Parameters.h:244
static void ImproperElem::getTupleInfo ( AtomSignature sig,
int *  count,
TupleSignature **  t 
)
inlinestatic

Definition at line 29 of file ComputeImpropers.h.

References AtomSignature::improperCnt, and AtomSignature::improperSigs.

29  {
30  *count = sig->improperCnt;
31  *t = sig->improperSigs;
32  }
TupleSignature * improperSigs
Definition: structures.h:336
int ImproperElem::hash ( void  ) const
inline

Definition at line 43 of file ComputeImpropers.h.

References atomID.

43  {
44  return 0x7FFFFFFF &((atomID[0]<<24) + (atomID[1]<<16) + (atomID[2]<<8) + atomID[3]);
45  }
AtomID atomID[size]
int ImproperElem::operator< ( const ImproperElem a) const
inline

Definition at line 86 of file ComputeImpropers.h.

References atomID.

86  {
87  return (atomID[0] < a.atomID[0] ||
88  (atomID[0] == a.atomID[0] &&
89  (atomID[1] < a.atomID[1] ||
90  (atomID[1] == a.atomID[1] &&
91  (atomID[2] < a.atomID[2] ||
92  (atomID[2] == a.atomID[2] &&
93  atomID[3] < a.atomID[3]
94  ))))));
95  }
AtomID atomID[size]
int ImproperElem::operator== ( const ImproperElem a) const
inline

Definition at line 81 of file ComputeImpropers.h.

References atomID.

81  {
82  return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
83  a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3]);
84  }
AtomID atomID[size]
void ImproperElem::submitReductionData ( BigReal data,
SubmitReduction reduction 
)
static

Member Data Documentation

AtomID ImproperElem::atomID[size]

Definition at line 21 of file ComputeImpropers.h.

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

int ImproperElem::localIndex[size]

Definition at line 22 of file ComputeImpropers.h.

Referenced by computeForce().

TuplePatchElem* ImproperElem::p[size]

Definition at line 23 of file ComputeImpropers.h.

Referenced by computeForce().

int ImproperElem::pressureProfileAtomTypes = 1
static

Definition at line 36 of file ComputeImpropers.h.

Referenced by computeForce().

BigReal ImproperElem::pressureProfileMin = 0
static

Definition at line 38 of file ComputeImpropers.h.

Referenced by computeForce().

int ImproperElem::pressureProfileSlabs = 0
static

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 35 of file ComputeImpropers.h.

Referenced by computeForce().

BigReal ImproperElem::pressureProfileThickness = 0
static

Definition at line 37 of file ComputeImpropers.h.

Referenced by computeForce().

Real ImproperElem::scale

Definition at line 24 of file ComputeImpropers.h.

Referenced by computeForce().

const ImproperValue* ImproperElem::value

Definition at line 41 of file ComputeImpropers.h.

Referenced by computeForce(), and ImproperElem().


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