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

Generated on Sat May 25 04:07:15 2013 for NAMD by  doxygen 1.3.9.1