CrosstermElem Class Reference

#include <ComputeCrossterms.h>

List of all members.

Public Types

enum  { size = 8 }
enum  {
  crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2,
  TENSOR = (virialIndex), reductionDataSize
}
enum  { reductionChecksumLabel = REDUCTION_CROSSTERM_CHECKSUM }

Public Member Functions

int hash () const
 CrosstermElem ()
 CrosstermElem (AtomID atom0, const TupleSignature *sig, const CrosstermValue *v)
 CrosstermElem (const Crossterm *a, const CrosstermValue *v)
 CrosstermElem (AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3, AtomID atom4, AtomID atom5, AtomID atom6, AtomID atom7)
 ~CrosstermElem ()
int operator== (const CrosstermElem &a) const
int operator< (const CrosstermElem &a) const

Static Public Member Functions

static void computeForce (CrosstermElem *, int, BigReal *, BigReal *)
static void getMoleculePointers (Molecule *, int *, int32 ***, Crossterm **)
static void getParameterPointers (Parameters *, const CrosstermValue **)
static void getTupleInfo (AtomSignature *sig, int *count, TupleSignature **t)
static void submitReductionData (BigReal *, SubmitReduction *)

Public Attributes

AtomID atomID [size]
int localIndex [size]
TuplePatchElemp [size]
Real scale
const CrosstermValuevalue

Static Public Attributes

static int pressureProfileSlabs = 0
static int pressureProfileAtomTypes = 1
static BigReal pressureProfileThickness = 0
static BigReal pressureProfileMin = 0

Detailed Description

Definition at line 17 of file ComputeCrossterms.h.


Member Enumeration Documentation

anonymous enum
Enumerator:
size 

Definition at line 20 of file ComputeCrossterms.h.

00020 { size = 8 };

anonymous enum
Enumerator:
crosstermEnergyIndex 
crosstermEnergyIndex_f 
crosstermEnergyIndex_ti_1 
crosstermEnergyIndex_ti_2 
TENSOR 
reductionDataSize 

Definition at line 46 of file ComputeCrossterms.h.

anonymous enum
Enumerator:
reductionChecksumLabel 

Definition at line 49 of file ComputeCrossterms.h.


Constructor & Destructor Documentation

CrosstermElem::CrosstermElem (  )  [inline]

Definition at line 52 of file ComputeCrossterms.h.

00052 { ; }

CrosstermElem::CrosstermElem ( AtomID  atom0,
const TupleSignature sig,
const CrosstermValue v 
) [inline]

Definition at line 54 of file ComputeCrossterms.h.

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

00054                                                                                  {
00055     atomID[0] = atom0;
00056     atomID[1] = atom0 + sig->offset[0];
00057     atomID[2] = atom0 + sig->offset[1];
00058     atomID[3] = atom0 + sig->offset[2];
00059     atomID[4] = atom0 + sig->offset[3];
00060     atomID[5] = atom0 + sig->offset[4];
00061     atomID[6] = atom0 + sig->offset[5];
00062     atomID[7] = atom0 + sig->offset[6];
00063     value = &v[sig->tupleParamType];
00064 
00065   }  

CrosstermElem::CrosstermElem ( const Crossterm a,
const CrosstermValue v 
) [inline]

Definition at line 67 of file ComputeCrossterms.h.

References crossterm::atom1, crossterm::atom2, crossterm::atom3, crossterm::atom4, crossterm::atom5, crossterm::atom6, crossterm::atom7, crossterm::atom8, atomID, crossterm::crossterm_type, and value.

00067                                                              {
00068     atomID[0] = a->atom1;
00069     atomID[1] = a->atom2;
00070     atomID[2] = a->atom3;
00071     atomID[3] = a->atom4;
00072     atomID[4] = a->atom5;
00073     atomID[5] = a->atom6;
00074     atomID[6] = a->atom7;
00075     atomID[7] = a->atom8;
00076     value = &v[a->crossterm_type];
00077   }

CrosstermElem::CrosstermElem ( AtomID  atom0,
AtomID  atom1,
AtomID  atom2,
AtomID  atom3,
AtomID  atom4,
AtomID  atom5,
AtomID  atom6,
AtomID  atom7 
) [inline]

Definition at line 79 of file ComputeCrossterms.h.

References atomID.

00080                                                                         {
00081     atomID[0] = atom0;
00082     atomID[1] = atom1;
00083     atomID[2] = atom2;
00084     atomID[3] = atom3;
00085     atomID[4] = atom4;
00086     atomID[5] = atom5;
00087     atomID[6] = atom6;
00088     atomID[7] = atom7;
00089   }

