CrosstermElem Class Reference

#include <ComputeCrossterms.h>

List of all members.

Public Types

 size = 8
 crosstermEnergyIndex
 crosstermEnergyIndex_f
 crosstermEnergyIndex_ti_1
 crosstermEnergyIndex_ti_2
 virialIndex
 reductionDataSize
 reductionChecksumLabel = REDUCTION_CROSSTERM_CHECKSUM
enum  { size = 8 }
enum  {
  crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2,
  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 
virialIndex 
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 ,
int  ,
BigReal ,
BigReal  
) [static]

Definition at line 53 of file ComputeCrossterms.C.

References A, TuplePatchElem::af, 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, f, Patch::flags, Molecule::get_fep_type(), 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->alchLambda2;
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 ,
int *  ,
int32 ***  ,
Crossterm **   
) [static]

Definition at line 38 of file ComputeCrossterms.C.

References Molecule::crossterms, Molecule::crosstermsByAtom, 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 ,
const CrosstermValue **   
) [static]

Definition at line 49 of file ComputeCrossterms.C.

References p.

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 ,
SubmitReduction  
) [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, REDUCTION_CROSSTERM_ENERGY, REDUCTION_VIRIAL_AMD_DIHE, REDUCTION_VIRIAL_NORMAL, and virialIndex.


Member Data Documentation

AtomID CrosstermElem::atomID[size]

Definition at line 21 of file ComputeCrossterms.h.

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

int CrosstermElem::localIndex[size]

Definition at line 22 of file ComputeCrossterms.h.

Referenced by computeForce().

TuplePatchElem* CrosstermElem::p[size]

Definition at line 23 of file ComputeCrossterms.h.

Referenced by computeForce(), and getParameterPointers().

int CrosstermElem::pressureProfileAtomTypes = 1 [static]

Definition at line 36 of file ComputeCrossterms.h.

Referenced by computeForce().

BigReal CrosstermElem::pressureProfileMin = 0 [static]

Definition at line 38 of file ComputeCrossterms.h.

Referenced by computeForce().

int CrosstermElem::pressureProfileSlabs = 0 [static]

Definition at line 35 of file ComputeCrossterms.h.

Referenced by computeForce().

BigReal CrosstermElem::pressureProfileThickness = 0 [static]

Definition at line 37 of file ComputeCrossterms.h.

Referenced by computeForce().

Real CrosstermElem::scale

Definition at line 24 of file ComputeCrossterms.h.

Referenced by computeForce().

const CrosstermValue* CrosstermElem::value

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 Thu Jun 21 01:17:20 2018 for NAMD by  doxygen 1.4.7