ComputeFullDirect.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 "ComputeFullDirect.h"
00013 #include "ComputeNonbondedUtil.h"
00014 #include "PatchMgr.h"
00015 #include "Molecule.h"
00016 #include "ReductionMgr.h"
00017 #include "Communicate.h"
00018 #include "Lattice.h"
00019 //#define DEBUGM
00020 #define MIN_DEBUG_LEVEL 3
00021 #include "Debug.h"
00022 #include "SimParameters.h"
00023 
00024 ComputeFullDirect::ComputeFullDirect(ComputeID c) : ComputeHomePatches(c)
00025 {
00026   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00027   SimParameters *simParams = Node::Object()->simParameters;
00028   useAvgPositions = 1;
00029 }
00030 
00031 ComputeFullDirect::~ComputeFullDirect()
00032 {
00033   delete reduction;
00034 }
00035 
00036 BigReal calc_fulldirect(BigReal *data1, BigReal *results1, int n1,
00037                         BigReal *data2, BigReal *results2, int n2,
00038                         int selfmode, Lattice *lattice, Tensor &virial)
00039 {
00040   if ( lattice->a_p() || lattice->b_p() || lattice->c_p() ) {
00041     #define FULLDIRECT_PERIODIC
00042     #include "ComputeFullDirectBase.h"
00043   } else {
00044     #undef FULLDIRECT_PERIODIC
00045     #include "ComputeFullDirectBase.h"
00046   }
00047 }
00048 
00049 void ComputeFullDirect::doWork()
00050 {
00051   int numLocalAtoms;
00052   BigReal *localData;
00053   BigReal *localResults;
00054   BigReal *newLocalResults;
00055   register BigReal *local_ptr;
00056   Lattice *lattice;
00057 
00058   int numWorkingPes = (PatchMap::Object())->numNodesWithPatches();
00059 
00060   if ( numWorkingPes > 1 ) NAMD_die("FullDirect not supported for parallel runs.");
00061 
00062   ResizeArrayIter<PatchElem> ap(patchList);
00063 
00064   // Skip computations if nothing to do.
00065   if ( ! patchList[0].p->flags.doFullElectrostatics )
00066   {
00067     for (ap = ap.begin(); ap != ap.end(); ap++) {
00068       CompAtom *x = (*ap).positionBox->open();
00069       Results *r = (*ap).forceBox->open();
00070       (*ap).positionBox->close(&x);
00071       (*ap).forceBox->close(&r);
00072     }
00073     reduction->submit();
00074     return;
00075   }
00076 
00077   // allocate storage
00078   numLocalAtoms = 0;
00079   for (ap = ap.begin(); ap != ap.end(); ap++) {
00080     numLocalAtoms += (*ap).p->getNumAtoms();
00081   }
00082 
00083   localData = new BigReal[4*numLocalAtoms];     // freed at end of this method
00084   localResults = new BigReal[3*numLocalAtoms];  // freed at end of this method
00085   newLocalResults = new BigReal[3*numLocalAtoms];  // freed at end of this method
00086 
00087   lattice = &((*(ap.begin())).p->lattice);
00088 
00089   // get positions and charges
00090   local_ptr = localData;
00091   for (ap = ap.begin(); ap != ap.end(); ap++) {
00092     CompAtom *x = (*ap).positionBox->open();
00093     if ( patchList[0].p->flags.doMolly ) {
00094       (*ap).positionBox->close(&x);
00095       x = (*ap).avgPositionBox->open();
00096     }
00097     int numAtoms = (*ap).p->getNumAtoms();
00098 
00099     for(int i=0; i<numAtoms; ++i)
00100     {
00101       *(local_ptr++) = x[i].position.x;
00102       *(local_ptr++) = x[i].position.y;
00103       *(local_ptr++) = x[i].position.z;
00104       *(local_ptr++) = x[i].charge;
00105     }
00106 
00107     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00108     else { (*ap).positionBox->close(&x); }
00109   } 
00110 
00111   // zero out forces
00112   local_ptr = localResults;
00113   for(int j=0; j<numLocalAtoms; ++j)
00114   {
00115     *(local_ptr++) = 0.;
00116     *(local_ptr++) = 0.;
00117     *(local_ptr++) = 0.;
00118   }
00119 
00120   // perform calculations
00121   BigReal electEnergy = 0;
00122   Tensor virial;
00123 
00124 #define PEMOD(N) (((N)+numWorkingPes)%numWorkingPes)
00125 
00126   int numStages = numWorkingPes / 2 + 2;
00127   int lastStage = numStages - 2;
00128   int sendDataPE = PEMOD(CkMyPe()+1);
00129   int recvDataPE = PEMOD(CkMyPe()-1);
00130   int sendResultsPE = PEMOD(CkMyPe()-1);
00131   int recvResultsPE = PEMOD(CkMyPe()+1);
00132   int numRemoteAtoms = numLocalAtoms;
00133   int oldNumRemoteAtoms = 0;
00134   BigReal *remoteData = 0;
00135   BigReal *remoteResults = 0;
00136   register BigReal *remote_ptr;
00137   register BigReal *end_ptr;
00138 
00139   MOStream *sendDataMsg=CkpvAccess(comm)->
00140                 newOutputStream(sendDataPE, FULLTAG, BUFSIZE);
00141   MIStream *recvDataMsg=CkpvAccess(comm)->
00142                 newInputStream(recvDataPE, FULLTAG);
00143 
00144   for ( int stage = 0; stage < numStages; ++stage )
00145   {
00146     // send remoteResults to sendResultsPE
00147     if ( stage > 1 )
00148     {
00149       DebugM(4,"send remoteResults to sendResultsPE " << sendResultsPE << "\n");
00150       MOStream *msg=CkpvAccess(comm)->
00151                 newOutputStream(sendResultsPE, FULLFORCETAG, BUFSIZE);
00152       msg->put(3*oldNumRemoteAtoms,remoteResults);
00153       delete [] remoteResults;
00154       msg->end();
00155       delete msg;
00156       sendResultsPE = PEMOD(sendResultsPE-1);
00157     }
00158 
00159     // send remoteData to sendDataPE
00160     if ( stage < lastStage )
00161     {
00162       DebugM(4,"send remoteData to sendDataPE " << sendDataPE << "\n");
00163       sendDataMsg->put(numRemoteAtoms);
00164       sendDataMsg->put(4*numRemoteAtoms,(stage?remoteData:localData));
00165       sendDataMsg->end();
00166     }
00167 
00168     // allocate new result storage
00169     if ( stage > 0 && stage <= lastStage )
00170     {
00171       DebugM(4,"allocate new result storage\n");
00172       remoteResults = new BigReal[3*numRemoteAtoms];
00173       remote_ptr = remoteResults;
00174       end_ptr = remoteResults + 3*numRemoteAtoms;
00175       for ( ; remote_ptr != end_ptr; ++remote_ptr ) *remote_ptr = 0.;
00176     }
00177 
00178     // do calculation
00179     if ( stage == 0 )
00180     {  // self interaction
00181       DebugM(4,"self interaction\n");
00182       electEnergy += calc_fulldirect(
00183         localData,localResults,numLocalAtoms,
00184         localData,localResults,numLocalAtoms,1,lattice,virial);
00185     }
00186     else if ( stage < lastStage ||
00187             ( stage == lastStage && ( numWorkingPes % 2 ) ) )
00188     {  // full other interaction
00189       DebugM(4,"full other interaction\n");
00190       electEnergy += calc_fulldirect(
00191         localData,localResults,numLocalAtoms,
00192         remoteData,remoteResults,numRemoteAtoms,0,lattice,virial);
00193     }
00194     else if ( stage == lastStage )
00195     {  // half other interaction
00196       DebugM(4,"half other interaction\n");
00197       if ( CkMyPe() < ( numWorkingPes / 2 ) )
00198         electEnergy += calc_fulldirect(
00199           localData,localResults,numLocalAtoms/2,
00200           remoteData,remoteResults,numRemoteAtoms,0,lattice,virial);
00201       else
00202         electEnergy += calc_fulldirect(
00203           localData,localResults,numLocalAtoms,
00204           remoteData + 4*(numRemoteAtoms/2),
00205           remoteResults + 3*(numRemoteAtoms/2),
00206           numRemoteAtoms - (numRemoteAtoms/2), 0,lattice,virial);
00207     }
00208 
00209     delete [] remoteData;  remoteData = 0;
00210     oldNumRemoteAtoms = numRemoteAtoms;
00211 
00212     // receive newLocalResults from recvResultsPE
00213     if ( stage > 1 )
00214     {
00215       DebugM(4,"receive newLocalResults from recvResultsPE "
00216                                                 << recvResultsPE << "\n");
00217       MIStream *msg=CkpvAccess(comm)->
00218                 newInputStream(recvResultsPE, FULLFORCETAG);
00219       msg->get(3*numLocalAtoms,newLocalResults);
00220       delete msg;
00221       recvResultsPE = PEMOD(recvResultsPE+1);
00222       remote_ptr = newLocalResults;
00223       local_ptr = localResults;
00224       end_ptr = localResults + 3*numLocalAtoms;
00225       for ( ; local_ptr != end_ptr; ++local_ptr, ++remote_ptr )
00226         *local_ptr += *remote_ptr;
00227     }
00228 
00229     // receive remoteData from recvDataPE
00230     if ( stage < lastStage )
00231     {
00232       DebugM(4,"receive remoteData from recvDataPE "
00233                                                 << recvDataPE << "\n");
00234       recvDataMsg->get(numRemoteAtoms);
00235       remoteData = new BigReal[4*numRemoteAtoms];
00236       recvDataMsg->get(4*numRemoteAtoms,remoteData);
00237     }
00238 
00239   }
00240 
00241   delete sendDataMsg;
00242   delete recvDataMsg;
00243 
00244   // send out reductions
00245   DebugM(4,"Full-electrostatics energy: " << electEnergy << "\n");
00246   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += electEnergy;
00247   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += virial.xx;
00248   reduction->item(REDUCTION_VIRIAL_SLOW_XY) += virial.xy;
00249   reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += virial.xz;
00250   reduction->item(REDUCTION_VIRIAL_SLOW_YX) += virial.yx;
00251   reduction->item(REDUCTION_VIRIAL_SLOW_YY) += virial.yy;
00252   reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += virial.yz;
00253   reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += virial.zx;
00254   reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += virial.zy;
00255   reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += virial.zz;
00256   reduction->submit();
00257 
00258   // add in forces
00259   local_ptr = localResults;
00260   for (ap = ap.begin(); ap != ap.end(); ap++) {
00261     Results *r = (*ap).forceBox->open();
00262     Force *f = r->f[Results::slow];
00263     int numAtoms = (*ap).p->getNumAtoms();
00264 
00265     for(int i=0; i<numAtoms; ++i)
00266     {
00267       f[i].x += *(local_ptr++);
00268       f[i].y += *(local_ptr++);
00269       f[i].z += *(local_ptr++);
00270     }
00271 
00272     (*ap).forceBox->close(&r);
00273   }
00274 
00275   // free storage
00276   delete [] localData;          // allocated at beginning of this method
00277   delete [] localResults;       // allocated at beginning of this method
00278   delete [] newLocalResults;    // allocated at beginning of this method
00279 }
00280 

Generated on Thu Sep 21 01:17:10 2017 for NAMD by  doxygen 1.4.7