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

#include <ComputeAngles.h>

Public Types

enum  { size = 3 }
 
enum  {
  angleEnergyIndex, angleEnergyIndex_f, angleEnergyIndex_ti_1, angleEnergyIndex_ti_2,
  TENSOR =(virialIndex), reductionDataSize
}
 
enum  { reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM }
 

Public Member Functions

int hash () const
 
 AngleElem ()
 
 AngleElem (AtomID atom0, const TupleSignature *sig, const AngleValue *v)
 
 AngleElem (const Angle *a, const AngleValue *v)
 
 AngleElem (AtomID atom0, AtomID atom1, AtomID atom2)
 
 ~AngleElem ()
 
int operator== (const AngleElem &a) const
 
int operator< (const AngleElem &a) const
 

Static Public Member Functions

static void computeForce (AngleElem *, int, BigReal *, BigReal *)
 
static void getMoleculePointers (Molecule *, int *, int32 ***, Angle **)
 
static void getParameterPointers (Parameters *, const AngleValue **)
 
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 AngleValuevalue
 

Static Public Attributes

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

Detailed Description

Definition at line 16 of file ComputeAngles.h.

Member Enumeration Documentation

anonymous enum
Enumerator
size 

Definition at line 19 of file ComputeAngles.h.

19 { size = 3 };
anonymous enum
Enumerator
angleEnergyIndex 
angleEnergyIndex_f 
angleEnergyIndex_ti_1 
angleEnergyIndex_ti_2 
TENSOR 
reductionDataSize 

Definition at line 46 of file ComputeAngles.h.

anonymous enum
Enumerator
reductionChecksumLabel 

Definition at line 48 of file ComputeAngles.h.

Constructor & Destructor Documentation

AngleElem::AngleElem ( )
inline

Definition at line 12 of file ComputeAngles.inl.

12 { ; }
AngleElem::AngleElem ( AtomID  atom0,
const TupleSignature sig,
const AngleValue v 
)
inline

Definition at line 14 of file ComputeAngles.inl.

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

14  {
15  atomID[0] = atom0;
16  atomID[1] = atom0 + sig->offset[0];
17  atomID[2] = atom0 + sig->offset[1];
18  value = &v[sig->tupleParamType];
19 }
AtomID atomID[size]
Definition: ComputeAngles.h:20
Index tupleParamType
Definition: structures.h:202
const AngleValue * value
Definition: ComputeAngles.h:40
AngleElem::AngleElem ( const Angle a,
const AngleValue v 
)
inline

Definition at line 21 of file ComputeAngles.inl.

References angle::angle_type, angle::atom1, angle::atom2, angle::atom3, atomID, and value.

22  {
23  atomID[0] = a->atom1;
24  atomID[1] = a->atom2;
25  atomID[2] = a->atom3;
26  value = &v[a->angle_type];
27  }
int32 atom2
Definition: structures.h:56
int32 atom1
Definition: structures.h:55
AtomID atomID[size]
Definition: ComputeAngles.h:20
int32 atom3
Definition: structures.h:57
Index angle_type
Definition: structures.h:58
const AngleValue * value
Definition: ComputeAngles.h:40
AngleElem::AngleElem ( AtomID  atom0,
AtomID  atom1,
AtomID  atom2 
)
inline

Definition at line 29 of file ComputeAngles.inl.

References atomID.

30  {
31  if (atom0 > atom2) { // Swap end atoms so lowest is first!
32  AtomID tmp = atom2; atom2 = atom0; atom0 = tmp;
33  }
34  atomID[0] = atom0;
35  atomID[1] = atom1;
36  atomID[2] = atom2;
37  }
int AtomID
Definition: NamdTypes.h:29
AtomID atomID[size]
Definition: ComputeAngles.h:20
AngleElem::~AngleElem ( )
inline

Definition at line 55 of file ComputeAngles.h.

55 { };

Member Function Documentation

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

Definition at line 47 of file ComputeAngles.C.

