ComputeDPMTA.C

Go to the documentation of this file.
00001 
00007 #include "common.h"
00008 #include "InfoStream.h"
00009 #include "Node.h"
00010 #include "SimParameters.h"
00011 #include "ComputeNonbondedUtil.h"
00012 #include "PatchMap.inl"
00013 #include "AtomMap.h"
00014 #include "ComputeDPMTA.h"
00015 #include "PatchMgr.h"
00016 #include "Molecule.h"
00017 #include "ReductionMgr.h"
00018 #include "Lattice.h"
00019 #include "Communicate.h"
00020 #include "InfoStream.h"
00021 #include "ProcessorPrivate.h"
00022 
00023 #define MIN_DEBUG_LEVEL 2
00024 // #define DEBUGM
00025 #include "Debug.h"
00026 
00027 #ifdef DPMTA
00028 
00029 #include <pvmc.h>
00030 
00031 // #define DUMP_DPMTA
00032 
00033 void ComputeDPMTA::get_FMA_cube(int resize)
00034 {
00035   Vector boxSize,boxCenter;     // used to see if things change
00036   PatchMap *patchMap = PatchMap::Object();
00037 
00038   if (usePBC == FALSE)
00039   {
00040     //  From these extremes, figure out how many patches we will
00041     //  have to have in each direction
00042     SimParameters *simParams = Node::Object()->simParameters;
00043     int dim_x = patchMap->gridsize_a();
00044     int dim_y = patchMap->gridsize_b();
00045     int dim_z = patchMap->gridsize_c();
00046 
00047     boxSize.x = dim_x*simParams->patchDimension;
00048     boxSize.y = dim_y*simParams->patchDimension;
00049     boxSize.z = dim_z*simParams->patchDimension;
00050     BigReal skirt = 2*simParams->patchDimension;
00051 
00052     boxCenter = patchMap->origin();
00053     boxCenter.x += boxSize.x/2.0;
00054     boxCenter.y += boxSize.y/2.0;
00055     boxCenter.z += boxSize.z/2.0;
00056 
00057     //  add the skirt of empty patches by adding 2 patches in every direction
00058     boxSize.x += skirt;
00059     boxSize.y += skirt;
00060     boxSize.z += skirt;
00061   }
00062   else
00063   {
00064     DebugM(2,"getting patch info for FMA box\n");
00065 
00066     // determine boxSize from the PBC lattice
00067     // lattice is the same on all patches, so choose first patch
00068     ResizeArrayIter<PatchElem> ap(patchList);
00069     DebugM(2,"getting first patch info for FMA box\n");
00070     ap = ap.begin();
00071     DebugM(2,"getting lattice from patch for FMA box\n");
00072     Lattice lattice = (*ap).p->lattice;
00073     if ( ! lattice.orthogonal() ) {
00074       NAMD_die("DPMTA (FMA) only supports orthogonal PBC's.");
00075     }
00076     DebugM(2,"getting patch dimension for FMA box\n");
00077     boxSize.x = lattice.a().x;
00078     boxSize.y = lattice.b().y;
00079     boxSize.z = lattice.c().z;
00080     DebugM(2,"boxSize is " << boxSize << "\n");
00081     boxCenter = lattice.origin();
00082   }
00083 
00084   // don't bother checking if the center has moved since it depends on the size.
00085   if (boxsize != boxSize)
00086   {
00087         DebugM(2,"resetting FMA box\n");
00088         // reset the size and center
00089         boxsize = boxSize;
00090         boxcenter = boxCenter;
00091 
00092         // reset DPMTA (only reset it after it has been initialized!)
00093         if (resize && usePBC)
00094         {
00095           PmtaVector center,v1,v2,v3;
00096           center.x = boxcenter.x;
00097           center.y = boxcenter.y;
00098           center.z = boxcenter.z;
00099           v1.x = boxsize.x;
00100           v2.y = boxsize.y;
00101           v3.z = boxsize.z;
00102           iout << iINFO << "DPMTA box resized:\n";
00103           iout << iINFO << "BOX DIMENSIONS = (" << v1.x << ","
00104                 << v2.y << "," << v3.z << ")\n";
00105           iout << iINFO << "BOX CENTER = (" << center.x << ","
00106                 << center.y << "," << center.z << ")\n";
00107           iout << endi;
00108           DebugM(2,"calling PMTAresize()\n");
00109           PMTAresize(&v1,&v2,&v3,&center);
00110           DebugM(2,"called PMTAresize()\n");
00111         }
00112   }
00113   DebugM(2,"cube center: " << boxcenter << " size=" << boxsize << "\n");
00114 }
00115 
00116 ComputeDPMTA::ComputeDPMTA(ComputeID c) : ComputeHomePatches(c)
00117 {
00118   useAvgPositions = 1;
00119 }
00120 
00121 void ComputeDPMTA::initialize()
00122 {
00123   ComputeHomePatches::initialize();
00124 
00125   DebugM(2,"ComputeDPMTA creating\n");
00126   // comm should always be initialized by this point...
00127   // In the (bug) case that it isn't, then initialize it.
00128   if (CkpvAccess(comm) == NULL)
00129   {
00130     NAMD_die("Communication protocol (Converse, PVM, etc.) not initialized.");
00131   }
00132 
00133   // **** NOTE: node 0 must initialized before any other nodes register.
00134 
00135   //  Set everything to 0
00136   totalAtoms = 0;
00137   fmaResults = NULL;
00138   ljResults = NULL;
00139   boxcenter = 1;        // reset the array (no divide by zero)
00140   boxsize = 1;  // reset the array (no divide by zero)
00141   usePBC = FALSE;       // assume not...
00142 
00143   // all nodes should init
00144   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00145 
00146   // Don't need any more initialization  -JCP
00147   ResizeArrayIter<PatchElem> ap(patchList);
00148   DebugM(2,"init() getting first patch info for FMA box\n");
00149   ap = ap.begin();
00150   DebugM(2,"init() getting lattice from patch for FMA box\n");
00151   initLattice.x = (*ap).p->lattice.a().x;
00152   initLattice.y = (*ap).p->lattice.b().y;
00153   initLattice.z = (*ap).p->lattice.c().z;
00154   DebugM(2,"init() initLattice is " << initLattice << "\n");
00155 
00156   //  NOTE that the theta value is hardwired to the value of 0.715
00157   //  as per the recommendation of the Duke developers
00158 
00159   //  NOTE 2: Theta is now an optional config parameter,
00160   //  but it defaults to 0.715
00161 
00162   // check for PBC
00163   usePBC = ( patchMap->periodic_a() ? 1 : 0 )
00164          + ( patchMap->periodic_b() ? 1 : 0 )
00165          + ( patchMap->periodic_c() ? 1 : 0 );
00166   if ((usePBC != 0) && (usePBC != 3))
00167   {
00168     NAMD_die("DPMTA (FMA) does not support 1 or 2 dimensional PBC's.");
00169   }
00170   DebugM(2,"Use PBC = " << usePBC << "\n");
00171   usePBC = (usePBC == 3);       // either PBC "3D" or no PBC
00172   if ( usePBC ) {
00173     iout << iWARN << "DPMTA (FMA) pressure tensor is incorrect.\n"
00174         << iWARN << "Do not use DPMTA with anisotropic pressure control.\n"
00175         << endi;
00176   }
00177 
00178   //  Get the size of the FMA cube
00179   DebugM(2,"DPMTA getting FMA cube\n");
00180   get_FMA_cube(FALSE);
00181   DebugM(2,"DPMTA got FMA cube\n");
00182 
00183   if (CkMyPe() != 0)
00184   {
00185     DebugM(2,"waiting for Init go-ahead\n");
00186     MIStream *msg2 = CkpvAccess(comm)->newInputStream(ANY, DPMTATAG);
00187     int dummy;
00188     msg2->get(dummy);
00189     delete msg2;
00190     slavetids=NULL;
00191     if (PMTAregister() < 0)
00192     {
00193         NAMD_die("PMTARegister failed!!");
00194     }
00195     DebugM(2,"DPMTA done PMTAinit.\n");
00196     return;
00197   }
00198   DebugM(2,"DPMTA configuring\n");
00199 
00200   // *****************************************
00201   // ONLY THE MASTER (NODE 0) NEEDS TO DO THIS:
00202 
00203   int numProcs = (PatchMap::Object())->numNodesWithPatches();
00204 
00205   slavetids = new int[numProcs];
00206   if (slavetids == NULL)
00207   {
00208     NAMD_die("Memory allocation failed in FMAInterface::FMAInterface");
00209   }
00210 
00211   // pvm_spawn is a dummy function under Converse.  Just the array is required.
00212   pvm_spawn(NULL,NULL,0,NULL,numProcs,slavetids);
00213   DebugM(2,"DPMTA slavetids allocated\n");
00214 
00215   // reduce function calling time
00216   SimParameters *simParams = Node::Object()->simParameters;
00217 
00218   //  initialize DPMTA
00219   PmtaInitData pmta_data;
00220   memset(&pmta_data,0,sizeof(pmta_data));
00221   pmta_data.nprocs = numProcs;
00222   pmta_data.nlevels = simParams->FMALevels;
00223   pmta_data.mp = simParams->FMAMp;
00224   pmta_data.mp_lj = 4;
00225   pmta_data.fft = simParams->FMAFFTOn;
00226   pmta_data.fftblock = simParams->FMAFFTBlock;
00227   pmta_data.pbc = usePBC;       // use Periodic boundary condition
00228   pmta_data.kterm = 0;
00229   pmta_data.theta = simParams->fmaTheta;
00230   pmta_data.v1.x = boxsize.x;
00231   pmta_data.v1.y = 0.;
00232   pmta_data.v1.z = 0.;
00233   pmta_data.v2.x = 0.;
00234   pmta_data.v2.y = boxsize.y;
00235   pmta_data.v2.z = 0.;
00236   pmta_data.v3.x = 0.;
00237   pmta_data.v3.y = 0.;
00238   pmta_data.v3.z = boxsize.z;
00239   pmta_data.cellctr.x = boxcenter.x;
00240   pmta_data.cellctr.y = boxcenter.y;
00241   pmta_data.cellctr.z = boxcenter.z;
00242   pmta_data.calling_num = pmta_data.nprocs;
00243   pmta_data.calling_tids = slavetids;
00244 
00245   iout << iINFO << "DPMTA parameters are:\n";
00246   iout << iINFO << "  LEVELS = " << pmta_data.nlevels << "\n";
00247   iout << iINFO << "  NUMBER OF MULTIPOLE TERMS = " << pmta_data.mp << "\n";
00248   iout << iINFO << "  FFT FLAG = " << pmta_data.fft << "\n";
00249   iout << iINFO << "  FFT BLOCKING FACTOR = " << pmta_data.fftblock << "\n";
00250   if ( usePBC ) iout << iINFO << "  SYSTEM IS PERIODIC\n" << endi;
00251   iout << iINFO << "  BOX DIMENSIONS = (" << pmta_data.v1.x << ","
00252         << pmta_data.v2.y << "," << pmta_data.v3.z << ")\n";
00253   iout << iINFO << "  BOX CENTER = (" << pmta_data.cellctr.x << ","
00254         << pmta_data.cellctr.y << "," << pmta_data.cellctr.z << ")\n";
00255   iout << endi;
00256 
00257   if ( usePBC )
00258   {
00259     pmta_data.cellctr.x = 0.;
00260     pmta_data.cellctr.y = 0.;
00261     pmta_data.cellctr.z = 0.;
00262   }
00263 
00264 #ifdef DUMP_DPMTA
00265   FILE *fp;
00266   fp = fopen("DUMP_DPMTA.init","w");
00267   fwrite(&pmta_data,sizeof(PmtaInitData),1,fp);
00268   fclose(fp);
00269 #endif
00270 
00271   DebugM(2,"DPMTA calling PMTAinit.\n");
00272   if (PMTAinit(&pmta_data,slavetids) >= 0)
00273   {
00274         iout << iINFO << "SUCCESSFULLY STARTED DPMTA\n" << endi;
00275   }
00276   else
00277   {
00278         NAMD_die("Unable to start DPMTA!");
00279   }
00280 
00281   // tell all nodes that it is OK to register
00282   MOStream *msg = CkpvAccess(comm)->newOutputStream(ALL, DPMTATAG, BUFSIZE);
00283   // don't actually put in data...  Nodes just need it as a flag.
00284   msg->put(TRUE);
00285   msg->end();
00286   delete msg;
00287   DebugM(2,"Init go-ahead\n");
00288   MIStream *msg1 = CkpvAccess(comm)->newInputStream(ANY, DPMTATAG);
00289   int dummy1;
00290   msg1->get(dummy1);
00291   delete msg1;
00292   DebugM(2,"got Init go-ahead\n");
00293 
00294   //  Register this master with the other DPMTA processes
00295   if (PMTAregister() < 0)
00296   {
00297         NAMD_die("PMTARegister failed!!");
00298   }
00299   DebugM(2,"DPMTA done PMTAinit.\n");
00300   DebugM(2,"DPMTA configured\n");
00301 }
00302 
00303 ComputeDPMTA::~ComputeDPMTA()
00304 {
00305   DebugM(2,"DPMTA exiting\n");
00306   //  If this is the master node, then call PMTAexit()
00307   if (CkMyPe() == 0)    PMTAexit();
00308 
00309   if (fmaResults)
00310         {
00311         free(fmaResults);
00312         fmaResults = NULL;
00313         }
00314   delete [] ljResults;
00315   delete [] slavetids;
00316   DebugM(2,"DPMTA exited\n");
00317 
00318   delete reduction;
00319 }
00320 
00321 
00322 void ComputeDPMTA::doWork()
00323 {
00324   ResizeArrayIter<PatchElem> ap(patchList);
00325   PmtaParticle *particle_list = NULL;
00326 
00327   // 0. only run when necessary
00328   // Skip computations if nothing to do.
00329   if (!patchList[0].p->flags.doFullElectrostatics)
00330   {
00331     for (ap = ap.begin(); ap != ap.end(); ap++) {
00332       CompAtom *x = (*ap).positionBox->open();
00333       Results *r = (*ap).forceBox->open();
00334       (*ap).forceBox->close(&r);
00335       (*ap).positionBox->close(&x);
00336     }
00337     reduction->submit();
00338     return;
00339   }
00340 
00341   // setup
00342   // 1. get totalAtoms
00343   for (totalAtoms=0, ap = ap.begin(); ap != ap.end(); ap++)
00344      totalAtoms += (*ap).p->getNumAtoms();
00345 
00346   Vector newLattice;
00347   Vector rescaleFactor;
00348   if (usePBC)
00349     {
00350     ap = ap.begin();
00351     Lattice lattice = (*ap).p->lattice;
00352     if ( ! lattice.orthogonal() ) {
00353       NAMD_die("DPMTA (FMA) only supports orthogonal PBC's.");
00354     }
00355     newLattice.x = lattice.a().x;
00356     newLattice.y = lattice.b().y;
00357     newLattice.z = lattice.c().z;
00358     rescaleFactor.x = initLattice.x / newLattice.x;
00359     rescaleFactor.y = initLattice.y / newLattice.y;
00360     rescaleFactor.z = initLattice.z / newLattice.z;
00361     DebugM(2,"Rescale factor = " << initLattice << "/" << newLattice
00362                 << " = " << rescaleFactor << "\n");
00363     DebugM(2,"boxcenter = " << boxcenter << "\n");
00364     }
00365   else
00366     {
00367     rescaleFactor.x = 1;
00368     rescaleFactor.y = 1;
00369     rescaleFactor.z = 1;
00370     }
00371 
00372   // 2. setup atom list
00373   int i,j;
00374   particle_list = (PmtaParticle *)calloc(totalAtoms,sizeof(PmtaParticle));
00375   fmaResults =    (PmtaPartInfo *)calloc(totalAtoms,sizeof(PmtaPartInfo));
00376   if (!particle_list || !fmaResults)
00377         {
00378         NAMD_die("DPMTA Failed to allocate memory.");
00379         }
00380 
00381   BigReal unitFactor = sqrt(COULOMB * ComputeNonbondedUtil::scaling
00382                                 * ComputeNonbondedUtil::dielectric_1);
00383   DebugM(2,"Charge unit factor = " << unitFactor << "\n");
00384   for (i=0, ap = ap.begin(); ap != ap.end(); ap++)
00385   {
00386     CompAtom *x = (*ap).positionBox->open();
00387     if ( patchList[0].p->flags.doMolly ) {
00388       (*ap).positionBox->close(&x);
00389       x = (*ap).avgPositionBox->open();
00390     }
00391 
00392     // store each atom in the particle_list
00393     Vector pos;
00394     for(j=0; j<(*ap).p->getNumAtoms(); j++)
00395     {
00396       // explicitly copy -- two different data structures
00397       if (usePBC)
00398         {
00399         particle_list[i].p.x = rescaleFactor.x * (x[j].position.x-boxcenter.x);
00400         particle_list[i].p.y = rescaleFactor.y * (x[j].position.y-boxcenter.y);
00401         particle_list[i].p.z = rescaleFactor.z * (x[j].position.z-boxcenter.z);
00402         }
00403       else
00404         {
00405         particle_list[i].p.x = x[j].position.x;
00406         particle_list[i].p.y = x[j].position.y;
00407         particle_list[i].p.z = x[j].position.z;
00408         }
00409       particle_list[i].q = x[j].charge * unitFactor;
00410       DebugM(1,"atom[" << i << "]=" << x[j] << " "
00411               << x[j].charge*unitFactor << "\n");
00412       i++;
00413       if (i > totalAtoms)
00414         {
00415         iout << iERRORF << iPE << " totalAtoms=" << totalAtoms
00416              << " but " << i << " atoms are seen!\n" << endi;
00417         NAMD_die("FMA: atom counts unequal!");
00418         }
00419     }
00420 
00421     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00422     else { (*ap).positionBox->close(&x); }
00423   } 
00424 
00425   if (i != totalAtoms)
00426   {
00427     iout << iERRORF << iPE << " totalAtoms=" << totalAtoms
00428          << " but " << i << " atoms are seen!\n" << endi;
00429     NAMD_die("FMA: atom counts unequal!");
00430   }
00431 
00432   DebugM(2,"DPMTA doWork() there are " << totalAtoms << " atoms in this node.\n");
00433 
00434 #ifdef DUMP_DPMTA
00435   FILE *fp;
00436   char dump_file[32];
00437   sprintf(dump_file,"DUMP_DPMTA.%d",(int)CkMyPe());
00438   fp = fopen(dump_file,"w");
00439   int32 n32 = i;
00440   fwrite(&n32,sizeof(int32),1,fp);
00441   fwrite(particle_list,sizeof(PmtaParticle),i,fp);
00442   fclose(fp);
00443 #endif
00444 
00445   // 3. (run DPMTA) compute the forces
00446   if ( PMTAforce(i, particle_list, fmaResults, NULL) < 0 )
00447     {
00448       NAMD_die("PMTAforce failed!!");
00449     }
00450   DebugM(2,"DPMTA forces done.  Now depositing.\n");
00451 
00452   // 4. deposit
00453   BigReal potential=0;
00454   for (i=0, ap = ap.begin(); ap != ap.end(); ap++)
00455   {
00456     Results *r = (*ap).forceBox->open();
00457     Force *f = r->f[Results::slow];
00458 
00459     // deposit here
00460     for(j=0; j<(*ap).p->getNumAtoms(); j++)
00461     {
00462       if (usePBC)
00463         {
00464         f[j].x += fmaResults[i].f.x * rescaleFactor.x * rescaleFactor.x;
00465         f[j].y += fmaResults[i].f.y * rescaleFactor.y * rescaleFactor.y;
00466         f[j].z += fmaResults[i].f.z * rescaleFactor.z * rescaleFactor.z;
00467         potential += fmaResults[i].v * rescaleFactor.x;
00468         }
00469       else
00470         {
00471         f[j].x += fmaResults[i].f.x;
00472         f[j].y += fmaResults[i].f.y;
00473         f[j].z += fmaResults[i].f.z;
00474         potential += fmaResults[i].v;
00475         }
00476       i++;
00477     }
00478 
00479     (*ap).forceBox->close(&r);
00480   }
00481 
00482   potential *= 0.5;
00483   DebugM(4,"Full-electrostatics energy: " << potential << "\n");
00484   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += potential;
00485   // DPMTA won't work correctly if scaled anisotropically anyway.  -JCP
00486   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += potential / 3.;
00487   reduction->item(REDUCTION_VIRIAL_SLOW_YY) += potential / 3.;
00488   reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += potential / 3.;
00489   reduction->submit();
00490 
00491   // 5. clean-up
00492   if (totalAtoms > 0)
00493   {
00494     free(particle_list);
00495     free(fmaResults);
00496   }
00497 
00498   DebugM(2,"DPMTA doWork() done\n");
00499 }
00500 
00501 #endif
00502 

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