NAMD: ImproperElem Class Reference

anonymous enum

Enumerator:

size

Definition at line 20 of file ComputeImpropers.h.

00020 { size = 4 };

anonymous enum

Enumerator:

improperEnergyIndex
improperEnergyIndex_f
improperEnergyIndex_ti_1
improperEnergyIndex_ti_2
virialIndex
reductionDataSize

Definition at line 46 of file ComputeImpropers.h.

00046        { improperEnergyIndex, improperEnergyIndex_f, improperEnergyIndex_ti_1, 
00047          improperEnergyIndex_ti_2, TENSOR(virialIndex), reductionDataSize };

anonymous enum

Enumerator:

reductionChecksumLabel

Definition at line 48 of file ComputeImpropers.h.

00048 { reductionChecksumLabel = REDUCTION_IMPROPER_CHECKSUM };

ImproperElem::ImproperElem ( ) [inline]

Definition at line 51 of file ComputeImpropers.h.

00051 { ; }

ImproperElem::ImproperElem	(	AtomID	atom0,
		const TupleSignature *	sig,
		const ImproperValue *	v
	)			`[inline]`

Definition at line 53 of file ComputeImpropers.h.

References atomID, TupleSignature::offset, TupleSignature::tupleParamType, and value.

00053                                                                                {
00054       atomID[0] = atom0;
00055       atomID[1] = atom0 + sig->offset[0];
00056       atomID[2] = atom0 + sig->offset[1];
00057       atomID[3] = atom0 + sig->offset[2];
00058       value = &v[sig->tupleParamType];
00059   }

ImproperElem::ImproperElem	(	const Improper *	a,
		const ImproperValue *	v
	)			`[inline]`

Definition at line 61 of file ComputeImpropers.h.

References improper::atom1, improper::atom2, improper::atom3, improper::atom4, atomID, improper::improper_type, and value.

00061                                                           {
00062     atomID[0] = a->atom1;
00063     atomID[1] = a->atom2;
00064     atomID[2] = a->atom3;
00065     atomID[3] = a->atom4;
00066     value = &v[a->improper_type];
00067   }

ImproperElem::ImproperElem	(	AtomID	atom0,
		AtomID	atom1,
		AtomID	atom2,
		AtomID	atom3
	)			`[inline]`

Definition at line 69 of file ComputeImpropers.h.

References atomID.

00069                                                                        {
00070     if (atom0 > atom3) {  // Swap end atoms so lowest is first!
00071       AtomID tmp = atom3; atom3 = atom0; atom0 = tmp; 
00072       tmp = atom1; atom1 = atom2; atom2 = tmp;
00073     }
00074     atomID[0] = atom0;
00075     atomID[1] = atom1;
00076     atomID[2] = atom2;
00077     atomID[3] = atom3;
00078   }

ImproperElem::~ImproperElem ( ) [inline]

Definition at line 79 of file ComputeImpropers.h.

00079 {};

void ImproperElem::computeForce	(	ImproperElem *	,
		int	,
		BigReal *	,
		BigReal *
	)			`[static]`

Definition at line 40 of file ComputeImpropers.C.

References A, atomID, B, cross(), DebugM, four_body_consts::delta, Lattice::delta(), TuplePatchElem::f, Patch::flags, Molecule::get_fep_bonded_type(), if(), improperEnergyIndex, improperEnergyIndex_f, improperEnergyIndex_ti_1, improperEnergyIndex_ti_2, four_body_consts::k, Patch::lattice, Vector::length(), localIndex, Node::molecule, ImproperValue::multiplicity, four_body_consts::n, Node::Object(), TuplePatchElem::p, p, CompAtom::partition, PI, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, scale, Node::simParameters, simParams, size, Flags::step, TWOPI, value, ImproperValue::values, Vector::x, TuplePatchElem::x, Vector::y, and Vector::z.

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

void ImproperElem::getMoleculePointers	(	Molecule *	,
		int *	,
		int32 ***	,
		Improper **
	)			`[static]`

Definition at line 25 of file ComputeImpropers.C.

References Molecule::impropers, Molecule::impropersByAtom, NAMD_die(), and Molecule::numImpropers.

00026 {
00027 #ifdef MEM_OPT_VERSION
00028   NAMD_die("Should not be called in ImproperElem::getMoleculePointers in memory optimized version!");
00029 #else
00030   *count = mol->numImpropers;
00031   *byatom = mol->impropersByAtom;
00032   *structarray = mol->impropers;
00033 #endif
00034 }

void ImproperElem::getParameterPointers	(	Parameters *	,
		const ImproperValue **
	)			`[static]`

Definition at line 36 of file ComputeImpropers.C.