References SimParameters::alchOn, angleEnergyIndex, angleEnergyIndex_f, angleEnergyIndex_ti_1, angleEnergyIndex_ti_2, atomID, DebugM, Lattice::delta(), TuplePatchElem::f, Patch::flags, Molecule::get_fep_bonded_type(), SimParameters::getBondLambda(), SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), AngleValue::k, AngleValue::k_ub, Patch::lattice, Vector::length(), localIndex, Node::molecule, NAMD_die(), AngleValue::normal, Node::Object(), p, TuplePatchElem::p, CompAtom::partition, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, AngleValue::r_ub, Vector::rlength(), scale, Node::simParameters, simParams, SimParameters::singleTopology, size, Flags::step, AngleValue::theta0, value, TuplePatchElem::x, Vector::x, Vector::y, and Vector::z.

48 {
49  const Lattice & lattice = tuples[0].p[0]->p->lattice;
50 
51  //fepb BKR
53  const int step = tuples[0].p[0]->p->flags.step;
54  const BigReal alchLambda = simParams->getCurrentLambda(step);
55  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
56  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
57  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
58  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
59  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
60  Molecule *const mol = Node::Object()->molecule;
61  //fepe
62 
63  for ( int ituple=0; ituple<ntuple; ++ituple ) {
64  const AngleElem &tup = tuples[ituple];
65  enum { size = 3 };
66  const AtomID (&atomID)[size](tup.atomID);
67  const int (&localIndex)[size](tup.localIndex);
68  TuplePatchElem * const(&p)[size](tup.p);
69  const Real (&scale)(tup.scale);
70  const AngleValue * const(&value)(tup.value);
71 
72  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
73  << localIndex[1] << " " << localIndex[2] << std::endl);
74 
75  const Position & pos1 = p[0]->x[localIndex[0]].position;
76  const Position & pos2 = p[1]->x[localIndex[1]].position;
77  Vector r12 = lattice.delta(pos1,pos2);
78  register BigReal d12inv = r12.rlength();
79  const Position & pos3 = p[2]->x[localIndex[2]].position;
80  Vector r32 = lattice.delta(pos3,pos2);
81  register BigReal d32inv = r32.rlength();
82  int normal = value->normal;
83 
84  BigReal cos_theta = (r12*r32)*(d12inv*d32inv);
85  // Make sure that the cosine value is acceptable. With roundoff, you
86  // can get values like 1.0+2e-16, which makes acos puke. So instead,
87  // just set these kinds of values to exactly 1.0
88  if (cos_theta > 1.0) cos_theta = 1.0;
89  else if (cos_theta < -1.0) cos_theta = -1.0;
90 
91  BigReal k = value->k * scale;
92  BigReal theta0 = value->theta0;
93 
94  // Get theta
95  BigReal theta = acos(cos_theta);
96 
97  // Compare it to the rest angle
98  BigReal diff;
99 
100  if (normal == 1) {
101  diff = theta - theta0;
102  } else {
103  diff = cos_theta - cos(theta0);
104  }
105 
106  // Add the energy from this angle to the total energy
107  BigReal energy = k *diff*diff;
108 
109  // Normalize vector r12 and r32
110  //BigReal d12inv = 1. / d12;
111  //BigReal d32inv = 1. / d32;
112 
113 
114  // Calculate constant factor 2k(theta-theta0)/sin(theta)
115  BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
116  if (normal != 1) {
117  diff *= (2.0* k);
118  } else if ( sin_theta < 1.e-6 ) {
119  // catch case where bonds are parallel
120  // this gets the force approximately right for theta0 of 0 or pi
121  // and at least avoids small division for other values
122  if ( diff < 0. ) diff = 2.0 * k;
123  else diff = -2.0 * k;
124  } else {
125  diff *= (-2.0* k) / sin_theta;
126  }
127  BigReal c1 = diff * d12inv;
128  BigReal c2 = diff * d32inv;
129 
130  // Calculate the actual forces
131  Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
132  Force force2 = force1;
133  Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
134  force2 += force3; force2 *= -1;
135 
136  // Check to see if we need to do the Urey-Bradley term
137  if (value->k_ub)
138  {
139  // Non-zero k_ub value, so calculate the harmonic
140  // potential between the 1-3 atoms
141 
142  if (normal != 1) {
143  NAMD_die("ERROR: Can't use cosAngles with Urey-Bradley angles");
144  }
145  BigReal k_ub = value->k_ub;
146  BigReal r_ub = value->r_ub;
147  Vector r13 = r12 - r32;
148  BigReal d13 = r13.length();
149  diff = d13- r_ub;
150 
151  energy += k_ub *diff*diff;
152 
153  diff *= -2.0*k_ub / d13;
154  r13 *= diff;
155 
156  force1 += r13;
157  force3 -= r13;
158  }
159 
160  //fepb - BKR scaling of alchemical bonded terms
161  // NB: TI derivative is the _unscaled_ energy.
162  if ( simParams->alchOn && !simParams->singleTopology) {
163  switch ( mol->get_fep_bonded_type(atomID, 3) ) {
164  case 1:
165  reduction[angleEnergyIndex_ti_1] += energy;
166  reduction[angleEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
167  energy *= bond_lambda_1;
168  force1 *= bond_lambda_1;
169  force2 *= bond_lambda_1;
170  force3 *= bond_lambda_1;
171  break;
172  case 2:
173  reduction[angleEnergyIndex_ti_2] += energy;
174  reduction[angleEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
175  energy *= bond_lambda_2;
176  force1 *= bond_lambda_2;
177  force2 *= bond_lambda_2;
178  force3 *= bond_lambda_2;
179  break;
180  }
181  }
182  //fepe
183 
184  p[0]->f[localIndex[0]].x += force1.x;
185  p[0]->f[localIndex[0]].y += force1.y;
186  p[0]->f[localIndex[0]].z += force1.z;
187 
188  p[1]->f[localIndex[1]].x += force2.x;
189  p[1]->f[localIndex[1]].y += force2.y;
190  p[1]->f[localIndex[1]].z += force2.z;
191 
192  p[2]->f[localIndex[2]].x += force3.x;
193  p[2]->f[localIndex[2]].y += force3.y;
194  p[2]->f[localIndex[2]].z += force3.z;
195 
196  DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
197 
198  reduction[angleEnergyIndex] += energy;
199  reduction[virialIndex_XX] += ( force1.x * r12.x + force3.x * r32.x );
200  reduction[virialIndex_XY] += ( force1.x * r12.y + force3.x * r32.y );
201  reduction[virialIndex_XZ] += ( force1.x * r12.z + force3.x * r32.z );
202  reduction[virialIndex_YX] += ( force1.y * r12.x + force3.y * r32.x );
203  reduction[virialIndex_YY] += ( force1.y * r12.y + force3.y * r32.y );
204  reduction[virialIndex_YZ] += ( force1.y * r12.z + force3.y * r32.z );
205  reduction[virialIndex_ZX] += ( force1.z * r12.x + force3.z * r32.x );
206  reduction[virialIndex_ZY] += ( force1.z * r12.y + force3.z * r32.y );
207  reduction[virialIndex_ZZ] += ( force1.z * r12.z + force3.z * r32.z );
208 
209  if (pressureProfileData) {
210  BigReal z1 = p[0]->x[localIndex[0]].position.z;
211  BigReal z2 = p[1]->x[localIndex[1]].position.z;
212  BigReal z3 = p[2]->x[localIndex[2]].position.z;
213  int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
214  int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
215  int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
219  int p1 = p[0]->x[localIndex[0]].partition;
220  int p2 = p[1]->x[localIndex[1]].partition;
221  int p3 = p[2]->x[localIndex[2]].partition;
222  int pn = pressureProfileAtomTypes;
224  p1, p2, pn,
225  force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
226  pressureProfileData);
228  p3, p2, pn,
229  force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
230  pressureProfileData);
231  }
232 
233  }
234 }
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
TuplePatchElem * p[size]
Definition: ComputeAngles.h:22
static BigReal pressureProfileThickness
Definition: ComputeAngles.h:36
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
Real r_ub
Definition: Parameters.h:94
float Real
Definition: common.h:107
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
static BigReal pressureProfileMin
Definition: ComputeAngles.h:37
int get_fep_bonded_type(const int *atomID, unsigned int order) const
Definition: Molecule.h:1385
BigReal length(void) const
Definition: Vector.h:169
Flags flags
Definition: Patch.h:127
void pp_clamp(int &n, int nslabs)
BigReal rlength(void)
Definition: Vector.h:177
int normal
Definition: Parameters.h:95
BigReal getBondLambda(const BigReal)
BigReal getCurrentLambda2(const int)
AtomID atomID[size]
Definition: ComputeAngles.h:20
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 localIndex[size]
Definition: ComputeAngles.h:21
#define simParams
Definition: Output.C:127
BigReal y
Definition: Vector.h:66
BigReal getCurrentLambda(const int)
static int pressureProfileSlabs
Definition: ComputeAngles.h:34
Real theta0
Definition: Parameters.h:92
const AngleValue * value
Definition: ComputeAngles.h:40
static int pressureProfileAtomTypes
Definition: ComputeAngles.h:35
Molecule * molecule
Definition: Node.h:176
Real k_ub
Definition: Parameters.h:93
double BigReal
Definition: common.h:112
int step
Definition: PatchTypes.h:16
void AngleElem::getMoleculePointers ( Molecule mol,
int *  count,
int32 ***  byatom,
Angle **  structarray 
)
static

