NAMD: ImproperElem Class Reference

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


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

ImproperElem 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