CrosstermElem::~CrosstermElem (  )  [inline]

Definition at line 90 of file ComputeCrossterms.h.

00090 {};


Member Function Documentation

void CrosstermElem::computeForce ( CrosstermElem tuples,
int  ntuple,
BigReal reduction,
BigReal pressureProfileData 
) [static]

Definition at line 53 of file ComputeCrossterms.C.

References A, TuplePatchElem::af, SimParameters::alchBondDecouple, SimParameters::alchOn, atomID, B, CrosstermValue::c, CMAP_PHI_0, CMAP_PSI_0, CMAP_SETUP_DIM, CMAP_SPACING, CMAP_TABLE_DIM, cross(), crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2, CrosstermData::d00, CrosstermData::d01, CrosstermData::d10, CrosstermData::d11, DebugM, Lattice::delta(), TuplePatchElem::f, Patch::flags, Molecule::get_fep_type(), SimParameters::getBondLambda(), SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), INDEX, j, Patch::lattice, Vector::length(), localIndex, Node::molecule, Node::Object(), order, TuplePatchElem::p, p, CompAtom::partition, PI, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, scale, Node::simParameters, simParams, size, Flags::step, value, Vector::x, TuplePatchElem::x, Vector::y, and Vector::z.

00055 {
00056  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00057  //fepb BKR
00058  SimParameters *const simParams = Node::Object()->simParameters;
00059  const int step = tuples[0].p[0]->p->flags.step;
00060  const BigReal alchLambda = simParams->getCurrentLambda(step);
00061  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
00062  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
00063  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
00064  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
00065  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
00066  Molecule *const mol = Node::Object()->molecule;
00067  //fepe
00068 
00069  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00070   const CrosstermElem &tup = tuples[ituple];
00071   enum { size = 8 };
00072   const AtomID (&atomID)[size](tup.atomID);
00073   const int    (&localIndex)[size](tup.localIndex);
00074   TuplePatchElem * const(&p)[size](tup.p);
00075   const Real (&scale)(tup.scale);
00076   const CrosstermValue * const(&value)(tup.value);
00077 
00078   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00079                << localIndex[1] << " " << localIndex[2] << " " <<
00080                localIndex[3] << std::endl);
00081 
00082   // Vector r12, r23, r34;      // vector between atoms
00083   Vector A,B,C,D,E,F;           // cross products
00084   BigReal rA,rB,rC,rD,rE,rF;    // length of vectors
00085   BigReal energy=0;     // energy from the angle
00086   BigReal phi,psi;              // angle between the plans
00087   double cos_phi,cos_psi;       // cos(phi)
00088   double sin_phi,sin_psi;       // sin(phi)
00089   Vector dcosdA,dcosdD; // Derivative d(cos(phi))/dA
00090   Vector dcosdB,dcosdE; // Derivative d(cos(phi))/dB
00091   Vector dsindC,dsindF; // Derivative d(sin(phi))/dC
00092   Vector dsindB,dsindE; // Derivative d(sin(phi))/dB
00093   BigReal U,U_phi,U_psi;        // energy constants
00094   BigReal diff;         // for periodicity
00095   Force f1(0,0,0),f2(0,0,0),f3(0,0,0);  // force components
00096   Force f4(0,0,0),f5(0,0,0),f6(0,0,0);  // force components
00097 
00098   //DebugM(3, "::computeForce() -- starting with crossterm type " << crosstermType << std::endl);
00099 
00100   //  Calculate the vectors between atoms
00101   const Position & pos0 = p[0]->x[localIndex[0]].position;
00102   const Position & pos1 = p[1]->x[localIndex[1]].position;
00103   const Position & pos2 = p[2]->x[localIndex[2]].position;
00104   const Position & pos3 = p[3]->x[localIndex[3]].position;
00105   const Position & pos4 = p[4]->x[localIndex[4]].position;
00106   const Position & pos5 = p[5]->x[localIndex[5]].position;
00107   const Position & pos6 = p[6]->x[localIndex[6]].position;
00108   const Position & pos7 = p[7]->x[localIndex[7]].position;
00109   const Vector r12 = lattice.delta(pos0,pos1);
00110   const Vector r23 = lattice.delta(pos1,pos2);
00111   const Vector r34 = lattice.delta(pos2,pos3);
00112   const Vector r56 = lattice.delta(pos4,pos5);
00113   const Vector r67 = lattice.delta(pos5,pos6);
00114   const Vector r78 = lattice.delta(pos6,pos7);
00115 
00116   //  Calculate the cross products
00117   A = cross(r12,r23);
00118   B = cross(r23,r34);
00119   C = cross(r23,A);
00120   D = cross(r56,r67);
00121   E = cross(r67,r78);
00122   F = cross(r67,D);
00123 
00124   //  Calculate the distances
00125   rA = A.length();
00126   rB = B.length();
00127   rC = C.length();
00128   rD = D.length();
00129   rE = E.length();
00130   rF = F.length();
00131 
00132   //  Calculate the sin and cos
00133   cos_phi = A*B/(rA*rB);
00134   sin_phi = C*B/(rC*rB);
00135   cos_psi = D*E/(rD*rE);
00136   sin_psi = F*E/(rF*rE);
00137 
00138   //  Normalize B
00139   rB = 1.0/rB;
00140   B *= rB;
00141   rE = 1.0/rE;
00142   E *= rE;
00143 
00144   phi= -atan2(sin_phi,cos_phi);
00145   psi= -atan2(sin_psi,cos_psi);
00146 
00147   if (fabs(sin_phi) > 0.1) {
00148     //  Normalize A
00149     rA = 1.0/rA;
00150     A *= rA;
00151     dcosdA = rA*(cos_phi*A-B);
00152     dcosdB = rB*(cos_phi*B-A);
00153   } else {
00154     //  Normalize C
00155     rC = 1.0/rC;
00156     C *= rC;
00157     dsindC = rC*(sin_phi*C-B);
00158     dsindB = rB*(sin_phi*B-C);
00159   }
00160 
00161   if (fabs(sin_psi) > 0.1) {
00162     //  Normalize A
00163     rD = 1.0/rD;
00164     D *= rD;
00165     dcosdD = rD*(cos_psi*D-E);
00166     dcosdE = rE*(cos_psi*E-D);
00167   } else {
00168     //  Normalize C
00169     rF = 1.0/rF;
00170     F *= rF;
00171     dsindF = rF*(sin_psi*F-E);
00172     dsindE = rE*(sin_psi*E-F);
00173   }
00174 
00175     //  Calculate the energy
00176 {
00177   const double h = CMAP_SPACING * PI / 180.0;
00178   const double h_1 = 1.0 / h;
00179 #ifdef FASTER
00180   const double six_h = 6.0 * h_1;
00181 #endif
00182   const double phi_0 = CMAP_PHI_0 * PI / 180.0;
00183   const double psi_0 = CMAP_PSI_0 * PI / 180.0;
00184 
00185   enum { D = CMAP_TABLE_DIM };
00186   double xa[2], xb[2], dxa[2], dxb[2];
00187   double ya[2], yb[2], dya[2], dyb[2];
00188   double t, dx_h, dy_h;
00189 #ifdef FASTER
00190   double s1, s2, s3, s4, s5;
00191 #endif
00192   double f = 0, fx = 0, fy = 0;
00193   const CrosstermData *table = &value->c[0][0];
00194   int i, j, ilo, jlo, ij;
00195 
00196   /* distance measured in grid points between angle and smallest value */
00197   dx_h = (phi - phi_0) * h_1;
00198   dy_h = (psi - psi_0) * h_1;
00199 
00200   /* find smallest numbered grid point in stencil */
00201   ilo = (int) floor(dx_h);
00202   if ( ilo < 0 ) ilo = 0; 
00203   if ( ilo >= CMAP_SETUP_DIM ) ilo = CMAP_SETUP_DIM - 1; 
00204   jlo = (int) floor(dy_h);
00205   if ( jlo < 0 ) jlo = 0; 
00206   if ( jlo >= CMAP_SETUP_DIM ) jlo = CMAP_SETUP_DIM - 1; 
00207 
00208 #if !defined(FASTER)
00209 
00210   /* find t for x-dimension and compute xa, xb, dxa, dxb */
00211   t = dx_h - (double) ilo;
00212   xa[0] = (1 - t) * (1 - t) * (1 + 2*t);
00213   xb[0] = h * t * (1 - t) * (1 - t);
00214   dxa[0] = -6 * t * (1 - t) * h_1;
00215   dxb[0] = (1 - t) * (1 - 3*t);
00216   t--;
00217   xa[1] = (1 + t) * (1 + t) * (1 - 2*t);
00218   xb[1] = h * t * (1 + t) * (1 + t);
00219   dxa[1] = -6 * t * (1 + t) * h_1;
00220   dxb[1] = (1 + t) * (1 + 3*t);
00221 
00222   /* find t for y-dimension and compute ya, yb, dya, dyb */
00223   t = dy_h - (double) jlo;
00224   ya[0] = (1 - t) * (1 - t) * (1 + 2*t);
00225   yb[0] = h * t * (1 - t) * (1 - t);
00226   dya[0] = -6 * t * (1 - t) * h_1;
00227   dyb[0] = (1 - t) * (1 - 3*t);
00228   t--;
00229   ya[1] = (1 + t) * (1 + t) * (1 - 2*t);
00230   yb[1] = h * t * (1 + t) * (1 + t);
00231   dya[1] = -6 * t * (1 + t) * h_1;
00232   dyb[1] = (1 + t) * (1 + 3*t);
00233 
00234 #else
00235 
00236   /* find t for x-dimension and compute xa, xb, dxa, dxb */
00237   t = dx_h - (double) ilo;
00238   s1 = 1-t;
00239   s2 = 2*t;
00240   s3 = 3*t;
00241   s4 = t*s1;
00242   s5 = h*s4;
00243   xa[0] = s1*s1*(1+s2);
00244   xa[1] = t*t*(3-s2);
00245   xb[0] = s5*s1;
00246   xb[1] = -s5*t;
00247   dxa[0] = -six_h*s4;
00248   dxa[1] = -dxa[0];
00249   dxb[0] = s1*(1-s3);
00250   dxb[1] = t*(-2+s3);
00251 
00252   /* find t for y-dimension and compute ya, yb, dya, dyb */
00253   t = dy_h - (double) jlo;
00254   s1 = 1-t;
00255   s2 = 2*t;
00256   s3 = 3*t;
00257   s4 = t*s1;
00258   s5 = h*s4;
00259   ya[0] = s1*s1*(1+s2);
00260   ya[1] = t*t*(3-s2);
00261   yb[0] = s5*s1;
00262   yb[1] = -s5*t;
00263   dya[0] = -six_h*s4;
00264   dya[1] = -dya[0];
00265   dyb[0] = s1*(1-s3);
00266   dyb[1] = t*(-2+s3);
00267 
00268 #endif
00269 
00270   for (i = 0;  i < 2;  i++) {
00271     for (j = 0;  j < 2;  j++) {
00272       ij = INDEX(D,i+ilo,j+jlo);
00273 
00274 #if !defined(FASTER)
00275 
00276       f += xa[i] * ya[j] * table[ij].d00
00277         + xb[i] * ya[j] * table[ij].d10
00278         + xa[i] * yb[j] * table[ij].d01
00279         + xb[i] * yb[j] * table[ij].d11;
00280 
00281       fx += dxa[i] * ya[j] * table[ij].d00
00282         + dxb[i] * ya[j] * table[ij].d10
00283         + dxa[i] * yb[j] * table[ij].d01
00284         + dxb[i] * yb[j] * table[ij].d11;
00285 
00286       fy += xa[i] * dya[j] * table[ij].d00
00287         + xb[i] * dya[j] * table[ij].d10
00288         + xa[i] * dyb[j] * table[ij].d01
00289         + xb[i] * dyb[j] * table[ij].d11;
00290 
00291 #else
00292 
00293       s1=ya[j]*table[ij].d00;
00294       s2=yb[j]*table[ij].d01;
00295       s3=ya[j]*table[ij].d10;
00296       s4=yb[j]*table[ij].d11;
00297 
00298       f+=xa[i]*(s1+s2)+xb[i]*(s3+s4);
00299       fx+=dxa[i]*(s1+s2)+dxb[i]*(s3+s4);
00300       fy+=xa[i]*(dya[j]*table[ij].d00+dyb[j]*table[ij].d01)
00301         +xb[i]*(dya[j]*table[ij].d10+dyb[j]*table[ij].d11);
00302 
00303 #endif
00304     }
00305   }
00306 
00307   /* return accumulated values */
00308   U = f * scale;
00309   U_phi = fx * scale;
00310   U_psi = fy * scale;
00311 
00312 /*
00313 CkPrintf("crossterm %d-%d-%d-%d %d-%d-%d-%d %lf %lf %d %d %lf %lf %lf\n",
00314     atomID[0], atomID[1], atomID[2], atomID[3],
00315     atomID[4], atomID[5], atomID[6], atomID[7],
00316     phi, psi, ilo, jlo, U, U_phi, U_psi);
00317 CkPrintf("%d %d-%d-%d-%d %d-%d-%d-%d\n", CkMyPe(),
00318    p[0]->patchID, p[1]->patchID, p[2]->patchID, p[3]->patchID,
00319    p[4]->patchID, p[5]->patchID, p[6]->patchID, p[7]->patchID);
00320 */
00321 
00322 }
00323 
00324     //  Add the energy from this crossterm to the total energy
00325     energy += U;
00326     
00327     //  Next, we want to calculate the forces.  In order
00328     //  to do that, we first need to figure out whether the
00329     //  sin or cos form will be more stable.  For this,
00330     //  just look at the value of phi
00331     if (fabs(sin_phi) > 0.1)
00332     {
00333       //  use the sin version to avoid 1/cos terms
00334       U_phi = U_phi/sin_phi;
00335 
00336       f1.x += U_phi*(r23.y*dcosdA.z - r23.z*dcosdA.y);
00337       f1.y += U_phi*(r23.z*dcosdA.x - r23.x*dcosdA.z);
00338       f1.z += U_phi*(r23.x*dcosdA.y - r23.y*dcosdA.x);
00339 
00340       f3.x += U_phi*(r23.z*dcosdB.y - r23.y*dcosdB.z);
00341       f3.y += U_phi*(r23.x*dcosdB.z - r23.z*dcosdB.x);
00342       f3.z += U_phi*(r23.y*dcosdB.x - r23.x*dcosdB.y);
00343 
00344       f2.x += U_phi*(r12.z*dcosdA.y - r12.y*dcosdA.z
00345                + r34.y*dcosdB.z - r34.z*dcosdB.y);
00346       f2.y += U_phi*(r12.x*dcosdA.z - r12.z*dcosdA.x
00347                + r34.z*dcosdB.x - r34.x*dcosdB.z);
00348       f2.z += U_phi*(r12.y*dcosdA.x - r12.x*dcosdA.y
00349              + r34.x*dcosdB.y - r34.y*dcosdB.x);
00350     }
00351     else
00352     {
00353       //  This angle is closer to 0 or 180 than it is to
00354       //  90, so use the cos version to avoid 1/sin terms
00355       U_phi = -U_phi/cos_phi;
00356 
00357       f1.x += U_phi*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00358                 - r23.x*r23.y*dsindC.y
00359                 - r23.x*r23.z*dsindC.z);
00360       f1.y += U_phi*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00361                 - r23.y*r23.z*dsindC.z
00362                 - r23.y*r23.x*dsindC.x);
00363       f1.z += U_phi*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00364                 - r23.z*r23.x*dsindC.x
00365                 - r23.z*r23.y*dsindC.y);
00366 
00367       f3 += cross(U_phi,dsindB,r23);
00368 
00369       f2.x += U_phi*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00370              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00371              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00372              +dsindB.z*r34.y - dsindB.y*r34.z);
00373       f2.y += U_phi*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00374              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00375              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00376              +dsindB.x*r34.z - dsindB.z*r34.x);
00377       f2.z += U_phi*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00378              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00379              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00380              +dsindB.y*r34.x - dsindB.x*r34.y);
00381     }
00382 
00383     if (fabs(sin_psi) > 0.1)
00384     {
00385       //  use the sin version to avoid 1/cos terms
00386       U_psi = U_psi/sin_psi;
00387 
00388       f4.x += U_psi*(r67.y*dcosdD.z - r67.z*dcosdD.y);
00389       f4.y += U_psi*(r67.z*dcosdD.x - r67.x*dcosdD.z);
00390       f4.z += U_psi*(r67.x*dcosdD.y - r67.y*dcosdD.x);
00391 
00392       f6.x += U_psi*(r67.z*dcosdE.y - r67.y*dcosdE.z);
00393       f6.y += U_psi*(r67.x*dcosdE.z - r67.z*dcosdE.x);
00394       f6.z += U_psi*(r67.y*dcosdE.x - r67.x*dcosdE.y);
00395 
00396       f5.x += U_psi*(r56.z*dcosdD.y - r56.y*dcosdD.z
00397                + r78.y*dcosdE.z - r78.z*dcosdE.y);
00398       f5.y += U_psi*(r56.x*dcosdD.z - r56.z*dcosdD.x
00399                + r78.z*dcosdE.x - r78.x*dcosdE.z);
00400       f5.z += U_psi*(r56.y*dcosdD.x - r56.x*dcosdD.y
00401              + r78.x*dcosdE.y - r78.y*dcosdE.x);
00402     }
00403     else
00404     {
00405       //  This angle is closer to 0 or 180 than it is to
00406       //  90, so use the cos version to avoid 1/sin terms
00407       U_psi = -U_psi/cos_psi;
00408 
00409       f4.x += U_psi*((r67.y*r67.y + r67.z*r67.z)*dsindF.x
00410                 - r67.x*r67.y*dsindF.y
00411                 - r67.x*r67.z*dsindF.z);
00412       f4.y += U_psi*((r67.z*r67.z + r67.x*r67.x)*dsindF.y
00413                 - r67.y*r67.z*dsindF.z
00414                 - r67.y*r67.x*dsindF.x);
00415       f4.z += U_psi*((r67.x*r67.x + r67.y*r67.y)*dsindF.z
00416                 - r67.z*r67.x*dsindF.x
00417                 - r67.z*r67.y*dsindF.y);
00418 
00419       f6 += cross(U_psi,dsindE,r67);
00420 
00421       f5.x += U_psi*(-(r67.y*r56.y + r67.z*r56.z)*dsindF.x
00422              +(2.0*r67.x*r56.y - r56.x*r67.y)*dsindF.y
00423              +(2.0*r67.x*r56.z - r56.x*r67.z)*dsindF.z
00424              +dsindE.z*r78.y - dsindE.y*r78.z);
00425       f5.y += U_psi*(-(r67.z*r56.z + r67.x*r56.x)*dsindF.y
00426              +(2.0*r67.y*r56.z - r56.y*r67.z)*dsindF.z
00427              +(2.0*r67.y*r56.x - r56.y*r67.x)*dsindF.x
00428              +dsindE.x*r78.z - dsindE.z*r78.x);
00429       f5.z += U_psi*(-(r67.x*r56.x + r67.y*r56.y)*dsindF.z
00430              +(2.0*r67.z*r56.x - r56.z*r67.x)*dsindF.x
00431              +(2.0*r67.z*r56.y - r56.z*r67.y)*dsindF.y
00432              +dsindE.y*r78.x - dsindE.x*r78.y);
00433     }
00434 
00435     //fepb - BKR scaling of alchemical bonded terms
00436     //       NB: TI derivative is the _unscaled_ energy.
00437     if ( simParams->alchOn ) {
00438       int typeSum1, typeSum2;
00439       typeSum1 = typeSum2 = 0;
00440       for (int i=0; i < 4; ++i ) {
00441         typeSum1 += (mol->get_fep_type(atomID[i]) == 2 ? -1 :\
00442             mol->get_fep_type(atomID[i]));
00443         typeSum2 += (mol->get_fep_type(atomID[i+4]) == 2 ? -1 :\
00444             mol->get_fep_type(atomID[i+4]));
00445       }
00446       int order = (simParams->alchBondDecouple ? 5 : 4);
00447       if ( (0 < typeSum1 && typeSum1 < order) ||
00448            (0 < typeSum2 && typeSum2 < order) ) {
00449         reduction[crosstermEnergyIndex_ti_1] += energy;
00450         reduction[crosstermEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
00451         energy *= bond_lambda_1;
00452         f1 *= bond_lambda_1;
00453         f2 *= bond_lambda_1;
00454         f3 *= bond_lambda_1;
00455         f4 *= bond_lambda_1;
00456         f5 *= bond_lambda_1;
00457         f6 *= bond_lambda_1;
00458       } else if ( (0 > typeSum1 && typeSum1 > -order) ||
00459                   (0 > typeSum2 && typeSum2 > -order) ) {
00460         reduction[crosstermEnergyIndex_ti_2] += energy;
00461         reduction[crosstermEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
00462         energy *= bond_lambda_2;
00463         f1 *= bond_lambda_2;
00464         f2 *= bond_lambda_2;
00465         f3 *= bond_lambda_2;
00466         f4 *= bond_lambda_2;
00467         f5 *= bond_lambda_2;
00468         f6 *= bond_lambda_2; 
00469       }
00470     }
00471     //fepe
00472 
00473   /* store the forces */
00474   p[0]->f[localIndex[0]] += f1;
00475   p[1]->f[localIndex[1]] += f2 - f1;
00476   p[2]->f[localIndex[2]] += f3 - f2;
00477   p[3]->f[localIndex[3]] += -f3;
00478   p[4]->f[localIndex[4]] += f4;
00479   p[5]->f[localIndex[5]] += f5 - f4;
00480   p[6]->f[localIndex[6]] += f6 - f5;
00481   p[7]->f[localIndex[7]] += -f6;
00482 
00483   if ( p[0]->af ) {
00484     p[0]->af[localIndex[0]] += f1;
00485     p[1]->af[localIndex[1]] += f2 - f1;
00486     p[2]->af[localIndex[2]] += f3 - f2;
00487     p[3]->af[localIndex[3]] += -f3;
00488     p[4]->af[localIndex[4]] += f4;
00489     p[5]->af[localIndex[5]] += f5 - f4;
00490     p[6]->af[localIndex[6]] += f6 - f5;
00491     p[7]->af[localIndex[7]] += -f6;
00492   }
00493 
00494   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00495   reduction[crosstermEnergyIndex] += energy;
00496   reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
00497   reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
00498   reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
00499   reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
00500   reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
00501   reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
00502   reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
00503   reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
00504   reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
00505 
00506   reduction[virialIndex_XX] += ( f4.x * r56.x + f5.x * r67.x + f6.x * r78.x );
00507   reduction[virialIndex_XY] += ( f4.x * r56.y + f5.x * r67.y + f6.x * r78.y );
00508   reduction[virialIndex_XZ] += ( f4.x * r56.z + f5.x * r67.z + f6.x * r78.z );
00509   reduction[virialIndex_YX] += ( f4.y * r56.x + f5.y * r67.x + f6.y * r78.x );
00510   reduction[virialIndex_YY] += ( f4.y * r56.y + f5.y * r67.y + f6.y * r78.y );
00511   reduction[virialIndex_YZ] += ( f4.y * r56.z + f5.y * r67.z + f6.y * r78.z );
00512   reduction[virialIndex_ZX] += ( f4.z * r56.x + f5.z * r67.x + f6.z * r78.x );
00513   reduction[virialIndex_ZY] += ( f4.z * r56.y + f5.z * r67.y + f6.z * r78.y );
00514   reduction[virialIndex_ZZ] += ( f4.z * r56.z + f5.z * r67.z + f6.z * r78.z );
00515 
00516   if (pressureProfileData) {
00517     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00518     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00519     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00520     BigReal z4 = p[3]->x[localIndex[3]].position.z;
00521     BigReal z5 = p[4]->x[localIndex[4]].position.z;
00522     BigReal z6 = p[5]->x[localIndex[5]].position.z;
00523     BigReal z7 = p[6]->x[localIndex[6]].position.z;
00524     BigReal z8 = p[7]->x[localIndex[7]].position.z;
00525     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00526     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00527     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00528     int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
00529     int n5 = (int)floor((z5-pressureProfileMin)/pressureProfileThickness);
00530     int n6 = (int)floor((z6-pressureProfileMin)/pressureProfileThickness);
00531     int n7 = (int)floor((z7-pressureProfileMin)/pressureProfileThickness);
00532     int n8 = (int)floor((z8-pressureProfileMin)/pressureProfileThickness);
00533     pp_clamp(n1, pressureProfileSlabs);
00534     pp_clamp(n2, pressureProfileSlabs);
00535     pp_clamp(n3, pressureProfileSlabs);
00536     pp_clamp(n4, pressureProfileSlabs);
00537     pp_clamp(n5, pressureProfileSlabs);
00538     pp_clamp(n6, pressureProfileSlabs);
00539     pp_clamp(n7, pressureProfileSlabs);
00540     pp_clamp(n8, pressureProfileSlabs);
00541     int p1 = p[0]->x[localIndex[0]].partition;
00542     int p2 = p[1]->x[localIndex[1]].partition;
00543     int p3 = p[2]->x[localIndex[2]].partition;
00544     int p4 = p[3]->x[localIndex[3]].partition;
00545     int p5 = p[4]->x[localIndex[4]].partition;
00546     int p6 = p[5]->x[localIndex[5]].partition;
00547     int p7 = p[6]->x[localIndex[6]].partition;
00548     int p8 = p[7]->x[localIndex[7]].partition;
00549     int pn = pressureProfileAtomTypes;
00550     pp_reduction(pressureProfileSlabs, n1, n2, p1, p2, pn,
00551                 f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
00552                 pressureProfileData);
00553     pp_reduction(pressureProfileSlabs, n2, n3, p2, p3, pn,
00554                 f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
00555                 pressureProfileData);
00556     pp_reduction(pressureProfileSlabs, n3, n4, p3, p4, pn,
00557                 f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
00558                 pressureProfileData);
00559     pp_reduction(pressureProfileSlabs, n5, n6, p5, p6, pn,
00560                 f4.x * r56.x, f4.y * r56.y, f4.z * r56.z,
00561                 pressureProfileData);
00562     pp_reduction(pressureProfileSlabs, n6, n7, p6, p7, pn,
00563                 f5.x * r67.x, f5.y * r67.y, f5.z * r67.z,
00564                 pressureProfileData);
00565     pp_reduction(pressureProfileSlabs, n7, n8, p7, p8, pn,
00566                 f6.x * r78.x, f6.y * r78.y, f6.z * r78.z,
00567                 pressureProfileData);
00568   }
00569 
00570  }
00571 }