Definition at line 32 of file ComputeAngles.C.

References NAMD_die(), and Molecule::numAngles.

33 {
34 #ifdef MEM_OPT_VERSION
35  NAMD_die("Should not be called in AngleElem::getMoleculePointers in memory optimized version!");
36 #else
37  *count = mol->numAngles;
38  *byatom = mol->anglesByAtom;
39  *structarray = mol->angles;
40 #endif
41 }
int numAngles
Definition: Molecule.h:560
void NAMD_die(const char *err_msg)
Definition: common.C:83
void AngleElem::getParameterPointers ( Parameters p,
const AngleValue **  v 
)
static

Definition at line 43 of file ComputeAngles.C.

References Parameters::angle_array.

43  {
44  *v = p->angle_array;
45 }
AngleValue * angle_array
Definition: Parameters.h:242
static void AngleElem::getTupleInfo ( AtomSignature sig,
int *  count,
TupleSignature **  t 
)
inlinestatic

Definition at line 28 of file ComputeAngles.h.

References AtomSignature::angleCnt, and AtomSignature::angleSigs.

28  {
29  *count = sig->angleCnt;
30  *t = sig->angleSigs;
31  }
TupleSignature * angleSigs
Definition: structures.h:334
int AngleElem::hash ( void  ) const
inline

