AngleElem Class Reference

#include <ComputeAngles.h>

List of all members.

Public Types

 size = 3
 angleEnergyIndex
 angleEnergyIndex_f
 angleEnergyIndex_ti_1
 angleEnergyIndex_ti_2
 virialIndex
 reductionDataSize
 reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM
enum  { size = 3 }
enum  {
  angleEnergyIndex, angleEnergyIndex_f, angleEnergyIndex_ti_1, angleEnergyIndex_ti_2,
  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.

00019 { size = 3 };

anonymous enum

Enumerator:
angleEnergyIndex 
angleEnergyIndex_f 
angleEnergyIndex_ti_1 
angleEnergyIndex_ti_2 
virialIndex 
reductionDataSize 

Definition at line 46 of file ComputeAngles.h.

00046        { angleEnergyIndex, angleEnergyIndex_f, angleEnergyIndex_ti_1, 
00047          angleEnergyIndex_ti_2, TENSOR(virialIndex), reductionDataSize };

anonymous enum

Enumerator:
reductionChecksumLabel 

Definition at line 48 of file ComputeAngles.h.

00048 { reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM };

Constructor & Destructor Documentation

AngleElem::AngleElem (  )  [inline]

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

Definition at line 12 of file ComputeAngles.inl.

00012 { ; }

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.

00014                                                                                        {
00015     atomID[0] = atom0;
00016     atomID[1] = atom0 + sig->offset[0];
00017     atomID[2] = atom0 + sig->offset[1];
00018     value = &v[sig->tupleParamType];
00019 }

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.

00022   {
00023     atomID[0] = a->atom1;
00024     atomID[1] = a->atom2;
00025     atomID[2] = a->atom3;
00026     value = &v[a->angle_type];
00027   }

AngleElem::AngleElem ( AtomID  atom0,
AtomID  atom1,
AtomID  atom2 
) [inline]

Definition at line 29 of file ComputeAngles.inl.

References atomID.

00030   {
00031     if (atom0 > atom2) {  // Swap end atoms so lowest is first!
00032       AtomID tmp = atom2; atom2 = atom0; atom0 = tmp; 
00033     }
00034     atomID[0] = atom0;
00035     atomID[1] = atom1;
00036     atomID[2] = atom2;
00037   }

AngleElem::~AngleElem (  )  [inline]

Definition at line 55 of file ComputeAngles.h.

00055 { };

Member Function Documentation

void AngleElem::computeForce ( AngleElem *  ,
int  ,
BigReal *  ,
BigReal *   
) [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(), TuplePatchElem::p, 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, Vector::x, TuplePatchElem::x, Vector::y, and Vector::z.

00048 {
00049  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00050 
00051  //fepb BKR
00052  SimParameters *const simParams = Node::Object()->simParameters;
00053  const int step = tuples[0].p[0]->p->flags.step;
00054  const BigReal alchLambda = simParams->getCurrentLambda(step);
00055  const BigReal alchLambda2 = simParams->alchLambda2;
00056  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
00057  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
00058  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
00059  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
00060  Molecule *const mol = Node::Object()->molecule;
00061  //fepe
00062 
00063  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00064   const AngleElem &tup = tuples[ituple];
00065   enum { size = 3 };
00066   const AtomID (&atomID)[size](tup.atomID);
00067   const int    (&localIndex)[size](tup.localIndex);
00068   TuplePatchElem * const(&p)[size](tup.p);
00069   const Real (&scale)(tup.scale);
00070   const AngleValue * const(&value)(tup.value);
00071 
00072   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00073                << localIndex[1] << " " << localIndex[2] << std::endl);
00074 
00075   const Position & pos1 = p[0]->x[localIndex[0]].position;
00076   const Position & pos2 = p[1]->x[localIndex[1]].position;
00077   Vector r12 = lattice.delta(pos1,pos2);
00078   register BigReal d12inv = r12.rlength();
00079   const Position & pos3 = p[2]->x[localIndex[2]].position;
00080   Vector r32 = lattice.delta(pos3,pos2);
00081   register BigReal d32inv = r32.rlength();
00082   int normal = value->normal;
00083 
00084   BigReal cos_theta = (r12*r32)*(d12inv*d32inv);
00085   //  Make sure that the cosine value is acceptable.  With roundoff, you
00086   //  can get values like 1.0+2e-16, which makes acos puke.  So instead,
00087   //  just set these kinds of values to exactly 1.0
00088   if (cos_theta > 1.0) cos_theta = 1.0;
00089   else if (cos_theta < -1.0) cos_theta = -1.0;
00090 
00091   BigReal k = value->k * scale;
00092   BigReal theta0 = value->theta0;
00093 
00094   //  Get theta
00095   BigReal theta = acos(cos_theta);
00096 
00097   //  Compare it to the rest angle
00098   BigReal diff;
00099 
00100   if (normal == 1) {
00101     diff = theta - theta0;
00102   } else {
00103     diff = cos_theta - cos(theta0);
00104   }
00105   
00106   //  Add the energy from this angle to the total energy
00107   BigReal energy = k *diff*diff;
00108 
00109   //  Normalize vector r12 and r32
00110   //BigReal d12inv = 1. / d12;
00111   //BigReal d32inv = 1. / d32;
00112 
00113 
00114   //  Calculate constant factor 2k(theta-theta0)/sin(theta)
00115   BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
00116   if (normal != 1) {
00117     diff *= (2.0* k);
00118   } else if ( sin_theta < 1.e-6 ) {
00119     // catch case where bonds are parallel
00120     // this gets the force approximately right for theta0 of 0 or pi
00121     // and at least avoids small division for other values
00122     if ( diff < 0. ) diff = 2.0 * k;
00123     else diff = -2.0 * k;
00124   } else { 
00125     diff *= (-2.0* k) / sin_theta; 
00126   }
00127   BigReal c1 = diff * d12inv;
00128   BigReal c2 = diff * d32inv;
00129 
00130   //  Calculate the actual forces
00131   Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
00132   Force force2 = force1;
00133   Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
00134   force2 += force3;  force2 *= -1;
00135 
00136   //  Check to see if we need to do the Urey-Bradley term
00137   if (value->k_ub)
00138   {
00139         //  Non-zero k_ub value, so calculate the harmonic
00140         //  potential between the 1-3 atoms
00141 
00142   if (normal != 1) {
00143     NAMD_die("ERROR: Can't use cosAngles with Urey-Bradley angles");
00144   }
00145         BigReal k_ub = value->k_ub;
00146         BigReal r_ub = value->r_ub;
00147         Vector r13 = r12 - r32;
00148         BigReal d13 = r13.length();
00149         diff = d13- r_ub;
00150 
00151         energy += k_ub *diff*diff;
00152 
00153         diff *= -2.0*k_ub / d13;
00154         r13 *= diff;
00155 
00156         force1 += r13;
00157         force3 -= r13;
00158   }
00159 
00160   //fepb - BKR scaling of alchemical bonded terms
00161   //       NB: TI derivative is the _unscaled_ energy.
00162   if ( simParams->alchOn ) {
00163     switch ( mol->get_fep_bonded_type(atomID, 3) ) {
00164     case 1:
00165       reduction[angleEnergyIndex_ti_1] += energy;
00166       reduction[angleEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy; 
00167       energy *= bond_lambda_1;
00168       force1 *= bond_lambda_1;
00169       force2 *= bond_lambda_1;
00170       force3 *= bond_lambda_1;
00171       break;
00172     case 2:
00173       reduction[angleEnergyIndex_ti_2] += energy;
00174       reduction[angleEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
00175       energy *= bond_lambda_2;
00176       force1 *= bond_lambda_2;
00177       force2 *= bond_lambda_2;
00178       force3 *= bond_lambda_2;
00179       break;
00180     }
00181   }
00182   //fepe
00183 
00184   p[0]->f[localIndex[0]].x += force1.x;
00185   p[0]->f[localIndex[0]].y += force1.y;
00186   p[0]->f[localIndex[0]].z += force1.z;
00187 
00188   p[1]->f[localIndex[1]].x += force2.x;
00189   p[1]->f[localIndex[1]].y += force2.y;
00190   p[1]->f[localIndex[1]].z += force2.z;
00191 
00192   p[2]->f[localIndex[2]].x += force3.x;
00193   p[2]->f[localIndex[2]].y += force3.y;
00194   p[2]->f[localIndex[2]].z += force3.z;
00195 
00196   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00197 
00198   reduction[angleEnergyIndex] += energy;
00199   reduction[virialIndex_XX] += ( force1.x * r12.x + force3.x * r32.x );
00200   reduction[virialIndex_XY] += ( force1.x * r12.y + force3.x * r32.y );
00201   reduction[virialIndex_XZ] += ( force1.x * r12.z + force3.x * r32.z );
00202   reduction[virialIndex_YX] += ( force1.y * r12.x + force3.y * r32.x );
00203   reduction[virialIndex_YY] += ( force1.y * r12.y + force3.y * r32.y );
00204   reduction[virialIndex_YZ] += ( force1.y * r12.z + force3.y * r32.z );
00205   reduction[virialIndex_ZX] += ( force1.z * r12.x + force3.z * r32.x );
00206   reduction[virialIndex_ZY] += ( force1.z * r12.y + force3.z * r32.y );
00207   reduction[virialIndex_ZZ] += ( force1.z * r12.z + force3.z * r32.z );
00208 
00209   if (pressureProfileData) {
00210     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00211     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00212     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00213     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00214     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00215     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00216     pp_clamp(n1, pressureProfileSlabs);
00217     pp_clamp(n2, pressureProfileSlabs);
00218     pp_clamp(n3, pressureProfileSlabs);
00219     int p1 = p[0]->x[localIndex[0]].partition;
00220     int p2 = p[1]->x[localIndex[1]].partition;
00221     int p3 = p[2]->x[localIndex[2]].partition;
00222     int pn = pressureProfileAtomTypes;
00223     pp_reduction(pressureProfileSlabs, n1, n2, 
00224                 p1, p2, pn,
00225                 force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
00226                 pressureProfileData);
00227     pp_reduction(pressureProfileSlabs, n3, n2, 
00228                 p3, p2, pn,
00229                 force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
00230                 pressureProfileData);
00231   }
00232 
00233  }
00234 }

void AngleElem::getMoleculePointers ( Molecule *  ,
int *  ,
int32 ***  ,
Angle **   
) [static]

Definition at line 32 of file ComputeAngles.C.

References Molecule::angles, Molecule::anglesByAtom, NAMD_die(), and Molecule::numAngles.

00033 {
00034 #ifdef MEM_OPT_VERSION
00035   NAMD_die("Should not be called in AngleElem::getMoleculePointers in memory optimized version!");
00036 #else
00037   *count = mol->numAngles;
00038   *byatom = mol->anglesByAtom;
00039   *structarray = mol->angles;
00040 #endif
00041 }

void AngleElem::getParameterPointers ( Parameters *  ,
const AngleValue **   
) [static]

Definition at line 43 of file ComputeAngles.C.

References p.

00043                                                                         {
00044   *v = p->angle_array;
00045 }

static void AngleElem::getTupleInfo ( AtomSignature *  sig,
int *  count,
TupleSignature **  t 
) [inline, static]

Definition at line 28 of file ComputeAngles.h.

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

00028                                                                                  {
00029         *count = sig->angleCnt;
00030         *t = sig->angleSigs;
00031     }

int AngleElem::hash ( void   )  const [inline]

Definition at line 42 of file ComputeAngles.h.

References atomID.

00042                    { 
00043     return 0x7FFFFFFF & ((atomID[0]<<22) + (atomID[1]<<11) + (atomID[2])); 
00044   }

int AngleElem::operator< ( const AngleElem &  a  )  const [inline]

Definition at line 45 of file ComputeAngles.inl.

References atomID.

00046   {
00047     return (atomID[0] < a.atomID[0] ||
00048             (atomID[0] == a.atomID[0] &&
00049             (atomID[1] < a.atomID[1] ||
00050             (atomID[1] == a.atomID[1] &&
00051              atomID[2] < a.atomID[2]) )));
00052   }

int AngleElem::operator== ( const AngleElem &  a  )  const [inline]

Definition at line 39 of file ComputeAngles.inl.

References atomID.

00040   {
00041     return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
00042         a.atomID[2] == atomID[2]);
00043   }

void AngleElem::submitReductionData ( BigReal *  ,
SubmitReduction *   
) [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, REDUCTION_BONDED_ENERGY_TI_2, REDUCTION_VIRIAL_NORMAL, and virialIndex.

00238 {
00239   reduction->item(REDUCTION_ANGLE_ENERGY) += data[angleEnergyIndex];
00240   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[angleEnergyIndex_f];
00241   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[angleEnergyIndex_ti_1];
00242   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[angleEnergyIndex_ti_2];
00243   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00244 }

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(), and getParameterPointers().

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 on Sat Sep 23 01:17:17 2017 for NAMD by  doxygen 1.4.7