void CrosstermElem::getMoleculePointers ( Molecule mol,
int *  count,
int32 ***  byatom,
Crossterm **  structarray 
) [static]

Definition at line 38 of file ComputeCrossterms.C.

References NAMD_die(), and Molecule::numCrossterms.

00039 {
00040 #ifdef MEM_OPT_VERSION
00041   NAMD_die("Should not be called in CrosstermElem::getMoleculePointers in memory optimized version!");
00042 #else
00043   *count = mol->numCrossterms;
00044   *byatom = mol->crosstermsByAtom;
00045   *structarray = mol->crossterms;
00046 #endif
00047 }

void CrosstermElem::getParameterPointers ( Parameters p,
const CrosstermValue **  v 
) [static]

Definition at line 49 of file ComputeCrossterms.C.

References Parameters::crossterm_array.

00049                                                                                 {
00050   *v = p->crossterm_array;
00051 }

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

Definition at line 29 of file ComputeCrossterms.h.

References AtomSignature::crosstermCnt, and AtomSignature::crosstermSigs.

00029                                                                                  {
00030         *count = sig->crosstermCnt;
00031         *t = sig->crosstermSigs;
00032     }

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

Definition at line 43 of file ComputeCrossterms.h.

References atomID.

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

int CrosstermElem::operator< ( const CrosstermElem a  )  const [inline]

