NAMD: DihedralElem Class Reference

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


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	loadTuplesForAtom (void , AtomID, Molecule )
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