Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

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

Generated on Fri Jul 25 04:07:16 2008 for NAMD by  doxygen 1.3.9.1