NAMD: DihedralElem Class Reference

00042 {
00043   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00044                << localIndex[1] << " " << localIndex[2] << std::endl);
00045 
00046   //  Calculate the vectors between atoms
00047   const Position & pos0 = p[0]->x[localIndex[0]].position;
00048   const Lattice & lattice = p[0]->p->lattice;
00049   const Position & pos1 = p[1]->x[localIndex[1]].position;
00050   const Vector r12 = lattice.delta(pos0,pos1);
00051   const Position & pos2 = p[2]->x[localIndex[2]].position;
00052   const Vector r23 = lattice.delta(pos1,pos2);
00053   const Position & pos3 = p[3]->x[localIndex[3]].position;
00054   const Vector r34 = lattice.delta(pos2,pos3);
00055 
00056   //  Calculate the cross products and distances
00057   Vector A = cross(r12,r23);
00058   register  BigReal rAinv = A.rlength();
00059   Vector B = cross(r23,r34);
00060   register  BigReal rBinv = B.rlength();
00061   Vector C = cross(r23,A);
00062   register  BigReal rCinv = C.rlength();
00063 
00064   //  Calculate the sin and cos
00065   BigReal cos_phi = (A*B)*(rAinv*rBinv);
00066   BigReal sin_phi = (C*B)*(rCinv*rBinv);
00067 
00068   BigReal phi= -atan2(sin_phi,cos_phi);
00069 
00070   BigReal K=0;          // energy
00071   BigReal K1=0;         // force
00072 
00073   // get the dihedral information
00074   int multiplicity = value->multiplicity;
00075 
00076   //  Loop through the multiple parameter sets for this
00077   //  bond.  We will only loop more than once if this
00078   //  has multiple parameter sets from Charmm22
00079   for (int mult_num=0; mult_num<multiplicity; mult_num++)
00080   {
00081     /* get angle information */
00082     Real k = value->values[mult_num].k * scale;
00083     Real delta = value->values[mult_num].delta;
00084     int n = value->values[mult_num].n;
00085 
00086     //  Calculate the energy
00087     if (n)
00088     {
00089       //  Periodicity is greater than 0, so use cos form
00090       K += k*(1+cos(n*phi - delta));
00091       K1 += -n*k*sin(n*phi - delta);
00092     }
00093     else
00094     {
00095       //  Periodicity is 0, so just use the harmonic form
00096       BigReal diff = phi-delta;
00097       if (diff < -PI)           diff += TWOPI;
00098       else if (diff > PI)       diff -= TWOPI;
00099 
00100       K += k*diff*diff;
00101       K1 += 2.0*k*diff;
00102     }
00103   } /* for multiplicity */
00104 
00105   Force f1,f2,f3;
00106 
00107   //  Normalize B
00108   //rB = 1.0/rB;
00109   B *= rBinv;
00110 
00111     //  Next, we want to calculate the forces.  In order
00112     //  to do that, we first need to figure out whether the
00113     //  sin or cos form will be more stable.  For this,
00114     //  just look at the value of phi
00115     if (fabs(sin_phi) > 0.1)
00116     {
00117       //  use the sin version to avoid 1/cos terms
00118 
00119       //  Normalize A
00120       A *= rAinv;
00121       Vector dcosdA;
00122       Vector dcosdB;
00123 
00124       dcosdA.x = rAinv*(cos_phi*A.x-B.x);
00125       dcosdA.y = rAinv*(cos_phi*A.y-B.y);
00126       dcosdA.z = rAinv*(cos_phi*A.z-B.z);
00127             
00128       dcosdB.x = rBinv*(cos_phi*B.x-A.x);
00129       dcosdB.y = rBinv*(cos_phi*B.y-A.y);
00130       dcosdB.z = rBinv*(cos_phi*B.z-A.z);
00131 
00132       K1 = K1/sin_phi;
00133 
00134       f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
00135       f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
00136       f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
00137                                               
00138       f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
00139       f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
00140       f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
00141                                               
00142       f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z
00143                + r34.y*dcosdB.z - r34.z*dcosdB.y);
00144       f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x
00145                + r34.z*dcosdB.x - r34.x*dcosdB.z);
00146       f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y
00147                + r34.x*dcosdB.y - r34.y*dcosdB.x);
00148     }
00149     else
00150     {
00151       //  This angle is closer to 0 or 180 than it is to
00152       //  90, so use the cos version to avoid 1/sin terms
00153 
00154       //  Normalize C
00155       //      rC = 1.0/rC;
00156       C *= rCinv;
00157       
00158       Vector dsindC;
00159       Vector dsindB;
00160 
00161       dsindC.x = rCinv*(sin_phi*C.x-B.x);
00162       dsindC.y = rCinv*(sin_phi*C.y-B.y);
00163       dsindC.z = rCinv*(sin_phi*C.z-B.z);
00164 
00165       dsindB.x = rBinv*(sin_phi*B.x-C.x);
00166       dsindB.y = rBinv*(sin_phi*B.y-C.y);
00167       dsindB.z = rBinv*(sin_phi*B.z-C.z);
00168 
00169       K1 = -K1/cos_phi;
00170 
00171       f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00172                 - r23.x*r23.y*dsindC.y
00173                 - r23.x*r23.z*dsindC.z);
00174       f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00175                 - r23.y*r23.z*dsindC.z
00176                 - r23.y*r23.x*dsindC.x);
00177       f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00178                 - r23.z*r23.x*dsindC.x
00179                 - r23.z*r23.y*dsindC.y);
00180 
00181       f3 = cross(K1,dsindB,r23);
00182 
00183       f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00184              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00185              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00186              +dsindB.z*r34.y - dsindB.y*r34.z);
00187       f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00188              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00189              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00190              +dsindB.x*r34.z - dsindB.z*r34.x);
00191       f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00192              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00193              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00194              +dsindB.y*r34.x - dsindB.x*r34.y);
00195     }
00196 
00197   /* store the forces */
00198   //  p[0]->f[localIndex[0]] += f1;
00199   //  p[1]->f[localIndex[1]] += f2 - f1;
00200   //  p[2]->f[localIndex[2]] += f3 - f2;
00201   //  p[3]->f[localIndex[3]] += -f3;
00202 
00203   p[0]->f[localIndex[0]].x += f1.x;
00204   p[0]->f[localIndex[0]].y += f1.y;
00205   p[0]->f[localIndex[0]].z += f1.z;
00206 
00207   p[1]->f[localIndex[1]].x += f2.x - f1.x;
00208   p[1]->f[localIndex[1]].y += f2.y - f1.y;
00209   p[1]->f[localIndex[1]].z += f2.z - f1.z;
00210 
00211   p[2]->f[localIndex[2]].x += f3.x - f2.x;
00212   p[2]->f[localIndex[2]].y += f3.y - f2.y;
00213   p[2]->f[localIndex[2]].z += f3.z - f2.z;
00214 
00215   p[3]->f[localIndex[3]].x += -f3.x;
00216   p[3]->f[localIndex[3]].y += -f3.y;
00217   p[3]->f[localIndex[3]].z += -f3.z;  
00218 
00219     /* store the force for dihedral-only accelMD */
00220   if ( p[0]->af ) {
00221     p[0]->af[localIndex[0]].x += f1.x;
00222     p[0]->af[localIndex[0]].y += f1.y;
00223     p[0]->af[localIndex[0]].z += f1.z;
00224 
00225     p[1]->af[localIndex[1]].x += f2.x - f1.x;
00226     p[1]->af[localIndex[1]].y += f2.y - f1.y;
00227     p[1]->af[localIndex[1]].z += f2.z - f1.z;
00228 
00229     p[2]->af[localIndex[2]].x += f3.x - f2.x;
00230     p[2]->af[localIndex[2]].y += f3.y - f2.y;
00231     p[2]->af[localIndex[2]].z += f3.z - f2.z;
00232 
00233     p[3]->af[localIndex[3]].x += -f3.x;
00234     p[3]->af[localIndex[3]].y += -f3.y;
00235     p[3]->af[localIndex[3]].z += -f3.z;
00236   }
00237 
00238   DebugM(3, "::computeForce() -- ending with delta energy " << K << std::endl);
00239   reduction[dihedralEnergyIndex] += K;
00240   reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
00241   reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
00242   reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
00243   reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
00244   reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
00245   reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
00246   reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
00247   reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
00248   reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
00249 
00250   if (pressureProfileData) {
00251     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00252     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00253     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00254     BigReal z4 = p[3]->x[localIndex[3]].position.z;
00255     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00256     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00257     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00258     int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
00259     pp_clamp(n1, pressureProfileSlabs);
00260     pp_clamp(n2, pressureProfileSlabs);
00261     pp_clamp(n3, pressureProfileSlabs);
00262     pp_clamp(n4, pressureProfileSlabs);
00263     int p1 = p[0]->x[localIndex[0]].partition;
00264     int p2 = p[1]->x[localIndex[1]].partition;
00265     int p3 = p[2]->x[localIndex[2]].partition;
00266     int p4 = p[3]->x[localIndex[3]].partition;
00267     int pn = pressureProfileAtomTypes;
00268     pp_reduction(pressureProfileSlabs, n1, n2,
00269                 p1, p2, pn,
00270                 f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
00271                 pressureProfileData);
00272     pp_reduction(pressureProfileSlabs, n2, n3,
00273                 p2, p3, pn,
00274                 f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
00275                 pressureProfileData);
00276     pp_reduction(pressureProfileSlabs, n3, n4,
00277                 p3, p4, pn,
00278                 f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
00279                 pressureProfileData);
00280   }
00281 }


Public Types
enum	{ size = 4 }
enum	{ dihedralEnergyIndex, virialIndex, reductionDataSize }
enum	{ reductionChecksumLabel = REDUCTION_DIHEDRAL_CHECKSUM }
Public Member Functions
void	computeForce (BigReal , BigReal )
int	hash () const
	DihedralElem ()
	DihedralElem (AtomID atom0, const TupleSignature sig, const DihedralValue v)
	DihedralElem (const Dihedral a, const DihedralValue v)
	DihedralElem (AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3)
	~DihedralElem ()
int	operator== (const DihedralElem &a) const
int	operator< (const DihedralElem &a) const
Static Public Member Functions
void	getMoleculePointers (Molecule , int , int32 *, Dihedral )
void	getParameterPointers (Parameters , const DihedralValue *)
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 DihedralValue *	value
Static Public Attributes
int	pressureProfileSlabs = 0
int	pressureProfileAtomTypes = 1
BigReal	pressureProfileThickness = 0
BigReal	pressureProfileMin = 0

DihedralElem 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