Definition at line 99 of file ComputeCrossterms.h.

References atomID.

00099                                               {
00100     return  (atomID[0] < a.atomID[0] ||
00101             (atomID[0] == a.atomID[0] &&
00102             (atomID[1] < a.atomID[1] ||
00103             (atomID[1] == a.atomID[1] &&
00104             (atomID[2] < a.atomID[2] ||
00105             (atomID[2] == a.atomID[2] &&
00106             (atomID[3] < a.atomID[3] ||
00107             (atomID[3] == a.atomID[3] &&
00108             (atomID[4] < a.atomID[4] ||
00109             (atomID[4] == a.atomID[4] &&
00110             (atomID[5] < a.atomID[5] ||
00111             (atomID[5] == a.atomID[5] &&
00112             (atomID[6] < a.atomID[6] ||
00113             (atomID[6] == a.atomID[6] &&
00114              atomID[7] < a.atomID[7] 
00115              ))))))))))))));
00116   }

int CrosstermElem::operator== ( const CrosstermElem a  )  const [inline]

Definition at line 92 of file ComputeCrossterms.h.

References atomID.

00092                                                {
00093     return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
00094         a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3] &&
00095         a.atomID[4] == atomID[4] && a.atomID[5] == atomID[5] &&
00096         a.atomID[6] == atomID[6] && a.atomID[7] == atomID[7]);
00097   }

