ComputeMsmSerial.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "Node.h"
00009 #include "PatchMap.h"
00010 #include "PatchMap.inl"
00011 #include "AtomMap.h"
00012 #include "ComputeMsmSerial.h"
00013 #include "ComputeMsmSerialMgr.decl.h"
00014 #include "PatchMgr.h"
00015 #include "Molecule.h"
00016 #include "ReductionMgr.h"
00017 #include "ComputeMgr.h"
00018 #include "ComputeMgr.decl.h"
00019 // #define DEBUGM
00020 #define MIN_DEBUG_LEVEL 3
00021 #include "Debug.h"
00022 #include "SimParameters.h"
00023 #include "WorkDistrib.h"
00024 #include "varsizemsg.h"
00025 #include "msm.h"
00026 #include <stdlib.h>
00027 #include <stdio.h>
00028 #include <errno.h>
00029 
00030 
00031 struct ComputeMsmSerialAtom {
00032   Position position;
00033   float charge;
00034   int id;
00035 };
00036 
00037 typedef Force MsmSerialForce;
00038 
00039 class MsmSerialCoordMsg : public CMessage_MsmSerialCoordMsg {
00040 public:
00041   int sourceNode;
00042   int numAtoms;
00043   Lattice lattice;
00044   ComputeMsmSerialAtom *coord;
00045 };
00046 
00047 class MsmSerialForceMsg : public CMessage_MsmSerialForceMsg {
00048 public:
00049   BigReal energy;
00050   BigReal virial[3][3];
00051   MsmSerialForce *force;
00052 };
00053 
00054 class ComputeMsmSerialMgr : public CBase_ComputeMsmSerialMgr {
00055 public:
00056   ComputeMsmSerialMgr();
00057   ~ComputeMsmSerialMgr();
00058 
00059   void setCompute(ComputeMsmSerial *c) { msmCompute = c; }
00060 
00061   void recvCoord(MsmSerialCoordMsg *);
00062   void recvForce(MsmSerialForceMsg *);
00063 
00064 private:
00065   CProxy_ComputeMsmSerialMgr msmProxy;
00066   ComputeMsmSerial *msmCompute;
00067 
00068   int numSources;
00069   int numArrived;
00070   MsmSerialCoordMsg **coordMsgs;
00071   int numAtoms;
00072   ComputeMsmSerialAtom *coord;
00073   MsmSerialForce *force;
00074   MsmSerialForceMsg *oldmsg;
00075 
00076   NL_Msm *msmsolver;
00077   double *msmcoord;
00078   double *msmforce;
00079 };
00080 
00081 ComputeMsmSerialMgr::ComputeMsmSerialMgr() :
00082   msmProxy(thisgroup), msmCompute(0), numSources(0), numArrived(0),
00083   coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0),
00084   msmsolver(0), msmcoord(0), msmforce(0)
00085 {
00086   CkpvAccess(BOCclass_group).computeMsmSerialMgr = thisgroup;
00087 }
00088 
00089 ComputeMsmSerialMgr::~ComputeMsmSerialMgr()
00090 {
00091   for (int i=0;  i < numSources;  i++)  delete coordMsgs[i];
00092   delete [] coordMsgs;
00093   delete [] coord;
00094   delete [] force;
00095   delete oldmsg;
00096   if (msmsolver) NL_msm_destroy(msmsolver);
00097   if (msmcoord) delete[] msmcoord;
00098   if (msmforce) delete[] msmforce;
00099 }
00100 
00101 ComputeMsmSerial::ComputeMsmSerial(ComputeID c) :
00102   ComputeHomePatches(c)
00103 {
00104   CProxy_ComputeMsmSerialMgr::ckLocalBranch(
00105         CkpvAccess(BOCclass_group).computeMsmSerialMgr)->setCompute(this);
00106   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00107 }
00108 
00109 ComputeMsmSerial::~ComputeMsmSerial()
00110 {
00111 }
00112 
00113 void ComputeMsmSerial::doWork()
00114 {
00115   ResizeArrayIter<PatchElem> ap(patchList);
00116 
00117   // Skip computations if nothing to do.
00118   if ( ! patchList[0].p->flags.doFullElectrostatics )
00119   {
00120     for (ap = ap.begin(); ap != ap.end(); ap++) {
00121       CompAtom *x = (*ap).positionBox->open();
00122       Results *r = (*ap).forceBox->open();
00123       (*ap).positionBox->close(&x);
00124       (*ap).forceBox->close(&r);
00125     }
00126     reduction->submit();
00127     return;
00128   }
00129 
00130   // allocate message
00131   int numLocalAtoms = 0;
00132   for (ap = ap.begin(); ap != ap.end(); ap++) {
00133     numLocalAtoms += (*ap).p->getNumAtoms();
00134   }
00135 
00136   MsmSerialCoordMsg *msg = new (numLocalAtoms, 0) MsmSerialCoordMsg;
00137   msg->sourceNode = CkMyPe();
00138   msg->numAtoms = numLocalAtoms;
00139   msg->lattice = patchList[0].p->flags.lattice;
00140   ComputeMsmSerialAtom *data_ptr = msg->coord;
00141 
00142   // get positions
00143   for (ap = ap.begin();  ap != ap.end();  ap++) {
00144     CompAtom *x = (*ap).positionBox->open();
00145     CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
00146     if ( patchList[0].p->flags.doMolly ) {
00147       (*ap).positionBox->close(&x);
00148       x = (*ap).avgPositionBox->open();
00149     }
00150     int numAtoms = (*ap).p->getNumAtoms();
00151 
00152     for(int i=0;  i < numAtoms;  i++)
00153     {
00154       data_ptr->position = x[i].position;
00155       data_ptr->charge = x[i].charge;
00156       data_ptr->id = xExt[i].id;
00157       ++data_ptr;
00158     }
00159 
00160     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00161     else { (*ap).positionBox->close(&x); }
00162   }
00163 
00164   CProxy_ComputeMsmSerialMgr msmProxy(
00165       CkpvAccess(BOCclass_group).computeMsmSerialMgr);
00166   msmProxy[0].recvCoord(msg);
00167 }
00168 
00169 
00170 static void rescale_nonperiodic_cell(
00171     Vector &u, Vector &v, Vector &w, Vector &c,
00172     Vector &ru, Vector &rv, Vector &rw,
00173     int isperiodic, int numatoms, const ComputeMsmSerialAtom *coord,
00174     double padding, double gridspacing) {
00175   const double NL_ZERO_TOLERANCE = 1e-4;
00176   double xlen = 1, inv_xlen = 1;  /* don't rescale periodic directions */
00177   double ylen = 1, inv_ylen = 1;
00178   double zlen = 1, inv_zlen = 1;
00179   double ulen_1 = u.rlength();
00180   double vlen_1 = v.rlength();
00181   double wlen_1 = w.rlength();
00182   double rpadx = padding * ulen_1;  /* padding distance in recip space */
00183   double rpady = padding * vlen_1;
00184   double rpadz = padding * wlen_1;
00185   double rhx = gridspacing * ulen_1;  /* grid spacing in recip space */
00186   double rhy = gridspacing * vlen_1;
00187   double rhz = gridspacing * wlen_1;
00188   Vector s, d;
00189   Vector rc = 0;  // don't move center along periodic directions
00190   double xmin, xmax, ymin, ymax, zmin, zmax;
00191   int is_periodic_x = (isperiodic & NL_MSM_PERIODIC_VEC1);
00192   int is_periodic_y = (isperiodic & NL_MSM_PERIODIC_VEC2);
00193   int is_periodic_z = (isperiodic & NL_MSM_PERIODIC_VEC3);
00194   int i;
00195 
00196   /* affine linear transformation of coordinates to reciprocal space,
00197    * where periodic vector directions map into [-0.5, 0.5) */
00198   //printf("*** center=%.4f %.4f %.4f\n", c.x, c.y, c.z);
00199   d = coord[0].position - c;
00200   s.x = ru * d;  // dot product
00201   s.y = rv * d;
00202   s.z = rw * d;
00203   xmin = xmax = s.x;
00204   ymin = ymax = s.y;
00205   zmin = zmax = s.z;
00206   for (i = 1;  i < numatoms;  i++) {
00207     d = coord[i].position - c;
00208     s.x = ru * d;  // dot product
00209     s.y = rv * d;
00210     s.z = rw * d;
00211     if      (s.x < xmin)  xmin = s.x;
00212     else if (s.x > xmax)  xmax = s.x;
00213     if      (s.y < ymin)  ymin = s.y;
00214     else if (s.y > ymax)  ymax = s.y;
00215     if      (s.z < zmin)  zmin = s.z;
00216     else if (s.z > zmax)  zmax = s.z;
00217   }
00218 #if 0
00219   printf("*** xmin=%.4f  xmax=%.4f\n", xmin, xmax);
00220   printf("*** ymin=%.4f  ymax=%.4f\n", ymin, ymax);
00221   printf("*** zmin=%.4f  zmax=%.4f\n", zmin, zmax);
00222 #endif
00223 
00224   if ( ! is_periodic_x) {
00225     xmax += rpadx;  /* pad the edges */
00226     xmin -= rpadx;
00227     if (rhx > 0) {  /* restrict center to rhx lattice points */
00228       double mupper = ceil(xmax / (2*rhx));
00229       double mlower = floor(xmin / (2*rhx));
00230       xmax = 2*rhx*mupper;
00231       xmin = 2*rhx*mlower;
00232     }
00233     rc.x = 0.5*(xmin + xmax);
00234     xlen = xmax - xmin;
00235     if (xlen < NL_ZERO_TOLERANCE) {
00236       xlen = 1;  /* leave scaling unchanged */
00237     }
00238     else {
00239       inv_xlen = 1/xlen;
00240     }
00241   }
00242 
00243   if ( ! is_periodic_y) {
00244     ymax += rpady;  /* pad the edges */
00245     ymin -= rpady;
00246     if (rhy > 0) {  /* restrict center to rhy lattice points */
00247       double mupper = ceil(ymax / (2*rhy));
00248       double mlower = floor(ymin / (2*rhy));
00249       ymax = 2*rhy*mupper;
00250       ymin = 2*rhy*mlower;
00251     }
00252     rc.y = 0.5*(ymin + ymax);
00253     ylen = ymax - ymin;
00254     if (ylen < NL_ZERO_TOLERANCE) {
00255       ylen = 1;  /* leave scaling unchanged */
00256     }
00257     else {
00258       inv_ylen = 1/ylen;
00259     }
00260   }
00261 
00262   if ( ! is_periodic_z) {
00263     zmax += rpadz;  /* pad the edges */
00264     zmin -= rpadz;
00265     if (rhz > 0) {  /* restrict center to rhz lattice points */
00266       double mupper = ceil(zmax / (2*rhz));
00267       double mlower = floor(zmin / (2*rhz));
00268       zmax = 2*rhz*mupper;
00269       zmin = 2*rhz*mlower;
00270     }
00271     rc.z = 0.5*(zmin + zmax);
00272     zlen = zmax - zmin;
00273     if (zlen < NL_ZERO_TOLERANCE) {
00274       zlen = 1;  /* leave scaling unchanged */
00275     }
00276     else {
00277       inv_zlen = 1/zlen;
00278     }
00279   }
00280 
00281 #if 0
00282   printf("xmin=%g  xmax=%g\n", xmin, xmax);
00283   printf("ymin=%g  ymax=%g\n", ymin, ymax);
00284   printf("zmin=%g  zmax=%g\n", zmin, zmax);
00285   printf("xlen=%g  ylen=%g  zlen=%g\n", xlen, ylen, zlen);
00286   printf("rc_x=%g  rc_y=%g  rc_z=%g\n", rc.x, rc.y, rc.z);
00287 #endif
00288 
00289   /* transform new center back to real space, then rescale basis vectors */
00290   c.x = (u.x*rc.x + v.x*rc.y + w.x*rc.z) + c.x;
00291   c.y = (u.y*rc.x + v.y*rc.y + w.y*rc.z) + c.y;
00292   c.z = (u.z*rc.x + v.z*rc.y + w.z*rc.z) + c.z;
00293 
00294 #if 0
00295   printf("c_x=%g  c_y=%g  c_z=%g\n", c.x, c.y, c.z);
00296 #endif
00297 
00298   u *= xlen;
00299   v *= ylen;
00300   w *= zlen;
00301 
00302   ru *= inv_xlen;
00303   rv *= inv_ylen;
00304   rw *= inv_zlen;
00305 }
00306 
00307 
00308 void ComputeMsmSerialMgr::recvCoord(MsmSerialCoordMsg *msg) {
00309   if ( ! numSources ) {
00310     numSources = (PatchMap::Object())->numNodesWithPatches();
00311     coordMsgs = new MsmSerialCoordMsg*[numSources];
00312     for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
00313     numArrived = 0;
00314     numAtoms = Node::Object()->molecule->numAtoms;
00315     coord = new ComputeMsmSerialAtom[numAtoms];
00316     force = new MsmSerialForce[numAtoms];
00317   }
00318 
00319   int i;
00320   for ( i=0; i < msg->numAtoms; ++i ) {
00321     coord[msg->coord[i].id] = msg->coord[i];
00322   }
00323 
00324   coordMsgs[numArrived] = msg;
00325   ++numArrived;
00326 
00327   if ( numArrived < numSources ) return;
00328   numArrived = 0;
00329 
00330   // ALL DATA ARRIVED --- CALCULATE FORCES
00331   Lattice lattice = msg->lattice;
00332   SimParameters *simParams = Node::Object()->simParameters;
00333 
00334   double energy = 0;
00335   double virial[3][3];
00336 
00337   int rc = 0;  // return code
00338 
00339   if ( ! msmsolver ) {
00340     //
00341     // setup MSM solver
00342     //
00343     msmsolver = NL_msm_create();
00344     if ( ! msmsolver ) NAMD_die("unable to create MSM solver");
00345     double dielectric = simParams->dielectric;
00346     double cutoff = simParams->cutoff;
00347     double gridspacing = simParams->MSMGridSpacing;
00348     double padding = simParams->MSMPadding;
00349     int approx = simParams->MSMApprox;
00350     int split = simParams->MSMSplit;
00351     int nlevels = simParams->MSMLevels;
00352     int msmflags = 0;
00353     msmflags |= (lattice.a_p() ? NL_MSM_PERIODIC_VEC1 : 0);
00354     msmflags |= (lattice.b_p() ? NL_MSM_PERIODIC_VEC2 : 0);
00355     msmflags |= (lattice.c_p() ? NL_MSM_PERIODIC_VEC3 : 0);
00356     msmflags |= NL_MSM_COMPUTE_LONG_RANGE;  // compute only long-range part
00357     //msmflags |= NL_MSM_COMPUTE_ALL;
00358     //printf("msmflags = %x\n", msmflags);
00359     rc = NL_msm_configure(msmsolver, gridspacing, approx, split, nlevels);
00360     if (rc) NAMD_die("unable to configure MSM solver");
00361     Vector u=lattice.a(), v=lattice.b(), w=lattice.c(), c=lattice.origin();
00362     Vector ru=lattice.a_r(), rv=lattice.b_r(), rw=lattice.c_r();
00363     if ((msmflags & NL_MSM_PERIODIC_ALL) != NL_MSM_PERIODIC_ALL) {
00364       // called only if there is some non-periodic boundary
00365       int isperiodic = (msmflags & NL_MSM_PERIODIC_ALL);
00366       //printf("calling rescale\n");
00367       rescale_nonperiodic_cell(u, v, w, c, ru, rv, rw,
00368           isperiodic, numAtoms, coord, padding, gridspacing);
00369     }
00370     double vec1[3], vec2[3], vec3[3], center[3];
00371     vec1[0] = u.x;
00372     vec1[1] = u.y;
00373     vec1[2] = u.z;
00374     vec2[0] = v.x;
00375     vec2[1] = v.y;
00376     vec2[2] = v.z;
00377     vec3[0] = w.x;
00378     vec3[1] = w.y;
00379     vec3[2] = w.z;
00380     center[0] = c.x;
00381     center[1] = c.y;
00382     center[2] = c.z;
00383 #if 0
00384     printf("dielectric = %g\n", dielectric);
00385     printf("vec1 = %g %g %g\n", vec1[0], vec1[1], vec1[2]);
00386     printf("vec2 = %g %g %g\n", vec2[0], vec2[1], vec2[2]);
00387     printf("vec3 = %g %g %g\n", vec3[0], vec3[1], vec3[2]);
00388     printf("center = %g %g %g\n", center[0], center[1], center[2]);
00389     printf("cutoff = %g\n", cutoff);
00390     printf("numatoms = %d\n", numAtoms);
00391 #endif
00392     rc = NL_msm_setup(msmsolver, cutoff, vec1, vec2, vec3, center, msmflags);
00393     if (rc) NAMD_die("unable to set up MSM solver");
00394     msmcoord = new double[4*numAtoms];
00395     msmforce = new double[3*numAtoms];
00396     if (msmcoord==0 || msmforce==0) NAMD_die("can't allocate MSM atom buffers");
00397     // scale charges - these won't change
00398     double celec = sqrt(COULOMB / simParams->dielectric);
00399     for (i = 0;  i < numAtoms;  i++) {
00400       msmcoord[4*i+3] = celec * coord[i].charge;
00401     }
00402   }
00403 
00404   // evaluate long-range MSM forces
00405   for (i = 0;  i < numAtoms;  i++) {
00406     msmcoord[4*i  ] = coord[i].position.x;
00407     msmcoord[4*i+1] = coord[i].position.y;
00408     msmcoord[4*i+2] = coord[i].position.z;
00409   }
00410   for (i = 0;  i < numAtoms;  i++) {
00411     msmforce[3*i  ] = 0;
00412     msmforce[3*i+1] = 0;
00413     msmforce[3*i+2] = 0;
00414   }
00415   rc = NL_msm_compute_force(msmsolver, msmforce, &energy, msmcoord, numAtoms);
00416   if (rc) NAMD_die("error evaluating MSM forces");
00417   for (i = 0;  i < numAtoms;  i++) {
00418     force[i].x = msmforce[3*i  ];
00419     force[i].y = msmforce[3*i+1];
00420     force[i].z = msmforce[3*i+2];
00421   }
00422   // MSM does not yet calculate virial
00423   for (int k=0;  k < 3;  k++) {
00424     for (int l=0;  l < 3;  l++) {
00425       virial[k][l] = 0;
00426     }
00427   }
00428 
00429   // distribute forces
00430   for (int j=0;  j < numSources;  j++) {
00431     MsmSerialCoordMsg *cmsg = coordMsgs[j];
00432     coordMsgs[j] = 0;
00433     MsmSerialForceMsg *fmsg = new (cmsg->numAtoms, 0) MsmSerialForceMsg;
00434 
00435     for (int i=0;  i < cmsg->numAtoms;  i++) {
00436       fmsg->force[i] = force[cmsg->coord[i].id];
00437     }
00438 
00439     if ( ! j ) {  // set virial and energy only for first message
00440       fmsg->energy = energy;
00441       for (int k=0;  k < 3;  k++) {
00442         for (int l=0;  l < 3;  l++) {
00443           fmsg->virial[k][l] = virial[k][l];
00444         }
00445       }
00446     }
00447     else {  // set other messages to zero, add into reduction only once
00448       fmsg->energy = 0;
00449       for (int k=0;  k < 3;  k++) {
00450         for (int l=0;  l < 3;  l++) {
00451           fmsg->virial[k][l] = 0;
00452         }
00453       }
00454     }
00455 
00456     msmProxy[cmsg->sourceNode].recvForce(fmsg);
00457     delete cmsg;
00458   }
00459 }
00460 
00461 void ComputeMsmSerialMgr::recvForce(MsmSerialForceMsg *msg) {
00462   msmCompute->saveResults(msg);
00463   delete oldmsg;
00464   oldmsg = msg;
00465 }
00466 
00467 void ComputeMsmSerial::saveResults(MsmSerialForceMsg *msg)
00468 {
00469   ResizeArrayIter<PatchElem> ap(patchList);
00470 
00471   MsmSerialForce *results_ptr = msg->force;
00472 
00473   // add in forces
00474   for (ap = ap.begin(); ap != ap.end(); ap++) {
00475     Results *r = (*ap).forceBox->open();
00476     Force *f = r->f[Results::slow];
00477     int numAtoms = (*ap).p->getNumAtoms();
00478 
00479     for(int i=0; i<numAtoms; ++i) {
00480       f[i].x += results_ptr->x;
00481       f[i].y += results_ptr->y;
00482       f[i].z += results_ptr->z;
00483       ++results_ptr;
00484     }
00485   
00486     (*ap).forceBox->close(&r);
00487   }
00488 
00489   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += msg->energy;
00490   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += msg->virial[0][0];
00491   reduction->item(REDUCTION_VIRIAL_SLOW_XY) += msg->virial[0][1];
00492   reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += msg->virial[0][2];
00493   reduction->item(REDUCTION_VIRIAL_SLOW_YX) += msg->virial[1][0];
00494   reduction->item(REDUCTION_VIRIAL_SLOW_YY) += msg->virial[1][1];
00495   reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += msg->virial[1][2];
00496   reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += msg->virial[2][0];
00497   reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += msg->virial[2][1];
00498   reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->virial[2][2];
00499   reduction->submit();
00500 }
00501 
00502 #include "ComputeMsmSerialMgr.def.h"
00503 

Generated on Tue Nov 21 01:17:11 2017 for NAMD by  doxygen 1.4.7