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

ComputeHomeTuples.h

Go to the documentation of this file.
00001 
00007 #ifndef COMPUTEHOMETUPLES_H
00008 #define COMPUTEHOMETUPLES_H
00009 
00010 #include "NamdTypes.h"
00011 #include "common.h"
00012 #include "structures.h"
00013 #include "Compute.h"
00014 #include "HomePatch.h"
00015 
00016 #include "Box.h"
00017 #include "OwnerBox.h"
00018 #include "UniqueSet.h"
00019 
00020 #include "Node.h"
00021 #include "SimParameters.h"
00022 #include "PatchMap.inl"
00023 #include "AtomMap.h"
00024 #include "ComputeHomeTuples.h"
00025 #include "PatchMgr.h"
00026 #include "HomePatchList.h"
00027 #include "Molecule.h"
00028 #include "Parameters.h"
00029 #include "ReductionMgr.h"
00030 #include "UniqueSet.h"
00031 #include "UniqueSetIter.h"
00032 #include "Priorities.h"
00033 
00034 class TuplePatchElem {
00035   public:
00036     PatchID patchID;
00037     Patch *p;
00038     Box<Patch,CompAtom> *positionBox;
00039     Box<Patch,CompAtom> *avgPositionBox;
00040     Box<Patch,Results> *forceBox;
00041     CompAtom *x;
00042     CompAtomExt *xExt;
00043     CompAtom *x_avg;
00044     Results *r;
00045     Force *f;
00046 
00047     int hash() const { return patchID; }
00048 
00049   TuplePatchElem(PatchID pid = -1) {
00050     patchID = pid;
00051     p = NULL;
00052     positionBox = NULL;
00053     avgPositionBox = NULL;
00054     forceBox = NULL;
00055     x = NULL;
00056     xExt = NULL;
00057     x_avg = NULL;
00058     r = NULL;
00059     f = NULL;
00060   }
00061 
00062   TuplePatchElem(Patch *p_param, ComputeID cid) {
00063     patchID = p_param->getPatchID();
00064     p = p_param;
00065     positionBox = p_param->registerPositionPickup(cid);
00066     avgPositionBox = p_param->registerAvgPositionPickup(cid);
00067     forceBox = p_param->registerForceDeposit(cid);
00068     x = NULL;
00069     xExt = NULL;
00070     x_avg = NULL;
00071     r = NULL;
00072     f = NULL;
00073   }
00074     
00075   ~TuplePatchElem() {};
00076 
00077   int operator==(const TuplePatchElem &elem) const {
00078     return (elem.patchID == patchID);
00079   }
00080 
00081   int operator<(const TuplePatchElem &elem) const {
00082     return (patchID < elem.patchID);
00083   }
00084 };
00085 
00086 typedef UniqueSet<TuplePatchElem> TuplePatchList;
00087 typedef UniqueSetIter<TuplePatchElem> TuplePatchListIter;
00088 
00089 class AtomMap;
00090 class ReductionMgr;
00091 
00092 template <class T, class S, class P> class ComputeHomeTuples : public Compute {
00093 
00094   protected:
00095   
00096     virtual void loadTuples(void) {
00097       int numTuples;
00098 
00099       #ifdef MEM_OPT_VERSION
00100       AtomSignature *allSigs;      
00101       #else
00102       int32 **tuplesByAtom;
00103       /* const (need to propagate const) */ S *tupleStructs;
00104       #endif
00105       
00106       const P *tupleValues;
00107       Node *node = Node::Object();
00108 
00109       #ifdef MEM_OPT_VERSION
00110       allSigs = node->molecule->atomSigPool;
00111       #else      
00112       T::getMoleculePointers(node->molecule,
00113                     &numTuples, &tuplesByAtom, &tupleStructs);      
00114       #endif
00115       
00116       T::getParameterPointers(node->parameters, &tupleValues);
00117 
00118       tupleList.resize(0);
00119 
00120       LocalID aid[T::size];
00121 
00122       const int lesOn = node->simParameters->lesOn;
00123       Real invLesFactor = lesOn ? 
00124                           1.0/node->simParameters->lesFactor :
00125                           1.0;
00126 
00127       // cycle through each patch and gather all tuples
00128       TuplePatchListIter ai(tuplePatchList);
00129     
00130       for ( ai = ai.begin(); ai != ai.end(); ai++ )
00131       {
00132         // CompAtom *atom = (*ai).x;
00133         Patch *patch = (*ai).p;
00134         int numAtoms = patch->getNumAtoms();
00135         CompAtomExt *atomExt = (*ai).xExt; //patch->getCompAtomExtInfo();
00136     
00137         // cycle through each atom in the patch and load up tuples
00138         for (int j=0; j < numAtoms; j++)
00139         {              
00140            /* cycle through each tuple */
00141            #ifdef MEM_OPT_VERSION
00142            AtomSignature *thisAtomSig = &allSigs[atomExt[j].sigId];
00143            TupleSignature *allTuples;
00144            T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
00145            for(int k=0; k<numTuples; k++) {
00146                T t(atomExt[j].id, &allTuples[k], tupleValues);
00147            #else
00148            /* get list of all tuples for the atom */
00149            int32 *curTuple = tuplesByAtom[atomExt[j].id];
00150            for( ; *curTuple != -1; ++curTuple) {             
00151              T t(&tupleStructs[*curTuple],tupleValues);
00152            #endif            
00153              register int i;
00154              aid[0] = atomMap->localID(t.atomID[0]);
00155              int homepatch = aid[0].pid;
00156              int samepatch = 1;
00157              int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
00158              for (i=1; i < T::size; i++) {
00159                  aid[i] = atomMap->localID(t.atomID[i]);
00160                  samepatch = samepatch && ( homepatch == aid[i].pid );
00161                  has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
00162              }
00163              if ( samepatch ) continue;
00164              t.scale = has_les ? invLesFactor : 1;
00165              for (i=1; i < T::size; i++) {
00166                  homepatch = patchMap->downstream(homepatch,aid[i].pid);
00167              }
00168              if ( homepatch != notUsed && isBasePatch[homepatch] ) {
00169                for (i=0; i < T::size; i++) {
00170                  TuplePatchElem *p;
00171                  t.p[i] = p = tuplePatchList.find(TuplePatchElem(aid[i].pid));
00172                  if ( ! p ) {
00173                #ifdef MEM_OPT_VERSION
00174                iout << iWARN << "Tuple with atoms ";
00175                #else
00176                    iout << iWARN << "Tuple " << *curTuple << " with atoms ";
00177                #endif
00178                    int erri;
00179                    for( erri = 0; erri < T::size; erri++ ) {
00180                      iout << t.atomID[erri] << "(" <<  aid[erri].pid << ") ";
00181                    }
00182                    iout << "missing patch " << aid[i].pid << "\n" << endi;
00183                    
00184                    NAMD_die("Patch needed for tuple is missing.\n");
00185                  }
00186                  t.localIndex[i] = aid[i].index;
00187                }
00188                tupleList.add(t);
00189              }
00190            }
00191         }
00192       }
00193     }
00194 
00195     int doLoadTuples;
00196   
00197   protected:
00198   
00199     ResizeArray<T> tupleList;
00200     TuplePatchList tuplePatchList;
00201   
00202     PatchMap *patchMap;
00203     AtomMap *atomMap;
00204     SubmitReduction *reduction;
00205     SubmitReduction *pressureProfileReduction;
00206     BigReal *pressureProfileData;
00207     int pressureProfileSlabs;
00208     char *isBasePatch;
00209   
00210     ComputeHomeTuples(ComputeID c) : Compute(c) {
00211       patchMap = PatchMap::Object();
00212       atomMap = AtomMap::Object();
00213       reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00214       
00215       SimParameters *params = Node::Object()->simParameters;
00216       if (params->pressureProfileOn) {
00217         pressureProfileSlabs = T::pressureProfileSlabs = 
00218           params->pressureProfileSlabs;
00219         int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
00220         pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00221           REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00222         int numAtomTypePairs = n*n;
00223         pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
00224       } else {
00225         pressureProfileReduction = NULL;
00226         pressureProfileData = NULL;
00227       }
00228       doLoadTuples = false;
00229       isBasePatch = 0;
00230     }
00231 
00232     ComputeHomeTuples(ComputeID c, PatchIDList pids) : Compute(c) {
00233       patchMap = PatchMap::Object();
00234       atomMap = AtomMap::Object();
00235       reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00236       SimParameters *params = Node::Object()->simParameters;
00237       if (params->pressureProfileOn) {
00238         pressureProfileSlabs = T::pressureProfileSlabs = 
00239           params->pressureProfileSlabs;
00240         int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
00241         pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00242           REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00243         int numAtomTypePairs = n*n;
00244         pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
00245       } else {
00246         pressureProfileReduction = NULL;
00247         pressureProfileData = NULL;
00248       }
00249       doLoadTuples = false;
00250       int nPatches = patchMap->numPatches();
00251       isBasePatch = new char[nPatches];
00252       int i;
00253       for (i=0; i<nPatches; ++i) { isBasePatch[i] = 0; }
00254       for (i=0; i<pids.size(); ++i) { isBasePatch[pids[i]] = 1; }
00255     }
00256 
00257   public:
00258   
00259     virtual ~ComputeHomeTuples() {
00260       delete reduction;
00261       delete [] isBasePatch;
00262       delete pressureProfileReduction;
00263       delete pressureProfileData;
00264     }
00265 
00266     //======================================================================
00267     // initialize() - Method is invoked only the first time
00268     // atom maps, patchmaps etc are ready and we are about to start computations
00269     //======================================================================
00270     virtual void initialize(void) {
00271     
00272       // Start with empty list
00273       tuplePatchList.clear();
00274     
00275       int nPatches = patchMap->numPatches();
00276       int pid;
00277       for (pid=0; pid<nPatches; ++pid) {
00278         if ( isBasePatch[pid] ) {
00279           Patch *patch = patchMap->patch(pid);
00280           tuplePatchList.add(TuplePatchElem(patch, cid));
00281         }
00282       }
00283     
00284       // Gather all proxy patches (neighbors, that is)
00285       PatchID neighbors[PatchMap::MaxOneOrTwoAway];
00286     
00287       for (pid=0; pid<nPatches; ++pid) if ( isBasePatch[pid] ) {
00288         int numNeighbors = patchMap->upstreamNeighbors(pid,neighbors);
00289         for ( int i = 0; i < numNeighbors; ++i ) {
00290           if ( ! tuplePatchList.find(TuplePatchElem(neighbors[i])) ) {
00291             Patch *patch = patchMap->patch(neighbors[i]);
00292             tuplePatchList.add(TuplePatchElem(patch, cid));
00293           }
00294         }
00295       }
00296       setNumPatches(tuplePatchList.size());
00297       doLoadTuples = true;
00298 
00299       basePriority = COMPUTE_PROXY_PRIORITY;  // no patch dependence
00300     }
00301 
00302     //======================================================================
00303     // atomUpdate() - Method is invoked after anytime that atoms have been
00304     // changed in patches used by this Compute object.
00305     //======================================================================
00306     void atomUpdate(void) {
00307       doLoadTuples = true;
00308     }
00309 
00310 //-------------------------------------------------------------------
00311 // Routine which is called by enqueued work msg.  It wraps
00312 // actualy Force computation with the apparatus needed
00313 // to get access to atom positions, return forces etc.
00314 //-------------------------------------------------------------------
00315     virtual void doWork(void) {
00316 
00317       // Open Boxes - register that we are using Positions
00318       // and will be depositing Forces.
00319       UniqueSetIter<TuplePatchElem> ap(tuplePatchList);
00320       for (ap = ap.begin(); ap != ap.end(); ap++) {
00321         ap->x = ap->positionBox->open();
00322         ap->xExt = ap->p->getCompAtomExtInfo();
00323         if ( ap->p->flags.doMolly ) ap->x_avg = ap->avgPositionBox->open();
00324         ap->r = ap->forceBox->open();
00325         ap->f = ap->r->f[Results::normal];
00326       } 
00327     
00328       BigReal reductionData[T::reductionDataSize];
00329       int tupleCount = 0;
00330       int numAtomTypes = T::pressureProfileAtomTypes;
00331       int numAtomTypePairs = numAtomTypes*numAtomTypes;
00332 
00333       if ( ! Node::Object()->simParameters->commOnly ) {
00334       if ( doLoadTuples ) {
00335         loadTuples();
00336         doLoadTuples = false;
00337       }
00338     
00339       for ( int i = 0; i < T::reductionDataSize; ++i ) reductionData[i] = 0;
00340       if (pressureProfileData) {
00341         memset(pressureProfileData, 0, 3*pressureProfileSlabs*numAtomTypePairs*sizeof(BigReal));
00342         // Silly variable hiding of the previous iterator
00343         UniqueSetIter<TuplePatchElem> newap(tuplePatchList);
00344         newap = newap.begin();
00345         const Lattice &lattice = newap->p->lattice;
00346         T::pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
00347         T::pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
00348       }
00349       // take triplet and pass with tuple info to force eval
00350       T *al = tupleList.begin();
00351       const int ntuple = tupleList.size();
00352       for (int i=0; i<ntuple; ++i) {
00353         al[i].computeForce(reductionData, pressureProfileData);
00354       }
00355       tupleCount += ntuple;
00356       }
00357  
00358       T::submitReductionData(reductionData,reduction);
00359       reduction->item(T::reductionChecksumLabel) += (BigReal)tupleCount;
00360       reduction->submit();
00361 
00362       if (pressureProfileReduction) {
00363         // For ease of calculation we stored interactions between types
00364         // i and j in (ni+j).  For efficiency now we coalesce the
00365         // cross interactions so that just i<=j are stored.
00366         const int arraysize = 3*pressureProfileSlabs;
00367         const BigReal *data = pressureProfileData;
00368         for (int i=0; i<numAtomTypes; i++) {
00369           for (int j=0; j<numAtomTypes; j++) {
00370             int ii=i;
00371             int jj=j;
00372             if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
00373             const int reductionOffset = 
00374               (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
00375             for (int k=0; k<arraysize; k++) {
00376               pressureProfileReduction->item(reductionOffset+k) += data[k];
00377             }
00378             data += arraysize;
00379           }
00380         }
00381         pressureProfileReduction->submit();
00382       }
00383     
00384       // Close boxes - i.e. signal we are done with Positions and
00385       // AtomProperties and that we are depositing Forces
00386       for (ap = ap.begin(); ap != ap.end(); ap++) {
00387         ap->positionBox->close(&(ap->x));
00388         if ( ap->p->flags.doMolly ) ap->avgPositionBox->close(&(ap->x_avg));
00389         ap->forceBox->close(&(ap->r));
00390       }
00391     }
00392 };
00393 
00394 
00395 #endif
00396 

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