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

anonymous enum
Enumerator
size 

Definition at line 19 of file ComputeAngles.h.

19 { size = 3 };

◆ anonymous enum

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

anonymous enum
Enumerator
reductionChecksumLabel 

Definition at line 48 of file ComputeAngles.h.

48 { reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM };

Constructor & Destructor Documentation

◆ AngleElem() [1/4]

AngleElem::AngleElem ( )
inline

Definition at line 12 of file ComputeAngles.inl.

12 { ; }

◆ AngleElem() [2/4]

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:204
AngleElem::value
const AngleValue * value
Definition: ComputeAngles.h:40

◆ AngleElem() [3/4]

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:58
angle::atom1
int32 atom1
Definition: structures.h:57
AngleElem::atomID
AtomID atomID[size]
Definition: ComputeAngles.h:20
angle::atom3
int32 atom3
Definition: structures.h:59
angle::angle_type
Index angle_type
Definition: structures.h:60
AngleElem::value
const AngleValue * value
Definition: ComputeAngles.h:40

◆ AngleElem() [4/4]

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  }
AngleElem::atomID
AtomID atomID[size]
Definition: ComputeAngles.h:20
AtomID
int32 AtomID
Definition: NamdTypes.h:35

◆ ~AngleElem()

AngleElem::~AngleElem ( )
inline

Definition at line 55 of file ComputeAngles.h.

55 { };

Member Function Documentation

◆ computeForce()

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

Definition at line 47 of file ComputeAngles.C.

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

◆ getMoleculePointers()

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:589
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:147

◆ getParameterPointers()

void AngleElem::getParameterPointers ( Parameters *  p,
const AngleValue **  v 
)
static

Definition at line 43 of file ComputeAngles.C.

References p.

43  {
44  *v = p->angle_array;
45 }
AngleElem::p
TuplePatchElem * p[size]
Definition: ComputeAngles.h:22

◆ getTupleInfo()

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:336

◆ hash()

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

◆ operator<()

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

◆ operator==()

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

◆ submitReductionData()

void AngleElem::submitReductionData ( BigReal *  data,
SubmitReduction *  reduction 
)
static

Definition at line 258 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.

259 {
260  reduction->item(REDUCTION_ANGLE_ENERGY) += data[angleEnergyIndex];
261  reduction->item(REDUCTION_BONDED_ENERGY_F) += data[angleEnergyIndex_f];
262  reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[angleEnergyIndex_ti_1];
263  reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[angleEnergyIndex_ti_2];
264  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
265 }
ADD_TENSOR
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:33
SubmitReduction::item
BigReal & item(int i)
Definition: ReductionMgr.h:313

Member Data Documentation

◆ atomID

AtomID AngleElem::atomID[size]

Definition at line 20 of file ComputeAngles.h.

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

◆ localIndex

int AngleElem::localIndex[size]

Definition at line 21 of file ComputeAngles.h.

Referenced by computeForce().

◆ p

TuplePatchElem* AngleElem::p[size]

Definition at line 22 of file ComputeAngles.h.

Referenced by computeForce(), and getParameterPointers().

◆ pressureProfileAtomTypes

int AngleElem::pressureProfileAtomTypes = 1
static

Definition at line 35 of file ComputeAngles.h.

Referenced by computeForce().

◆ pressureProfileMin

BigReal AngleElem::pressureProfileMin = 0
static

Definition at line 37 of file ComputeAngles.h.

Referenced by computeForce().

◆ pressureProfileSlabs

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().

◆ pressureProfileThickness

BigReal AngleElem::pressureProfileThickness = 0
static

Definition at line 36 of file ComputeAngles.h.

Referenced by computeForce().

◆ scale

Real AngleElem::scale

Definition at line 23 of file ComputeAngles.h.

Referenced by computeForce().

◆ value

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.14