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]
 
TuplePatchElem * p [size]
 
Real scale
 
const AngleValue * value
 

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.

46  { angleEnergyIndex, angleEnergyIndex_f, angleEnergyIndex_ti_1,
47  angleEnergyIndex_ti_2, TENSOR(virialIndex), reductionDataSize };
anonymous enum
Enumerator
reductionChecksumLabel 

Definition at line 48 of file ComputeAngles.h.

48 { reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM };

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 }
AngleElem::atomID
AtomID atomID[size]
Definition: ComputeAngles.h:20
TupleSignature::tupleParamType
Index tupleParamType
Definition: structures.h:202
AngleElem::value
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  }
angle::atom2
int32 atom2
Definition: structures.h:56
angle::atom1
int32 atom1
Definition: structures.h:55
AngleElem::atomID
AtomID atomID[size]
Definition: ComputeAngles.h:20
angle::atom3
int32 atom3
Definition: structures.h:57
angle::angle_type
Index angle_type
Definition: structures.h:58
AngleElem::value
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  }
AtomID
int AtomID
Definition: NamdTypes.h:29
AngleElem::atomID
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
52  SimParameters *const simParams = Node::Object()->simParameters;
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);
216  pp_clamp(n1, pressureProfileSlabs);
217  pp_clamp(n2, pressureProfileSlabs);
218  pp_clamp(n3, pressureProfileSlabs);
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;
223  pp_reduction(pressureProfileSlabs, n1, n2,
224  p1, p2, pn,
225  force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
226  pressureProfileData);
227  pp_reduction(pressureProfileSlabs, n3, n2,
228  p3, p2, pn,
229  force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
230  pressureProfileData);
231  }
232 
233  }
234 }
Node::Object
static Node * Object()
Definition: Node.h:86
CompAtom::partition
unsigned char partition
Definition: NamdTypes.h:56
AngleElem::p
TuplePatchElem * p[size]
Definition: ComputeAngles.h:22
AngleElem::pressureProfileThickness
static BigReal pressureProfileThickness
Definition: ComputeAngles.h:36
AtomID
int AtomID
Definition: NamdTypes.h:29
Patch::lattice
Lattice & lattice
Definition: Patch.h:126
pp_reduction
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
Definition: PressureProfile.h:22
Vector
Definition: Vector.h:64
Node::simParameters
SimParameters * simParameters
Definition: Node.h:178
AngleValue::r_ub
Real r_ub
Definition: Parameters.h:96
Real
float Real
Definition: common.h:109
DebugM
#define DebugM(x, y)
Definition: Debug.h:59
Vector::z
BigReal z
Definition: Vector.h:66
CompAtom::position
Position position
Definition: NamdTypes.h:53
AngleElem::pressureProfileMin
static BigReal pressureProfileMin
Definition: ComputeAngles.h:37
Molecule::get_fep_bonded_type
int get_fep_bonded_type(const int *atomID, unsigned int order) const
Definition: Molecule.h:1389
Vector::length
BigReal length(void) const
Definition: Vector.h:169
Patch::flags
Flags flags
Definition: Patch.h:127
pp_clamp
void pp_clamp(int &n, int nslabs)
Definition: PressureProfile.h:15
Vector::rlength
BigReal rlength(void)
Definition: Vector.h:177
AngleValue::normal
int normal
Definition: Parameters.h:97
SimParameters::getBondLambda
BigReal getBondLambda(const BigReal)
Definition: SimParameters.C:7267
SimParameters::getCurrentLambda2
BigReal getCurrentLambda2(const int)
Definition: SimParameters.C:7147
AngleElem::atomID
AtomID atomID[size]
Definition: ComputeAngles.h:20
Lattice::delta
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
Vector::x
BigReal x
Definition: Vector.h:66
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
AngleElem::localIndex
int localIndex[size]
Definition: ComputeAngles.h:21
simParams
#define simParams
Definition: Output.C:127
Vector::y
BigReal y
Definition: Vector.h:66
SimParameters::getCurrentLambda
BigReal getCurrentLambda(const int)
Definition: SimParameters.C:7188
AngleElem::pressureProfileSlabs
static int pressureProfileSlabs
Definition: ComputeAngles.h:34
AngleValue::theta0
Real theta0
Definition: Parameters.h:94
AngleElem::value
const AngleValue * value
Definition: ComputeAngles.h:40
AngleElem::pressureProfileAtomTypes
static int pressureProfileAtomTypes
Definition: ComputeAngles.h:35
Node::molecule
Molecule * molecule
Definition: Node.h:176
AngleValue::k_ub
Real k_ub
Definition: Parameters.h:95
BigReal
double BigReal
Definition: common.h:114
Flags::step
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 }
Molecule::numAngles
int numAngles
Definition: Molecule.h:561
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:85
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 }
Parameters::angle_array
AngleValue * angle_array
Definition: Parameters.h:244
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  }
AtomSignature::angleSigs
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  }
AngleElem::atomID
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  }
AngleElem::atomID
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  }
AngleElem::atomID
AtomID atomID[size]
Definition: ComputeAngles.h:20
void AngleElem::submitReductionData ( BigReal *  data,
SubmitReduction *  reduction 
)
static

Definition at line 237 of file ComputeAngles.C.

References ADD_TENSOR, angleEnergyIndex, angleEnergyIndex_f, angleEnergyIndex_ti_1, angleEnergyIndex_ti_2, SubmitReduction::item(), REDUCTION_ANGLE_ENERGY, REDUCTION_BONDED_ENERGY_F, REDUCTION_BONDED_ENERGY_TI_1, and REDUCTION_BONDED_ENERGY_TI_2.

238 {
239  reduction->item(REDUCTION_ANGLE_ENERGY) += data[angleEnergyIndex];
240  reduction->item(REDUCTION_BONDED_ENERGY_F) += data[angleEnergyIndex_f];
241  reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[angleEnergyIndex_ti_1];
242  reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[angleEnergyIndex_ti_2];
243  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
244 }
ADD_TENSOR
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:32
SubmitReduction::item
BigReal & item(int i)
Definition: ReductionMgr.h:312

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:
Generated by   doxygen 1.8.5