ComputeFmmSerial.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 "ComputeFmmSerial.h"
00013 #include "ComputeFmmSerialMgr.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 <stdlib.h>
00026 #include <stdio.h>
00027 #include <errno.h>
00028 
00029 
00030 #ifdef FMM_SOLVER
00031 
00032 // Calculate FMM:
00033 //
00034 // ufmmlap library from J. Huang - http://fastmultipole.org/Main/FMMSuite/
00035 //   (ufmmlap = Uniform FMM Laplace Solver)
00036 //
00037 // The positions must translated and scaled to be within the unit
00038 // cube centered at the origin, i.e., [-1/2,1/2]^3.  The potentials 
00039 // and forces must then also be scaled.
00040 //
00041 extern "C" void fmmlap_uni_(
00042     int *natoms,        // input: number of particles
00043     double zat[][3],    // input: particle positions, length natoms
00044     double charge[],    // input: particle charges, length natoms
00045     double pot[],       // output: particle potentials, length natoms
00046     double field[][3],  // output: particle forces, length natoms
00047     int *nlev,          // input: number of oct-tree levels
00048     int *ier            // output: error code
00049     );
00050 
00051 #endif // FMM_SOLVER
00052 
00053 
00054 struct ComputeFmmSerialAtom {
00055   Position position;
00056   float charge;
00057   int id;
00058 };
00059 
00060 typedef Force FmmSerialForce;
00061 
00062 class FmmSerialCoordMsg : public CMessage_FmmSerialCoordMsg {
00063 public:
00064   int sourceNode;
00065   int numAtoms;
00066   Lattice lattice;
00067   ComputeFmmSerialAtom *coord;
00068 };
00069 
00070 class FmmSerialForceMsg : public CMessage_FmmSerialForceMsg {
00071 public:
00072   BigReal energy;
00073   BigReal virial[3][3];
00074   FmmSerialForce *force;
00075 };
00076 
00077 class ComputeFmmSerialMgr : public CBase_ComputeFmmSerialMgr {
00078 public:
00079   ComputeFmmSerialMgr();
00080   ~ComputeFmmSerialMgr();
00081 
00082   void setCompute(ComputeFmmSerial *c) { fmmCompute = c; }
00083 
00084   void recvCoord(FmmSerialCoordMsg *);
00085   void recvForce(FmmSerialForceMsg *);
00086 
00087 private:
00088   CProxy_ComputeFmmSerialMgr fmmProxy;
00089   ComputeFmmSerial *fmmCompute;
00090 
00091   int numSources;
00092   int numArrived;
00093   FmmSerialCoordMsg **coordMsgs;
00094   int numAtoms;
00095   ComputeFmmSerialAtom *coord;
00096   FmmSerialForce *force;
00097   FmmSerialForceMsg *oldmsg;
00098 
00099   int fmmsolver;
00100   int nlevels;
00101   double scaling;
00102   Vector center;
00103   double *chargebuffer;
00104   double *epotbuffer;
00105   double (*posbuffer)[3];
00106   double (*forcebuffer)[3];
00107 };
00108 
00109 ComputeFmmSerialMgr::ComputeFmmSerialMgr() :
00110   fmmProxy(thisgroup), fmmCompute(0), numSources(0), numArrived(0),
00111   coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0),
00112   fmmsolver(0), nlevels(0), scaling(0), center(0),
00113   chargebuffer(0), epotbuffer(0), posbuffer(0), forcebuffer(0)
00114 {
00115   CkpvAccess(BOCclass_group).computeFmmSerialMgr = thisgroup;
00116 }
00117 
00118 ComputeFmmSerialMgr::~ComputeFmmSerialMgr()
00119 {
00120   for (int i=0;  i < numSources;  i++)  delete coordMsgs[i];
00121   delete [] coordMsgs;
00122   delete [] coord;
00123   delete [] force;
00124   delete oldmsg;
00125   if (chargebuffer) delete[] chargebuffer;
00126   if (epotbuffer)   delete[] epotbuffer;
00127   if (posbuffer)    delete[] posbuffer;
00128   if (forcebuffer)  delete[] forcebuffer;
00129 }
00130 
00131 ComputeFmmSerial::ComputeFmmSerial(ComputeID c) :
00132   ComputeHomePatches(c)
00133 {
00134   CProxy_ComputeFmmSerialMgr::ckLocalBranch(
00135         CkpvAccess(BOCclass_group).computeFmmSerialMgr)->setCompute(this);
00136   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00137 }
00138 
00139 ComputeFmmSerial::~ComputeFmmSerial()
00140 {
00141 }
00142 
00143 void ComputeFmmSerial::doWork()
00144 {
00145   ResizeArrayIter<PatchElem> ap(patchList);
00146 
00147   // Skip computations if nothing to do.
00148   if ( ! patchList[0].p->flags.doFullElectrostatics )
00149   {
00150     for (ap = ap.begin(); ap != ap.end(); ap++) {
00151       CompAtom *x = (*ap).positionBox->open();
00152       Results *r = (*ap).forceBox->open();
00153       (*ap).positionBox->close(&x);
00154       (*ap).forceBox->close(&r);
00155     }
00156     reduction->submit();
00157     return;
00158   }
00159 
00160   // allocate message
00161   int numLocalAtoms = 0;
00162   for (ap = ap.begin(); ap != ap.end(); ap++) {
00163     numLocalAtoms += (*ap).p->getNumAtoms();
00164   }
00165 
00166   FmmSerialCoordMsg *msg = new (numLocalAtoms, 0) FmmSerialCoordMsg;
00167   msg->sourceNode = CkMyPe();
00168   msg->numAtoms = numLocalAtoms;
00169   msg->lattice = patchList[0].p->flags.lattice;
00170   ComputeFmmSerialAtom *data_ptr = msg->coord;
00171 
00172   // get positions
00173   for (ap = ap.begin();  ap != ap.end();  ap++) {
00174     CompAtom *x = (*ap).positionBox->open();
00175     CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
00176     if ( patchList[0].p->flags.doMolly ) {
00177       (*ap).positionBox->close(&x);
00178       x = (*ap).avgPositionBox->open();
00179     }
00180     int numAtoms = (*ap).p->getNumAtoms();
00181 
00182     for(int i=0;  i < numAtoms;  i++)
00183     {
00184       data_ptr->position = x[i].position;
00185       data_ptr->charge = x[i].charge;
00186       data_ptr->id = xExt[i].id;
00187       ++data_ptr;
00188     }
00189 
00190     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00191     else { (*ap).positionBox->close(&x); }
00192   }
00193 
00194   CProxy_ComputeFmmSerialMgr fmmProxy(
00195       CkpvAccess(BOCclass_group).computeFmmSerialMgr);
00196   fmmProxy[0].recvCoord(msg);
00197 }
00198 
00199 
00200 void ComputeFmmSerialMgr::recvCoord(FmmSerialCoordMsg *msg) {
00201   int i;
00202   if ( ! numSources ) {
00203     numSources = (PatchMap::Object())->numNodesWithPatches();
00204     coordMsgs = new FmmSerialCoordMsg*[numSources];
00205     for ( i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
00206     numArrived = 0;
00207     numAtoms = Node::Object()->molecule->numAtoms;
00208     coord = new ComputeFmmSerialAtom[numAtoms];
00209     force = new FmmSerialForce[numAtoms];
00210   }
00211 
00212   for ( i=0; i < msg->numAtoms; ++i ) {
00213     coord[msg->coord[i].id] = msg->coord[i];
00214   }
00215 
00216   coordMsgs[numArrived] = msg;
00217   ++numArrived;
00218 
00219   if ( numArrived < numSources ) return;
00220   numArrived = 0;
00221 
00222   // ALL DATA ARRIVED --- CALCULATE FORCES
00223   Lattice lattice = msg->lattice;
00224   SimParameters *simParams = Node::Object()->simParameters;
00225 
00226   double energy = 0;
00227   double virial[3][3] = { 0 };  // FMM does not calculate virial
00228 
00229   int rc = 0;  // return code
00230 
00231   if ( ! fmmsolver ) {
00232     //
00233     // setup FMM solver
00234     //
00235 
00236     // check boundary conditions
00237     if (lattice.a_p() || lattice.b_p() || lattice.c_p()) {
00238       NAMD_die("FMM solver requires non-periodic boundaries");
00239     }
00240 
00241     // setup number of levels
00242     if (simParams->FMMLevels > 0) {
00243       // explicitly set in config file
00244       nlevels = simParams->FMMLevels;
00245     }
00246     else {
00247       // otherwise estimate number of levels as round(log_8(numAtoms))
00248       nlevels = (int) floor(log((double)numAtoms) / log(8.) + 0.5);
00249     }
00250 
00251     // find bounding cube length
00252     if (numAtoms <= 0) {
00253       NAMD_die("setting up FMM with no atoms");
00254     }
00255     Vector min, max;
00256     min = max = coord[0].position;
00257     for (i = 1;  i < numAtoms;  i++) {
00258       Vector r = coord[i].position;
00259       if      (r.x < min.x)  min.x = r.x;
00260       else if (r.x > max.x)  max.x = r.x;
00261       if      (r.y < min.y)  min.y = r.y;
00262       else if (r.y > max.y)  max.y = r.y;
00263       if      (r.z < min.z)  min.z = r.z;
00264       else if (r.z > max.z)  max.z = r.z;
00265     }
00266     double padding = simParams->FMMPadding;
00267     if (padding <= 0) padding = 0.01;  // pad by at least a delta margin
00268     min -= padding;
00269     max += padding;
00270     double len = max.x - min.x;
00271     if (len < max.y - min.y)  len = max.y - min.y;
00272     if (len < max.z - min.z)  len = max.z - min.z;
00273     scaling = 1.0 / len;  // scale coordinates by length of the cube
00274     center = 0.5*(min + max);  // center of cube
00275     iout << iINFO << "FMM scaling length set to " << len << " A\n" << endi;
00276     iout << iINFO << "FMM center set to " << center << "\n" << endi;
00277 
00278     // allocate buffer space
00279     chargebuffer = new double[numAtoms];     // double *
00280     epotbuffer   = new double[numAtoms];     // double *
00281     posbuffer    = new double[numAtoms][3];  // double(*)[3]
00282     forcebuffer  = new double[numAtoms][3];  // double(*)[3]
00283     if (chargebuffer == 0 || epotbuffer == 0 ||
00284         posbuffer == 0 || forcebuffer == 0) {
00285       NAMD_die("can't allocate buffer space for FMM data");
00286     }
00287 
00288     // scale the charges - these won't change
00289     double celec = sqrt(COULOMB / simParams->dielectric);
00290     for (i = 0;  i < numAtoms;  i++) {
00291       chargebuffer[i] = celec * coord[i].charge;
00292     }
00293     fmmsolver = 1;  // is setup
00294   }
00295 
00296   // translate and scale positions into [-1/2,1/2]^3
00297   for (i = 0;  i < numAtoms;  i++) {
00298     posbuffer[i][0] = scaling*(coord[i].position.x - center.x);
00299     posbuffer[i][1] = scaling*(coord[i].position.y - center.y);
00300     posbuffer[i][2] = scaling*(coord[i].position.z - center.z);
00301   }
00302 
00303   // call the FMM solver
00304   int errcode = 0;
00305 #ifdef FMM_SOLVER
00306   fmmlap_uni_(&numAtoms, posbuffer, chargebuffer, epotbuffer, forcebuffer,
00307       &nlevels, &errcode);
00308 #else
00309   NAMD_die("Must link to FMM library to use FMM\n");
00310 #endif
00311 
00312   // scale force and potentials
00313   for (i = 0;  i < numAtoms;  i++) {
00314     double qs = chargebuffer[i] * scaling;
00315     force[i].x = qs * scaling * forcebuffer[i][0];
00316     force[i].y = qs * scaling * forcebuffer[i][1];
00317     force[i].z = qs * scaling * forcebuffer[i][2];
00318     energy += qs * epotbuffer[i];
00319   }
00320   energy *= 0.5;
00321 
00322   // distribute forces
00323   for (int j=0;  j < numSources;  j++) {
00324     FmmSerialCoordMsg *cmsg = coordMsgs[j];
00325     coordMsgs[j] = 0;
00326     FmmSerialForceMsg *fmsg = new (cmsg->numAtoms, 0) FmmSerialForceMsg;
00327 
00328     for (int i=0;  i < cmsg->numAtoms;  i++) {
00329       fmsg->force[i] = force[cmsg->coord[i].id];
00330     }
00331 
00332     if ( ! j ) {  // set virial and energy only for first message
00333       fmsg->energy = energy;
00334       for (int k=0;  k < 3;  k++) {
00335         for (int l=0;  l < 3;  l++) {
00336           fmsg->virial[k][l] = virial[k][l];
00337         }
00338       }
00339     }
00340     else {  // set other messages to zero, add into reduction only once
00341       fmsg->energy = 0;
00342       for (int k=0;  k < 3;  k++) {
00343         for (int l=0;  l < 3;  l++) {
00344           fmsg->virial[k][l] = 0;
00345         }
00346       }
00347     }
00348 
00349     fmmProxy[cmsg->sourceNode].recvForce(fmsg);
00350     delete cmsg;
00351   }
00352 }
00353 
00354 void ComputeFmmSerialMgr::recvForce(FmmSerialForceMsg *msg) {
00355   fmmCompute->saveResults(msg);
00356   delete oldmsg;
00357   oldmsg = msg;
00358 }
00359 
00360 void ComputeFmmSerial::saveResults(FmmSerialForceMsg *msg)
00361 {
00362   ResizeArrayIter<PatchElem> ap(patchList);
00363 
00364   FmmSerialForce *results_ptr = msg->force;
00365 
00366   // add in forces
00367   for (ap = ap.begin(); ap != ap.end(); ap++) {
00368     Results *r = (*ap).forceBox->open();
00369     Force *f = r->f[Results::slow];
00370     int numAtoms = (*ap).p->getNumAtoms();
00371 
00372     for(int i=0; i<numAtoms; ++i) {
00373       f[i].x += results_ptr->x;
00374       f[i].y += results_ptr->y;
00375       f[i].z += results_ptr->z;
00376       ++results_ptr;
00377     }
00378   
00379     (*ap).forceBox->close(&r);
00380   }
00381 
00382   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += msg->energy;
00383   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += msg->virial[0][0];
00384   reduction->item(REDUCTION_VIRIAL_SLOW_XY) += msg->virial[0][1];
00385   reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += msg->virial[0][2];
00386   reduction->item(REDUCTION_VIRIAL_SLOW_YX) += msg->virial[1][0];
00387   reduction->item(REDUCTION_VIRIAL_SLOW_YY) += msg->virial[1][1];
00388   reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += msg->virial[1][2];
00389   reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += msg->virial[2][0];
00390   reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += msg->virial[2][1];
00391   reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->virial[2][2];
00392   reduction->submit();
00393 }
00394 
00395 #include "ComputeFmmSerialMgr.def.h"

Generated on Tue Sep 19 01:17:10 2017 for NAMD by  doxygen 1.4.7