References p.

00036                                                                               {
00037   *v = p->improper_array;
00038 }

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

Definition at line 29 of file ComputeImpropers.h.

References AtomSignature::improperCnt, and AtomSignature::improperSigs.

00029                                                                                  {
00030         *count = sig->improperCnt;
00031         *t = sig->improperSigs;
00032     }

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

Definition at line 43 of file ComputeImpropers.h.

References atomID.

00043                    {
00044     return 0x7FFFFFFF &((atomID[0]<<24) + (atomID[1]<<16) + (atomID[2]<<8) + atomID[3]);
00045   }

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

Definition at line 86 of file ComputeImpropers.h.

References atomID.

00086                                              {
00087     return  (atomID[0] < a.atomID[0] ||
00088             (atomID[0] == a.atomID[0] &&
00089             (atomID[1] < a.atomID[1] ||
00090             (atomID[1] == a.atomID[1] &&
00091             (atomID[2] < a.atomID[2] ||
00092             (atomID[2] == a.atomID[2] &&
00093              atomID[3] < a.atomID[3] 
00094              ))))));
00095   }

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

Definition at line 81 of file ComputeImpropers.h.

References atomID.

00081                                               {
00082     return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
00083         a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3]);
00084   }

void ImproperElem::submitReductionData	(	BigReal *	,
		SubmitReduction *
	)			`[static]`

Definition at line 303 of file ComputeImpropers.C.

References ADD_TENSOR, improperEnergyIndex, improperEnergyIndex_f, improperEnergyIndex_ti_1, improperEnergyIndex_ti_2, SubmitReduction::item(), REDUCTION_BONDED_ENERGY_F, REDUCTION_BONDED_ENERGY_TI_1, REDUCTION_BONDED_ENERGY_TI_2, REDUCTION_IMPROPER_ENERGY, REDUCTION_VIRIAL_NORMAL, and virialIndex.

00304 {
00305   reduction->item(REDUCTION_IMPROPER_ENERGY) += data[improperEnergyIndex];
00306   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[improperEnergyIndex_f];
00307   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[improperEnergyIndex_ti_1];
00308   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[improperEnergyIndex_ti_2];
00309   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00310 }

AtomID ImproperElem::atomID[size]

Definition at line 21 of file ComputeImpropers.h.

Referenced by computeForce(), hash(), ImproperElem(), operator<(), and operator==().

int ImproperElem::localIndex[size]

Definition at line 22 of file ComputeImpropers.h.

Referenced by computeForce().

TuplePatchElem* ImproperElem::p[size]

Definition at line 23 of file ComputeImpropers.h.

Referenced by computeForce(), and getParameterPointers().

int ImproperElem::pressureProfileAtomTypes = 1 [static]

Definition at line 36 of file ComputeImpropers.h.

Referenced by computeForce().

BigReal ImproperElem::pressureProfileMin = 0 [static]

Definition at line 38 of file ComputeImpropers.h.

Referenced by computeForce().

int ImproperElem::pressureProfileSlabs = 0 [static]

Definition at line 35 of file ComputeImpropers.h.

Referenced by computeForce().

BigReal ImproperElem::pressureProfileThickness = 0 [static]

Definition at line 37 of file ComputeImpropers.h.

Referenced by computeForce().

Real ImproperElem::scale

Definition at line 24 of file ComputeImpropers.h.

Referenced by computeForce().

const ImproperValue* ImproperElem::value

Definition at line 41 of file ComputeImpropers.h.

Referenced by computeForce(), and ImproperElem().

ImproperElem Class Reference

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Static Public Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Types
	size = 4
	improperEnergyIndex
	improperEnergyIndex_f
	improperEnergyIndex_ti_1
	improperEnergyIndex_ti_2
	virialIndex
	reductionDataSize
	reductionChecksumLabel = REDUCTION_IMPROPER_CHECKSUM
enum	{ size = 4 }
enum	{ improperEnergyIndex, improperEnergyIndex_f, improperEnergyIndex_ti_1, improperEnergyIndex_ti_2, virialIndex, reductionDataSize }
enum	{ reductionChecksumLabel = REDUCTION_IMPROPER_CHECKSUM }
Public Member Functions
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
static void	computeForce (ImproperElem , int, BigReal , BigReal *)
static void	getMoleculePointers (Molecule , int , int32 *, Improper )
static void	getParameterPointers (Parameters , const ImproperValue *)
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 ImproperValue *	value
Static Public Attributes
static int	pressureProfileSlabs = 0
static int	pressureProfileAtomTypes = 1
static BigReal	pressureProfileThickness = 0
static BigReal	pressureProfileMin = 0