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: 2009/10/09 20:09:05 $
00011  * $Revision: 1.1197 $
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 "Lattice.h"
00033 #include "main.decl.h"
00034 #include "main.h"
00035 #include "Node.h"
00036 #include "PatchMgr.h"
00037 #include "PatchMap.inl"
00038 #include "NamdTypes.h"
00039 #include "PDB.h"
00040 #include "SimParameters.h"
00041 #include "Molecule.h"
00042 #include "NamdOneTools.h"
00043 #include "Compute.h"
00044 #include "ComputeMap.h"
00045 #include "RecBisection.h"
00046 #include "Random.h"
00047 #include "varsizemsg.h"
00048 #include "ProxyMgr.h"
00049 #include "Priorities.h"
00050 
00051 //#define DEBUGM
00052 #define MIN_DEBUG_LEVEL 2
00053 #include "Debug.h"
00054 
00055 
00056 class ComputeMapChangeMsg : public CMessage_ComputeMapChangeMsg
00057 {
00058 public:
00059 
00060   int numNewNodes;
00061   int *newNodes;
00062 
00063 //  VARSIZE_DECL(ComputeMapChangeMsg);
00064 };
00065 
00066 /*
00067 VARSIZE_MSG(ComputeMapChangeMsg,
00068   VARSIZE_ARRAY(newNodes);
00069 )
00070 */
00071 
00072 
00073 //======================================================================
00074 // Public functions
00075 //----------------------------------------------------------------------
00076 WorkDistrib::WorkDistrib()
00077 {
00078   CkpvAccess(BOCclass_group).workDistrib = thisgroup;
00079   mapsArrived = false;
00080   awaitingMaps = false;
00081 }
00082 
00083 //----------------------------------------------------------------------
00084 WorkDistrib::~WorkDistrib(void)
00085 { }
00086 
00087 
00088 //----------------------------------------------------------------------
00089 void WorkDistrib::saveComputeMapChanges(int ep, CkGroupID chareID)
00090 {
00091   saveComputeMapReturnEP = ep;
00092   saveComputeMapReturnChareID = chareID;
00093   saveComputeMapCount = CkNumPes();
00094 
00095   ComputeMap *computeMap = ComputeMap::Object();
00096 
00097   int i;
00098   int nc = computeMap->numComputes();
00099   for (i=0; i<nc; i++) {
00100     DebugM(3, "ComputeMap (" << i << ") node = " << computeMap->node(i) << " newNode = " << computeMap->newNode(i) << "\n");
00101   }
00102   
00103   ComputeMapChangeMsg *mapMsg = new (nc, 0) ComputeMapChangeMsg ;
00104 
00105   mapMsg->numNewNodes = nc;
00106   for(i=0; i<nc; i++)
00107     mapMsg->newNodes[i] = computeMap->newNode(i);
00108 
00109   CProxy_WorkDistrib(thisgroup).recvComputeMapChanges(mapMsg);
00110 }
00111 
00112 void WorkDistrib::recvComputeMapChanges(ComputeMapChangeMsg *msg) {
00113   
00114   ComputeMap *computeMap = ComputeMap::Object();
00115   int i;
00116   for(i=0; i<computeMap->numComputes(); i++)
00117     computeMap->setNewNode(i,msg->newNodes[i]);
00118 
00119   delete msg;
00120 
00121   CProxy_WorkDistrib workProxy(thisgroup);
00122   workProxy[0].doneSaveComputeMap();
00123 
00124   DebugM(2, "ComputeMap after send!\n");
00125   for (i=0; i<computeMap->numComputes(); i++) {
00126     DebugM(2, "ComputeMap (" << i << ") node = " << computeMap->node(i) << " newNode = " << computeMap->newNode(i) << " type=" << computeMap->type(i) << "\n");
00127   }
00128   DebugM(2, "===================================================\n");
00129 }
00130 
00131 void WorkDistrib::doneSaveComputeMap() {
00132   if (!--saveComputeMapCount) { 
00133     CkSendMsgBranch(saveComputeMapReturnEP, CkAllocMsg(0,0,0), 0, saveComputeMapReturnChareID);
00134   }
00135 }
00136 
00137 
00138 //only called on node 0
00139 int *WorkDistrib::caclNumAtomsInEachPatch(){
00140     StringList *current;
00141     int i;
00142     CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00143     Node *node = nd.ckLocalBranch();
00144     PatchMap *patchMap = PatchMap::Object();
00145     CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00146     PatchMgr *patchMgr = pm.ckLocalBranch();
00147     SimParameters *params = node->simParameters;
00148     Molecule *molecule = node->molecule;
00149     PDB *pdb = node->pdb;
00150 
00151     int numPatches = patchMap->numPatches();
00152     int numAtoms = pdb->num_atoms();    
00153 
00154     vector<int> *eachPatchAtomList = patchMap->getTmpPatchAtomsList();
00155 
00156     int *patchAtomCnt = new int[numPatches];    
00157     memset(patchAtomCnt, 0, sizeof(int)*numPatches);
00158 
00159     const Lattice lattice = params->lattice;
00160 
00161     Position eachAtomPos;
00162     if (params->splitPatch == SPLIT_PATCH_HYDROGEN)
00163       {
00164       // split atoms into patched based on helix-group and position
00165       int aid, pid=0;
00166       for(i=0; i < numAtoms; i++)
00167         {        
00168         // Assign atoms to patches without splitting hydrogen groups.
00169         // We know that the hydrogenGroup array is sorted with group parents
00170         // listed first.  Thus, only change the pid if an atom is a group parent.
00171         aid = molecule->hydrogenGroup[i].atomID;
00172         pdb->get_position_for_atom(&eachAtomPos, aid);
00173         if (molecule->hydrogenGroup[i].isGP)            
00174             pid = patchMap->assignToPatch(eachAtomPos,lattice);
00175         // else: don't change pid        
00176         patchAtomCnt[pid]++;
00177         eachPatchAtomList[pid].push_back(aid);
00178         }
00179       }
00180     else
00181       {
00182       // split atoms into patched based on position
00183       for(i=0; i < numAtoms; i++)
00184         {
00185         pdb->get_position_for_atom(&eachAtomPos, i);
00186         int pid = patchMap->assignToPatch(eachAtomPos,lattice);        
00187         patchAtomCnt[pid]++;
00188         eachPatchAtomList[pid].push_back(i);
00189         }
00190       }   
00191 
00192     return patchAtomCnt;    
00193 }
00194 
00195 //----------------------------------------------------------------------
00196 // This should only be called on node 0.
00197 //----------------------------------------------------------------------
00198 FullAtomList *WorkDistrib::createAtomLists(void)
00199 {
00200   int i;
00201   StringList *current;  //  Pointer used to retrieve configuration items
00202   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00203   Node *node = nd.ckLocalBranch();
00204   PatchMap *patchMap = PatchMap::Object();
00205   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00206   PatchMgr *patchMgr = pm.ckLocalBranch();
00207   SimParameters *params = node->simParameters;
00208   Molecule *molecule = node->molecule;
00209   PDB *pdb = node->pdb;
00210 
00211   int numPatches = patchMap->numPatches();
00212   int numAtoms = pdb->num_atoms();
00213 
00214   Vector *positions = new Position[numAtoms];
00215   pdb->get_all_positions(positions);
00216 
00217   Vector *velocities = new Velocity[numAtoms];
00218 
00219   if ( params->initialTemp < 0.0 ) {
00220     Bool binvels=FALSE;
00221 
00222     //  Reading the veolcities from a PDB
00223     current = node->configList->find("velocities");
00224 
00225     if (current == NULL) {
00226       current = node->configList->find("binvelocities");
00227       binvels = TRUE;
00228     }
00229 
00230     if (!binvels) {
00231       velocities_from_PDB(current->data, velocities, numAtoms);
00232     }
00233     else {
00234       velocities_from_binfile(current->data, velocities, numAtoms);
00235     }
00236   }
00237   else {
00238     // Random velocities for a given temperature
00239     random_velocities(params->initialTemp, molecule, velocities, numAtoms);
00240   }
00241 
00242   //  If COMMotion == no, remove center of mass motion
00243   if (!(params->comMove)) {
00244     remove_com_motion(velocities, molecule, numAtoms);
00245   }
00246 
00247   FullAtomList *atoms = new FullAtomList[numPatches];
00248 
00249   const Lattice lattice = params->lattice;
00250 
00251   if (params->splitPatch == SPLIT_PATCH_HYDROGEN)
00252     {
00253     // split atoms into patched based on helix-group and position
00254     int aid, pid=0;
00255     for(i=0; i < numAtoms; i++)
00256       {
00257       if ( ! ( i % 1000 ) )
00258         {
00259         DebugM(3,"Assigned " << i << " atoms to patches so far.\n");
00260         }
00261       // Assign atoms to patches without splitting hydrogen groups.
00262       // We know that the hydrogenGroup array is sorted with group parents
00263       // listed first.  Thus, only change the pid if an atom is a group parent.
00264       aid = molecule->hydrogenGroup[i].atomID;
00265       if (molecule->hydrogenGroup[i].isGP)
00266         pid = patchMap->assignToPatch(positions[aid],lattice);
00267       // else: don't change pid
00268       FullAtom a;
00269       a.id = aid;
00270       a.position = positions[aid];
00271       a.velocity = velocities[aid];
00272       #ifdef MEM_OPT_VERSION
00273       a.sigId = molecule->getAtomSigId(aid);
00274       a.exclId = molecule->getAtomExclSigId(aid);
00275       #endif
00276       a.vdwType = molecule->atomvdwtype(aid);
00277       atoms[pid].add(a);
00278       }
00279     }
00280   else
00281     {
00282     // split atoms into patched based on position
00283     for(i=0; i < numAtoms; i++)
00284       {
00285       if ( ! ( i % 1000 ) )
00286         {
00287         DebugM(3,"Assigned " << i << " atoms to patches so far.\n");
00288         }
00289       int pid = patchMap->assignToPatch(positions[i],lattice);
00290       FullAtom a;
00291       a.id = i;
00292       a.position = positions[i];
00293       a.velocity = velocities[i];
00294       #ifdef MEM_OPT_VERSION
00295       a.sigId = molecule->getAtomSigId(i);
00296       a.exclId = molecule->getAtomExclSigId(i);
00297       #endif
00298       a.vdwType = molecule->atomvdwtype(i);
00299       atoms[pid].add(a);
00300       }
00301     }
00302 
00303   delete [] positions;
00304   delete [] velocities;
00305 
00306   for(i=0; i < numPatches; i++)
00307   {
00308     ScaledPosition center(0.5*(patchMap->min_a(i)+patchMap->max_a(i)),
00309                           0.5*(patchMap->min_b(i)+patchMap->max_b(i)),
00310                           0.5*(patchMap->min_c(i)+patchMap->max_c(i)));
00311 
00312     int n = atoms[i].size();
00313     FullAtom *a = atoms[i].begin();
00314     int j;
00315 //Modifications for alchemical fep
00316     Bool alchFepOn = params->alchFepOn;
00317     Bool alchThermIntOn = params->alchThermIntOn;
00318 //fepe
00319     Bool lesOn = params->lesOn;
00320   
00321     Bool pairInteractionOn = params->pairInteractionOn;
00322 
00323     Bool pressureProfileTypes = (params->pressureProfileAtomTypes > 1);
00324 
00325     Transform mother_transform;
00326     for(j=0; j < n; j++)
00327     {
00328       int aid = a[j].id;
00329 
00330       if (params->splitPatch == SPLIT_PATCH_HYDROGEN) {
00331         if ( molecule->is_hydrogenGroupParent(aid) ) {
00332           a[j].hydrogenGroupSize = molecule->get_groupSize(aid);
00333         } else {
00334           a[j].hydrogenGroupSize = 0;
00335         }
00336       } else {
00337         a[j].hydrogenGroupSize = 1;
00338       }
00339 
00340       a[j].nonbondedGroupSize = 0;  // must be set based on coordinates
00341 
00342       a[j].atomFixed = molecule->is_atom_fixed(aid) ? 1 : 0;
00343       a[j].fixedPosition = a[j].position;
00344 
00345       if ( a[j].hydrogenGroupSize ) {
00346         a[j].position = lattice.nearest(
00347                 a[j].position, center, &(a[j].transform));
00348         mother_transform = a[j].transform;
00349       } else {
00350         a[j].position = lattice.apply_transform(a[j].position,mother_transform);
00351         a[j].transform = mother_transform;
00352       }
00353 
00354       a[j].mass = molecule->atommass(aid);
00355       a[j].charge = molecule->atomcharge(aid);
00356 
00357 //Modifications for alchemical fep
00358       if ( alchFepOn || alchThermIntOn || lesOn || pairInteractionOn || pressureProfileTypes) {
00359         a[j].partition = molecule->get_fep_type(aid);
00360       } 
00361       else {
00362         a[j].partition = 0;
00363       }
00364 //fepe
00365 
00366     }
00367 
00368     int size, allfixed, k;
00369     for(j=0; j < n; j+=size) {
00370       size = a[j].hydrogenGroupSize;
00371       if ( ! size ) {
00372         NAMD_bug("Mother atom with hydrogenGroupSize of 0!");
00373       }
00374       allfixed = 1;
00375       for ( k = 0; k < size; ++k ) {
00376         allfixed = ( allfixed && (a[j+k].atomFixed) );
00377       }
00378       for ( k = 0; k < size; ++k ) {
00379         a[j+k].groupFixed = allfixed ? 1 : 0;
00380       }
00381     }
00382 
00383     if ( params->outputPatchDetails ) {
00384       int patchId = i;
00385       int numAtomsInPatch = n;
00386       int numFixedAtomsInPatch = 0;
00387       int numAtomsInFixedGroupsInPatch = 0;
00388       for(j=0; j < n; j++) {
00389         numFixedAtomsInPatch += ( a[j].atomFixed ? 1 : 0 );
00390         numAtomsInFixedGroupsInPatch += ( a[j].groupFixed ? 1 : 0 );
00391       }
00392       iout << "PATCH_DETAILS:"
00393            << " patch " << patchId
00394            << " atoms " << numAtomsInPatch
00395            << " fixed_atoms " << numFixedAtomsInPatch
00396            << " fixed_groups " << numAtomsInFixedGroupsInPatch
00397            << "\n" << endi;
00398     }
00399   }
00400 
00401   return atoms;
00402 
00403 }
00404 
00405 //This should only be called on node 0
00406 //This differs from creatHomePatches in: create home patches
00407 //without populating them with actual atoms' data. The only field
00408 //set is the number of atoms each patch contains.
00409 void WorkDistrib::preCreateHomePatches(){
00410 
00411   PatchMap *patchMap = PatchMap::Object();
00412   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00413   PatchMgr *patchMgr = pm.ckLocalBranch();
00414 
00415   int numPatches = patchMap->numPatches();
00416 
00417   patchMap->initTmpPatchAtomsList();
00418 
00419   int *patchAtomCnt = caclNumAtomsInEachPatch();
00420 
00421   int maxAtoms = -1;
00422   int maxPatch = -1;
00423   for(int i=0; i < numPatches; i++) {
00424     int numAtoms = patchAtomCnt[i];
00425     if ( numAtoms > maxAtoms ) { maxAtoms = numAtoms; maxPatch = i; }
00426   }
00427   iout << iINFO << "LARGEST PATCH (" << maxPatch <<
00428         ") HAS " << maxAtoms << " ATOMS\n" << endi;
00429 
00430   for(int i=0; i < numPatches; i++)
00431   {
00432     if ( ! ( i % 100 ) )
00433     {
00434       DebugM(3,"Pre-created " << i << " patches so far.\n");
00435     }
00436 
00437     patchMgr->preCreateHomePatch(i,patchAtomCnt[i]);
00438   }
00439 
00440   delete [] patchAtomCnt;
00441 }
00442 
00443 
00444 //----------------------------------------------------------------------
00445 // This should only be called on node 0.
00446 //----------------------------------------------------------------------
00447 void WorkDistrib::createHomePatches(void)
00448 {
00449   int i;
00450   PatchMap *patchMap = PatchMap::Object();
00451   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00452   PatchMgr *patchMgr = pm.ckLocalBranch();
00453 
00454   int numPatches = patchMap->numPatches();
00455 
00456   FullAtomList *atoms = createAtomLists();
00457     
00458 #ifdef MEM_OPT_VERSION
00459 /*  CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00460   Node *node = nd.ckLocalBranch();
00461   node->molecule->delEachAtomSigs();
00462   node->molecule->delMassChargeSpace();
00463 */
00464 #endif
00465 
00466   int maxAtoms = -1;
00467   int maxPatch = -1;
00468   for(i=0; i < numPatches; i++) {
00469     int numAtoms = atoms[i].size();
00470     if ( numAtoms > maxAtoms ) { maxAtoms = numAtoms; maxPatch = i; }
00471   }
00472   iout << iINFO << "LARGEST PATCH (" << maxPatch <<
00473         ") HAS " << maxAtoms << " ATOMS\n" << endi;
00474 
00475   for(i=0; i < numPatches; i++)
00476   {
00477     if ( ! ( i % 100 ) )
00478     {
00479       DebugM(3,"Created " << i << " patches so far.\n");
00480     }
00481 
00482     patchMgr->createHomePatch(i,atoms[i]);
00483   }
00484 
00485   delete [] atoms;
00486 }
00487 
00488 void WorkDistrib::fillOnePatchAtoms(int patchId, FullAtomList *onePatchAtoms, Vector *velocities){
00489     // ref BOC
00490     CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00491     Node *node = nd.ckLocalBranch();
00492     PDB *pdb = node->pdb;
00493     SimParameters *params = node->simParameters;
00494     Molecule *molecule = node->molecule;
00495 
00496     const Lattice lattice = params->lattice;
00497 
00498     CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00499     PatchMgr *patchMgr = pm.ckLocalBranch();
00500     // ref singleton
00501     PatchMap *patchMap = PatchMap::Object();
00502 
00503     vector<int> *eachPatchAtomsList = patchMap->getTmpPatchAtomsList();
00504     vector<int> *thisPatchAtomsList = &eachPatchAtomsList[patchId];
00505 
00506     for(int i=0; i<thisPatchAtomsList->size(); i++){
00507         int aid = thisPatchAtomsList->at(i);
00508         FullAtom a;
00509         a.id = aid;
00510         Position pos;
00511         pdb->get_position_for_atom(&pos, aid);
00512         a.position = pos;
00513         a.velocity = velocities[aid];
00514         #ifdef MEM_OPT_VERSION
00515         a.sigId = molecule->getAtomSigId(aid);
00516         a.exclId = molecule->getAtomExclSigId(aid);
00517         #endif
00518         a.vdwType = molecule->atomvdwtype(aid);
00519         onePatchAtoms->add(a);
00520     }
00521 
00522     ScaledPosition center(0.5*(patchMap->min_a(patchId)+patchMap->max_a(patchId)),
00523                           0.5*(patchMap->min_b(patchId)+patchMap->max_b(patchId)),
00524                           0.5*(patchMap->min_c(patchId)+patchMap->max_c(patchId)));
00525     int n = onePatchAtoms->size();
00526     FullAtom *a = onePatchAtoms->begin();
00527     int j;
00528 //Modifications for alchemical fep
00529     Bool alchFepOn = params->alchFepOn;
00530     Bool alchThermIntOn = params->alchThermIntOn;
00531 //fepe
00532     Bool lesOn = params->lesOn;
00533   
00534     Bool pairInteractionOn = params->pairInteractionOn;
00535 
00536     Bool pressureProfileTypes = (params->pressureProfileAtomTypes > 1);
00537 
00538     Transform mother_transform;
00539     for(j=0; j < n; j++)
00540     {
00541       int aid = a[j].id;
00542 
00543       if (params->splitPatch == SPLIT_PATCH_HYDROGEN) {
00544         if ( molecule->is_hydrogenGroupParent(aid) ) {
00545           a[j].hydrogenGroupSize = molecule->get_groupSize(aid);
00546         } else {
00547           a[j].hydrogenGroupSize = 0;
00548         }
00549       } else {
00550         a[j].hydrogenGroupSize = 1;
00551       }
00552 
00553       a[j].nonbondedGroupSize = 0;  // must be set based on coordinates
00554 
00555       a[j].atomFixed = molecule->is_atom_fixed(aid) ? 1 : 0;
00556       a[j].fixedPosition = a[j].position;
00557 
00558       if ( a[j].hydrogenGroupSize ) {
00559         a[j].position = lattice.nearest(
00560                 a[j].position, center, &(a[j].transform));
00561         mother_transform = a[j].transform;
00562       } else {
00563         a[j].position = lattice.apply_transform(a[j].position,mother_transform);
00564         a[j].transform = mother_transform;
00565       }
00566 
00567       a[j].mass = molecule->atommass(aid);
00568       a[j].charge = molecule->atomcharge(aid);
00569 
00570 //Modifications for alchemical fep
00571       if ( alchFepOn || alchThermIntOn || lesOn || pairInteractionOn || pressureProfileTypes) {
00572         a[j].partition = molecule->get_fep_type(aid);
00573       } 
00574       else {
00575         a[j].partition = 0;
00576       }
00577 //fepe
00578     }
00579 
00580     int size, allfixed, k;
00581     for(j=0; j < n; j+=size) {
00582       size = a[j].hydrogenGroupSize;
00583       if ( ! size ) {
00584         NAMD_bug("Mother atom with hydrogenGroupSize of 0!");
00585       }
00586       allfixed = 1;
00587       for ( k = 0; k < size; ++k ) {
00588         allfixed = ( allfixed && (a[j+k].atomFixed) );
00589       }
00590       for ( k = 0; k < size; ++k ) {
00591         a[j+k].groupFixed = allfixed ? 1 : 0;
00592       }
00593     }
00594 
00595     if(params->fixedAtomsOn){
00596         int fixedCnt=0;
00597         for(j=0; j<n; j++)
00598             fixedCnt += (a[j].atomFixed ? 1:0);
00599         patchMgr->setHomePatchFixedAtomNum(patchId, fixedCnt);
00600     }    
00601 
00602     if ( params->outputPatchDetails ) {    
00603       int numAtomsInPatch = n;
00604       int numFixedAtomsInPatch = 0;
00605       int numAtomsInFixedGroupsInPatch = 0;
00606       for(j=0; j < n; j++) {
00607         numFixedAtomsInPatch += ( a[j].atomFixed ? 1 : 0 );
00608         numAtomsInFixedGroupsInPatch += ( a[j].groupFixed ? 1 : 0 );
00609       }
00610       iout << "PATCH_DETAILS:"
00611            << " patch " << patchId
00612            << " atoms " << numAtomsInPatch
00613            << " fixed_atoms " << numFixedAtomsInPatch
00614            << " fixed_groups " << numAtomsInFixedGroupsInPatch
00615            << "\n" << endi;
00616     }
00617 }
00618 
00619 //should be called only on node 0
00620 void WorkDistrib::initAndSendHomePatch(){
00621   StringList *current;
00622   // ref BOC
00623   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00624   Node *node = nd.ckLocalBranch();
00625   Molecule *molecule = node->molecule;
00626   PDB *pdb = node->pdb;
00627   SimParameters *params = node->simParameters;
00628 
00629   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00630   PatchMgr *patchMgr = pm.ckLocalBranch();
00631   // ref singleton
00632   PatchMap *patchMap = PatchMap::Object();
00633 
00634   int numAtoms = pdb->num_atoms();
00635 
00636   //1. create atoms' velocities
00637   Vector *velocities = new Velocity[numAtoms];
00638   if ( params->initialTemp < 0.0 ) {
00639     Bool binvels=FALSE;
00640 
00641     //  Reading the veolcities from a PDB
00642     current = node->configList->find("velocities");
00643 
00644     if (current == NULL) {
00645       current = node->configList->find("binvelocities");
00646       binvels = TRUE;
00647     }
00648 
00649     if (!binvels) {
00650       velocities_from_PDB(current->data, velocities, numAtoms);
00651     }
00652     else {
00653       velocities_from_binfile(current->data, velocities, numAtoms);
00654     }
00655   }
00656   else {
00657     // Random velocities for a given temperature
00658     random_velocities(params->initialTemp, molecule, velocities, numAtoms);
00659   }
00660 
00661   //  If COMMotion == no, remove center of mass motion
00662   if (!(params->comMove)) {
00663     remove_com_motion(velocities, molecule, numAtoms);
00664   }
00665 
00666 
00667   for(int i=0; i<patchMap->numPatches(); i++){
00668       //2. fill each patch with actual atom data
00669       //the work flow should looks like the createAtomList
00670       FullAtomList *onePatchAtoms = new FullAtomList;
00671       fillOnePatchAtoms(i, onePatchAtoms, velocities);
00672       
00673       patchMgr->fillHomePatchAtomList(i, onePatchAtoms);
00674 
00675       delete onePatchAtoms;
00676 
00677       //distribute this home patch
00678       if(patchMap->node(i) != node->myid()){
00679           //need to move to other nodes
00680           patchMgr->sendOneHomePatch(i, patchMap->node(i));
00681       }
00682       
00683   }  
00684 
00685   delete [] velocities;
00686   patchMap->delTmpPatchAtomsList();
00687 }
00688 
00689 void WorkDistrib::distributeHomePatches() {
00690   // ref BOC
00691   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00692   Node *node = nd.ckLocalBranch();
00693   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00694   PatchMgr *patchMgr = pm.ckLocalBranch();
00695   // ref singleton
00696   PatchMap *patchMap = PatchMap::Object();
00697 
00698   // Move patches to the proper node
00699   for(int i=0;i < patchMap->numPatches(); i++)
00700   {
00701     if (patchMap->node(i) != node->myid() )
00702     {
00703       DebugM(3,"patchMgr->movePatch("
00704         << i << "," << patchMap->node(i) << ")\n");
00705       patchMgr->movePatch(i,patchMap->node(i));
00706     }
00707   }
00708   patchMgr->sendMovePatches();
00709 }
00710 
00711 void WorkDistrib::reinitAtoms() {
00712 
00713   PatchMap *patchMap = PatchMap::Object();
00714   CProxy_PatchMgr pm(CkpvAccess(BOCclass_group).patchMgr);
00715   PatchMgr *patchMgr = pm.ckLocalBranch();
00716 
00717   int numPatches = patchMap->numPatches();
00718 
00719   FullAtomList *atoms = createAtomLists();
00720 
00721   for(int i=0; i < numPatches; i++) {
00722     patchMgr->sendAtoms(i,atoms[i]);
00723   }
00724 
00725   delete [] atoms;
00726 
00727 }
00728 
00729 
00730 //----------------------------------------------------------------------
00731 
00732 class MapDistribMsg: public CMessage_MapDistribMsg {
00733   public:
00734     char *patchMapData;
00735     char *computeMapData;
00736 
00737 //  VARSIZE_DECL(MapDistribMsg);
00738 };
00739 
00740 /*
00741 VARSIZE_MSG(MapDistribMsg,
00742   VARSIZE_ARRAY(patchMapData);
00743   VARSIZE_ARRAY(computeMapData);
00744 )
00745 */
00746 
00747 void WorkDistrib::sendMaps(void)
00748 {
00749   if ( CkNumPes() == 1 ) {
00750     mapsArrived = true;
00751     return;
00752   }
00753 
00754   //Automatically enable spanning tree
00755   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00756   Node *node = nd.ckLocalBranch();
00757   SimParameters *params = node->simParameters;
00758   if(PatchMap::Object()->numPatches() <= CkNumPes()/4 && params->isSendSpanningTreeUnset())
00759     ProxyMgr::Object()->setSendSpanning();
00760 
00761   int sizes[2];
00762   sizes[0] = PatchMap::Object()->packSize();
00763   sizes[1] = ComputeMap::Object()->packSize();
00764 
00765   MapDistribMsg *mapMsg = new (sizes[0], sizes[1], 0) MapDistribMsg;
00766 
00767   PatchMap::Object()->pack(mapMsg->patchMapData);
00768   ComputeMap::Object()->pack(mapMsg->computeMapData);
00769 
00770   CProxy_WorkDistrib workProxy(thisgroup);
00771   workProxy[0].saveMaps(mapMsg);
00772 }
00773 
00774 // saveMaps() is called when the map message is received
00775 void WorkDistrib::saveMaps(MapDistribMsg *msg)
00776 {
00777   // Use a resend to forward messages before processing.  Otherwise the
00778   // map distribution is slow on many CPUs.  We need to use a tree
00779   // rather than a broadcast because some implementations of broadcast
00780   // generate a copy of the message on the sender for each recipient.
00781   // This is because MPI doesn't allow re-use of an outstanding buffer.
00782 
00783   if ( mapsArrived && CkMyPe() ) {
00784     PatchMap::Object()->unpack(msg->patchMapData);
00785     ComputeMap::Object()->unpack(msg->computeMapData);
00786 
00787     //Automatically enable spanning tree
00788     CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00789     Node *node = nd.ckLocalBranch();
00790     SimParameters *params = node->simParameters;
00791     if(PatchMap::Object()->numPatches() <= CkNumPes()/4 && params->isSendSpanningTreeUnset())
00792       ProxyMgr::Object()->setSendSpanning();
00793   }
00794   if ( mapsArrived ) {
00795     delete msg;
00796     return;
00797   }
00798 
00799   mapsArrived = true;
00800 
00801   int pids[3];
00802   int basePe = 2 * CkMyPe() + 1;
00803   int npid = 0;
00804   if ( (basePe+npid) < CkNumPes() ) { pids[npid] = basePe + npid; ++npid; }
00805   if ( (basePe+npid) < CkNumPes() ) { pids[npid] = basePe + npid; ++npid; }
00806   pids[npid] = CkMyPe(); ++npid;  // always send the message to ourselves
00807   CProxy_WorkDistrib(thisgroup).saveMaps(msg,npid,pids);
00808 }
00809 
00810 
00811 //----------------------------------------------------------------------
00812 void WorkDistrib::patchMapInit(void)
00813 {
00814   PatchMap *patchMap = PatchMap::Object();
00815   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
00816   Node *node = nd.ckLocalBranch();
00817   SimParameters *params = node->simParameters;
00818   Lattice lattice = params->lattice;
00819 
00820   BigReal patchSize = params->patchDimension;
00821 
00822   ScaledPosition xmin, xmax;
00823 
00824   double maxNumPatches = 1000000;  // need to adjust fractional values
00825   if ( params->minAtomsPerPatch > 0 )
00826     maxNumPatches = node->pdb->num_atoms() / params->minAtomsPerPatch;
00827 
00828   DebugM(3,"Mapping patches\n");
00829   if ( lattice.a_p() && lattice.b_p() && lattice.c_p() ) {
00830     xmin = 0.;  xmax = 0.;
00831   } else {
00832   // Need to use full box for FMA to match NAMD 1.X results.
00833   if ( params->FMAOn ) {
00834     node->pdb->find_extremes(&(xmin.x),&(xmax.x),lattice.a_r());
00835     node->pdb->find_extremes(&(xmin.y),&(xmax.y),lattice.b_r());
00836     node->pdb->find_extremes(&(xmin.z),&(xmax.z),lattice.c_r());
00837   // Otherwise, this allows a small number of stray atoms.
00838   } else {
00839     node->pdb->find_extremes(&(xmin.x),&(xmax.x),lattice.a_r(),0.9);
00840     node->pdb->find_extremes(&(xmin.y),&(xmax.y),lattice.b_r(),0.9);
00841     node->pdb->find_extremes(&(xmin.z),&(xmax.z),lattice.c_r(),0.9);
00842   }
00843   }
00844 
00845   BigReal origin_shift;
00846   origin_shift = lattice.a_r() * lattice.origin();
00847   xmin.x -= origin_shift;
00848   xmax.x -= origin_shift;
00849   origin_shift = lattice.b_r() * lattice.origin();
00850   xmin.y -= origin_shift;
00851   xmax.y -= origin_shift;
00852   origin_shift = lattice.c_r() * lattice.origin();
00853   xmin.z -= origin_shift;
00854   xmax.z -= origin_shift;
00855 
00856   // SimParameters default is -1 for unset
00857   int twoAwayX = params->twoAwayX;
00858   int twoAwayY = params->twoAwayY;
00859   int twoAwayZ = params->twoAwayZ;
00860 
00861 #ifdef NAMD_CUDA
00862   // for CUDA be sure there are more patches than pes
00863 
00864   int numPatches = patchMap->sizeGrid(
00865         xmin,xmax,lattice,patchSize,maxNumPatches,
00866         twoAwayX>0 ? 2 : 1, twoAwayY>0 ? 2 : 1, twoAwayZ>0 ? 2 : 1);
00867   if ( numPatches < CkNumPes() && twoAwayX < 0 ) {
00868     twoAwayX = 1;
00869     numPatches = patchMap->sizeGrid(
00870         xmin,xmax,lattice,patchSize,maxNumPatches,
00871         twoAwayX>0 ? 2 : 1, twoAwayY>0 ? 2 : 1, twoAwayZ>0 ? 2 : 1);
00872   }
00873   if ( numPatches < CkNumPes() && twoAwayY < 0 ) {
00874     twoAwayY = 1;
00875     numPatches = patchMap->sizeGrid(
00876         xmin,xmax,lattice,patchSize,maxNumPatches,
00877         twoAwayX>0 ? 2 : 1, twoAwayY>0 ? 2 : 1, twoAwayZ>0 ? 2 : 1);
00878   }
00879   if ( numPatches < CkNumPes() && twoAwayZ < 0 ) {
00880     twoAwayZ = 1;
00881     numPatches = patchMap->sizeGrid(
00882         xmin,xmax,lattice,patchSize,maxNumPatches,
00883         twoAwayX>0 ? 2 : 1, twoAwayY>0 ? 2 : 1, twoAwayZ>0 ? 2 : 1);
00884   }
00885   if ( numPatches < CkNumPes() ) {
00886     NAMD_die("CUDA-enabled NAMD requires more patches than processes.");
00887   }
00888 
00889   patchMap->makePatches(xmin,xmax,lattice,patchSize,maxNumPatches,
00890         twoAwayX>0 ? 2 : 1, twoAwayY>0 ? 2 : 1, twoAwayZ>0 ? 2 : 1);
00891 
00892 #else
00893 
00894   int numPatches = patchMap->sizeGrid(
00895         xmin,xmax,lattice,patchSize,maxNumPatches,
00896         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00897   if ( numPatches > (CkNumPes() - 1) && twoAwayZ < 0 ) {
00898     twoAwayZ = 0;
00899     numPatches = patchMap->sizeGrid(
00900         xmin,xmax,lattice,patchSize,maxNumPatches,
00901         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00902   }
00903   if ( numPatches > (CkNumPes() - 1) && twoAwayY < 0 ) {
00904     twoAwayY = 0;
00905     numPatches = patchMap->sizeGrid(
00906         xmin,xmax,lattice,patchSize,maxNumPatches,
00907         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00908   }
00909   if ( numPatches > (CkNumPes() - 1) && twoAwayX < 0 ) {
00910     twoAwayX = 0;
00911     numPatches = patchMap->sizeGrid(
00912         xmin,xmax,lattice,patchSize,maxNumPatches,
00913         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00914   }
00915 
00916   patchMap->makePatches(xmin,xmax,lattice,patchSize,maxNumPatches,
00917         twoAwayX ? 2 : 1, twoAwayY ? 2 : 1, twoAwayZ ? 2 : 1);
00918 
00919 #endif
00920 
00921 }
00922 
00923 
00924 //----------------------------------------------------------------------
00925 void WorkDistrib::assignNodeToPatch()
00926 {
00927   PatchMap *patchMap = PatchMap::Object();
00928   int nNodes = Node::Object()->numNodes();
00929 
00930 #if USE_TOPOMAP 
00931   TopoManager tmgr;
00932   int numPes = tmgr.getDimNX() * tmgr.getDimNY() * tmgr.getDimNZ();
00933   if (numPes > patchMap->numPatches() && (assignPatchesTopoGridRecBisection() > 0)) {
00934     CkPrintf ("Blue Gene/L topology partitioner finished successfully \n");
00935   }
00936   else
00937 #endif
00938     if (nNodes > patchMap->numPatches())
00939       assignPatchesBitReversal();
00940     else
00941 #if (CMK_BLUEGENEP | CMK_BLUEGENEL)
00942       assignPatchesRecursiveBisection();
00943 #else
00944       assignPatchesSpaceFillingCurve();
00945 #endif
00946       // assignPatchesRoundRobin();
00947       // assignPatchesToLowestLoadNode();
00948   
00949   int *nAtoms = new int[nNodes];
00950   int numAtoms=0;
00951   int i;
00952   for(i=0; i < nNodes; i++)
00953     nAtoms[i] = 0;
00954 
00955   for(i=0; i < patchMap->numPatches(); i++)
00956   {
00957     //    iout << iINFO << "Patch " << i << " has " 
00958     //   << patchMap->patch(i)->getNumAtoms() << " atoms and "
00959     //   << patchMap->patch(i)->getNumAtoms() * 
00960     //            patchMap->patch(i)->getNumAtoms() 
00961     //   << " pairs.\n" << endi;
00962 
00963     if (patchMap->patch(i)) {
00964       numAtoms += patchMap->patch(i)->getNumAtoms();
00965       nAtoms[patchMap->node(i)] += patchMap->patch(i)->getNumAtoms();
00966     }
00967   }
00968 
00969 //  for(i=0; i < nNodes; i++)
00970 //    iout << iINFO 
00971 //       << nAtoms[i] << " atoms assigned to node " << i << "\n" << endi;
00972   if ( numAtoms != Node::Object()->molecule->numAtoms ) {
00973     NAMD_die("Incorrect atom count in WorkDistrib::assignNodeToPatch\n");
00974   }
00975 
00976   delete [] nAtoms;
00977  
00978   //  PatchMap::Object()->printPatchMap();
00979 }
00980 
00981 //----------------------------------------------------------------------
00982 // void WorkDistrib::assignPatchesSlices() 
00983 // {
00984 //   int pid; 
00985 //   int assignedNode = 0;
00986 //   PatchMap *patchMap = PatchMap::Object();
00987 //   Node *node = CLocalBranch(Node, CkpvAccess(BOCclass_group).node);
00988 
00989 //   int *numAtoms = new int[node->numNodes()];
00990 //   for (int i=0; i<node->numNodes(); i++) {
00991 //     numAtoms[i] = 0;
00992 //   }
00993 
00994 //   // Assign patch to node with least atoms assigned.
00995 //   for(pid=0; pid < patchMap->numPatches(); pid++) {
00996 //     assignedNode = 0;
00997 //     for (i=1; i < node->numNodes(); i++) {
00998 //       if (numAtoms[i] < numAtoms[assignedNode]) assignedNode = i;
00999 //     }
01000 //     patchMap->assignNode(pid, assignedNode);
01001 //     numAtoms[assignedNode] += patchMap->patch(pid)->getNumAtoms();
01002 
01003 //     /*
01004 //     iout << iINFO << "Patch (" << pid << ") has " 
01005 //       << patchMap->patch(pid)->getNumAtoms() 
01006 //       << " atoms:  Assigned to Node(" << assignedNode << ")\n" 
01007 //       << endi;
01008 //     */
01009 //   }
01010 
01011 //   delete[] numAtoms;
01012 // }
01013 
01014 //----------------------------------------------------------------------
01015 void WorkDistrib::assignPatchesToLowestLoadNode() 
01016 {
01017   int pid; 
01018   int assignedNode = 0;
01019   PatchMap *patchMap = PatchMap::Object();
01020   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01021   Node *node = nd.ckLocalBranch();
01022   SimParameters *simParams = node->simParameters;
01023 
01024   int *load = new int[node->numNodes()];
01025   int *assignedNodes = new int[patchMap->numPatches()];
01026   for (int i=0; i<node->numNodes(); i++) {
01027     load[i] = 0;
01028   }
01029 
01030   int defaultNode = 0;
01031   if ( simParams->noPatchesOnZero && node->numNodes() > 1 ){
01032     defaultNode = 1;
01033     if( simParams->noPatchesOnOne && node->numNodes() > 2)
01034       defaultNode = 2;
01035   }
01036   // Assign patch to node with least atoms assigned.
01037   for(pid=0; pid < patchMap->numPatches(); pid++) {
01038     assignedNode = defaultNode;
01039     for (int i=assignedNode + 1; i < node->numNodes(); i++) {
01040       if (load[i] < load[assignedNode]) assignedNode = i;
01041     }
01042     assignedNodes[pid] = assignedNode;
01043     load[assignedNode] += patchMap->patch(pid)->getNumAtoms() + 1;
01044   }
01045 
01046   delete[] load;
01047   sortNodesAndAssign(assignedNodes);
01048   delete[] assignedNodes;
01049 }
01050 
01051 //----------------------------------------------------------------------
01052 void WorkDistrib::assignPatchesBitReversal() 
01053 {
01054   int pid; 
01055   PatchMap *patchMap = PatchMap::Object();
01056   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01057   Node *node = nd.ckLocalBranch();
01058 
01059   int ncpus = node->numNodes();
01060   int npatches = patchMap->numPatches();
01061   if ( ncpus <= npatches )
01062     NAMD_bug("WorkDistrib::assignPatchesBitReversal called improperly");
01063 
01064   // find next highest power of two
01065   int npow2 = 1;  int nbits = 0;
01066   while ( npow2 < ncpus ) { npow2 *= 2; nbits += 1; }
01067 
01068   // build bit reversal sequence
01069   SortableResizeArray<int> seq(ncpus);
01070   // avoid using node 0 (reverse of 0 is 0 so start at 1)
01071   int i = 1;
01072   for ( int icpu=0; icpu<(ncpus-1); ++icpu ) {
01073     int ri;
01074     for ( ri = ncpus; ri >= ncpus; ++i ) {
01075       ri = 0;
01076       int pow2 = 1;
01077       int rpow2 = npow2 / 2;
01078       for ( int j=0; j<nbits; ++j ) {
01079         ri += rpow2 * ( ( i / pow2 ) % 2 );
01080         pow2 *= 2;  rpow2 /= 2;
01081       }
01082     }
01083     seq[icpu] = ri;
01084   }
01085 
01086   // extract and sort patch locations
01087   sortNodesAndAssign(seq.begin());
01088   if ( ncpus > 2*npatches ) sortNodesAndAssign(seq.begin()+npatches, 1);
01089 }
01090 
01091 //----------------------------------------------------------------------
01092 struct nodesort {
01093   int node;
01094   int a_total;
01095   int b_total;
01096   int c_total;
01097   int npatches;
01098   nodesort() : node(-1),a_total(0),b_total(0),c_total(0),npatches(0) { ; }
01099   int operator==(const nodesort &o) const {
01100     float a1 = ((float)a_total)/((float)npatches);
01101     float a2 = ((float)o.a_total)/((float)o.npatches);
01102     float b1 = ((float)b_total)/((float)npatches);
01103     float b2 = ((float)o.b_total)/((float)o.npatches);
01104     float c1 = ((float)c_total)/((float)npatches);
01105     float c2 = ((float)o.c_total)/((float)o.npatches);
01106     return ((a1 == a2) && (b1 == b2) && (c1 == c2));
01107   }
01108   int operator<(const nodesort &o) const {
01109     float a1 = ((float)a_total)/((float)npatches);
01110     float a2 = ((float)o.a_total)/((float)o.npatches);
01111     float b1 = ((float)b_total)/((float)npatches);
01112     float b2 = ((float)o.b_total)/((float)o.npatches);
01113     float c1 = ((float)c_total)/((float)npatches);
01114     float c2 = ((float)o.c_total)/((float)o.npatches);
01115     return ( (a1 < a2) || ((a1 == a2) && (b1 < b2)) ||
01116                 ((a1 == a2) && (b1 == b2) && (c1 < c2)) );
01117   }
01118 };
01119 
01120 
01121 void WorkDistrib::sortNodesAndAssign(int *assignedNode, int baseNodes) {
01122   // if baseNodes is zero (default) then set both nodes and basenodes
01123   // if baseNodes is nonzero then this is a second call to set basenodes only
01124   int i, pid; 
01125   PatchMap *patchMap = PatchMap::Object();
01126   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01127   Node *node = nd.ckLocalBranch();
01128   int nnodes = node->numNodes();
01129   int npatches = patchMap->numPatches();
01130 
01131   ResizeArray<nodesort> allnodes(nnodes);
01132   for ( i=0; i < nnodes; ++i ) {
01133     allnodes[i].node = i;
01134   }
01135   for ( pid=0; pid<npatches; ++pid ) {
01136     // iout << pid << " " << assignedNode[pid] << "\n" << endi;
01137     allnodes[assignedNode[pid]].npatches++;
01138     allnodes[assignedNode[pid]].a_total += patchMap->index_a(pid);
01139     allnodes[assignedNode[pid]].b_total += patchMap->index_b(pid);
01140     allnodes[assignedNode[pid]].c_total += patchMap->index_c(pid);
01141   }
01142   SortableResizeArray<nodesort> usednodes(nnodes);
01143   usednodes.resize(0);
01144   for ( i=0; i < nnodes; ++i ) {
01145     if ( allnodes[i].npatches ) usednodes.add(allnodes[i]);
01146   }
01147   usednodes.sort();
01148   int nused = usednodes.size();
01149   int i2 = nused/2;
01150   for ( i=0; i < nnodes; ++i ) {
01151     if ( allnodes[i].npatches ) allnodes[usednodes[i2++].node].node = i;
01152     if ( i2 == nused ) i2 = 0;
01153   }
01154 
01155   for ( pid=0; pid<npatches; ++pid ) {
01156     // iout << pid << " " <<  allnodes[assignedNode[pid]].node << "\n" << endi;
01157     if ( ! baseNodes ) {
01158       patchMap->assignNode(pid, allnodes[assignedNode[pid]].node);
01159     }
01160     patchMap->assignBaseNode(pid, allnodes[assignedNode[pid]].node);
01161   }
01162 }
01163 
01164 //----------------------------------------------------------------------
01165 void WorkDistrib::assignPatchesRoundRobin() 
01166 {
01167   int pid; 
01168   PatchMap *patchMap = PatchMap::Object();
01169   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01170   Node *node = nd.ckLocalBranch();
01171   int *assignedNode = new int[patchMap->numPatches()];
01172 
01173   for(pid=0; pid < patchMap->numPatches(); pid++) {
01174     assignedNode[pid] = pid % node->numNodes();
01175   }
01176 
01177   sortNodesAndAssign(assignedNode);
01178   delete [] assignedNode;
01179 }
01180 
01181 //----------------------------------------------------------------------
01182 void WorkDistrib::assignPatchesRecursiveBisection() 
01183 {
01184   PatchMap *patchMap = PatchMap::Object();
01185   int *assignedNode = new int[patchMap->numPatches()];
01186   int numNodes = Node::Object()->numNodes();
01187   SimParameters *simParams = Node::Object()->simParameters;
01188   int usedNodes = numNodes;
01189   int unusedNodes = 0;
01190   if ( simParams->noPatchesOnZero && numNodes > 1 ){
01191     usedNodes -= 1;
01192     if(simParams->noPatchesOnOne && numNodes > 2)
01193       usedNodes -= 1;
01194   }  
01195   unusedNodes = numNodes - usedNodes;
01196   RecBisection recBisec(usedNodes,PatchMap::Object());
01197   if ( recBisec.partition(assignedNode) ) {
01198     if ( unusedNodes !=0 ) {
01199       for ( int i=0; i<patchMap->numPatches(); ++i ) {
01200         assignedNode[i] += unusedNodes;
01201       }
01202     }
01203     sortNodesAndAssign(assignedNode);
01204     delete [] assignedNode;
01205   } else {
01206     //free the array here since a same array will be allocated
01207     //in assignPatchesToLowestLoadNode function, thus reducting
01208     //temporary memory usage
01209     delete [] assignedNode; 
01210     
01211     iout << iWARN 
01212          << "WorkDistrib: Recursive bisection fails,"
01213          << "invoking least-load algorithm\n";
01214     assignPatchesToLowestLoadNode();
01215   }
01216 }
01217 
01218 //----------------------------------------------------------------------
01219 void WorkDistrib::assignPatchesSpaceFillingCurve() 
01220 {
01221   PatchMap *patchMap = PatchMap::Object();
01222   int *assignedNode = new int[patchMap->numPatches()];
01223   int numNodes = Node::Object()->numNodes();
01224   SimParameters *simParams = Node::Object()->simParameters;
01225   int usedNodes = numNodes;
01226   int unusedNodes = 0;
01227   if ( simParams->noPatchesOnZero && numNodes > 1 ){
01228     usedNodes -= 1;
01229     if(simParams->noPatchesOnOne && numNodes > 2)
01230       usedNodes -= 1;
01231   }  
01232   unusedNodes = numNodes - usedNodes;
01233 
01234   int numPatches = patchMap->numPatches();
01235   if ( numPatches < usedNodes )
01236     NAMD_bug("WorkDistrib::assignPatchesSpaceFillingCurve() called with more nodes than patches");
01237 
01238   ResizeArray<double> patchLoads(numPatches);
01239   SortableResizeArray<double> sortedLoads(numPatches);
01240   for ( int i=0; i<numPatches; ++i ) {
01241     double load = patchMap->patch(i)->getNumAtoms() + 10;
01242     patchLoads[i] = load;
01243     sortedLoads[i] = load;
01244   }
01245   sortedLoads.sort();
01246 
01247   // limit maxPatchLoad to adjusted average load per node
01248   double sumLoad = 0;
01249   double maxPatchLoad = 1;
01250   for ( int i=0; i<numPatches; ++i ) {
01251     double load = sortedLoads[i];
01252     double total = sumLoad + (numPatches-i) * load;
01253     if ( usedNodes * load > total ) break;
01254     sumLoad += load;
01255     maxPatchLoad = load;
01256   }
01257   double totalLoad = 0;
01258   for ( int i=0; i<numPatches; ++i ) {
01259     if ( patchLoads[i] > maxPatchLoad ) patchLoads[i] = maxPatchLoad;
01260     totalLoad += patchLoads[i];
01261   }
01262   if ( usedNodes * maxPatchLoad > totalLoad )
01263     NAMD_bug("algorithm failure in WorkDistrib::assignPatchesSpaceFillingCurve()");
01264 
01265   // walk through patches in space-filling curve
01266   sumLoad = 0;
01267   int node = 0;
01268   int adim = patchMap->gridsize_a();
01269   int bdim = patchMap->gridsize_b();
01270   int cdim = patchMap->gridsize_c();
01271   int b = 0;
01272   int c = 0;
01273   int binc = 1;
01274   int cinc = 1;
01275   for ( int a = 0; a < adim; ++a ) {
01276     while ( b >= 0 && b < bdim ) {
01277       while ( c >= 0 && c < cdim ) {
01278         int pid = patchMap->pid(a,b,c);
01279         assignedNode[pid] = node;
01280         sumLoad += patchLoads[pid];
01281         double targetLoad = (double)(node+1) / (double)usedNodes;
01282         targetLoad *= totalLoad;
01283         if ( node+1 < usedNodes && sumLoad >= targetLoad ) ++node;
01284         c += cinc;
01285       }
01286       cinc *= -1;  c += cinc;
01287       b += binc;
01288     }
01289     binc *= -1;  b += binc;
01290   }
01291 
01292   for ( int i=0; i<patchMap->numPatches(); ++i ) {
01293     assignedNode[i] += unusedNodes;
01294   }
01295   sortNodesAndAssign(assignedNode);
01296   delete [] assignedNode; 
01297 }
01298 
01299 //----------------------------------------------------------------------
01300 void WorkDistrib::mapComputes(void)
01301 {
01302   PatchMap *patchMap = PatchMap::Object();
01303   ComputeMap *computeMap = ComputeMap::Object();
01304   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01305   Node *node = nd.ckLocalBranch();
01306 
01307   DebugM(3,"Mapping computes\n");
01308 
01309   computeMap->allocateCids();
01310 
01311   // Handle full electrostatics
01312   if ( node->simParameters->fullDirectOn )
01313     mapComputeHomePatches(computeFullDirectType);
01314   if ( node->simParameters->FMAOn )
01315 #ifdef DPMTA
01316     mapComputeHomePatches(computeDPMTAType);
01317 #else
01318     NAMD_die("This binary does not include DPMTA (FMA).");
01319 #endif
01320   if ( node->simParameters->PMEOn ) {
01321 #ifdef DPME
01322     if ( node->simParameters->useDPME )
01323       mapComputeHomePatches(computeDPMEType);
01324     else {
01325       if (node->simParameters->useOptPME) {
01326         mapComputeHomePatches(optPmeType);
01327         if ( node->simParameters->pressureProfileEwaldOn )
01328           mapComputeHomePatches(computeEwaldType);
01329       }
01330       else {
01331         mapComputeHomePatches(computePmeType);
01332         if ( node->simParameters->pressureProfileEwaldOn )
01333           mapComputeHomePatches(computeEwaldType);
01334       }
01335     }
01336 #else
01337     if (node->simParameters->useOptPME) {
01338       mapComputeHomePatches(optPmeType);
01339       if ( node->simParameters->pressureProfileEwaldOn )
01340         mapComputeHomePatches(computeEwaldType);
01341     }
01342     else {      
01343       mapComputeHomePatches(computePmeType);
01344       if ( node->simParameters->pressureProfileEwaldOn )
01345         mapComputeHomePatches(computeEwaldType);
01346     }
01347 #endif
01348   }
01349 
01350   if ( node->simParameters->globalForcesOn ) {
01351     DebugM(2,"adding ComputeGlobal\n");
01352     mapComputeHomePatches(computeGlobalType);
01353   }
01354 
01355   if ( node->simParameters->extForcesOn )
01356     mapComputeHomePatches(computeExtType);
01357 
01358 #ifdef NAMD_CUDA
01359   mapComputeNode(computeNonbondedCUDAType);
01360   mapComputeHomeTuples(computeExclsType);
01361   mapComputePatch(computeSelfExclsType);
01362 #endif
01363 
01364   mapComputeNonbonded();
01365 
01366   // If we're doing true pair interactions, no need for bonded terms.
01367   // But if we're doing within-group interactions, we do need them.
01368   if ( !node->simParameters->pairInteractionOn || 
01369       node->simParameters->pairInteractionSelf) { 
01370     mapComputeHomeTuples(computeBondsType);
01371     mapComputeHomeTuples(computeAnglesType);
01372     mapComputeHomeTuples(computeDihedralsType);
01373     mapComputeHomeTuples(computeImpropersType);
01374     mapComputeHomeTuples(computeCrosstermsType);
01375     mapComputePatch(computeSelfBondsType);
01376     mapComputePatch(computeSelfAnglesType);
01377     mapComputePatch(computeSelfDihedralsType);
01378     mapComputePatch(computeSelfImpropersType);
01379     mapComputePatch(computeSelfCrosstermsType);
01380   }
01381 
01382 
01383   if ( node->simParameters->eFieldOn )
01384     mapComputePatch(computeEFieldType);
01385   /* BEGIN gf */
01386   if ( node->simParameters->mgridforceOn )
01387     mapComputeHomePatches(computeGridForceType);
01388   /* END gf */
01389   if ( node->simParameters->stirOn )
01390     mapComputePatch(computeStirType);
01391   if ( node->simParameters->sphericalBCOn )
01392     mapComputePatch(computeSphericalBCType);
01393   if ( node->simParameters->cylindricalBCOn )
01394     mapComputePatch(computeCylindricalBCType);
01395   if ( node->simParameters->tclBCOn ) {
01396     mapComputeHomePatches(computeTclBCType);
01397   }
01398   if ( node->simParameters->constraintsOn )
01399     mapComputePatch(computeRestraintsType);
01400   if ( node->simParameters->consForceOn )
01401     mapComputePatch(computeConsForceType);
01402   if ( node->simParameters->consTorqueOn )
01403     mapComputePatch(computeConsTorqueType);
01404 }
01405 
01406 //----------------------------------------------------------------------
01407 void WorkDistrib::mapComputeHomeTuples(ComputeType type)
01408 {
01409   PatchMap *patchMap = PatchMap::Object();
01410   ComputeMap *computeMap = ComputeMap::Object();
01411   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01412   Node *node = nd.ckLocalBranch();
01413 
01414   int numNodes = node->numNodes();
01415   int numPatches = patchMap->numPatches();
01416   ComputeID cid;
01417   PatchIDList basepids;
01418 
01419   for(int i=0; i<numNodes; i++) {
01420     patchMap->basePatchIDList(i,basepids);
01421     if ( basepids.size() ) {
01422       cid=computeMap->storeCompute(i,basepids.size(),type);
01423       for(int j=0; j<basepids.size(); ++j) {
01424         patchMap->newCid(basepids[j],cid);
01425         computeMap->newPid(cid,basepids[j]);
01426       }
01427     }
01428   }
01429 }
01430 
01431 //----------------------------------------------------------------------
01432 void WorkDistrib::mapComputeHomePatches(ComputeType type)
01433 {
01434   PatchMap *patchMap = PatchMap::Object();
01435   ComputeMap *computeMap = ComputeMap::Object();
01436   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01437   Node *node = nd.ckLocalBranch();
01438 
01439   int numNodes = node->numNodes();
01440   int numPatches = patchMap->numPatches();
01441   ComputeID *cid = new ComputeID[numNodes];
01442 
01443   for(int i=0; i<numNodes; i++) {
01444     if ( patchMap->numPatchesOnNode(i) ) {
01445       cid[i]=computeMap->storeCompute(i,patchMap->numPatchesOnNode(i),type);
01446     }
01447   }
01448 
01449   PatchID j;
01450 
01451   for(j=0;j<numPatches;j++)
01452   {
01453     patchMap->newCid(j,cid[patchMap->node(j)]);
01454     computeMap->newPid(cid[patchMap->node(j)],j);
01455   }
01456 
01457   delete [] cid;
01458 }
01459 
01460 //----------------------------------------------------------------------
01461 void WorkDistrib::mapComputePatch(ComputeType type)
01462 {
01463   PatchMap *patchMap = PatchMap::Object();
01464   ComputeMap *computeMap = ComputeMap::Object();
01465 
01466   PatchID i;
01467   ComputeID cid;
01468 
01469   for(i=0; i<patchMap->numPatches(); i++)
01470   {
01471     cid=computeMap->storeCompute(patchMap->node(i),1,type);
01472     computeMap->newPid(cid,i);
01473     patchMap->newCid(i,cid);
01474   }
01475 
01476 }
01477 
01478 
01479 //----------------------------------------------------------------------
01480 void WorkDistrib::mapComputeNode(ComputeType type)
01481 {
01482   PatchMap *patchMap = PatchMap::Object();
01483   ComputeMap *computeMap = ComputeMap::Object();
01484 
01485   PatchID i;
01486   ComputeID cid;
01487 
01488   for(int i=0; i<CkNumPes(); i++) {
01489     computeMap->storeCompute(i,0,type);
01490   }
01491 
01492 }
01493 
01494 
01495 //----------------------------------------------------------------------
01496 void WorkDistrib::mapComputeNonbonded(void)
01497 {
01498   // For each patch, create 1 electrostatic object for self-interaction.
01499   // Then create 1 for each 1-away and 2-away neighbor which has a larger
01500   // pid.
01501 
01502   PatchMap *patchMap = PatchMap::Object();
01503   ComputeMap *computeMap = ComputeMap::Object();
01504   CProxy_Node nd(CkpvAccess(BOCclass_group).node);
01505   Node *node = nd.ckLocalBranch();
01506 
01507   PatchID oneAway[PatchMap::MaxOneOrTwoAway];
01508   PatchID oneAwayTrans[PatchMap::MaxOneOrTwoAway];
01509 
01510   PatchID i;
01511   ComputeID cid;
01512   int numNeighbors;
01513   int j;
01514 
01515   for(i=0; i<patchMap->numPatches(); i++) // do the self 
01516   {
01517     int64 numAtoms = patchMap->patch(i)->getNumAtoms();  // avoid overflow
01518     int64 numFixed = patchMap->patch(i)->getNumFixedAtoms();  // avoid overflow
01519     int numPartitions = 0;
01520     int divide = node->simParameters->numAtomsSelf;
01521     if (divide > 0) {
01522       numPartitions = (int) ( 0.5 +
01523         (numAtoms*numAtoms-numFixed*numFixed) / (double)(2*divide*divide) );
01524     }
01525     if (numPartitions < 1) numPartitions = 1;
01526     if ( numPartitions > node->simParameters->maxSelfPart )
01527                         numPartitions = node->simParameters->maxSelfPart;
01528     // self-interaction
01529     DebugM(4,"Mapping " << numPartitions << " ComputeNonbondedSelf objects for patch " << i << "\n");
01530 //    iout <<"Self numPartitions = " <<numPartitions <<" numAtoms " <<numAtoms <<std::endl;
01531     for(int partition=0; partition < numPartitions; partition++)
01532     {
01533       cid=computeMap->storeCompute(patchMap->node(i),1,
01534                                    computeNonbondedSelfType,
01535                                    partition,numPartitions);
01536       computeMap->newPid(cid,i);
01537       patchMap->newCid(i,cid);
01538     }
01539   }
01540 
01541   for(int p1=0; p1 <patchMap->numPatches(); p1++) // do the pairs
01542   {
01543     // this only returns half of neighbors, which is what we want
01544     numNeighbors=patchMap->oneOrTwoAwayNeighbors(p1,oneAway,oneAwayTrans);
01545     for(j=0;j<numNeighbors;j++)
01546     {
01547         int p2 = oneAway[j];
01548         int64 numAtoms1 = patchMap->patch(p1)->getNumAtoms();
01549         int64 numAtoms2 = patchMap->patch(p2)->getNumAtoms();
01550         int64 numFixed1 = patchMap->patch(p1)->getNumFixedAtoms();
01551         int64 numFixed2 = patchMap->patch(p2)->getNumFixedAtoms();
01552         const int nax = patchMap->numaway_a();  // 1 or 2
01553         const int nay = patchMap->numaway_b();  // 1 or 2
01554         const int naz = patchMap->numaway_c();  // 1 or 2
01555         const int ia1 = patchMap->index_a(p1);
01556         const int ia2 = patchMap->index_a(p2);
01557         const int ib1 = patchMap->index_b(p1);
01558         const int ib2 = patchMap->index_b(p2);
01559         const int ic1 = patchMap->index_c(p1);
01560         const int ic2 = patchMap->index_c(p2);
01561         int distance = 3;
01562         if ( ia1 == ia2 ) --distance;
01563         else if ( ia1 == ia2 + nax - 1 ) --distance;
01564         else if ( ia1 + nax - 1 == ia2 ) --distance;
01565         if ( ib1 == ib2 ) --distance;
01566         else if ( ib1 == ib2 + nay - 1 ) --distance;
01567         else if ( ib1 + nay - 1 == ib2 ) --distance;
01568         if ( ic1 == ic2 ) --distance;
01569         else if ( ic1 == ic2 + naz - 1 ) --distance;
01570         else if ( ic1 + naz - 1 == ic2 ) --distance;
01571         int divide = 0;
01572         if ( distance == 0 ) {
01573           divide = node->simParameters->numAtomsSelf2;
01574         } else if (distance == 1) {
01575           divide = node->simParameters->numAtomsPair;
01576         } else {
01577           divide = node->simParameters->numAtomsPair2;
01578         }
01579         int numPartitions = 0;
01580         if (divide > 0) {
01581           numPartitions = (int) ( 0.5 +
01582             (numAtoms1*numAtoms2-numFixed1*numFixed2)/(double)(divide*divide) );
01583         }
01584         if ( numPartitions < 1 ) numPartitions = 1;
01585         if ( numPartitions > node->simParameters->maxPairPart )
01586                         numPartitions = node->simParameters->maxPairPart;
01587 //      if ( numPartitions > 1 ) iout << "Mapping " << numPartitions << " ComputeNonbondedPair objects for patches " << p1 << "(" << numAtoms1 << ") and " << p2 << "(" << numAtoms2 << ")\n" << endi;
01588         for(int partition=0; partition < numPartitions; partition++)
01589         {
01590           cid=computeMap->storeCompute(
01591                 patchMap->basenode(patchMap->downstream2(p1,p2)),
01592                 2,computeNonbondedPairType,partition,numPartitions);
01593           computeMap->newPid(cid,p1);
01594           computeMap->newPid(cid,p2,oneAwayTrans[j]);
01595           patchMap->newCid(p1,cid);
01596           patchMap->newCid(p2,cid);
01597         }
01598     }
01599   }
01600 
01601 }
01602 
01603 //----------------------------------------------------------------------
01604 void WorkDistrib::messageEnqueueWork(Compute *compute) {
01605   LocalWorkMsg *msg = compute->localWorkMsg;
01606   int seq = compute->sequence();
01607 
01608   if ( seq < 0 ) {
01609     NAMD_bug("compute->sequence() < 0 in WorkDistrib::messageEnqueueWork");
01610   } else {
01611     SET_PRIORITY(msg,seq,compute->priority());
01612   }
01613 
01614   msg->compute = compute; // pointer is valid since send is to local Pe
01615   int type = compute->type();
01616 
01617   CProxy_WorkDistrib wdProxy(CkpvAccess(BOCclass_group).workDistrib);
01618   switch ( type ) {
01619   case computeExclsType:
01620   case computeSelfExclsType:
01621     wdProxy[CkMyPe()].enqueueExcls(msg);
01622     break;
01623   case computeBondsType:
01624   case computeSelfBondsType:
01625     wdProxy[CkMyPe()].enqueueBonds(msg);
01626     break;
01627   case computeAnglesType:
01628   case computeSelfAnglesType:
01629     wdProxy[CkMyPe()].enqueueAngles(msg);
01630     break;
01631   case computeDihedralsType:
01632   case computeSelfDihedralsType:
01633     wdProxy[CkMyPe()].enqueueDihedrals(msg);
01634     break;
01635   case computeImpropersType:
01636   case computeSelfImpropersType:
01637     wdProxy[CkMyPe()].enqueueImpropers(msg);
01638     break;
01639   case computeCrosstermsType:
01640   case computeSelfCrosstermsType:
01641     wdProxy[CkMyPe()].enqueueCrossterms(msg);
01642     break;
01643   case computeNonbondedSelfType:
01644     switch ( seq % 2 ) {
01645     case 0:
01646       wdProxy[CkMyPe()].enqueueSelfA(msg);
01647       break;
01648     case 1:
01649       wdProxy[CkMyPe()].enqueueSelfB(msg);
01650       break;
01651     default:
01652       NAMD_bug("WorkDistrib::messageEnqueueSelf case statement error!");
01653     }
01654     break;
01655   case computeNonbondedPairType:
01656     switch ( seq % 2 ) {
01657     case 0:
01658       wdProxy[CkMyPe()].enqueueWorkA(msg);
01659       break;
01660     case 1:
01661       wdProxy[CkMyPe()].enqueueWorkB(msg);
01662       break;
01663     case 2:
01664       wdProxy[CkMyPe()].enqueueWorkC(msg);
01665       break;
01666     default:
01667       NAMD_bug("WorkDistrib::messageEnqueueWork case statement error!");
01668     }
01669     break;
01670   case computeNonbondedCUDAType:
01671 #ifdef NAMD_CUDA
01672     // CkPrintf("CUDA %d %d %x\n", CkMyPe(), seq, compute->priority());
01673     wdProxy[CkMyPe()].enqueueCUDA(msg);
01674 #else
01675     msg->compute->doWork();
01676 #endif
01677     break;
01678   case computePmeType:
01679     // CkPrintf("PME %d %d %x\n", CkMyPe(), seq, compute->priority());
01680 #ifdef NAMD_CUDA
01681     wdProxy[CkMyPe()].enqueuePme(msg);
01682 #else
01683     msg->compute->doWork();
01684 #endif
01685     break;
01686   case optPmeType:
01687     // CkPrintf("PME %d %d %x\n", CkMyPe(), seq, compute->priority());
01688 #ifdef NAMD_CUDA
01689     wdProxy[CkMyPe()].enqueuePme(msg);
01690 #else
01691     msg->compute->doWork();
01692 #endif
01693     break;
01694   default:
01695     wdProxy[CkMyPe()].enqueueWork(msg);
01696   }
01697 }
01698 
01699 void WorkDistrib::enqueueWork(LocalWorkMsg *msg) {
01700   msg->compute->doWork();
01701   if ( msg->compute->localWorkMsg != msg )
01702     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01703 }
01704 
01705 void WorkDistrib::enqueueExcls(LocalWorkMsg *msg) {
01706   msg->compute->doWork();
01707   if ( msg->compute->localWorkMsg != msg )
01708     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01709 }
01710 
01711 void WorkDistrib::enqueueBonds(LocalWorkMsg *msg) {
01712   msg->compute->doWork();
01713   if ( msg->compute->localWorkMsg != msg )
01714     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01715 }
01716 
01717 void WorkDistrib::enqueueAngles(LocalWorkMsg *msg) {
01718   msg->compute->doWork();
01719   if ( msg->compute->localWorkMsg != msg )
01720     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01721 }
01722 
01723 void WorkDistrib::enqueueDihedrals(LocalWorkMsg *msg) {
01724   msg->compute->doWork();
01725   if ( msg->compute->localWorkMsg != msg )
01726     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01727 }
01728 
01729 void WorkDistrib::enqueueImpropers(LocalWorkMsg *msg) {
01730   msg->compute->doWork();
01731   if ( msg->compute->localWorkMsg != msg )
01732     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01733 }
01734 
01735 void WorkDistrib::enqueueCrossterms(LocalWorkMsg *msg) {
01736   msg->compute->doWork();
01737   if ( msg->compute->localWorkMsg != msg )
01738     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01739 }
01740 
01741 void WorkDistrib::enqueuePme(LocalWorkMsg *msg) {
01742   msg->compute->doWork();
01743   if ( msg->compute->localWorkMsg != msg )
01744     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01745 }
01746 
01747 void WorkDistrib::enqueueSelfA(LocalWorkMsg *msg) {
01748   msg->compute->doWork();
01749   if ( msg->compute->localWorkMsg != msg )
01750     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01751 }
01752 
01753 void WorkDistrib::enqueueSelfB(LocalWorkMsg *msg) {
01754   msg->compute->doWork();
01755   if ( msg->compute->localWorkMsg != msg )
01756     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01757 }
01758 
01759 void WorkDistrib::enqueueWorkA(LocalWorkMsg *msg) {
01760   msg->compute->doWork();
01761   if ( msg->compute->localWorkMsg != msg )
01762     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01763 }
01764 
01765 void WorkDistrib::enqueueWorkB(LocalWorkMsg *msg) {
01766   msg->compute->doWork();
01767   if ( msg->compute->localWorkMsg != msg )
01768     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01769 }
01770 
01771 void WorkDistrib::enqueueWorkC(LocalWorkMsg *msg) {
01772   msg->compute->doWork();
01773   if ( msg->compute->localWorkMsg != msg )
01774     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01775 }
01776 
01777 void WorkDistrib::enqueueCUDA(LocalWorkMsg *msg) {
01778   msg->compute->doWork();
01779   if ( msg->compute->localWorkMsg != msg )
01780     NAMD_bug("WorkDistrib LocalWorkMsg recycling failed!");
01781 }
01782 
01783 
01784 //**********************************************************************
01785 //
01786 //                      FUNCTION velocities_from_PDB
01787 //
01788 //   INPUTS:
01789 //      v - Array of vectors to populate
01790 //      filename - name of the PDB filename to read in
01791 //
01792 //      This function reads in a set of initial velocities from a
01793 //      PDB file.  It places the velocities into the array of Vectors
01794 //      passed to it.
01795 //
01796 //***********************************************************************/
01797 
01798 void WorkDistrib::velocities_from_PDB(char *filename, 
01799                                       Vector *v, int totalAtoms)
01800 {
01801   PDB *v_pdb;           //  PDB info from velocity PDB
01802   int i;
01803 
01804   //  Read the PDB
01805   v_pdb = new PDB(filename);
01806   if ( v_pdb == NULL )
01807   {
01808     NAMD_die("memory allocation failed in Node::velocities_from_PDB");
01809   }
01810 
01811   //  Make sure the number of velocities read in matches
01812   //  the number of atoms we have
01813   if (v_pdb->num_atoms() != totalAtoms)
01814   {
01815     char err_msg[129];
01816 
01817     sprintf(err_msg, "FOUND %d COORDINATES IN VELOCITY PDB!!",
01818             v_pdb->num_atoms());
01819 
01820     NAMD_die(err_msg);
01821   }
01822 
01823   //  Get the entire list of atom info and loop through
01824   //  them assigning the velocity vector for each one
01825   v_pdb->get_all_positions(v);
01826 
01827   for (i=0; i<totalAtoms; i++)
01828   {
01829     v[i].x *= PDBVELINVFACTOR;
01830     v[i].y *= PDBVELINVFACTOR;
01831     v[i].z *= PDBVELINVFACTOR;
01832   }
01833 
01834   delete v_pdb;
01835 }
01836 //              END OF FUNCTION velocities_from_PDB
01837 
01838 //**********************************************************************
01839 //
01840 //                      FUNCTION velocities_from_binfile
01841 //
01842 //    INPUTS:
01843 //      fname - File name to write velocities to
01844 //      n - Number of atoms in system
01845 //      vels - Array of velocity vectors
01846 //                                      
01847 //      This function writes out the velocities in binary format.  This is
01848 //     done to preserve accuracy between restarts of namd.
01849 //
01850 //**********************************************************************
01851 
01852 void WorkDistrib::velocities_from_binfile(char *fname, Vector *vels, int n)
01853 {
01854   read_binary_file(fname,vels,n);
01855 }
01856 //               END OF FUNCTION velocities_from_binfile
01857 
01858 //**********************************************************************
01859 //
01860 //                      FUNCTION random_velocities
01861 //
01862 //   INPUTS:
01863 //      v - array of vectors to populate
01864 //      Temp - Temperature to acheive
01865 //
01866 //      This function assigns a random velocity distribution to a
01867 //   simulation to achieve a desired initial temperature.  The method
01868 //   used here was stolen from the program X-PLOR.
01869 //
01870 //**********************************************************************
01871 
01872 void WorkDistrib::random_velocities(BigReal Temp,Molecule *structure,
01873                                     Vector *v, int totalAtoms)
01874 {
01875   int i, j;             //  Loop counter
01876   BigReal kbT;          //  Boltzman constant * Temp
01877   BigReal randnum;      //  Random number from -6.0 to 6.0
01878   BigReal kbToverM;     //  sqrt(Kb*Temp/Mass)
01879   SimParameters *simParams = Node::Object()->simParameters;
01880   Bool lesOn = simParams->lesOn;
01881   Random vel_random(simParams->randomSeed);
01882 
01883   int lesReduceTemp = lesOn && simParams->lesReduceTemp;
01884   BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01885 
01886   kbT = Temp*BOLTZMAN;
01887 
01888   //  Loop through all the atoms and assign velocities in
01889   //  the x, y and z directions for each one
01890   for (i=0; i<totalAtoms; i++)
01891   {
01892     if (structure->atommass(i) <= 0.) {
01893       kbToverM = 0.;
01894     } else {
01895       kbToverM = sqrt(kbT *
01896         ( lesOn && structure->get_fep_type(i) ? tempFactor : 1.0 ) /
01897                           structure->atommass(i) );
01898     }
01899 
01900     //  The following comment was stolen from X-PLOR where
01901     //  the following section of code was adapted from.
01902     
01903     //  This section generates a Gaussian random
01904     //  deviate of 0.0 mean and standard deviation RFD for
01905     //  each of the three spatial dimensions.
01906     //  The algorithm is a "sum of uniform deviates algorithm"
01907     //  which may be found in Abramowitz and Stegun,
01908     //  "Handbook of Mathematical Functions", pg 952.
01909     for (randnum=0.0, j=0; j<12; j++)
01910     {
01911       randnum += vel_random.uniform();
01912     }
01913 
01914     randnum -= 6.0;
01915 
01916     v[i].x = randnum*kbToverM;
01917 
01918     for (randnum=0.0, j=0; j<12; j++)
01919     {
01920       randnum += vel_random.uniform();
01921     }
01922 
01923     randnum -= 6.0;
01924 
01925     v[i].y = randnum*kbToverM;
01926 
01927     for (randnum=0.0, j=0; j<12; j++)
01928     {
01929       randnum += vel_random.uniform();
01930     }
01931 
01932     randnum -= 6.0;
01933     
01934     v[i].z = randnum*kbToverM;
01935   }
01936 }
01937 /*                      END OF FUNCTION random_velocities               */
01938 
01939 //**********************************************************************
01940 //
01941 //                      FUNCTION remove_com_motion
01942 //
01943 //   INPUTS:
01944 //      vel - Array of initial velocity vectors
01945 //
01946 //      This function removes the center of mass motion from a molecule.
01947 //
01948 //**********************************************************************
01949 
01950 void WorkDistrib::remove_com_motion(Vector *vel, Molecule *structure, int n)
01951 {
01952   Vector mv(0,0,0);             //  Sum of (mv)_i
01953   BigReal totalMass=0;  //  Total mass of system
01954   int i;                        //  Loop counter
01955 
01956   //  Loop through and compute the net momentum
01957   for (i=0; i<n; i++)
01958   {
01959     BigReal mass = structure->atommass(i);
01960     mv += mass * vel[i];
01961     totalMass += mass;
01962   }
01963 
01964   mv /= totalMass;
01965 
01966   iout << iINFO << "REMOVING COM VELOCITY "
01967         << ( PDBVELFACTOR * mv ) << "\n" << endi;
01968 
01969   for (i=0; i<n; i++) { vel[i] -= mv; }
01970 
01971 }
01972 /*                      END OF FUNCTION remove_com_motion               */
01973 
01974 #if USE_TOPOMAP 
01975 
01976 //Specifically designed for BGL and other 3d Tori architectures
01977 //Partition Torus and Patch grid together using recursive bisection.
01978 int WorkDistrib::assignPatchesTopoGridRecBisection() {
01979   
01980   PatchMap *patchMap = PatchMap::Object();
01981   int *assignedNode = new int[patchMap->numPatches()];
01982   int numNodes = Node::Object()->numNodes();
01983   SimParameters *simParams = Node::Object()->simParameters;
01984   int usedNodes = numNodes;
01985   
01986   if ( simParams->noPatchesOnZero && numNodes > 1 ) usedNodes -= 1;
01987   RecBisection recBisec(patchMap->numPatches(), PatchMap::Object());
01988   
01989   int xsize = 0, ysize = 0, zsize = 0;
01990   
01991   // Right now assumes a T*** (e.g. TXYZ) mapping
01992   TopoManager tmgr;
01993   xsize = tmgr.getDimNX();
01994   ysize = tmgr.getDimNY();
01995   zsize = tmgr.getDimNZ();
01996   
01997   //Fix to not assign patches to processor 0
01998   int rc = recBisec.partitionProcGrid(xsize, ysize, zsize, assignedNode);
01999  
02000   delete [] assignedNode;
02001 
02002   return rc;
02003 }
02004 #endif
02005 
02006 #include "WorkDistrib.def.h"
02007 

Generated on Mon Nov 23 04:59:26 2009 for NAMD by  doxygen 1.3.9.1