Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

WorkDistrib.C

Go to the documentation of this file.
00001 
00007 /*****************************************************************************
00008  * $Source: /home/cvs/namd/cvsroot/namd2/src/WorkDistrib.C,v $
00009  * $Author: bhatele $
00010  * $Date: 2008/08/27 03:10:26 $
00011  * $Revision: 1.1183 $
00012  *****************************************************************************/
00013 
00020 #include <stdio.h>
00021 
00022 #include "InfoStream.h"
00023 #include "ProcessorPrivate.h"
00024 #include "BOCgroup.h"
00025 #include "WorkDistrib.decl.h"
00026 #include "WorkDistrib.h"
00027 
00028 #ifdef USE_COMM_LIB
00029 #include "ComlibManager.h"
00030 #endif
00031 
00032 #include "main.decl.h"
00033 #include "main.h"
00034 #include "Node.h"
00035 #include "PatchMgr.h"
00036 #include "PatchMap.inl"
00037 #include "NamdTypes.h"
00038 #include "PDB.h"
00039 #include "SimParameters.h"
00040 #include "Molecule.h"
00041 #include "NamdOneTools.h"
00042 #include "Compute.h"
00043 #include "ComputeMap.h"
00044 #include "RecBisection.h"
00045 #include "Random.h"
00046 #include "varsizemsg.h"
00047 #include "ProxyMgr.h"
00048 #include "Priorities.h"
00049 
00050 //#define DEBUGM
00051 #define MIN_DEBUG_LEVEL 2
00052 #include "Debug.h"
00053 
00054 
00055 class ComputeMapChangeMsg : public CMessage_ComputeMapChangeMsg
00056 {
00057 public:
00058 
00059   int numNewNodes;
00060   int *newNodes;
00061 
00062 //  VARSIZE_DECL(ComputeMapChangeMsg);
00063 };
00064 
00065 /*
00066 VARSIZE_MSG(ComputeMapChangeMsg,
00067   VARSIZE_ARRAY(newNodes);
00068 )
00069 */
00070 
00071 
00072 //======================================================================
00073 // Public functions
00074 //----------------------------------------------------------------------
00075 WorkDistrib::WorkDistrib()
00076 {
00077   CkpvAccess(BOCclass_group).workDistrib = thisgroup;
00078   mapsArrived = false;
00079   awaitingMaps = false;
00080 }
00081 
00082 //----------------------------------------------------------------------
00083 WorkDistrib::~WorkDistrib(void)
00084 { }
00085 
00086 
00087 //----------------------------------------------------------------------
00088 void WorkDistrib::saveComputeMapChanges(int ep, CkGroupID chareID)
00089 {
00090   saveComputeMapReturnEP = ep;
00091   saveComputeMapReturnChareID = chareID;
00092   saveComputeMapCount = CkNumPes();
00093 
00094   ComputeMap *computeMap = ComputeMap::Object();
00095 
00096   int i;
00097   int nc = computeMap->numComputes();
00098   for (i=0; i<nc; i++) {
00099     DebugM(3, "ComputeMap (" << i << ") node = " << computeMap->node(i) << " newNode = " << computeMap->newNode(i) << "\n");
00100   }
00101   
00102   ComputeMapChangeMsg *mapMsg = new (nc, 0) ComputeMapChangeMsg ;
00103 
00104   mapMsg->numNewNodes = nc;
00105   for(i=0; i<nc; i++)
00106     mapMsg->newNodes[i] = computeMap->newNode(i);
00107 
00108   CProxy_WorkDistrib(thisgroup).recvComputeMapChanges(mapMsg);
00109 }
00110 
00111 void WorkDistrib::recvComputeMapChanges(ComputeMapChangeMsg *msg) {
00112   
00113   ComputeMap *computeMap = ComputeMap::Object();
00114   int i;
00115   for(i=0; i<computeMap->numComputes(); i++)
00116     computeMap->setNewNode(i,msg->newNodes[i]);
00117 
00118   delete msg;
00119 
00120 #if CHARM_VERSION > 050402
00121   CProxy_WorkDistrib workProxy(thisgroup);
00122   workProxy[0].doneSaveComputeMap();
00123 #else
00124   CProxy_WorkDistrib(thisgroup).doneSaveComputeMap(0);
00125 #endif
00126 
00127   DebugM(2, "ComputeMap after send!\n");
00128   for (i=0; i<computeMap->numComputes(); i++) {
00129     DebugM(2, "ComputeMap (" << i << ") node = " << computeMap->node(i) << " newNode = " << computeMap->newNode(i) << " type=" << computeMap->type(i) << "\n");
00130   }
00131   DebugM(2, "===================================================\n");
00132 }
00133 
00134 void WorkDistrib::doneSaveComputeMap() {
00135   if (!--saveComputeMapCount) { 
00136     CkSendMsgBranch(saveComputeMapReturnEP, CkAllocMsg(0,0,0), 0, saveComputeMapReturnChareID);
00137   }
00138 }
00139 
00140 
00141 //only called on node 0
00142 int *WorkDistrib::caclNumAtomsInEachPatch(){
00143     StringList *current;
00144     int i;
00145     CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00146     Node *node = nd.ckLocalBranch();
00147     PatchMap *patchMap = PatchMap::Object();
00148     CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00149     PatchMgr *patchMgr = pm.ckLocalBranch();
00150     SimParameters *params = node->simParameters;
00151     Molecule *molecule = node->molecule;
00152     PDB *pdb = node->pdb;
00153 
00154     int numPatches = patchMap->numPatches();
00155     int numAtoms = pdb->num_atoms();    
00156 
00157     vector<int> *eachPatchAtomList = patchMap->getTmpPatchAtomsList();
00158 
00159     int *patchAtomCnt = new int[numPatches];    
00160     memset(patchAtomCnt, 0, sizeof(int)*numPatches);
00161 
00162     const Lattice lattice = params->lattice;
00163 
00164     Position eachAtomPos;
00165     if (params->splitPatch == SPLIT_PATCH_HYDROGEN)
00166       {
00167       // split atoms into patched based on helix-group and position
00168       int aid, pid=0;
00169       for(i=0; i < numAtoms; i++)
00170         {        
00171         // Assign atoms to patches without splitting hydrogen groups.
00172         // We know that the hydrogenGroup array is sorted with group parents
00173         // listed first.  Thus, only change the pid if an atom is a group parent.
00174         aid = molecule->hydrogenGroup[i].atomID;
00175         pdb->get_position_for_atom(&eachAtomPos, aid);
00176         if (molecule->hydrogenGroup[i].isGP)            
00177             pid = patchMap->assignToPatch(eachAtomPos,lattice);
00178         // else: don't change pid        
00179         patchAtomCnt[pid]++;
00180         eachPatchAtomList[pid].push_back(aid);
00181         }
00182       }
00183     else
00184       {
00185       // split atoms into patched based on position
00186       for(i=0; i < numAtoms; i++)
00187         {
00188         pdb->get_position_for_atom(&eachAtomPos, i);
00189         int pid = patchMap->assignToPatch(eachAtomPos,lattice);        
00190         patchAtomCnt[pid]++;
00191         eachPatchAtomList[pid].push_back(i);
00192         }
00193       }   
00194 
00195     return patchAtomCnt;    
00196 }
00197 
00198 //----------------------------------------------------------------------
00199 // This should only be called on node 0.
00200 //----------------------------------------------------------------------
00201 FullAtomList *WorkDistrib::createAtomLists(void)
00202 {
00203   int i;
00204   StringList *current;  //  Pointer used to retrieve configuration items
00205   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00206   Node *node = nd.ckLocalBranch();
00207   PatchMap *patchMap = PatchMap::Object();
00208   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00209   PatchMgr *patchMgr = pm.ckLocalBranch();
00210   SimParameters *params = node->simParameters;
00211   Molecule *molecule = node->molecule;
00212   PDB *pdb = node->pdb;
00213 
00214   int numPatches = patchMap->numPatches();
00215   int numAtoms = pdb->num_atoms();
00216 
00217   Vector *positions = new Position[numAtoms];
00218   pdb->get_all_positions(positions);
00219 
00220   Vector *velocities = new Velocity[numAtoms];
00221 
00222   if ( params->initialTemp < 0.0 ) {
00223     Bool binvels=FALSE;
00224 
00225     //  Reading the veolcities from a PDB
00226     current = node->configList->find("velocities");
00227 
00228     if (current == NULL) {
00229       current = node->configList->find("binvelocities");
00230       binvels = TRUE;
00231     }
00232 
00233     if (!binvels) {
00234       velocities_from_PDB(current->data, velocities, numAtoms);
00235     }
00236     else {
00237       velocities_from_binfile(current->data, velocities, numAtoms);
00238     }
00239   }
00240   else {
00241     // Random velocities for a given temperature
00242     random_velocities(params->initialTemp, molecule, velocities, numAtoms);
00243   }
00244 
00245   //  If COMMotion == no, remove center of mass motion
00246   if (!(params->comMove)) {
00247     remove_com_motion(velocities, molecule, numAtoms);
00248   }
00249 
00250   FullAtomList *atoms = new FullAtomList[numPatches];
00251 
00252   const Lattice lattice = params->lattice;
00253 
00254   if (params->splitPatch == SPLIT_PATCH_HYDROGEN)
00255     {
00256     // split atoms into patched based on helix-group and position
00257     int aid, pid=0;
00258     for(i=0; i < numAtoms; i++)
00259       {
00260       if ( ! ( i % 1000 ) )
00261         {
00262         DebugM(3,"Assigned " << i << " atoms to patches so far.\n");
00263         }
00264       // Assign atoms to patches without splitting hydrogen groups.
00265       // We know that the hydrogenGroup array is sorted with group parents
00266       // listed first.  Thus, only change the pid if an atom is a group parent.
00267       aid = molecule->hydrogenGroup[i].atomID;
00268       if (molecule->hydrogenGroup[i].isGP)
00269         pid = patchMap->assignToPatch(positions[aid],lattice);
00270       // else: don't change pid
00271       FullAtom a;
00272       a.id = aid;
00273       a.position = positions[aid];
00274       a.velocity = velocities[aid];
00275       #ifdef MEM_OPT_VERSION
00276       a.sigId = molecule->getAtomSigId(aid);
00277       a.exclId = molecule->getAtomExclSigId(aid);
00278       a.vdwType = molecule->atomvdwtype(aid);
00279       #endif
00280       atoms[pid].add(a);
00281       }
00282     }
00283   else
00284     {
00285     // split atoms into patched based on position
00286     for(i=0; i < numAtoms; i++)
00287       {
00288       if ( ! ( i % 1000 ) )
00289         {
00290         DebugM(3,"Assigned " << i << " atoms to patches so far.\n");
00291         }
00292       int pid = patchMap->assignToPatch(positions[i],lattice);
00293       FullAtom a;
00294       a.id = i;
00295       a.position = positions[i];
00296       a.velocity = velocities[i];
00297       #ifdef MEM_OPT_VERSION
00298       a.sigId = molecule->getAtomSigId(i);
00299       a.exclId = molecule->getAtomExclSigId(i);
00300       a.vdwType = molecule->atomvdwtype(i);
00301       #endif
00302       atoms[pid].add(a);
00303       }
00304     }
00305 
00306   delete [] positions;
00307   delete [] velocities;
00308 
00309   for(i=0; i < numPatches; i++)
00310   {
00311     ScaledPosition center(0.5*(patchMap->min_a(i)+patchMap->max_a(i)),
00312                           0.5*(patchMap->min_b(i)+patchMap->max_b(i)),
00313                           0.5*(patchMap->min_c(i)+patchMap->max_c(i)));
00314 
00315     int n = atoms[i].size();
00316     FullAtom *a = atoms[i].begin();
00317     int j;
00318 //Modifications for alchemical fep
00319 //SD & CC, CNRS - LCTN, Nancy
00320     Bool fepOn = params->fepOn;
00321     Bool thermInt = params->thermInt;
00322 //fepe
00323     Bool lesOn = params->lesOn;
00324   
00325     Bool pairInteractionOn = params->pairInteractionOn;
00326 
00327     Bool pressureProfileTypes = (params->pressureProfileAtomTypes > 1);
00328 
00329     Transform mother_transform;
00330     for(j=0; j < n; j++)
00331     {
00332       int aid = a[j].id;
00333 
00334       if (params->splitPatch == SPLIT_PATCH_HYDROGEN) {
00335         if ( molecule->is_hydrogenGroupParent(aid) ) {
00336           a[j].hydrogenGroupSize = molecule->get_groupSize(aid);
00337         } else {
00338           a[j].hydrogenGroupSize = 0;
00339         }
00340       } else {
00341         a[j].hydrogenGroupSize = 1;
00342       }
00343 
00344       a[j].nonbondedGroupIsAtom = 0;
00345 
00346       a[j].atomFixed = molecule->is_atom_fixed(aid) ? 1 : 0;
00347       a[j].fixedPosition = a[j].position;
00348 
00349       if ( a[j].hydrogenGroupSize ) {
00350         a[j].position = lattice.nearest(
00351                 a[j].position, center, &(a[j].transform));
00352         mother_transform = a[j].transform;
00353       } else {
00354         a[j].position = lattice.apply_transform(a[j].position,mother_transform);
00355         a[j].transform = mother_transform;
00356       }
00357 
00358       a[j].mass = molecule->atommass(aid);
00359       a[j].charge = molecule->atomcharge(aid);
00360 
00361 //Modifications for alchemical fep
00362 //SD & CC, CNRS - LCTN, Nancy
00363       if ( fepOn || thermInt || lesOn || pairInteractionOn || pressureProfileTypes) {
00364         a[j].partition = molecule->get_fep_type(aid);
00365       } 
00366       else {
00367         a[j].partition = 0;
00368       }
00369 //fepe
00370 
00371     }
00372 
00373     int size, allfixed, k;
00374     for(j=0; j < n; j+=size) {
00375       size = a[j].hydrogenGroupSize;
00376       if ( ! size ) {
00377         NAMD_bug("Mother atom with hydrogenGroupSize of 0!");
00378       }
00379       allfixed = 1;
00380       for ( k = 0; k < size; ++k ) {
00381         allfixed = ( allfixed && (a[j+k].atomFixed) );
00382       }
00383       for ( k = 0; k < size; ++k ) {
00384         a[j+k].groupFixed = allfixed ? 1 : 0;
00385       }
00386     }
00387 
00388     if ( params->outputPatchDetails ) {
00389       int patchId = i;
00390       int numAtomsInPatch = n;
00391       int numFixedAtomsInPatch = 0;
00392       int numAtomsInFixedGroupsInPatch = 0;
00393       for(j=0; j < n; j++) {
00394         numFixedAtomsInPatch += ( a[j].atomFixed ? 1 : 0 );
00395         numAtomsInFixedGroupsInPatch += ( a[j].groupFixed ? 1 : 0 );
00396       }
00397       iout << "PATCH_DETAILS:"
00398            << " patch " << patchId
00399            << " atoms " << numAtomsInPatch
00400            << " fixed_atoms " << numFixedAtomsInPatch
00401            << " fixed_groups " << numAtomsInFixedGroupsInPatch
00402            << "\n" << endi;
00403     }
00404   }
00405 
00406   return atoms;
00407 
00408 }
00409 
00410 //This should only be called on node 0
00411 //This differs from creatHomePatches in: create home patches
00412 //without populating them with actual atoms' data. The only field
00413 //set is the number of atoms each patch contains.
00414 void WorkDistrib::preCreateHomePatches(){
00415 
00416   PatchMap *patchMap = PatchMap::Object();
00417   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00418   PatchMgr *patchMgr = pm.ckLocalBranch();
00419 
00420   int numPatches = patchMap->numPatches();
00421 
00422   patchMap->initTmpPatchAtomsList();
00423 
00424   int *patchAtomCnt = caclNumAtomsInEachPatch();
00425 
00426   int maxAtoms = -1;
00427   int maxPatch = -1;
00428   for(int i=0; i < numPatches; i++) {
00429     int numAtoms = patchAtomCnt[i];
00430     if ( numAtoms > maxAtoms ) { maxAtoms = numAtoms; maxPatch = i; }
00431   }
00432   iout << iINFO << "LARGEST PATCH (" << maxPatch <<
00433         ") HAS " << maxAtoms << " ATOMS\n" << endi;
00434 
00435   for(int i=0; i < numPatches; i++)
00436   {
00437     if ( ! ( i % 100 ) )
00438     {
00439       DebugM(3,"Pre-created " << i << " patches so far.\n");
00440     }
00441 
00442     patchMgr->preCreateHomePatch(i,patchAtomCnt[i]);
00443   }
00444 
00445   delete [] patchAtomCnt;
00446 }
00447 
00448 
00449 //----------------------------------------------------------------------
00450 // This should only be called on node 0.
00451 //----------------------------------------------------------------------
00452 void WorkDistrib::createHomePatches(void)
00453 {
00454   int i;
00455   PatchMap *patchMap = PatchMap::Object();
00456   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00457   PatchMgr *patchMgr = pm.ckLocalBranch();
00458 
00459   int numPatches = patchMap->numPatches();
00460 
00461   FullAtomList *atoms = createAtomLists();
00462     
00463 #ifdef MEM_OPT_VERSION
00464 /*  CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00465   Node *node = nd.ckLocalBranch();
00466   node->molecule->delEachAtomSigs();
00467   node->molecule->delMassChargeSpace();
00468 */
00469 #endif
00470 
00471   int maxAtoms = -1;
00472   int maxPatch = -1;
00473   for(i=0; i < numPatches; i++) {
00474     int numAtoms = atoms[i].size();
00475     if ( numAtoms > maxAtoms ) { maxAtoms = numAtoms; maxPatch = i; }
00476   }
00477   iout << iINFO << "LARGEST PATCH (" << maxPatch <<
00478         ") HAS " << maxAtoms << " ATOMS\n" << endi;
00479 
00480   for(i=0; i < numPatches; i++)
00481   {
00482     if ( ! ( i % 100 ) )
00483     {
00484       DebugM(3,"Created " << i << " patches so far.\n");
00485     }
00486 
00487     patchMgr->createHomePatch(i,atoms[i]);
00488   }
00489 
00490   delete [] atoms;
00491 }
00492 
00493 void WorkDistrib::fillOnePatchAtoms(int patchId, FullAtomList *onePatchAtoms, Vector *velocities){
00494     // ref BOC
00495     CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00496     Node *node = nd.ckLocalBranch();
00497     PDB *pdb = node->pdb;
00498     SimParameters *params = node->simParameters;
00499     Molecule *molecule = node->molecule;
00500 
00501     const Lattice lattice = params->lattice;
00502 
00503     CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00504     PatchMgr *patchMgr = pm.ckLocalBranch();
00505     // ref singleton
00506     PatchMap *patchMap = PatchMap::Object();
00507 
00508     vector<int> *eachPatchAtomsList = patchMap->getTmpPatchAtomsList();
00509     vector<int> *thisPatchAtomsList = &eachPatchAtomsList[patchId];
00510 
00511     for(int i=0; i<thisPatchAtomsList->size(); i++){
00512         int aid = thisPatchAtomsList->at(i);
00513         FullAtom a;
00514         a.id = aid;
00515         Position pos;
00516         pdb->get_position_for_atom(&pos, aid);
00517         a.position = pos;
00518         a.velocity = velocities[aid];
00519         #ifdef MEM_OPT_VERSION
00520         a.sigId = molecule->getAtomSigId(aid);
00521         a.exclId = molecule->getAtomExclSigId(aid);
00522         a.vdwType = molecule->atomvdwtype(aid);
00523         #endif
00524         onePatchAtoms->add(a);
00525     }
00526 
00527     ScaledPosition center(0.5*(patchMap->min_a(patchId)+patchMap->max_a(patchId)),
00528                           0.5*(patchMap->min_b(patchId)+patchMap->max_b(patchId)),
00529                           0.5*(patchMap->min_c(patchId)+patchMap->max_c(patchId)));
00530     int n = onePatchAtoms->size();
00531     FullAtom *a = onePatchAtoms->begin();
00532     int j;
00533 //Modifications for alchemical fep
00534 //SD & CC, CNRS - LCTN, Nancy
00535     Bool fepOn = params->fepOn;
00536     Bool thermInt = params->thermInt;
00537 //fepe
00538     Bool lesOn = params->lesOn;
00539   
00540     Bool pairInteractionOn = params->pairInteractionOn;
00541 
00542     Bool pressureProfileTypes = (params->pressureProfileAtomTypes > 1);
00543 
00544     Transform mother_transform;
00545     for(j=0; j < n; j++)
00546     {
00547       int aid = a[j].id;
00548 
00549       if (params->splitPatch == SPLIT_PATCH_HYDROGEN) {
00550         if ( molecule->is_hydrogenGroupParent(aid) ) {
00551           a[j].hydrogenGroupSize = molecule->get_groupSize(aid);
00552         } else {
00553           a[j].hydrogenGroupSize = 0;
00554         }
00555       } else {
00556         a[j].hydrogenGroupSize = 1;
00557       }
00558 
00559       a[j].nonbondedGroupIsAtom = 0;
00560 
00561       a[j].atomFixed = molecule->is_atom_fixed(aid) ? 1 : 0;
00562       a[j].fixedPosition = a[j].position;
00563 
00564       if ( a[j].hydrogenGroupSize ) {
00565         a[j].position = lattice.nearest(
00566                 a[j].position, center, &(a[j].transform));
00567         mother_transform = a[j].transform;
00568       } else {
00569         a[j].position = lattice.apply_transform(a[j].position,mother_transform);
00570         a[j].transform = mother_transform;
00571       }
00572 
00573       a[j].mass = molecule->atommass(aid);
00574       a[j].charge = molecule->atomcharge(aid);
00575 
00576 //Modifications for alchemical fep
00577 //SD & CC, CNRS - LCTN, Nancy
00578       if ( fepOn || thermInt || lesOn || pairInteractionOn || pressureProfileTypes) {
00579         a[j].partition = molecule->get_fep_type(aid);
00580       } 
00581       else {
00582         a[j].partition = 0;
00583       }
00584 //fepe
00585     }
00586 
00587     int size, allfixed, k;
00588     for(j=0; j < n; j+=size) {
00589       size = a[j].hydrogenGroupSize;
00590       if ( ! size ) {
00591         NAMD_bug("Mother atom with hydrogenGroupSize of 0!");
00592       }
00593       allfixed = 1;
00594       for ( k = 0; k < size; ++k ) {
00595         allfixed = ( allfixed && (a[j+k].atomFixed) );
00596       }
00597       for ( k = 0; k < size; ++k ) {
00598         a[j+k].groupFixed = allfixed ? 1 : 0;
00599       }
00600     }
00601 
00602     if(params->fixedAtomsOn){
00603         int fixedCnt=0;
00604         for(j=0; j<n; j++)
00605             fixedCnt += (a[j].atomFixed ? 1:0);
00606         patchMgr->setHomePatchFixedAtomNum(patchId, fixedCnt);
00607     }    
00608 
00609     if ( params->outputPatchDetails ) {    
00610       int numAtomsInPatch = n;
00611       int numFixedAtomsInPatch = 0;
00612       int numAtomsInFixedGroupsInPatch = 0;
00613       for(j=0; j < n; j++) {
00614         numFixedAtomsInPatch += ( a[j].atomFixed ? 1 : 0 );
00615         numAtomsInFixedGroupsInPatch += ( a[j].groupFixed ? 1 : 0 );
00616       }
00617       iout << "PATCH_DETAILS:"
00618            << " patch " << patchId
00619            << " atoms " << numAtomsInPatch
00620            << " fixed_atoms " << numFixedAtomsInPatch
00621            << " fixed_groups " << numAtomsInFixedGroupsInPatch
00622            << "\n" << endi;
00623     }
00624 }
00625 
00626 //should be called only on node 0
00627 void WorkDistrib::initAndSendHomePatch(){
00628   StringList *current;
00629   // ref BOC
00630   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00631   Node *node = nd.ckLocalBranch();
00632   Molecule *molecule = node->molecule;
00633   PDB *pdb = node->pdb;
00634   SimParameters *params = node->simParameters;
00635 
00636   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00637   PatchMgr *patchMgr = pm.ckLocalBranch();
00638   // ref singleton
00639   PatchMap *patchMap = PatchMap::Object();
00640 
00641   int numAtoms = pdb->num_atoms();
00642 
00643   //1. create atoms' velocities
00644   Vector *velocities = new Velocity[numAtoms];
00645   if ( params->initialTemp < 0.0 ) {
00646     Bool binvels=FALSE;
00647 
00648     //  Reading the veolcities from a PDB
00649     current = node->configList->find("velocities");
00650 
00651     if (current == NULL) {
00652       current = node->configList->find("binvelocities");
00653       binvels = TRUE;
00654     }
00655 
00656     if (!binvels) {
00657       velocities_from_PDB(current->data, velocities, numAtoms);
00658     }
00659     else {
00660       velocities_from_binfile(current->data, velocities, numAtoms);
00661     }
00662   }
00663   else {
00664     // Random velocities for a given temperature
00665     random_velocities(params->initialTemp, molecule, velocities, numAtoms);
00666   }
00667 
00668   //  If COMMotion == no, remove center of mass motion
00669   if (!(params->comMove)) {
00670     remove_com_motion(velocities, molecule, numAtoms);
00671   }
00672 
00673 
00674   for(int i=0; i<patchMap->numPatches(); i++){
00675       //2. fill each patch with actual atom data
00676       //the work flow should looks like the createAtomList
00677       FullAtomList *onePatchAtoms = new FullAtomList;
00678       fillOnePatchAtoms(i, onePatchAtoms, velocities);
00679       
00680       patchMgr->fillHomePatchAtomList(i, onePatchAtoms);
00681 
00682       delete onePatchAtoms;
00683 
00684       //distribute this home patch
00685       if(patchMap->node(i) != node->myid()){
00686           //need to move to other nodes
00687           patchMgr->sendOneHomePatch(i, patchMap->node(i));
00688       }
00689       
00690   }  
00691 
00692   delete [] velocities;
00693   patchMap->delTmpPatchAtomsList();
00694 }
00695 
00696 void WorkDistrib::distributeHomePatches() {
00697   // ref BOC
00698   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00699   Node *node = nd.ckLocalBranch();
00700   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00701   PatchMgr *patchMgr = pm.ckLocalBranch();
00702   // ref singleton
00703   PatchMap *patchMap = PatchMap::Object();
00704 
00705   // Move patches to the proper node
00706   for(int i=0;i < patchMap->numPatches(); i++)
00707   {
00708     if (patchMap->node(i) != node->myid() )
00709     {
00710       DebugM(3,"patchMgr->movePatch("
00711         << i << "," << patchMap->node(i) << ")\n");
00712       patchMgr->movePatch(i,patchMap->node(i));
00713     }
00714   }
00715   patchMgr->sendMovePatches();
00716 }
00717 
00718 void WorkDistrib::reinitAtoms() {
00719 
00720   PatchMap *patchMap = PatchMap::Object();
00721   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00722   PatchMgr *patchMgr = pm.ckLocalBranch();
00723 
00724   int numPatches = patchMap->numPatches();
00725 
00726   FullAtomList *atoms = createAtomLists();
00727 
00728   for(int i=0; i < numPatches; i++) {
00729     patchMgr->sendAtoms(i,atoms[i]);
00730   }
00731 
00732   delete [] atoms;
00733 
00734 }
00735 
00736 
00737 //----------------------------------------------------------------------
00738 
00739 class MapDistribMsg: public CMessage_MapDistribMsg {
00740   public:
00741     char *patchMapData;
00742     char *computeMapData;
00743 
00744 //  VARSIZE_DECL(MapDistribMsg);
00745 };
00746 
00747 /*
00748 VARSIZE_MSG(MapDistribMsg,
00749   VARSIZE_ARRAY(patchMapData);
00750   VARSIZE_ARRAY(computeMapData);
00751 )
00752 */
00753 
00754 void WorkDistrib::sendMaps(void)
00755 {
00756   if ( CkNumPes() == 1 ) {
00757     mapsArrived = true;
00758     return;
00759   }
00760 
00761   //Automatically enable spanning tree
00762   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00763   Node *node = nd.ckLocalBranch();
00764   SimParameters *params = node->simParameters;
00765   if(PatchMap::Object()->numPatches() <= CkNumPes()/4 && params->isSendProxySTAuto())
00766     ProxyMgr::Object()->setSendSpanning();
00767 
00768   int sizes[2];
00769   sizes[0] = PatchMap::Object()->packSize();
00770   sizes[1] = ComputeMap::Object()->packSize();
00771 
00772   MapDistribMsg *mapMsg = new (sizes[0], sizes[1], 0) MapDistribMsg;
00773 
00774   PatchMap::Object()->pack(mapMsg->patchMapData);
00775   ComputeMap::Object()->pack(mapMsg->computeMapData);
00776 
00777 #if CHARM_VERSION > 050402
00778   CProxy_WorkDistrib workProxy(thisgroup);
00779   workProxy[0].saveMaps(mapMsg);
00780 #else
00781   CProxy_WorkDistrib(thisgroup).saveMaps(mapMsg,0);
00782 #endif
00783 }
00784 
00785 // saveMaps() is called when the map message is received
00786 void WorkDistrib::saveMaps(MapDistribMsg *msg)
00787 {
00788   // Use a resend to forward messages before processing.  Otherwise the
00789   // map distribution is slow on many CPUs.  We need to use a tree
00790   // rather than a broadcast because some implementations of broadcast
00791   // generate a copy of the message on the sender for each recipient.
00792   // This is because MPI doesn't allow re-use of an outstanding buffer.
00793 
00794   if ( mapsArrived && CkMyPe() ) {
00795     PatchMap::Object()->unpack(msg->patchMapData);
00796     ComputeMap::Object()->unpack(msg->computeMapData);
00797 
00798     //Automatically enable spanning tree
00799     CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00800     Node *node = nd.ckLocalBranch();
00801     SimParameters *params = node->simParameters;
00802     if(PatchMap::Object()->numPatches() <= CkNumPes()/4 && params->isSendProxySTAuto())
00803       ProxyMgr::Object()->setSendSpanning();
00804   }
00805   if ( mapsArrived ) {
00806     delete msg;
00807     return;
00808   }
00809 
00810   mapsArrived = true;
00811 
00812   int pids[3];
00813   int basePe = 2 * CkMyPe() + 1;
00814   int npid = 0;
00815   if ( (basePe+npid) < CkNumPes() ) { pids[npid] = basePe + npid; ++npid; }
00816   if ( (basePe+npid) < CkNumPes() ) { pids[npid] = basePe + npid; ++npid; }
00817   pids[npid] = CkMyPe(); ++npid;  // always send the message to ourselves
00818   CProxy_WorkDistrib(thisgroup).saveMaps(msg,npid,pids);
00819 }
00820 
00821 
00822 //----------------------------------------------------------------------
00823 void WorkDistrib::patchMapInit(void)
00824 {
00825   PatchMap *patchMap = PatchMap::Object();
00826   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00827   Node *node = nd.ckLocalBranch();
00828   SimParameters *params = node->simParameters;
00829   Lattice lattice = params->lattice;
00830 
00831   BigReal patchSize = params->patchDimension;
00832 
00833   ScaledPosition xmin, xmax;
00834   ScaledPosition sysDim, sysMin;
00835 
00836   double maxNumPatches = 1000000;  // need to adjust fractional values
00837   if ( params->minAtomsPerPatch > 0 )
00838     maxNumPatches = node->pdb->num_atoms() / params->minAtomsPerPatch;
00839 
00840   DebugM(3,"Mapping patches\n");
00841   // Need to use full box for FMA to match NAMD 1.X results.
00842   if ( params->FMAOn ) {
00843     node->pdb->find_extremes(&(xmin.x),&(xmax.x),lattice.a_r());
00844     node->pdb->find_extremes(&(xmin.y),&(xmax.y),lattice.b_r());
00845     node->pdb->find_extremes(&(xmin.z),&(xmax.z),lattice.c_r());
00846   // Otherwise, this allows a small number of stray atoms.
00847   } else {
00848     node->pdb->find_extremes(&(xmin.x),&(xmax.x),lattice.a_r(),0.9);
00849     node->pdb->find_extremes(&(xmin.y),&(xmax.y),lattice.b_r(),0.9);
00850     node->pdb->find_extremes(&(xmin.z),&(xmax.z),lattice.c_r(),0.9);
00851   }
00852 
00853   BigReal origin_shift;
00854   origin_shift = lattice.a_r() * lattice.origin();
00855   xmin.x -= origin_shift;
00856   xmax.x -= origin_shift;
00857   origin_shift = lattice.b_r() * lattice.origin();
00858   xmin.y -= origin_shift;
00859   xmax.y -= origin_shift;
00860   origin_shift = lattice.c_r() * lattice.origin();
00861   xmin.z -= origin_shift;
00862   xmax.z -= origin_shift;
00863 
00864   // SimParameters default is -1 for unset
00865   int twoAwayX = params->twoAwayX;
00866   int twoAwayY = params->twoAwayY;
00867   int twoAwayZ = params->twoAwayZ;
00868 
00869   int numPatches = patchMap->sizeGrid(
00870         xmin,xmax,lattice,patchSize,maxNumPatches,
00871         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00872   if ( numPatches > (CkNumPes() - 1) && twoAwayZ < 0 ) {
00873     twoAwayZ = 0;
00874     numPatches = patchMap->sizeGrid(
00875         xmin,xmax,lattice,patchSize,maxNumPatches,
00876         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00877   }
00878   if ( numPatches > (CkNumPes() - 1) && twoAwayY < 0 ) {
00879     twoAwayY = 0;
00880     numPatches = patchMap->sizeGrid(
00881         xmin,xmax,lattice,patchSize,maxNumPatches,
00882         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00883   }
00884   if ( numPatches > (CkNumPes() - 1) && twoAwayX < 0 ) {
00885     twoAwayX = 0;
00886     numPatches = patchMap->sizeGrid(
00887         xmin,xmax,lattice,patchSize,maxNumPatches,
00888         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00889   }
00890 
00891   patchMap->makePatches(xmin,xmax,lattice,patchSize,maxNumPatches,
00892         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00893 
00894 }
00895 
00896 
00897 //----------------------------------------------------------------------
00898 void WorkDistrib::assignNodeToPatch()
00899 {
00900   int method=1;
00901 
00902   PatchMap *patchMap = PatchMap::Object();
00903   int nNodes = Node::Object()->numNodes();
00904 
00905 #if USE_TOPOMAP 
00906   TopoManager tmgr;
00907   int numPes = tmgr.getDimX() * tmgr.getDimY() * tmgr.getDimZ();
00908   if (numPes > patchMap->numPatches() && (assignPatchesTopoGridRecBisection() > 0)) {
00909     CkPrintf ("Blue Gene/L topology partitioner finished successfully \n");
00910   }
00911   else
00912 #endif
00913     if (nNodes > patchMap->numPatches())
00914       assignPatchesBitReversal();
00915   // else if (nNodes == patchMap->numPatches())
00916   //   assignPatchesRoundRobin();
00917     else if (method==1)
00918       assignPatchesRecursiveBisection();
00919     else
00920       assignPatchesToLowestLoadNode();
00921   
00922   int *nAtoms = new int[nNodes];
00923   int numAtoms=0;
00924   int i;
00925   for(i=0; i < nNodes; i++)
00926     nAtoms[i] = 0;
00927 
00928   for(i=0; i < patchMap->numPatches(); i++)
00929   {
00930     //    iout << iINFO << "Patch " << i << " has " 
00931     //   << patchMap->patch(i)->getNumAtoms() << " atoms and "
00932     //   << patchMap->patch(i)->getNumAtoms() * 
00933     //            patchMap->patch(i)->getNumAtoms() 
00934     //   << " pairs.\n" << endi;
00935 
00936     if (patchMap->patch(i)) {
00937       numAtoms += patchMap->patch(i)->getNumAtoms();
00938       nAtoms[patchMap->node(i)] += patchMap->patch(i)->getNumAtoms();
00939     }
00940   }
00941 
00942 //  for(i=0; i < nNodes; i++)
00943 //    iout << iINFO 
00944 //       << nAtoms[i] << " atoms assigned to node " << i << "\n" << endi;
00945   if ( numAtoms != Node::Object()->molecule->numAtoms ) {
00946     NAMD_die("Incorrect atom count in WorkDistrib::assignNodeToPatch\n");
00947   }
00948 
00949   delete [] nAtoms;
00950  
00951   //  PatchMap::Object()->printPatchMap();
00952 }
00953 
00954 //----------------------------------------------------------------------
00955 // void WorkDistrib::assignPatchesSlices() 
00956 // {
00957 //   int pid; 
00958 //   int assignedNode = 0;
00959 //   PatchMap *patchMap = PatchMap::Object();
00960 //   Node *node = CLocalBranch(Node, CkpvAccess(BOCclass_group).node);
00961 
00962 //   int *numAtoms = new int[node->numNodes()];
00963 //   for (int i=0; i<node->numNodes(); i++) {
00964 //     numAtoms[i] = 0;
00965 //   }
00966 
00967 //   // Assign patch to node with least atoms assigned.
00968 //   for(pid=0; pid < patchMap->numPatches(); pid++) {
00969 //     assignedNode = 0;
00970 //     for (i=1; i < node->numNodes(); i++) {
00971 //       if (numAtoms[i] < numAtoms[assignedNode]) assignedNode = i;
00972 //     }
00973 //     patchMap->assignNode(pid, assignedNode);
00974 //     numAtoms[assignedNode] += patchMap->patch(pid)->getNumAtoms();
00975 
00976 //     /*
00977 //     iout << iINFO << "Patch (" << pid << ") has " 
00978 //       << patchMap->patch(pid)->getNumAtoms() 
00979 //       << " atoms:  Assigned to Node(" << assignedNode << ")\n" 
00980 //       << endi;
00981 //     */
00982 //   }
00983 
00984 //   delete[] numAtoms;
00985 // }
00986 
00987 //----------------------------------------------------------------------
00988 void WorkDistrib::assignPatchesToLowestLoadNode() 
00989 {
00990   int pid; 
00991   int assignedNode = 0;
00992   PatchMap *patchMap = PatchMap::Object();
00993   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00994   Node *node = nd.ckLocalBranch();
00995   SimParameters *simParams = node->simParameters;
00996 
00997   int *load = new int[node->numNodes()];
00998   int *assignedNodes = new int[patchMap->numPatches()];
00999   for (int i=0; i<node->numNodes(); i++) {
01000     load[i] = 0;
01001   }
01002 
01003   int defaultNode = 0;
01004   if ( simParams->noPatchesOnZero && node->numNodes() > 1 ){
01005     defaultNode = 1;
01006     if( simParams->noPatchesOnOne && node->numNodes() > 2)
01007       defaultNode = 2;
01008   }
01009   // Assign patch to node with least atoms assigned.
01010   for(pid=0; pid < patchMap->numPatches(); pid++) {
01011     assignedNode = defaultNode;
01012     for (int i=assignedNode + 1; i < node->numNodes(); i++) {
01013       if (load[i] < load[assignedNode]) assignedNode = i;
01014     }
01015     assignedNodes[pid] = assignedNode;
01016     load[assignedNode] += patchMap->patch(pid)->getNumAtoms() + 1;
01017   }
01018 
01019   delete[] load;
01020   sortNodesAndAssign(assignedNodes);
01021   delete[] assignedNodes;
01022 }
01023 
01024 //----------------------------------------------------------------------
01025 void WorkDistrib::assignPatchesBitReversal() 
01026 {
01027   int pid; 
01028   PatchMap *patchMap = PatchMap::Object();
01029   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01030   Node *node = nd.ckLocalBranch();
01031 
01032   int ncpus = node->numNodes();
01033   int npatches = patchMap->numPatches();
01034   if ( ncpus <= npatches )
01035     NAMD_bug("WorkDistrib::assignPatchesBitReversal called improperly");
01036 
01037   // find next highest power of two
01038   int npow2 = 1;  int nbits = 0;
01039   while ( npow2 < ncpus ) { npow2 *= 2; nbits += 1; }
01040 
01041   // build bit reversal sequence
01042   SortableResizeArray<int> seq(ncpus);
01043   // avoid using node 0 (reverse of 0 is 0 so start at 1)
01044   int i = 1;
01045   for ( int icpu=0; icpu<(ncpus-1); ++icpu ) {
01046     int ri;
01047     for ( ri = ncpus; ri >= ncpus; ++i ) {
01048       ri = 0;
01049       int pow2 = 1;
01050       int rpow2 = npow2 / 2;
01051       for ( int j=0; j<nbits; ++j ) {
01052         ri += rpow2 * ( ( i / pow2 ) % 2 );
01053         pow2 *= 2;  rpow2 /= 2;
01054       }
01055     }
01056     seq[icpu] = ri;
01057   }
01058 
01059   // extract and sort patch locations
01060   sortNodesAndAssign(seq.begin());
01061   if ( ncpus > 2*npatches ) sortNodesAndAssign(seq.begin()+npatches, 1);
01062 }
01063 
01064 //----------------------------------------------------------------------
01065 struct nodesort {
01066   int node;
01067   int a_total;
01068   int b_total;
01069   int c_total;
01070   int npatches;
01071   nodesort() : node(-1),a_total(0),b_total(0),c_total(0),npatches(0) { ; }
01072   int operator==(const nodesort &o) const {
01073     float a1 = ((float)a_total)/((float)npatches);
01074     float a2 = ((float)o.a_total)/((float)o.npatches);
01075     float b1 = ((float)b_total)/((float)npatches);
01076     float b2 = ((float)o.b_total)/((float)o.npatches);
01077     float c1 = ((float)c_total)/((float)npatches);
01078     float c2 = ((float)o.c_total)/((float)o.npatches);
01079     return ((a1 == a2) && (b1 == b2) && (c1 == c2));
01080   }
01081   int operator<(const nodesort &o) const {
01082     float a1 = ((float)a_total)/((float)npatches);
01083     float a2 = ((float)o.a_total)/((float)o.npatches);
01084     float b1 = ((float)b_total)/((float)npatches);
01085     float b2 = ((float)o.b_total)/((float)o.npatches);
01086     float c1 = ((float)c_total)/((float)npatches);
01087     float c2 = ((float)o.c_total)/((float)o.npatches);
01088     return ( (a1 < a2) || ((a1 == a2) && (b1 < b2)) ||
01089                 ((a1 == a2) && (b1 == b2) && (c1 < c2)) );
01090   }
01091 };
01092 
01093 
01094 void WorkDistrib::sortNodesAndAssign(int *assignedNode, int baseNodes) {
01095   // if baseNodes is zero (default) then set both nodes and basenodes
01096   // if baseNodes is nonzero then this is a second call to set basenodes only
01097   int i, pid; 
01098   PatchMap *patchMap = PatchMap::Object();
01099   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01100   Node *node = nd.ckLocalBranch();
01101   int nnodes = node->numNodes();
01102   int npatches = patchMap->numPatches();
01103 
01104   ResizeArray<nodesort> allnodes(nnodes);
01105   for ( i=0; i < nnodes; ++i ) {
01106     allnodes[i].node = i;
01107   }
01108   for ( pid=0; pid<npatches; ++pid ) {
01109     // iout << pid << " " << assignedNode[pid] << "\n" << endi;
01110     allnodes[assignedNode[pid]].npatches++;
01111     allnodes[assignedNode[pid]].a_total += patchMap->index_a(pid);
01112     allnodes[assignedNode[pid]].b_total += patchMap->index_b(pid);
01113     allnodes[assignedNode[pid]].c_total += patchMap->index_c(pid);
01114   }
01115   SortableResizeArray<nodesort> usednodes(nnodes);
01116   usednodes.resize(0);
01117   for ( i=0; i < nnodes; ++i ) {
01118     if ( allnodes[i].npatches ) usednodes.add(allnodes[i]);
01119   }
01120   usednodes.sort();
01121   int nused = usednodes.size();
01122   int i2 = nused/2;
01123   for ( i=0; i < nnodes; ++i ) {
01124     if ( allnodes[i].npatches ) allnodes[usednodes[i2++].node].node = i;
01125     if ( i2 == nused ) i2 = 0;
01126   }
01127 
01128   for ( pid=0; pid<npatches; ++pid ) {
01129     // iout << pid << " " <<  allnodes[assignedNode[pid]].node << "\n" << endi;
01130     if ( ! baseNodes ) {
01131       patchMap->assignNode(pid, allnodes[assignedNode[pid]].node);
01132     }
01133     patchMap->assignBaseNode(pid, allnodes[assignedNode[pid]].node);
01134   }
01135 }
01136 
01137 //----------------------------------------------------------------------
01138 void WorkDistrib::assignPatchesRoundRobin() 
01139 {
01140   int pid; 
01141   PatchMap *patchMap = PatchMap::Object();
01142   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01143   Node *node = nd.ckLocalBranch();
01144   int *assignedNode = new int[patchMap->numPatches()];
01145 
01146   for(pid=0; pid < patchMap->numPatches(); pid++) {
01147     assignedNode[pid] = pid % node->numNodes();
01148   }
01149 
01150   sortNodesAndAssign(assignedNode);
01151   delete [] assignedNode;
01152 }
01153 
01154 //----------------------------------------------------------------------
01155 void WorkDistrib::assignPatchesRecursiveBisection() 
01156 {
01157   PatchMap *patchMap = PatchMap::Object();
01158   int *assignedNode = new int[patchMap->numPatches()];
01159   int numNodes = Node::Object()->numNodes();
01160   SimParameters *simParams = Node::Object()->simParameters;
01161   int usedNodes = numNodes;
01162   int unusedNodes = 0;
01163   if ( simParams->noPatchesOnZero && numNodes > 1 ){
01164     usedNodes -= 1;
01165     if(simParams->noPatchesOnOne && numNodes > 2)
01166       usedNodes -= 1;
01167   }  
01168   unusedNodes = numNodes - usedNodes;
01169   RecBisection recBisec(usedNodes,PatchMap::Object());
01170   if ( recBisec.partition(assignedNode) ) {
01171     if ( unusedNodes !=0 ) {
01172       for ( int i=0; i<patchMap->numPatches(); ++i ) {
01173         assignedNode[i] += unusedNodes;
01174       }
01175     }
01176     sortNodesAndAssign(assignedNode);
01177     delete [] assignedNode;
01178   } else {
01179     //free the array here since a same array will be allocated
01180     //in assignPatchesToLowestLoadNode function, thus reducting
01181     //temporary memory usage
01182     delete [] assignedNode; 
01183     
01184     iout << iWARN 
01185          << "WorkDistrib: Recursive bisection fails,"
01186          << "invoking least-load algorithm\n";
01187     assignPatchesToLowestLoadNode();
01188   }
01189 }
01190 
01191 //----------------------------------------------------------------------
01192 void WorkDistrib::mapComputes(void)
01193 {
01194   PatchMap *patchMap = PatchMap::Object();
01195   ComputeMap *computeMap = ComputeMap::Object();
01196   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01197   Node *node = nd.ckLocalBranch();
01198 
01199   DebugM(3,"Mapping computes\n");
01200 
01201   computeMap->allocateCids();
01202 
01203   // Handle full electrostatics
01204   if ( node->simParameters->fullDirectOn )
01205     mapComputeHomePatches(computeFullDirectType);
01206   if ( node->simParameters->FMAOn )
01207 #ifdef DPMTA
01208     mapComputeHomePatches(computeDPMTAType);
01209 #else
01210     NAMD_die("This binary does not include DPMTA (FMA).");
01211 #endif
01212   if ( node->simParameters->PMEOn ) {
01213 #ifdef DPME
01214     if ( node->simParameters->useDPME )
01215       mapComputeHomePatches(computeDPMEType);
01216     else {
01217       mapComputeHomePatches(computePmeType);
01218       if ( node->simParameters->pressureProfileEwaldOn )
01219         mapComputeHomePatches(computeEwaldType);
01220     }
01221 #else
01222     mapComputeHomePatches(computePmeType);
01223     if ( node->simParameters->pressureProfileEwaldOn )
01224       mapComputeHomePatches(computeEwaldType);
01225 #endif
01226   }
01227 
01228   if ( node->simParameters->globalForcesOn ) {
01229     DebugM(2,"adding ComputeGlobal\n");
01230     mapComputeHomePatches(computeGlobalType);
01231   }
01232 
01233   if ( node->simParameters->extForcesOn )
01234     mapComputeHomePatches(computeExtType);
01235 
01236   mapComputeNonbonded();
01237 
01238   // If we're doing true pair interactions, no need for bonded terms.
01239   // But if we're doing within-group interactions, we do need them.
01240   if ( !node->simParameters->pairInteractionOn || 
01241       node->simParameters->pairInteractionSelf) { 
01242     mapComputeHomeTuples(computeBondsType);
01243     mapComputeHomeTuples(computeAnglesType);
01244     mapComputeHomeTuples(computeDihedralsType);
01245     mapComputeHomeTuples(computeImpropersType);
01246     mapComputeHomeTuples(computeCrosstermsType);
01247     mapComputePatch(computeSelfBondsType);
01248     mapComputePatch(computeSelfAnglesType);
01249     mapComputePatch(computeSelfDihedralsType);
01250     mapComputePatch(computeSelfImpropersType);
01251     mapComputePatch(computeSelfCrosstermsType);
01252   }
01253 
01254 
01255   if ( node->simParameters->eFieldOn )
01256     mapComputePatch(computeEFieldType);
01257   /* BEGIN gf */
01258   if ( node->simParameters->gridforceOn )
01259     mapComputePatch(computeGridForceType);
01260   /* END gf */
01261   if ( node->simParameters->stirOn )
01262     mapComputePatch(computeStirType);
01263   if ( node->simParameters->sphericalBCOn )
01264     mapComputePatch(computeSphericalBCType);
01265   if ( node->simParameters->cylindricalBCOn )
01266     mapComputePatch(computeCylindricalBCType);
01267   if ( node->simParameters->tclBCOn ) {
01268     mapComputeHomePatches(computeTclBCType);
01269   }
01270   if ( node->simParameters->constraintsOn )
01271     mapComputePatch(computeRestraintsType);
01272   if ( node->simParameters->consForceOn )
01273     mapComputePatch(computeConsForceType);
01274   if ( node->simParameters->consTorqueOn )
01275     mapComputePatch(computeConsTorqueType);
01276 }
01277 
01278 //----------------------------------------------------------------------
01279 void WorkDistrib::mapComputeHomeTuples(ComputeType type)
01280 {
01281   PatchMap *patchMap = PatchMap::Object();
01282   ComputeMap *computeMap = ComputeMap::Object();
01283   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01284   Node *node = nd.ckLocalBranch();
01285 
01286   int numNodes = node->numNodes();
01287   int numPatches = patchMap->numPatches();
01288   ComputeID cid;
01289   PatchIDList basepids;
01290 
01291   for(int i=0; i<numNodes; i++) {
01292     patchMap->basePatchIDList(i,basepids);
01293     if ( basepids.size() ) {
01294       cid=computeMap->storeCompute(i,basepids.size(),type);
01295       for(int j=0; j<basepids.size(); ++j) {
01296         patchMap->newCid(basepids[j],cid);
01297         computeMap->newPid(cid,basepids[j]);
01298       }
01299     }
01300   }
01301 }
01302 
01303 //----------------------------------------------------------------------
01304 void WorkDistrib::mapComputeHomePatches(ComputeType type)
01305 {
01306   PatchMap *patchMap = PatchMap::Object();
01307   ComputeMap *computeMap = ComputeMap::Object();
01308   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01309   Node *node = nd.ckLocalBranch();
01310 
01311   int numNodes = node->numNodes();
01312   int numPatches = patchMap->numPatches();
01313   ComputeID *cid = new ComputeID[numNodes];
01314 
01315   for(int i=0; i<numNodes; i++) {
01316     if ( patchMap->numPatchesOnNode(i) ) {
01317       cid[i]=computeMap->storeCompute(i,patchMap->numPatchesOnNode(i),type);
01318     }
01319   }
01320 
01321   PatchID j;
01322 
01323   for(j=0;j<numPatches;j++)
01324   {
01325     patchMap->newCid(j,cid[patchMap->node(j)]);
01326     computeMap->newPid(cid[patchMap->node(j)],j);
01327   }
01328 
01329   delete [] cid;
01330 }
01331 
01332 //----------------------------------------------------------------------
01333 void WorkDistrib::mapComputePatch(ComputeType type)
01334 {
01335   PatchMap *patchMap = PatchMap::Object();
01336   ComputeMap *computeMap = ComputeMap::Object();
01337 
01338   PatchID i;
01339   ComputeID cid;
01340 
01341   for(i=0; i<patchMap->numPatches(); i++)
01342   {
01343     cid=computeMap->storeCompute(patchMap->node(i),1,type);
01344     computeMap->newPid(cid,i);
01345     patchMap->newCid(i,cid);
01346   }
01347 
01348 }
01349 
01350 
01351 //----------------------------------------------------------------------
01352 void WorkDistrib::mapComputeNonbonded(void)
01353 {
01354   // For each patch, create 1 electrostatic object for self-interaction.
01355   // Then create 1 for each 1-away and 2-away neighbor which has a larger
01356   // pid.
01357 
01358   PatchMap *patchMap = PatchMap::Object();
01359   ComputeMap *computeMap = ComputeMap::Object();
01360   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01361   Node *node = nd.ckLocalBranch();
01362 
01363   PatchID oneAway[PatchMap::MaxOneOrTwoAway];
01364   PatchID oneAwayTrans[PatchMap::MaxOneOrTwoAway];
01365 
01366   PatchID i;
01367   ComputeID cid;
01368   int numNeighbors;
01369   int j;
01370 
01371   for(i=0; i<patchMap->numPatches(); i++) // do the self 
01372   {
01373     int64 numAtoms = patchMap->patch(i)->getNumAtoms();  // avoid overflow
01374     int64 numFixed = patchMap->patch(i)->getNumFixedAtoms();  // avoid overflow
01375     int numPartitions = 0;
01376     int divide = node->simParameters->numAtomsSelf;
01377     if (divide > 0) {
01378       numPartitions = (int) ( 0.5 +
01379         (numAtoms*numAtoms-numFixed*numFixed) / (double)(2*divide*divide) );
01380     }
01381     if (numPartitions < 1) numPartitions = 1;
01382     if ( numPartitions > node->simParameters->maxSelfPart )
01383                         numPartitions = node->simParameters->maxSelfPart;
01384     // self-interaction
01385     DebugM(4,"Mapping " << numPartitions << " ComputeNonbondedSelf objects for patch " << i << "\n");
01386 //    iout <<"Self numPartitions = " <<numPartitions <<" numAtoms " <<numAtoms <<std::endl;
01387     for(int partition=0; partition < numPartitions; partition++)
01388     {
01389       cid=computeMap->storeCompute(patchMap->node(i),1,
01390                                    computeNonbondedSelfType,
01391                                    partition,numPartitions);
01392       computeMap->newPid(cid,i);
01393       patchMap->newCid(i,cid);
01394     }
01395   }
01396 
01397   for(int p1=0; p1 <patchMap->numPatches(); p1++) // do the pairs
01398   {
01399     // this only returns half of neighbors, which is what we want
01400     numNeighbors=patchMap->oneOrTwoAwayNeighbors(p1,oneAway,oneAwayTrans);
01401     for(j=0;j<numNeighbors;j++)
01402     {
01403         int p2 = oneAway[j];
01404         int64 numAtoms1 = patchMap->patch(p1)->getNumAtoms();
01405         int64 numAtoms2 = patchMap->patch(p2)->getNumAtoms();
01406         int64 numFixed1 = patchMap->patch(p1)->getNumFixedAtoms();
01407         int64 numFixed2 = patchMap->patch(p2)->getNumFixedAtoms();
01408         const int nax = patchMap->numaway_a();  // 1 or 2
01409         const int nay = patchMap->numaway_b();  // 1 or 2
01410         const int naz = patchMap->numaway_c();  // 1 or 2
01411         const int ia1 = patchMap->index_a(p1);
01412         const int ia2 = patchMap->index_a(p2);
01413         const int ib1 = patchMap->index_b(p1);
01414         const int ib2 = patchMap->index_b(p2);
01415         const int ic1 = patchMap->index_c(p1);
01416         const int ic2 = patchMap->index_c(p2);
01417         int distance = 3;
01418         if ( ia1 == ia2 ) --distance;
01419         else if ( ia1 == ia2 + nax - 1 ) --distance;
01420         else if ( ia1 + nax - 1 == ia2 ) --distance;
01421         if ( ib1 == ib2 ) --distance;
01422         else if ( ib1 == ib2 + nay - 1 ) --distance;
01423         else if ( ib1 + nay - 1 == ib2 ) --distance;
01424         if ( ic1 == ic2 ) --distance;
01425         else if ( ic1 == ic2 + naz - 1 ) --distance;
01426         else if ( ic1 + naz - 1 == ic2 ) --distance;
01427         int divide = 0;
01428         if ( distance == 0 ) {
01429           divide = node->simParameters->numAtomsSelf2;
01430         } else if (distance == 1) {
01431           divide = node->simParameters->numAtomsPair;
01432         } else {
01433           divide = node->simParameters->numAtomsPair2;
01434         }
01435         int numPartitions = 0;
01436         if (divide > 0) {
01437           numPartitions = (int) ( 0.5 +
01438             (numAtoms1*numAtoms2-numFixed1*numFixed2)/(double)(divide*divide) );
01439         }
01440         if ( numPartitions < 1 ) numPartitions = 1;
01441         if ( numPartitions > node->simParameters->maxPairPart )
01442                         numPartitions = node->simParameters->maxPairPart;
01443 //      if ( numPartitions > 1 ) iout << "Mapping " << numPartitions << " ComputeNonbondedPair objects for patches " << p1 << "(" << numAtoms1 << ") and " << p2 << "(" << numAtoms2 << ")\n" << endi;
01444         for(int partition=0; partition < numPartitions; partition++)
01445         {
01446           cid=computeMap->storeCompute(
01447                 patchMap->basenode(patchMap->downstream2(p1,p2)),
01448                 2,computeNonbondedPairType,partition,numPartitions);
01449           computeMap->newPid(cid,p1);
01450           computeMap->newPid(cid,p2,oneAwayTrans[j]);
01451           patchMap->newCid(p1,cid);
01452           patchMap->newCid(p2,cid);
01453         }
01454     }
01455   }
01456 
01457 }
01458 
01459 //----------------------------------------------------------------------
01460 void WorkDistrib::messageEnqueueWork(Compute *compute) {
01461   LocalWorkMsg *msg = compute->localWorkMsg;
01462   int seq = compute->sequence();
01463 
01464   if ( seq < 0 ) {
01465     NAMD_bug("compute->sequence() < 0 in WorkDistrib::messageEnqueueWork");
01466   } else {
01467     SET_PRIORITY(msg,seq,compute->priority());
01468   }
01469 
01470   msg->compute = compute; // pointer is valid since send is to local Pe
01471   int type = compute->type();
01472 
01473   CProxy_WorkDistrib wdProxy(CkpvAccess(BOCclass_group).workDistrib);
01474   switch ( type ) {
01475   case computeBondsType:
01476   case computeSelfBondsType:
01477 #if CHARM_VERSION > 050402
01478     wdProxy[CkMyPe()].enqueueBonds(msg);
01479 #else
01480     wdProxy.enqueueBonds(msg,CkMyPe());
01481 #endif
01482     break;
01483   case computeAnglesType:
01484   case computeSelfAnglesType:
01485 #if CHARM_VERSION > 050402
01486     wdProxy[CkMyPe()].enqueueAngles(msg);
01487 #else
01488     wdProxy.enqueueAngles(msg,CkMyPe());
01489 #endif
01490     break;
01491   case computeDihedralsType:
01492   case computeSelfDihedralsType:
01493 #if CHARM_VERSION > 050402
01494     wdProxy[CkMyPe()].enqueueDihedrals(msg);
01495 #else
01496     wdProxy.enqueueDihedrals(msg,CkMyPe());
01497 #endif
01498     break;
01499   case computeImpropersType:
01500   case computeSelfImpropersType:
01501 #if CHARM_VERSION > 050402
01502     wdProxy[CkMyPe()].enqueueImpropers(msg);
01503 #else
01504     wdProxy.enqueueImpropers(msg,CkMyPe());
01505 #endif
01506     break;
01507   case computeCrosstermsType:
01508   case computeSelfCrosstermsType:
01509 #if CHARM_VERSION > 050402
01510     wdProxy[CkMyPe()].enqueueCrossterms(msg);
01511 #else
01512     wdProxy.enqueueCrossterms(msg,CkMyPe());
01513 #endif
01514     break;
01515   case computeNonbondedSelfType:
01516     switch ( seq % 2 ) {
01517     case 0:
01518 #if CHARM_VERSION > 050402
01519       wdProxy[CkMyPe()].enqueueSelfA(msg);
01520 #else
01521       wdProxy.enqueueSelfA(msg,CkMyPe());
01522 #endif
01523       break;
01524     case 1:
01525 #if CHARM_VERSION > 050402
01526       wdProxy[CkMyPe()].enqueueSelfB(msg);
01527 #else
01528       wdProxy.enqueueSelfB(msg,CkMyPe());
01529 #endif
01530       break;
01531     default:
01532       NAMD_bug("WorkDistrib::messageEnqueueSelf case statement error!");
01533     }
01534     break;
01535   case computeNonbondedPairType:
01536     switch ( seq % 2 ) {
01537     case 0:
01538 #if CHARM_VERSION > 050402
01539       wdProxy[CkMyPe()].enqueueWorkA(msg);
01540 #else 
01541       wdProxy.enqueueWorkA(msg,CkMyPe());
01542 #endif
01543       break;
01544     case 1:
01545 #if CHARM_VERSION > 050402
01546       wdProxy[CkMyPe()].enqueueWorkB(msg);
01547 #else
01548       wdProxy.enqueueWorkB(msg,CkMyPe());
01549 #endif
01550       break;
01551     case 2:
01552 #if CHARM_VERSION > 050402
01553       wdProxy[CkMyPe()].enqueueWorkC(msg);
01554 #else
01555       wdProxy.enqueueWorkC(msg,CkMyPe());
01556 #endif
01557       break;
01558     default:
01559       NAMD_bug("WorkDistrib::messageEnqueueWork case statement error!");
01560     }
01561     break;
01562   case computePmeType:
01563 #if 0
01564 #if CHARM_VERSION > 050402
01565     wdProxy[CkMyPe()].enqueuePme(msg);
01566 #else
01567     wdProxy.enqueuePme(msg,CkMyPe());
01568 #endif
01569 #else
01570     msg->compute->doWork();
01571 #endif
01572     break;
01573   default:
01574 #if CHARM_VERSION > 050402
01575     wdProxy[CkMyPe()].enqueueWork(msg);
01576 #else
01577     wdProxy.enqueueWork(msg,CkMyPe());
01578 #endif
01579   }
01580 }
01581 
01582 void WorkDistrib::enqueueWork(LocalWorkMsg *msg) {
01583   msg->compute->doWork();
01584   if ( msg->compute->localWorkMsg != msg )
01585     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01586 }
01587 
01588 void WorkDistrib::enqueueBonds(LocalWorkMsg *msg) {
01589   msg->compute->doWork();
01590   if ( msg->compute->localWorkMsg != msg )
01591     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01592 }
01593 
01594 void WorkDistrib::enqueueA