void CrosstermElem::submitReductionData ( BigReal data,
SubmitReduction reduction 
) [static]

Definition at line 573 of file ComputeCrossterms.C.

References ADD_TENSOR, crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2, SubmitReduction::item(), REDUCTION_BONDED_ENERGY_F, REDUCTION_BONDED_ENERGY_TI_1, REDUCTION_BONDED_ENERGY_TI_2, and REDUCTION_CROSSTERM_ENERGY.

00574 {
00575   reduction->item(REDUCTION_CROSSTERM_ENERGY) += data[crosstermEnergyIndex];
00576   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[crosstermEnergyIndex_f];
00577   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[crosstermEnergyIndex_ti_1];
00578   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[crosstermEnergyIndex_ti_2];
00579   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00580   ADD_TENSOR(reduction,REDUCTION_VIRIAL_AMD_DIHE,data,virialIndex);
00581 }


Member Data Documentation

Definition at line 21 of file ComputeCrossterms.h.

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

Definition at line 22 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 23 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 36 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 38 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 35 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 37 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 24 of file ComputeCrossterms.h.

Referenced by computeForce().

Definition at line 41 of file ComputeCrossterms.h.

Referenced by computeForce(), and CrosstermElem().


The documentation for this class was generated from the following files:

Generated on 11 Nov 2019 for NAMD by  doxygen 1.6.1