NAMD: AngleElem Class Reference

00073 {
00074   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00075                << localIndex[1] << " " << localIndex[2] << std::endl);
00076 
00077   const Position & pos1 = p[0]->x[localIndex[0]].position;
00078   const Lattice & lattice = p[0]->p->lattice;
00079   const Position & pos2 = p[1]->x[localIndex[1]].position;
00080   Vector r12 = lattice.delta(pos1,pos2);
00081   register BigReal d12inv = r12.rlength();
00082   const Position & pos3 = p[2]->x[localIndex[2]].position;
00083   Vector r32 = lattice.delta(pos3,pos2);
00084   register BigReal d32inv = r32.rlength();
00085 
00086   BigReal cos_theta = (r12*r32)*(d12inv*d32inv);
00087   //  Make sure that the cosine value is acceptable.  With roundoff, you
00088   //  can get values like 1.0+2e-16, which makes acos puke.  So instead,
00089   //  just set these kinds of values to exactly 1.0
00090   if (cos_theta > 1.0) cos_theta = 1.0;
00091   else if (cos_theta < -1.0) cos_theta = -1.0;
00092 
00093   BigReal k = value->k * scale;
00094   BigReal theta0 = value->theta0;
00095 
00096   //  Get theta
00097   BigReal theta = acos(cos_theta);
00098 
00099   //  Compare it to the rest angle
00100   BigReal diff = theta - theta0;
00101 
00102   //  Add the energy from this angle to the total energy
00103   BigReal energy = k *diff*diff;
00104 
00105   //  Normalize vector r12 and r32
00106   //BigReal d12inv = 1. / d12;
00107   //BigReal d32inv = 1. / d32;
00108 
00109   //  Calculate constant factor 2k(theta-theta0)/sin(theta)
00110   BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
00111   if ( sin_theta < 1.e-6 ) {
00112     // catch case where bonds are parallel
00113     // this gets the force approximately right for theta0 of 0 or pi
00114     // and at least avoids small division for other values
00115     if ( diff < 0. ) diff = 2.0 * k;
00116     else diff = -2.0 * k;
00117   } else { diff *= (-2.0* k) / sin_theta; }
00118   BigReal c1 = diff * d12inv;
00119   BigReal c2 = diff * d32inv;
00120 
00121   //  Calculate the actual forces
00122   Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
00123   Force force2 = force1;
00124   Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
00125   force2 += force3;  force2 *= -1;
00126 
00127   //  Check to see if we need to do the Urey-Bradley term
00128   if (value->k_ub)
00129   {
00130         //  Non-zero k_ub value, so calculate the harmonic
00131         //  potential between the 1-3 atoms
00132         BigReal k_ub = value->k_ub;
00133         BigReal r_ub = value->r_ub;
00134         Vector r13 = r12 - r32;
00135         BigReal d13 = r13.length();
00136         diff = d13- r_ub;
00137 
00138         energy += k_ub *diff*diff;
00139 
00140         diff *= -2.0*k_ub / d13;
00141         r13 *= diff;
00142 
00143         force1 += r13;
00144         force3 -= r13;
00145   }
00146   
00147   p[0]->f[localIndex[0]].x += force1.x;
00148   p[0]->f[localIndex[0]].y += force1.y;
00149   p[0]->f[localIndex[0]].z += force1.z;
00150 
00151   p[1]->f[localIndex[1]].x += force2.x;
00152   p[1]->f[localIndex[1]].y += force2.y;
00153   p[1]->f[localIndex[1]].z += force2.z;
00154 
00155   p[2]->f[localIndex[2]].x += force3.x;
00156   p[2]->f[localIndex[2]].y += force3.y;
00157   p[2]->f[localIndex[2]].z += force3.z;
00158 
00159   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00160 
00161   reduction[angleEnergyIndex] += energy;
00162   reduction[virialIndex_XX] += ( force1.x * r12.x + force3.x * r32.x );
00163   reduction[virialIndex_XY] += ( force1.x * r12.y + force3.x * r32.y );
00164   reduction[virialIndex_XZ] += ( force1.x * r12.z + force3.x * r32.z );
00165   reduction[virialIndex_YX] += ( force1.y * r12.x + force3.y * r32.x );
00166   reduction[virialIndex_YY] += ( force1.y * r12.y + force3.y * r32.y );
00167   reduction[virialIndex_YZ] += ( force1.y * r12.z + force3.y * r32.z );
00168   reduction[virialIndex_ZX] += ( force1.z * r12.x + force3.z * r32.x );
00169   reduction[virialIndex_ZY] += ( force1.z * r12.y + force3.z * r32.y );
00170   reduction[virialIndex_ZZ] += ( force1.z * r12.z + force3.z * r32.z );
00171 
00172   if (pressureProfileData) {
00173     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00174     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00175     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00176     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00177     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00178     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00179     pp_clamp(n1, pressureProfileSlabs);
00180     pp_clamp(n2, pressureProfileSlabs);
00181     pp_clamp(n3, pressureProfileSlabs);
00182     int p1 = p[0]->x[localIndex[0]].partition;
00183     int p2 = p[1]->x[localIndex[1]].partition;
00184     int p3 = p[2]->x[localIndex[2]].partition;
00185     int pn = pressureProfileAtomTypes;
00186     pp_reduction(pressureProfileSlabs, n1, n2, 
00187                 p1, p2, pn,
00188                 force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
00189                 pressureProfileData);
00190     pp_reduction(pressureProfileSlabs, n3, n2, 
00191                 p3, p2, pn,
00192                 force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
00193                 pressureProfileData);
00194   }
00195 }


Public Types
enum	{ size = 3 }
enum	{ angleEnergyIndex, virialIndex, reductionDataSize }
enum	{ reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM }
Public Member Functions
void	computeForce (BigReal , BigReal )
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
void	loadTuplesForAtom (void , AtomID, Molecule )
void	getMoleculePointers (Molecule , int , int32 *, Angle )
void	getParameterPointers (Parameters , const AngleValue *)
void	getTupleInfo (AtomSignature sig, int count, TupleSignature **t)
void	submitReductionData (BigReal , SubmitReduction )
Public Attributes
AtomID	atomID [size]
int	localIndex [size]
TuplePatchElem *	p [size]
Real	scale
const AngleValue *	value
Static Public Attributes
int	pressureProfileSlabs = 0
int	pressureProfileAtomTypes = 1
BigReal	pressureProfileThickness = 0
BigReal	pressureProfileMin = 0

AngleElem Class Reference

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Static Public Attributes

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation