ComputeDPME.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 "ComputeDPME.h"
00013 #include "ComputeDPMEMsgs.h"
00014 #include "ComputeNonbondedUtil.h"
00015 #include "PatchMgr.h"
00016 #include "Molecule.h"
00017 #include "ReductionMgr.h"
00018 #include "SimParameters.h"
00019 #include "ComputeMgr.h"
00020 #include "ComputeMgr.decl.h"
00021 // #define DEBUGM
00022 #define MIN_DEBUG_LEVEL 3
00023 #include "Debug.h"
00024 
00025 #ifdef DPME
00026 #include "dpme2.h"
00027 
00028 class ComputeDPMEMaster {
00029 private:
00030   friend class ComputeDPME;
00031   ComputeDPME *host;
00032   ComputeDPMEMaster(ComputeDPME *);
00033   ~ComputeDPMEMaster();
00034   void recvData(ComputeDPMEDataMsg *);
00035   ResizeArray<int> homeNode;
00036   ResizeArray<int> endForNode;
00037   int numWorkingPes;
00038   int numLocalAtoms;
00039   Pme2Particle *localData;
00040   SubmitReduction *reduction;
00041   int runcount;
00042 };
00043 
00044 ComputeDPME::ComputeDPME(ComputeID c, ComputeMgr *m) :
00045   ComputeHomePatches(c), comm(m)
00046 {
00047   DebugM(4,"ComputeDPME created.\n");
00048   useAvgPositions = 1;
00049 
00050   int numWorkingPes = (PatchMap::Object())->numNodesWithPatches();
00051 
00052   masterNode = numWorkingPes - 1;
00053   if ( CkMyPe() == masterNode ) {
00054     master = new ComputeDPMEMaster(this);
00055     master->numWorkingPes = numWorkingPes;
00056   }
00057   else master = 0;
00058 }
00059 
00060 ComputeDPME::~ComputeDPME()
00061 {
00062   delete master;
00063 }
00064 
00065 // These are needed to fix up argument mismatches in DPME.
00066 
00067 extern "C" {
00068   extern int cfftf(int *, double *, double *);
00069   extern int cfftb(int *, double *, double *);
00070   extern int cffti1(int *, double *, int *);
00071   extern int cfftf1(int *n, double *c, double *ch, double *wa, int *ifac);
00072   extern int cfftb1(int *n, double *c, double *ch, double *wa, int *ifac);
00073 }
00074 
00075 int cfftf(int *n, doublecomplex *c, double *wsave) {
00076   // Casting (doublecomplex*) to (double*) is probably OK.
00077   return cfftf(n, (double *)c, wsave);
00078 }
00079 
00080 int cfftb(int *n, doublecomplex *c, double *wsave) {
00081   // Casting (doublecomplex*) to (double*) is probably OK.
00082   return cfftb(n, (double *)c, wsave);
00083 }
00084 
00085 int cffti1(int *n, double *wa, double *ifac) {
00086   // Casting (double*) to (int*) is dangerous if sizes differ!!!
00087   return cffti1(n, wa, (int *)ifac);
00088 }
00089 
00090 int cfftf1(int *n, double *c, double *ch, double *wa, double *ifac) {
00091   // Casting (double*) to (int*) is dangerous if sizes differ!!!
00092   return cfftf1(n, c, ch, wa, (int *)ifac);
00093 }
00094 
00095 int cfftb1(int *n, double *c, double *ch, double *wa, double *ifac) {
00096   // Casting (double*) to (int*) is dangerous if sizes differ!!!
00097   return cfftb1(n, c, ch, wa, (int *)ifac);
00098 }
00099 
00100 void ComputeDPME::doWork()
00101 {
00102   DebugM(4,"Entering ComputeDPME::doWork().\n");
00103 
00104   Pme2Particle *localData;
00105 
00106   ResizeArrayIter<PatchElem> ap(patchList);
00107 
00108   // Skip computations if nothing to do.
00109   if ( ! patchList[0].p->flags.doFullElectrostatics )
00110   {
00111     for (ap = ap.begin(); ap != ap.end(); ap++) {
00112       CompAtom *x = (*ap).positionBox->open();
00113       Results *r = (*ap).forceBox->open();
00114       (*ap).positionBox->close(&x);
00115       (*ap).forceBox->close(&r);
00116     }
00117     if ( master ) {
00118       master->reduction->submit();
00119     }
00120     return;
00121   }
00122 
00123   // allocate storage
00124   numLocalAtoms = 0;
00125   for (ap = ap.begin(); ap != ap.end(); ap++) {
00126     numLocalAtoms += (*ap).p->getNumAtoms();
00127   }
00128 
00129   Lattice lattice = patchList[0].p->flags.lattice;
00130 
00131   localData = new Pme2Particle[numLocalAtoms];  // given to message
00132 
00133   // get positions and charges
00134   Pme2Particle * data_ptr = localData;
00135   const BigReal coulomb_sqrt = sqrt( COULOMB * ComputeNonbondedUtil::scaling
00136                                 * ComputeNonbondedUtil::dielectric_1 );
00137   for (ap = ap.begin(); ap != ap.end(); ap++) {
00138     CompAtom *x = (*ap).positionBox->open();
00139     if ( patchList[0].p->flags.doMolly ) {
00140       (*ap).positionBox->close(&x);
00141       x = (*ap).avgPositionBox->open();
00142     }
00143     int numAtoms = (*ap).p->getNumAtoms();
00144 
00145     for(int i=0; i<numAtoms; ++i)
00146     {
00147       Vector tmp = lattice.delta(x[i].position);
00148       data_ptr->x = tmp.x;
00149       data_ptr->y = tmp.y;
00150       data_ptr->z = tmp.z;
00151       data_ptr->cg = coulomb_sqrt * x[i].charge;
00152       data_ptr->id = x[i].id;
00153       ++data_ptr;
00154     }
00155 
00156     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00157     else { (*ap).positionBox->close(&x); }
00158   }
00159 
00160   // send data to master
00161   ComputeDPMEDataMsg *msg = new ComputeDPMEDataMsg;
00162   msg->node = CkMyPe();
00163   msg->numParticles = numLocalAtoms;
00164   msg->particles = localData;
00165   comm->sendComputeDPMEData(msg);
00166 }
00167 
00168 void ComputeDPME::recvData(ComputeDPMEDataMsg *msg)
00169 {
00170   if ( master ) {
00171     master->recvData(msg);
00172   }
00173   else NAMD_die("ComputeDPME::master is NULL!");
00174 }
00175 
00176 ComputeDPMEMaster::ComputeDPMEMaster(ComputeDPME *h) :
00177   host(h), numLocalAtoms(0), runcount(0)
00178 {
00179   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00180   Molecule * molecule = Node::Object()->molecule;
00181   localData = new Pme2Particle[molecule->numAtoms];
00182 }
00183 
00184 ComputeDPMEMaster::~ComputeDPMEMaster()
00185 {
00186   delete reduction;
00187   delete [] localData;
00188 }
00189 
00190 void ComputeDPMEMaster::recvData(ComputeDPMEDataMsg *msg)
00191 { 
00192   DebugM(4,"ComputeDPMEMaster::recvData() " << msg->numParticles
00193         << " particles from node " << msg->node << "\n");
00194 
00195   {
00196     homeNode.add(msg->node);
00197     Pme2Particle *data_ptr = localData + numLocalAtoms;
00198     for ( int j = 0; j < msg->numParticles; ++j, ++data_ptr ) {
00199       *data_ptr = msg->particles[j];
00200     }
00201     numLocalAtoms += msg->numParticles;
00202     endForNode.add(numLocalAtoms);
00203     delete msg;
00204   }
00205 
00206   if ( homeNode.size() < numWorkingPes ) return;  // messages outstanding
00207 
00208   DebugM(4,"ComputeDPMEMaster::recvData() running serial code.\n");
00209 
00210   // single processor version
00211 
00212   Lattice lattice = host->getFlags()->lattice;
00213   SimParameters * simParams = Node::Object()->simParameters;
00214   int i;
00215 
00216   AtomInfo atom_info;
00217   atom_info.numatoms = numLocalAtoms;
00218   atom_info.atompnt = 0;  // not used
00219   atom_info.freepnt = 0;  // not used
00220   atom_info.nlocal = numLocalAtoms;
00221   atom_info.nother = 0;
00222 
00223   if ( ! lattice.orthogonal() ) {
00224     NAMD_die("DPME only supports orthogonal PBC's.");
00225   }
00226 
00227   BoxInfo box_info;
00228   box_info.box.x = lattice.a().x;
00229   box_info.box.y = lattice.b().y;
00230   box_info.box.z = lattice.c().z;
00231   box_info.box.origin = 0.;  // why only one number?
00232   box_info.prd.x = box_info.box.x;
00233   box_info.prd.y = box_info.box.y;
00234   box_info.prd.z = box_info.box.z;
00235   box_info.prd2.x = 0.5 * box_info.box.x;
00236   box_info.prd2.y = 0.5 * box_info.box.y;
00237   box_info.prd2.z = 0.5 * box_info.box.z;
00238   box_info.cutoff = simParams->cutoff;
00239   box_info.rs = simParams->cutoff;
00240   box_info.mc2.x = 2. * ( box_info.prd.x - box_info.cutoff );
00241   box_info.mc2.y = 2. * ( box_info.prd.y - box_info.cutoff );
00242   box_info.mc2.z = 2. * ( box_info.prd.z - box_info.cutoff );
00243   box_info.ewaldcof = ComputeNonbondedUtil::ewaldcof;
00244   box_info.dtol = simParams->PMETolerance;
00245   for (i = 0; i <= 8; i++) {
00246     box_info.recip[i] = 0.; /* assume orthogonal box */
00247   }
00248   box_info.recip[0] = 1.0/box_info.box.x;
00249   box_info.recip[4] = 1.0/box_info.box.y;
00250   box_info.recip[8] = 1.0/box_info.box.z;
00251 
00252   GridInfo grid_info;
00253   grid_info.order = simParams->PMEInterpOrder;
00254   grid_info.nfftgrd.x = simParams->PMEGridSizeX;
00255   grid_info.nfftgrd.y = simParams->PMEGridSizeY;
00256   grid_info.nfftgrd.z = simParams->PMEGridSizeZ;
00257   grid_info.nfft = 0;
00258   grid_info.volume = lattice.volume();
00259 
00260   PeInfo pe_info;  // hopefully this isn't used for anything
00261   pe_info.myproc.node = 0;
00262   pe_info.myproc.nprocs = 1;
00263   pe_info.myproc.ncube = 0;
00264   pe_info.inst_node[0] = 0;
00265   pe_info.igrid = 0;
00266 
00267   PmeVector *localResults;
00268   double recip_vir[6];
00269   int time_count = 0;
00270   int tsteps = 1;
00271   double mytime = 0.;
00272 
00273   // perform calculations
00274   BigReal electEnergy = 0;
00275 
00276   // calculate self energy
00277   Pme2Particle *data_ptr = localData;
00278   for(i=0; i<numLocalAtoms; ++i)
00279   {
00280     electEnergy += data_ptr->cg * data_ptr->cg;
00281     ++data_ptr;
00282   }
00283   electEnergy *= -1. * box_info.ewaldcof / SQRT_PI;
00284 
00285   DebugM(4,"Ewald self energy: " << electEnergy << "\n");
00286 
00287   DebugM(4,"Calling dpme_eval_recip().\n");
00288 
00289   double pme_start_time = 0;
00290   if ( runcount == 1 ) pme_start_time = CmiTimer();
00291 
00292   electEnergy += dpme_eval_recip( atom_info, localData - 1, &localResults,
00293                         recip_vir, grid_info, box_info, pe_info,
00294                         time_count, tsteps, &mytime );
00295 
00296   if ( runcount == 1 ) {
00297     iout << iINFO << "PME reciprocal sum CPU time per evaluation: "
00298          << (CmiTimer() - pme_start_time) << "\n" << endi;
00299   }
00300 
00301   DebugM(4,"Returned from dpme_eval_recip().\n");
00302 
00303   // REVERSE SIGN OF VIRIAL RETURNED BY DPME
00304   for(i=0; i<6; ++i) recip_vir[i] *= -1.;
00305 
00306   // send out reductions
00307   DebugM(4,"Timestep : " << host->getFlags()->step << "\n");
00308   DebugM(4,"Reciprocal sum energy: " << electEnergy << "\n");
00309   DebugM(4,"Reciprocal sum virial: " << recip_vir[0] << " " <<
00310         recip_vir[1] << " " << recip_vir[2] << " " << recip_vir[3] << " " <<
00311         recip_vir[4] << " " << recip_vir[5] << "\n");
00312   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += electEnergy;
00313   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += (BigReal)(recip_vir[0]);
00314   reduction->item(REDUCTION_VIRIAL_SLOW_XY) += (BigReal)(recip_vir[1]);
00315   reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += (BigReal)(recip_vir[2]);
00316   reduction->item(REDUCTION_VIRIAL_SLOW_YX) += (BigReal)(recip_vir[1]);
00317   reduction->item(REDUCTION_VIRIAL_SLOW_YY) += (BigReal)(recip_vir[3]);
00318   reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += (BigReal)(recip_vir[4]);
00319   reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += (BigReal)(recip_vir[2]);
00320   reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += (BigReal)(recip_vir[4]);
00321   reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += (BigReal)(recip_vir[5]);
00322   reduction->submit();
00323 
00324   PmeVector *results_ptr = localResults + 1;
00325 
00326   numLocalAtoms = 0;
00327   for ( i = 0; i < homeNode.size(); ++i ) {
00328     ComputeDPMEResultsMsg *msg = new ComputeDPMEResultsMsg;
00329     msg->node = homeNode[i];
00330     msg->numParticles = endForNode[i] - numLocalAtoms;
00331     msg->forces = new PmeVector[msg->numParticles];
00332     for ( int j = 0; j < msg->numParticles; ++j, ++results_ptr ) {
00333       msg->forces[j] = *results_ptr;
00334     }
00335     numLocalAtoms = endForNode[i];
00336     host->comm->sendComputeDPMEResults(msg,homeNode[i]);
00337   }
00338 
00339   // reset
00340   runcount += 1;
00341   numLocalAtoms = 0;
00342   homeNode.resize(0);
00343   endForNode.resize(0);
00344 
00345 }
00346 
00347 void ComputeDPME::recvResults(ComputeDPMEResultsMsg *msg)
00348 {
00349   if ( CkMyPe() != msg->node ) {
00350     NAMD_die("ComputeDPME results sent to wrong node!\n");
00351     return;
00352   }
00353   if ( numLocalAtoms != msg->numParticles ) {
00354     NAMD_die("ComputeDPME sent wrong number of results!\n");
00355     return;
00356   }
00357 
00358   PmeVector *results_ptr = msg->forces;
00359   ResizeArrayIter<PatchElem> ap(patchList);
00360 
00361   // add in forces
00362   for (ap = ap.begin(); ap != ap.end(); ap++) {
00363     Results *r = (*ap).forceBox->open();
00364     Force *f = r->f[Results::slow];
00365     int numAtoms = (*ap).p->getNumAtoms();
00366 
00367     for(int i=0; i<numAtoms; ++i)
00368     {
00369       f[i].x += results_ptr->x;
00370       f[i].y += results_ptr->y;
00371       f[i].z += results_ptr->z;
00372       ++results_ptr;
00373     }
00374 
00375     (*ap).forceBox->close(&r);
00376   }
00377 
00378   delete msg;
00379 
00380 }
00381 
00382 #endif  // DPME
00383 

Generated on Fri Sep 22 01:17:11 2017 for NAMD by  doxygen 1.4.7