Definition at line 42 of file ComputeAngles.h.

References atomID.

42  {
43  return 0x7FFFFFFF & ((atomID[0]<<22) + (atomID[1]<<11) + (atomID[2]));
44  }
AtomID atomID[size]
Definition: ComputeAngles.h:20
int AngleElem::operator< ( const AngleElem a) const
inline

Definition at line 45 of file ComputeAngles.inl.

References atomID.

46  {
47  return (atomID[0] < a.atomID[0] ||
48  (atomID[0] == a.atomID[0] &&
49  (atomID[1] < a.atomID[1] ||
50  (atomID[1] == a.atomID[1] &&
51  atomID[2] < a.atomID[2]) )));
52  }
AtomID atomID[size]
Definition: ComputeAngles.h:20
int AngleElem::operator== ( const AngleElem a) const
inline

Definition at line 39 of file ComputeAngles.inl.

References atomID.

40  {
41  return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
42  a.atomID[2] == atomID[2]);
43  }
AtomID atomID[size]
Definition: ComputeAngles.h:20
void AngleElem::submitReductionData ( BigReal data,
SubmitReduction reduction 
)
static

Member Data Documentation

AtomID AngleElem::atomID[size]

Definition at line 20 of file ComputeAngles.h.

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

int AngleElem::localIndex[size]

Definition at line 21 of file ComputeAngles.h.

Referenced by computeForce().

TuplePatchElem* AngleElem::p[size]

Definition at line 22 of file ComputeAngles.h.

Referenced by computeForce().

int AngleElem::pressureProfileAtomTypes = 1
static

Definition at line 35 of file ComputeAngles.h.

Referenced by computeForce().

BigReal AngleElem::pressureProfileMin = 0
static

Definition at line 37 of file ComputeAngles.h.

Referenced by computeForce().

int AngleElem::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 34 of file ComputeAngles.h.

Referenced by computeForce().

BigReal AngleElem::pressureProfileThickness = 0
static

Definition at line 36 of file ComputeAngles.h.

Referenced by computeForce().

Real AngleElem::scale

Definition at line 23 of file ComputeAngles.h.

Referenced by computeForce().

const AngleValue* AngleElem::value

Definition at line 40 of file ComputeAngles.h.

Referenced by AngleElem(), and computeForce().


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