ComputeCrossterms.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeCrossterms.h"
00009 #include "Molecule.h"
00010 #include "Parameters.h"
00011 #include "Node.h"
00012 #include "ReductionMgr.h"
00013 #include "Lattice.h"
00014 #include "PressureProfile.h"
00015 #include "Debug.h"
00016 
00017 #define FASTER
00018 
00019 enum {
00020   CMAP_TABLE_DIM = 25,
00021   CMAP_SETUP_DIM = 24,
00022   CMAP_SPACING = 15,
00023   CMAP_PHI_0 = -180,
00024   CMAP_PSI_0 = -180
00025 };
00026 
00027 #define INDEX(ncols,i,j)  ((i)*ncols + (j))
00028 
00029 
00030 // static initialization
00031 int CrosstermElem::pressureProfileSlabs = 0;
00032 int CrosstermElem::pressureProfileAtomTypes = 1;
00033 BigReal CrosstermElem::pressureProfileThickness = 0;
00034 BigReal CrosstermElem::pressureProfileMin = 0;
00035 
00036 
00037 void CrosstermElem::getMoleculePointers
00038     (Molecule* mol, int* count, int32*** byatom, Crossterm** structarray)
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 }
00048 
00049 void CrosstermElem::getParameterPointers(Parameters *p, const CrosstermValue **v) {
00050   *v = p->crossterm_array;
00051 }
00052 
00053 void CrosstermElem::computeForce(CrosstermElem *tuples, int ntuple, BigReal *reduction,
00054                                 BigReal *pressureProfileData)
00055 {
00056  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00057 
00058  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00059   const CrosstermElem &tup = tuples[ituple];
00060   enum { size = 8 };
00061   const AtomID (&atomID)[size](tup.atomID);
00062   const int    (&localIndex)[size](tup.localIndex);
00063   TuplePatchElem * const(&p)[size](tup.p);
00064   const Real (&scale)(tup.scale);
00065   const CrosstermValue * const(&value)(tup.value);
00066 
00067   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00068                << localIndex[1] << " " << localIndex[2] << " " <<
00069                localIndex[3] << std::endl);
00070 
00071   // Vector r12, r23, r34;      // vector between atoms
00072   Vector A,B,C,D,E,F;           // cross products
00073   BigReal rA,rB,rC,rD,rE,rF;    // length of vectors
00074   BigReal energy=0;     // energy from the angle
00075   BigReal phi,psi;              // angle between the plans
00076   double cos_phi,cos_psi;       // cos(phi)
00077   double sin_phi,sin_psi;       // sin(phi)
00078   Vector dcosdA,dcosdD; // Derivative d(cos(phi))/dA
00079   Vector dcosdB,dcosdE; // Derivative d(cos(phi))/dB
00080   Vector dsindC,dsindF; // Derivative d(sin(phi))/dC
00081   Vector dsindB,dsindE; // Derivative d(sin(phi))/dB
00082   BigReal U,U_phi,U_psi;        // energy constants
00083   BigReal diff;         // for periodicity
00084   Force f1(0,0,0),f2(0,0,0),f3(0,0,0);  // force components
00085   Force f4(0,0,0),f5(0,0,0),f6(0,0,0);  // force components
00086 
00087   //DebugM(3, "::computeForce() -- starting with crossterm type " << crosstermType << std::endl);
00088 
00089   //  Calculate the vectors between atoms
00090   const Position & pos0 = p[0]->x[localIndex[0]].position;
00091   const Position & pos1 = p[1]->x[localIndex[1]].position;
00092   const Position & pos2 = p[2]->x[localIndex[2]].position;
00093   const Position & pos3 = p[3]->x[localIndex[3]].position;
00094   const Position & pos4 = p[4]->x[localIndex[4]].position;
00095   const Position & pos5 = p[5]->x[localIndex[5]].position;
00096   const Position & pos6 = p[6]->x[localIndex[6]].position;
00097   const Position & pos7 = p[7]->x[localIndex[7]].position;
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   const Vector r56 = lattice.delta(pos4,pos5);
00102   const Vector r67 = lattice.delta(pos5,pos6);
00103   const Vector r78 = lattice.delta(pos6,pos7);
00104 
00105   //  Calculate the cross products
00106   A = cross(r12,r23);
00107   B = cross(r23,r34);
00108   C = cross(r23,A);
00109   D = cross(r56,r67);
00110   E = cross(r67,r78);
00111   F = cross(r67,D);
00112 
00113   //  Calculate the distances
00114   rA = A.length();
00115   rB = B.length();
00116   rC = C.length();
00117   rD = D.length();
00118   rE = E.length();
00119   rF = F.length();
00120 
00121   //  Calculate the sin and cos
00122   cos_phi = A*B/(rA*rB);
00123   sin_phi = C*B/(rC*rB);
00124   cos_psi = D*E/(rD*rE);
00125   sin_psi = F*E/(rF*rE);
00126 
00127   //  Normalize B
00128   rB = 1.0/rB;
00129   B *= rB;
00130   rE = 1.0/rE;
00131   E *= rE;
00132 
00133   phi= -atan2(sin_phi,cos_phi);
00134   psi= -atan2(sin_psi,cos_psi);
00135 
00136   if (fabs(sin_phi) > 0.1) {
00137     //  Normalize A
00138     rA = 1.0/rA;
00139     A *= rA;
00140     dcosdA = rA*(cos_phi*A-B);
00141     dcosdB = rB*(cos_phi*B-A);
00142   } else {
00143     //  Normalize C
00144     rC = 1.0/rC;
00145     C *= rC;
00146     dsindC = rC*(sin_phi*C-B);
00147     dsindB = rB*(sin_phi*B-C);
00148   }
00149 
00150   if (fabs(sin_psi) > 0.1) {
00151     //  Normalize A
00152     rD = 1.0/rD;
00153     D *= rD;
00154     dcosdD = rD*(cos_psi*D-E);
00155     dcosdE = rE*(cos_psi*E-D);
00156   } else {
00157     //  Normalize C
00158     rF = 1.0/rF;
00159     F *= rF;
00160     dsindF = rF*(sin_psi*F-E);
00161     dsindE = rE*(sin_psi*E-F);
00162   }
00163 
00164     //  Calculate the energy
00165 {
00166   const double h = CMAP_SPACING * PI / 180.0;
00167   const double h_1 = 1.0 / h;
00168 #ifdef FASTER
00169   const double six_h = 6.0 * h_1;
00170 #endif
00171   const double phi_0 = CMAP_PHI_0 * PI / 180.0;
00172   const double psi_0 = CMAP_PSI_0 * PI / 180.0;
00173 
00174   enum { D = CMAP_TABLE_DIM };
00175   double xa[2], xb[2], dxa[2], dxb[2];
00176   double ya[2], yb[2], dya[2], dyb[2];
00177   double t, dx_h, dy_h;
00178 #ifdef FASTER
00179   double s1, s2, s3, s4, s5;
00180 #endif
00181   double f = 0, fx = 0, fy = 0;
00182   const CrosstermData *table = &value->c[0][0];
00183   int i, j, ilo, jlo, ij;
00184 
00185   /* distance measured in grid points between angle and smallest value */
00186   dx_h = (phi - phi_0) * h_1;
00187   dy_h = (psi - psi_0) * h_1;
00188 
00189   /* find smallest numbered grid point in stencil */
00190   ilo = (int) floor(dx_h);
00191   if ( ilo < 0 ) ilo = 0; 
00192   if ( ilo >= CMAP_SETUP_DIM ) ilo = CMAP_SETUP_DIM - 1; 
00193   jlo = (int) floor(dy_h);
00194   if ( jlo < 0 ) jlo = 0; 
00195   if ( jlo >= CMAP_SETUP_DIM ) jlo = CMAP_SETUP_DIM - 1; 
00196 
00197 #if !defined(FASTER)
00198 
00199   /* find t for x-dimension and compute xa, xb, dxa, dxb */
00200   t = dx_h - (double) ilo;
00201   xa[0] = (1 - t) * (1 - t) * (1 + 2*t);
00202   xb[0] = h * t * (1 - t) * (1 - t);
00203   dxa[0] = -6 * t * (1 - t) * h_1;
00204   dxb[0] = (1 - t) * (1 - 3*t);
00205   t--;
00206   xa[1] = (1 + t) * (1 + t) * (1 - 2*t);
00207   xb[1] = h * t * (1 + t) * (1 + t);
00208   dxa[1] = -6 * t * (1 + t) * h_1;
00209   dxb[1] = (1 + t) * (1 + 3*t);
00210 
00211   /* find t for y-dimension and compute ya, yb, dya, dyb */
00212   t = dy_h - (double) jlo;
00213   ya[0] = (1 - t) * (1 - t) * (1 + 2*t);
00214   yb[0] = h * t * (1 - t) * (1 - t);
00215   dya[0] = -6 * t * (1 - t) * h_1;
00216   dyb[0] = (1 - t) * (1 - 3*t);
00217   t--;
00218   ya[1] = (1 + t) * (1 + t) * (1 - 2*t);
00219   yb[1] = h * t * (1 + t) * (1 + t);
00220   dya[1] = -6 * t * (1 + t) * h_1;
00221   dyb[1] = (1 + t) * (1 + 3*t);
00222 
00223 #else
00224 
00225   /* find t for x-dimension and compute xa, xb, dxa, dxb */
00226   t = dx_h - (double) ilo;
00227   s1 = 1-t;
00228   s2 = 2*t;
00229   s3 = 3*t;
00230   s4 = t*s1;
00231   s5 = h*s4;
00232   xa[0] = s1*s1*(1+s2);
00233   xa[1] = t*t*(3-s2);
00234   xb[0] = s5*s1;
00235   xb[1] = -s5*t;
00236   dxa[0] = -six_h*s4;
00237   dxa[1] = -dxa[0];
00238   dxb[0] = s1*(1-s3);
00239   dxb[1] = t*(-2+s3);
00240 
00241   /* find t for y-dimension and compute ya, yb, dya, dyb */
00242   t = dy_h - (double) jlo;
00243   s1 = 1-t;
00244   s2 = 2*t;
00245   s3 = 3*t;
00246   s4 = t*s1;
00247   s5 = h*s4;
00248   ya[0] = s1*s1*(1+s2);
00249   ya[1] = t*t*(3-s2);
00250   yb[0] = s5*s1;
00251   yb[1] = -s5*t;
00252   dya[0] = -six_h*s4;
00253   dya[1] = -dya[0];
00254   dyb[0] = s1*(1-s3);
00255   dyb[1] = t*(-2+s3);
00256 
00257 #endif
00258 
00259   for (i = 0;  i < 2;  i++) {
00260     for (j = 0;  j < 2;  j++) {
00261       ij = INDEX(D,i+ilo,j+jlo);
00262 
00263 #if !defined(FASTER)
00264 
00265       f += xa[i] * ya[j] * table[ij].d00
00266         + xb[i] * ya[j] * table[ij].d10
00267         + xa[i] * yb[j] * table[ij].d01
00268         + xb[i] * yb[j] * table[ij].d11;
00269 
00270       fx += dxa[i] * ya[j] * table[ij].d00
00271         + dxb[i] * ya[j] * table[ij].d10
00272         + dxa[i] * yb[j] * table[ij].d01
00273         + dxb[i] * yb[j] * table[ij].d11;
00274 
00275       fy += xa[i] * dya[j] * table[ij].d00
00276         + xb[i] * dya[j] * table[ij].d10
00277         + xa[i] * dyb[j] * table[ij].d01
00278         + xb[i] * dyb[j] * table[ij].d11;
00279 
00280 #else
00281 
00282       s1=ya[j]*table[ij].d00;
00283       s2=yb[j]*table[ij].d01;
00284       s3=ya[j]*table[ij].d10;
00285       s4=yb[j]*table[ij].d11;
00286 
00287       f+=xa[i]*(s1+s2)+xb[i]*(s3+s4);
00288       fx+=dxa[i]*(s1+s2)+dxb[i]*(s3+s4);
00289       fy+=xa[i]*(dya[j]*table[ij].d00+dyb[j]*table[ij].d01)
00290         +xb[i]*(dya[j]*table[ij].d10+dyb[j]*table[ij].d11);
00291 
00292 #endif
00293     }
00294   }
00295 
00296   /* return accumulated values */
00297   U = f * scale;
00298   U_phi = fx * scale;
00299   U_psi = fy * scale;
00300 
00301 /*
00302 CkPrintf("crossterm %d-%d-%d-%d %d-%d-%d-%d %lf %lf %d %d %lf %lf %lf\n",
00303     atomID[0], atomID[1], atomID[2], atomID[3],
00304     atomID[4], atomID[5], atomID[6], atomID[7],
00305     phi, psi, ilo, jlo, U, U_phi, U_psi);
00306 CkPrintf("%d %d-%d-%d-%d %d-%d-%d-%d\n", CkMyPe(),
00307    p[0]->patchID, p[1]->patchID, p[2]->patchID, p[3]->patchID,
00308    p[4]->patchID, p[5]->patchID, p[6]->patchID, p[7]->patchID);
00309 */
00310 
00311 }
00312 
00313     //  Add the energy from this crossterm to the total energy
00314     energy += U;
00315 
00316     //  Next, we want to calculate the forces.  In order
00317     //  to do that, we first need to figure out whether the
00318     //  sin or cos form will be more stable.  For this,
00319     //  just look at the value of phi
00320     if (fabs(sin_phi) > 0.1)
00321     {
00322       //  use the sin version to avoid 1/cos terms
00323       U_phi = U_phi/sin_phi;
00324 
00325       f1.x += U_phi*(r23.y*dcosdA.z - r23.z*dcosdA.y);
00326       f1.y += U_phi*(r23.z*dcosdA.x - r23.x*dcosdA.z);
00327       f1.z += U_phi*(r23.x*dcosdA.y - r23.y*dcosdA.x);
00328 
00329       f3.x += U_phi*(r23.z*dcosdB.y - r23.y*dcosdB.z);
00330       f3.y += U_phi*(r23.x*dcosdB.z - r23.z*dcosdB.x);
00331       f3.z += U_phi*(r23.y*dcosdB.x - r23.x*dcosdB.y);
00332 
00333       f2.x += U_phi*(r12.z*dcosdA.y - r12.y*dcosdA.z
00334                + r34.y*dcosdB.z - r34.z*dcosdB.y);
00335       f2.y += U_phi*(r12.x*dcosdA.z - r12.z*dcosdA.x
00336                + r34.z*dcosdB.x - r34.x*dcosdB.z);
00337       f2.z += U_phi*(r12.y*dcosdA.x - r12.x*dcosdA.y
00338              + r34.x*dcosdB.y - r34.y*dcosdB.x);
00339     }
00340     else
00341     {
00342       //  This angle is closer to 0 or 180 than it is to
00343       //  90, so use the cos version to avoid 1/sin terms
00344       U_phi = -U_phi/cos_phi;
00345 
00346       f1.x += U_phi*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00347                 - r23.x*r23.y*dsindC.y
00348                 - r23.x*r23.z*dsindC.z);
00349       f1.y += U_phi*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00350                 - r23.y*r23.z*dsindC.z
00351                 - r23.y*r23.x*dsindC.x);
00352       f1.z += U_phi*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00353                 - r23.z*r23.x*dsindC.x
00354                 - r23.z*r23.y*dsindC.y);
00355 
00356       f3 += cross(U_phi,dsindB,r23);
00357 
00358       f2.x += U_phi*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00359              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00360              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00361              +dsindB.z*r34.y - dsindB.y*r34.z);
00362       f2.y += U_phi*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00363              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00364              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00365              +dsindB.x*r34.z - dsindB.z*r34.x);
00366       f2.z += U_phi*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00367              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00368              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00369              +dsindB.y*r34.x - dsindB.x*r34.y);
00370     }
00371 
00372     if (fabs(sin_psi) > 0.1)
00373     {
00374       //  use the sin version to avoid 1/cos terms
00375       U_psi = U_psi/sin_psi;
00376 
00377       f4.x += U_psi*(r67.y*dcosdD.z - r67.z*dcosdD.y);
00378       f4.y += U_psi*(r67.z*dcosdD.x - r67.x*dcosdD.z);
00379       f4.z += U_psi*(r67.x*dcosdD.y - r67.y*dcosdD.x);
00380 
00381       f6.x += U_psi*(r67.z*dcosdE.y - r67.y*dcosdE.z);
00382       f6.y += U_psi*(r67.x*dcosdE.z - r67.z*dcosdE.x);
00383       f6.z += U_psi*(r67.y*dcosdE.x - r67.x*dcosdE.y);
00384 
00385       f5.x += U_psi*(r56.z*dcosdD.y - r56.y*dcosdD.z
00386                + r78.y*dcosdE.z - r78.z*dcosdE.y);
00387       f5.y += U_psi*(r56.x*dcosdD.z - r56.z*dcosdD.x
00388                + r78.z*dcosdE.x - r78.x*dcosdE.z);
00389       f5.z += U_psi*(r56.y*dcosdD.x - r56.x*dcosdD.y
00390              + r78.x*dcosdE.y - r78.y*dcosdE.x);
00391     }
00392     else
00393     {
00394       //  This angle is closer to 0 or 180 than it is to
00395       //  90, so use the cos version to avoid 1/sin terms
00396       U_psi = -U_psi/cos_psi;
00397 
00398       f4.x += U_psi*((r67.y*r67.y + r67.z*r67.z)*dsindF.x
00399                 - r67.x*r67.y*dsindF.y
00400                 - r67.x*r67.z*dsindF.z);
00401       f4.y += U_psi*((r67.z*r67.z + r67.x*r67.x)*dsindF.y
00402                 - r67.y*r67.z*dsindF.z
00403                 - r67.y*r67.x*dsindF.x);
00404       f4.z += U_psi*((r67.x*r67.x + r67.y*r67.y)*dsindF.z
00405                 - r67.z*r67.x*dsindF.x
00406                 - r67.z*r67.y*dsindF.y);
00407 
00408       f6 += cross(U_psi,dsindE,r67);
00409 
00410       f5.x += U_psi*(-(r67.y*r56.y + r67.z*r56.z)*dsindF.x
00411              +(2.0*r67.x*r56.y - r56.x*r67.y)*dsindF.y
00412              +(2.0*r67.x*r56.z - r56.x*r67.z)*dsindF.z
00413              +dsindE.z*r78.y - dsindE.y*r78.z);
00414       f5.y += U_psi*(-(r67.z*r56.z + r67.x*r56.x)*dsindF.y
00415              +(2.0*r67.y*r56.z - r56.y*r67.z)*dsindF.z
00416              +(2.0*r67.y*r56.x - r56.y*r67.x)*dsindF.x
00417              +dsindE.x*r78.z - dsindE.z*r78.x);
00418       f5.z += U_psi*(-(r67.x*r56.x + r67.y*r56.y)*dsindF.z
00419              +(2.0*r67.z*r56.x - r56.z*r67.x)*dsindF.x
00420              +(2.0*r67.z*r56.y - r56.z*r67.y)*dsindF.y
00421              +dsindE.y*r78.x - dsindE.x*r78.y);
00422     }
00423 
00424   /* store the forces */
00425   p[0]->f[localIndex[0]] += f1;
00426   p[1]->f[localIndex[1]] += f2 - f1;
00427   p[2]->f[localIndex[2]] += f3 - f2;
00428   p[3]->f[localIndex[3]] += -f3;
00429   p[4]->f[localIndex[4]] += f4;
00430   p[5]->f[localIndex[5]] += f5 - f4;
00431   p[6]->f[localIndex[6]] += f6 - f5;
00432   p[7]->f[localIndex[7]] += -f6;
00433 
00434   if ( p[0]->af ) {
00435     p[0]->af[localIndex[0]] += f1;
00436     p[1]->af[localIndex[1]] += f2 - f1;
00437     p[2]->af[localIndex[2]] += f3 - f2;
00438     p[3]->af[localIndex[3]] += -f3;
00439     p[4]->af[localIndex[4]] += f4;
00440     p[5]->af[localIndex[5]] += f5 - f4;
00441     p[6]->af[localIndex[6]] += f6 - f5;
00442     p[7]->af[localIndex[7]] += -f6;
00443   }
00444 
00445   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00446   reduction[crosstermEnergyIndex] += energy;
00447   reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
00448   reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
00449   reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
00450   reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
00451   reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
00452   reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
00453   reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
00454   reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
00455   reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
00456 
00457   reduction[virialIndex_XX] += ( f4.x * r56.x + f5.x * r67.x + f6.x * r78.x );
00458   reduction[virialIndex_XY] += ( f4.x * r56.y + f5.x * r67.y + f6.x * r78.y );
00459   reduction[virialIndex_XZ] += ( f4.x * r56.z + f5.x * r67.z + f6.x * r78.z );
00460   reduction[virialIndex_YX] += ( f4.y * r56.x + f5.y * r67.x + f6.y * r78.x );
00461   reduction[virialIndex_YY] += ( f4.y * r56.y + f5.y * r67.y + f6.y * r78.y );
00462   reduction[virialIndex_YZ] += ( f4.y * r56.z + f5.y * r67.z + f6.y * r78.z );
00463   reduction[virialIndex_ZX] += ( f4.z * r56.x + f5.z * r67.x + f6.z * r78.x );
00464   reduction[virialIndex_ZY] += ( f4.z * r56.y + f5.z * r67.y + f6.z * r78.y );
00465   reduction[virialIndex_ZZ] += ( f4.z * r56.z + f5.z * r67.z + f6.z * r78.z );
00466 
00467   if (pressureProfileData) {
00468     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00469     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00470     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00471     BigReal z4 = p[3]->x[localIndex[3]].position.z;
00472     BigReal z5 = p[4]->x[localIndex[4]].position.z;
00473     BigReal z6 = p[5]->x[localIndex[5]].position.z;
00474     BigReal z7 = p[6]->x[localIndex[6]].position.z;
00475     BigReal z8 = p[7]->x[localIndex[7]].position.z;
00476     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00477     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00478     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00479     int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
00480     int n5 = (int)floor((z5-pressureProfileMin)/pressureProfileThickness);
00481     int n6 = (int)floor((z6-pressureProfileMin)/pressureProfileThickness);
00482     int n7 = (int)floor((z7-pressureProfileMin)/pressureProfileThickness);
00483     int n8 = (int)floor((z8-pressureProfileMin)/pressureProfileThickness);
00484     pp_clamp(n1, pressureProfileSlabs);
00485     pp_clamp(n2, pressureProfileSlabs);
00486     pp_clamp(n3, pressureProfileSlabs);
00487     pp_clamp(n4, pressureProfileSlabs);
00488     pp_clamp(n5, pressureProfileSlabs);
00489     pp_clamp(n6, pressureProfileSlabs);
00490     pp_clamp(n7, pressureProfileSlabs);
00491     pp_clamp(n8, pressureProfileSlabs);
00492     int p1 = p[0]->x[localIndex[0]].partition;
00493     int p2 = p[1]->x[localIndex[1]].partition;
00494     int p3 = p[2]->x[localIndex[2]].partition;
00495     int p4 = p[3]->x[localIndex[3]].partition;
00496     int p5 = p[4]->x[localIndex[4]].partition;
00497     int p6 = p[5]->x[localIndex[5]].partition;
00498     int p7 = p[6]->x[localIndex[6]].partition;
00499     int p8 = p[7]->x[localIndex[7]].partition;
00500     int pn = pressureProfileAtomTypes;
00501     pp_reduction(pressureProfileSlabs, n1, n2, p1, p2, pn,
00502                 f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
00503                 pressureProfileData);
00504     pp_reduction(pressureProfileSlabs, n2, n3, p2, p3, pn,
00505                 f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
00506                 pressureProfileData);
00507     pp_reduction(pressureProfileSlabs, n3, n4, p3, p4, pn,
00508                 f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
00509                 pressureProfileData);
00510     pp_reduction(pressureProfileSlabs, n5, n6, p5, p6, pn,
00511                 f4.x * r56.x, f4.y * r56.y, f4.z * r56.z,
00512                 pressureProfileData);
00513     pp_reduction(pressureProfileSlabs, n6, n7, p6, p7, pn,
00514                 f5.x * r67.x, f5.y * r67.y, f5.z * r67.z,
00515                 pressureProfileData);
00516     pp_reduction(pressureProfileSlabs, n7, n8, p7, p8, pn,
00517                 f6.x * r78.x, f6.y * r78.y, f6.z * r78.z,
00518                 pressureProfileData);
00519   }
00520 
00521  }
00522 }
00523 
00524 void CrosstermElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00525 {
00526   reduction->item(REDUCTION_CROSSTERM_ENERGY) += data[crosstermEnergyIndex];
00527   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00528   ADD_TENSOR(reduction,REDUCTION_VIRIAL_AMD_DIHE,data,virialIndex);
00529 }
00530 
00531 
00532 /******************************************************************************/
00533 
00534 static void lu_decomp_nopivot(double *m, int n);
00535 static void forward_back_sub(double *b, double *m, int n);
00536 
00537 void crossterm_setup(CrosstermData *table)
00538 {
00539   enum { D = CMAP_TABLE_DIM };
00540   enum { N = CMAP_SETUP_DIM };
00541   const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
00542   const double h_2 = h_1 * h_1;
00543   const double tr_h = 3.0 * h_1;
00544   int i, j;
00545   int ij;
00546   int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
00547   int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
00548 
00549   /* allocate spline coefficient matrix */
00550   double* const m = new double[N*N];
00551   memset(m,0,N*N*sizeof(double));
00552 
00553   /* initialize spline coefficient matrix */
00554   m[0] = 4;
00555   for (i = 1;  i < N;  i++) {
00556     m[INDEX(N,i-1,i)] = 1;
00557     m[INDEX(N,i,i-1)] = 1;
00558     m[INDEX(N,i,i)] = 4;
00559   }
00560   /* periodic boundary conditions for spline */
00561   m[INDEX(N,0,N-1)] = 1;
00562   m[INDEX(N,N-1,0)] = 1;
00563 
00564   /* compute LU-decomposition for this matrix */
00565   lu_decomp_nopivot(m, N);
00566 
00567   /* allocate vector for solving spline derivatives */
00568   double* const v = new double[N];
00569   memset(v,0,N*sizeof(double));
00570 
00571     /* march through rows of table */
00572     for (i = 0;  i < N;  i++) {
00573 
00574       /* setup RHS vector for solving spline derivatives */
00575       v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
00576       for (j = 1;  j < N;  j++) {
00577         v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
00578       }
00579 
00580       /* solve system, returned into vector */
00581       forward_back_sub(v, m, N);
00582 
00583       /* store values as derivatives wrt differenced table values */
00584       for (j = 0;  j < N;  j++) {
00585         table[INDEX(D,i,j)].d01 = v[j];
00586       }
00587       table[INDEX(D,i,N)].d01 = v[0];
00588     }
00589     for (j = 0;  j <= N;  j++) {
00590       table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
00591     }
00592 
00593     /* march through columns of table */
00594     for (j = 0;  j < N;  j++) {
00595 
00596       /* setup RHS vector for solving spline derivatives */
00597       v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
00598       for (i = 1;  i < N;  i++) {
00599         v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
00600       }
00601 
00602       /* solve system, returned into vector */
00603       forward_back_sub(v, m, N);
00604 
00605       /* store values as derivatives wrt differenced table values */
00606       for (i = 0;  i < N;  i++) {
00607         table[INDEX(D,i,j)].d10 = v[i];
00608       }
00609       table[INDEX(D,N,j)].d10 = v[0];
00610     }
00611     for (i = 0;  i <= N;  i++) {
00612       table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
00613     }
00614 
00615     /* march back through rows of table
00616      *
00617      * This is CHARMM's approach for calculating mixed partial derivatives,
00618      * by splining the first derivative values.
00619      *
00620      * Here we spline the dx values along y to calculate dxdy derivatives.
00621      *
00622      * Test cases show error with CHARMM is within 1e-5 roundoff error.
00623      */
00624     for (i = 0;  i < N;  i++) {
00625 
00626       /* setup RHS vector for solving dxy derivatives from dx */
00627       v[0] = tr_h * (table[INDEX(D,i,1)].d10 - table[INDEX(D,i,N-1)].d10);
00628       for (j = 1;  j < N;  j++) {
00629         v[j] = tr_h * (table[INDEX(D,i,j+1)].d10 - table[INDEX(D,i,j-1)].d10);
00630       }
00631 
00632       /* solve system, returned into vector */
00633       forward_back_sub(v, m, N);
00634 
00635       /* store values as dxy derivatives wrt differenced table values */
00636       for (j = 0;  j < N;  j++) {
00637         table[INDEX(D,i,j)].d11 = v[j];
00638       }
00639       table[INDEX(D,i,N)].d11 = v[0];
00640     }
00641     for (j = 0;  j <= N;  j++) {
00642       table[INDEX(D,N,j)].d11 = table[INDEX(D,0,j)].d11;
00643     }
00644 
00645   /* done with temp storage */
00646   delete [] m;
00647   delete [] v;
00648 
00649 }
00650 
00651 
00652 void lu_decomp_nopivot(double *m, int n)
00653 {
00654   double l_ik;
00655   int i, j, k;
00656 
00657   for (k = 0;  k < n-1;  k++) {
00658     for (i = k+1;  i < n;  i++) {
00659       l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
00660       for (j = k;  j < n;  j++) {
00661         m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
00662       }
00663       m[INDEX(n,i,k)] = l_ik;
00664     }
00665   }
00666 }
00667 
00668 
00669 void forward_back_sub(double *b, double *m, int n)
00670 {
00671   int i, j;
00672 
00673   /* in place forward elimination (solves Ly=b) using lower triangle of m */
00674   for (j = 0;  j < n-1;  j++) {
00675     for (i = j+1;  i < n;  i++) {
00676       b[i] -= m[INDEX(n,i,j)] * b[j];
00677     }
00678   }
00679   /* in place back substitution (solves Ux=y) using upper triangle of m */
00680   for (j = n-1;  j >= 0;  j--) {
00681     b[j] /= m[INDEX(n,j,j)];
00682     for (i = j-1;  i >= 0;  i--) {
00683       b[i] -= m[INDEX(n,i,j)] * b[j];
00684     }
00685   }
00686 }
00687 
00688 

Generated on Sat Nov 18 01:17:12 2017 for NAMD by  doxygen 1.4.7