ComputeHomeTuples.h

Go to the documentation of this file.
00001 
00007 #ifndef COMPUTEHOMETUPLES_H
00008 #define COMPUTEHOMETUPLES_H
00009 
00010 #ifdef USE_HOMETUPLES
00011 #include <vector>
00012 #endif
00013 #include "NamdTypes.h"
00014 #include "common.h"
00015 #include "structures.h"
00016 #include "Compute.h"
00017 #include "HomePatch.h"
00018 
00019 #include "Box.h"
00020 #include "OwnerBox.h"
00021 #include "UniqueSet.h"
00022 
00023 #include "Node.h"
00024 #include "SimParameters.h"
00025 #include "PatchMap.inl"
00026 #include "AtomMap.h"
00027 #include "ComputeHomeTuples.h"
00028 #include "PatchMgr.h"
00029 #include "ProxyMgr.h"
00030 #include "HomePatchList.h"
00031 #include "Molecule.h"
00032 #include "Parameters.h"
00033 #include "ReductionMgr.h"
00034 #include "UniqueSet.h"
00035 #include "UniqueSetIter.h"
00036 #include "Priorities.h"
00037 #include "LdbCoordinator.h"
00038 
00039 class TuplePatchElem {
00040   public:
00041     PatchID patchID;
00042     Patch *p;
00043     Box<Patch,CompAtom> *positionBox;
00044     Box<Patch,CompAtom> *avgPositionBox;
00045     Box<Patch,Results> *forceBox;
00046     CompAtom *x;
00047     CompAtomExt *xExt;
00048     CompAtom *x_avg;
00049     Results *r;
00050     Force *f;
00051     Force *af;
00052 
00053     int hash() const { return patchID; }
00054 
00055   TuplePatchElem(PatchID pid = -1) {
00056     patchID = pid;
00057     p = NULL;
00058     positionBox = NULL;
00059     avgPositionBox = NULL;
00060     forceBox = NULL;
00061     x = NULL;
00062     xExt = NULL;
00063     x_avg = NULL;
00064     r = NULL;
00065     f = NULL;
00066     af = NULL;
00067   }
00068 
00069   TuplePatchElem(Patch *p_param, Compute *cid) {
00070     patchID = p_param->getPatchID();
00071     p = p_param;
00072     positionBox = p_param->registerPositionPickup(cid);
00073     avgPositionBox = p_param->registerAvgPositionPickup(cid);
00074     forceBox = p_param->registerForceDeposit(cid);
00075     x = NULL;
00076     xExt = NULL;
00077     x_avg = NULL;
00078     r = NULL;
00079     f = NULL;
00080     af = NULL;
00081   }
00082     
00083   ~TuplePatchElem() {};
00084 
00085   int operator==(const TuplePatchElem &elem) const {
00086     return (elem.patchID == patchID);
00087   }
00088 
00089   int operator<(const TuplePatchElem &elem) const {
00090     return (patchID < elem.patchID);
00091   }
00092 };
00093 
00094 typedef UniqueSet<TuplePatchElem> TuplePatchList;
00095 typedef UniqueSetIter<TuplePatchElem> TuplePatchListIter;
00096 
00097 class AtomMap;
00098 class ReductionMgr;
00099 
00100 #ifdef MEM_OPT_VERSION
00101 template <class T> struct ElemTraits {
00102   typedef AtomSignature signature;
00103   static signature* get_sig_pointer(Molecule *mol) { return mol->atomSigPool; }
00104   static int get_sig_id(const CompAtomExt &a) { return a.sigId; }
00105 };
00106 
00107 template <> struct ElemTraits <ExclElem> {
00108   typedef ExclusionSignature signature;
00109   static signature* get_sig_pointer(Molecule *mol) { return mol->exclSigPool; }
00110   static int get_sig_id(const CompAtomExt &a) { return a.exclId; }
00111 };
00112 #endif
00113 
00114 #ifdef USE_HOMETUPLES
00115 //
00116 // Simple base class for HomeTuples and SelfTuples that stores the type of the tuple
00117 //
00118 class Tuples {
00119 private:
00120   int type;
00121 protected:
00122   Tuples(int type) : type(type) {}
00123 public:
00124   // Tuple types
00125   enum {BOND=0, ANGLE, DIHEDRAL, IMPROPER, EXCLUSION, CROSSTERM, NUM_TUPLE_TYPES};
00126 
00127   int getType() {return type;}
00128   virtual void submitTupleCount(SubmitReduction *reduction, int tupleCount)=0;
00129   // virtual void copyTupleData(void* tupleData)=0;
00130   virtual int getNumTuples()=0;
00131   virtual void* getTupleList()=0;
00132   virtual void loadTuples(TuplePatchList& tuplePatchList, const char* isBasePatch, AtomMap *atomMap,
00133       const std::vector<int>& pids = std::vector<int>())=0;
00134 };
00135 
00136 //
00137 // HomeTuples class. These are created and stored in ComputeBondedCUDA::registerCompute()
00138 // e.g.: new HomeTuples<BondElem, Bond, BondValue>(BOND)
00139 //
00140 template <class T, class S, class P> class HomeTuples : public Tuples {
00141   protected:
00142     std::vector<T> tupleList;
00143 
00144   public:
00145 
00146     HomeTuples(int type=-1) : Tuples(type) {}
00147 
00148 #if __cplusplus < 201103L
00149 #define final
00150 #endif
00151 
00152     virtual void* getTupleList() final {
00153       return (void*)tupleList.data();
00154     }
00155 
00156     virtual void submitTupleCount(SubmitReduction *reduction, int tupleCount) final {
00157       reduction->item(T::reductionChecksumLabel) += (BigReal)tupleCount;
00158     }
00159 
00160     // virtual void copyTupleData(void* tupleData) final {
00161       // for (int i=0;i < tupleList.size();i++) {
00162       //   tupleData[i] = 
00163       // }
00164       // T::loadTupleData(tupleData);
00165     // }
00166 
00167     virtual int getNumTuples() final {
00168       return tupleList.size();
00169     }
00170 
00171     virtual void loadTuples(TuplePatchList& tuplePatchList, const char* isBasePatch, AtomMap *atomMap,
00172       const std::vector<int>& pids = std::vector<int>()) {
00173 
00174       if (isBasePatch == NULL) {
00175         NAMD_bug("NULL isBasePatch detected in HomeTuples::loadTuples()");
00176       }
00177 
00178       int numTuples;
00179 
00180 #ifdef MEM_OPT_VERSION
00181       typename ElemTraits<T>::signature *allSigs;      
00182 #else
00183       int32 **tuplesByAtom;
00184       /* const (need to propagate const) */ S *tupleStructs;
00185 #endif
00186       
00187       const P *tupleValues;
00188       Node *node = Node::Object();
00189       PatchMap *patchMap = PatchMap::Object();
00190       // AtomMap *atomMap = AtomMap::Object();
00191 
00192 #ifdef MEM_OPT_VERSION
00193       allSigs = ElemTraits<T>::get_sig_pointer(node->molecule);
00194 #else      
00195       T::getMoleculePointers(node->molecule,
00196         &numTuples, &tuplesByAtom, &tupleStructs);      
00197 #endif
00198       
00199       T::getParameterPointers(node->parameters, &tupleValues);
00200 
00201       tupleList.clear();
00202 
00203       LocalID aid[T::size];
00204 
00205       const int lesOn = node->simParameters->lesOn;
00206       const int soluteScalingOn = node->simParameters->soluteScalingOn;
00207       Real invLesFactor = lesOn ? 
00208                           1.0/node->simParameters->lesFactor :
00209                           1.0;
00210       const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
00211       const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
00212 
00213       // cycle through each patch and gather all tuples
00214       TuplePatchListIter ai(tuplePatchList);
00215       if (pids.size() == 0) ai = ai.begin();
00216 
00217       int numPid = (pids.size() == 0) ? tuplePatchList.size() : pids.size();
00218 
00219       for (int ipid=0;ipid < numPid;ipid++) {
00220         // Patch *patch;
00221         int numAtoms;
00222         CompAtomExt *atomExt;
00223         // Take next patch
00224         if (pids.size() == 0) {
00225           Patch* patch = (*ai).p;
00226           numAtoms = patch->getNumAtoms();
00227           atomExt = (*ai).xExt;
00228           ai++;
00229         } else {
00230           TuplePatchElem *tpe = tuplePatchList.find(TuplePatchElem(pids[ipid]));
00231           Patch* patch = tpe->p;
00232           numAtoms = patch->getNumAtoms();
00233           atomExt = tpe->xExt;          
00234         }
00235    
00236         // cycle through each atom in the patch and load up tuples
00237         for (int j=0; j < numAtoms; j++)
00238         {
00239           /* cycle through each tuple */
00240 #ifdef MEM_OPT_VERSION
00241           typename ElemTraits<T>::signature *thisAtomSig =
00242                    &allSigs[ElemTraits<T>::get_sig_id(atomExt[j])];
00243           TupleSignature *allTuples;
00244           T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
00245           for(int k=0; k<numTuples; k++) {
00246             T t(atomExt[j].id, &allTuples[k], tupleValues);
00247 #else
00248           /* get list of all tuples for the atom */
00249           int32 *curTuple = tuplesByAtom[atomExt[j].id];
00250           for( ; *curTuple != -1; ++curTuple) {
00251             T t(&tupleStructs[*curTuple],tupleValues);
00252 #endif            
00253             register int i;
00254             aid[0] = atomMap->localID(t.atomID[0]);
00255             int homepatch = aid[0].pid;
00256             int samepatch = 1;
00257             int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
00258             int has_ss = soluteScalingOn && node->molecule->get_ss_type(t.atomID[0]);
00259             for (i=1; i < T::size; i++) {
00260               aid[i] = atomMap->localID(t.atomID[i]);
00261               samepatch = samepatch && ( homepatch == aid[i].pid );
00262               has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
00263               has_ss |= soluteScalingOn && node->molecule->get_ss_type(t.atomID[i]);
00264             }
00265             if (T::size < 4 && !soluteScalingAll) has_ss = false;
00266             if ( samepatch ) continue;
00267             t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
00268             for (i=1; i < T::size; i++) {
00269               homepatch = patchMap->downstream(homepatch,aid[i].pid);
00270             }
00271             if ( homepatch != notUsed && isBasePatch[homepatch] ) {
00272               TuplePatchElem *p;
00273               for (i=0; i < T::size; i++) {
00274                 t.p[i] = p = tuplePatchList.find(TuplePatchElem(aid[i].pid));
00275                 if ( ! p ) {
00276 #ifdef MEM_OPT_VERSION
00277                   iout << iWARN << "Tuple with atoms ";
00278 #else
00279                   iout << iWARN << "Tuple " << *curTuple << " with atoms ";
00280 #endif
00281                   int erri;
00282                   for( erri = 0; erri < T::size; erri++ ) {
00283                     iout << t.atomID[erri] << "(" <<  aid[erri].pid << ") ";
00284                   }
00285                   iout << "missing patch " << aid[i].pid << "\n" << endi;
00286                   break;
00287                 }
00288                 t.localIndex[i] = aid[i].index;
00289               }
00290               if ( ! p ) continue;
00291 #ifdef MEM_OPT_VERSION
00292               //avoid adding Tuples whose atoms are all fixed
00293               if(node->simParameters->fixedAtomsOn && !node->simParameters->fixedAtomsForces) {
00294                 int allfixed = 1;
00295                 for(i=0; i<T::size; i++){
00296                   CompAtomExt *one = &(t.p[i]->xExt[aid[i].index]);
00297                   allfixed = allfixed & one->atomFixed;
00298                 }
00299                 if(!allfixed) tupleList.push_back(t);
00300               }else{
00301                 tupleList.push_back(t);
00302               }
00303 #else
00304               tupleList.push_back(t);
00305 #endif               
00306             }
00307           }
00308         }
00309       }
00310     }
00311 
00312 };
00313 #endif
00314 
00315 template <class T, class S, class P> class ComputeHomeTuples : public Compute {
00316 
00317   protected:
00318   
00319 #ifndef USE_HOMETUPLES
00320     virtual void loadTuples(void) {
00321       int numTuples;
00322 
00323       #ifdef MEM_OPT_VERSION
00324       typename ElemTraits<T>::signature *allSigs;      
00325       #else
00326       int32 **tuplesByAtom;
00327       /* const (need to propagate const) */ S *tupleStructs;
00328       #endif
00329       
00330       const P *tupleValues;
00331       Node *node = Node::Object();
00332 
00333       #ifdef MEM_OPT_VERSION
00334       allSigs = ElemTraits<T>::get_sig_pointer(node->molecule);
00335       #else      
00336       T::getMoleculePointers(node->molecule,
00337                     &numTuples, &tuplesByAtom, &tupleStructs);      
00338       #endif
00339       
00340       T::getParameterPointers(node->parameters, &tupleValues);
00341 
00342       tupleList.resize(0);
00343 
00344       LocalID aid[T::size];
00345 
00346       const int lesOn = node->simParameters->lesOn;
00347       const int soluteScalingOn = node->simParameters->soluteScalingOn;
00348       Real invLesFactor = lesOn ? 
00349                           1.0/node->simParameters->lesFactor :
00350                           1.0;
00351       const Real soluteScalingFactor = node->simParameters->soluteScalingFactor;
00352       const Bool soluteScalingAll = node->simParameters->soluteScalingAll;
00353 
00354       // cycle through each patch and gather all tuples
00355       TuplePatchListIter ai(tuplePatchList);
00356     
00357       for ( ai = ai.begin(); ai != ai.end(); ai++ )
00358       {
00359         // CompAtom *atom = (*ai).x;
00360         Patch *patch = (*ai).p;
00361         int numAtoms = patch->getNumAtoms();
00362         CompAtomExt *atomExt = (*ai).xExt; //patch->getCompAtomExtInfo();
00363     
00364         // cycle through each atom in the patch and load up tuples
00365         for (int j=0; j < numAtoms; j++)
00366         {              
00367            /* cycle through each tuple */
00368            #ifdef MEM_OPT_VERSION
00369            typename ElemTraits<T>::signature *thisAtomSig =
00370                    &allSigs[ElemTraits<T>::get_sig_id(atomExt[j])];
00371            TupleSignature *allTuples;
00372            T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
00373            for(int k=0; k<numTuples; k++) {
00374                T t(atomExt[j].id, &allTuples[k], tupleValues);
00375            #else
00376            /* get list of all tuples for the atom */
00377            int32 *curTuple = tuplesByAtom[atomExt[j].id];
00378            for( ; *curTuple != -1; ++curTuple) {             
00379              T t(&tupleStructs[*curTuple],tupleValues);
00380            #endif            
00381              register int i;
00382              aid[0] = atomMap->localID(t.atomID[0]);
00383              int homepatch = aid[0].pid;
00384              int samepatch = 1;
00385              int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
00386              int has_ss = soluteScalingOn && node->molecule->get_ss_type(t.atomID[0]);
00387              for (i=1; i < T::size; i++) {
00388                  aid[i] = atomMap->localID(t.atomID[i]);
00389                  samepatch = samepatch && ( homepatch == aid[i].pid );
00390                  has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
00391                  has_ss |= soluteScalingOn && node->molecule->get_ss_type(t.atomID[i]);
00392              }
00393              if (T::size < 4 && !soluteScalingAll) has_ss = false;
00394              if ( samepatch ) continue;
00395              t.scale = (!has_les && !has_ss) ? 1.0 : ( has_les ? invLesFactor : soluteScalingFactor );
00396              for (i=1; i < T::size; i++) {
00397                  homepatch = patchMap->downstream(homepatch,aid[i].pid);
00398              }
00399              if ( homepatch != notUsed && isBasePatch[homepatch] ) {
00400                TuplePatchElem *p;
00401                for (i=0; i < T::size; i++) {
00402                  t.p[i] = p = tuplePatchList.find(TuplePatchElem(aid[i].pid));
00403                  if ( ! p ) {
00404                      #ifdef MEM_OPT_VERSION
00405                      iout << iWARN << "Tuple with atoms ";
00406                      #else
00407                    iout << iWARN << "Tuple " << *curTuple << " with atoms ";
00408                      #endif
00409                    int erri;
00410                    for( erri = 0; erri < T::size; erri++ ) {
00411                      iout << t.atomID[erri] << "(" <<  aid[erri].pid << ") ";
00412                    }
00413                    iout << "missing patch " << aid[i].pid << "\n" << endi;
00414                    break;
00415                  }
00416                  t.localIndex[i] = aid[i].index;
00417                }
00418                if ( ! p ) continue;
00419              #ifdef MEM_OPT_VERSION
00420                //avoid adding Tuples whose atoms are all fixed
00421                if(node->simParameters->fixedAtomsOn &&
00422                   !node->simParameters->fixedAtomsForces) {
00423                  int allfixed = 1;
00424                  for(i=0; i<T::size; i++){
00425                    CompAtomExt *one = &(t.p[i]->xExt[aid[i].index]);
00426                    allfixed = allfixed & one->atomFixed;
00427                  }
00428                  if(!allfixed) tupleList.add(t);
00429                }else{
00430                  tupleList.add(t);
00431                }
00432              #else
00433                tupleList.add(t);
00434              #endif               
00435              }
00436            }
00437         }
00438       }
00439     }
00440 #endif
00441 
00442     int doLoadTuples;
00443   
00444   protected:
00445   
00446 #ifdef USE_HOMETUPLES
00447     HomeTuples<T, S, P>* tuples;
00448     TuplePatchList tuplePatchList;
00449 #else
00450     ResizeArray<T> tupleList;
00451     TuplePatchList tuplePatchList;
00452 #endif
00453 
00454     PatchMap *patchMap;
00455     AtomMap *atomMap;
00456     SubmitReduction *reduction;
00457     int accelMDdoDihe;
00458     SubmitReduction *pressureProfileReduction;
00459     BigReal *pressureProfileData;
00460     int pressureProfileSlabs;
00461     char *isBasePatch;
00462   
00463     ComputeHomeTuples(ComputeID c) : Compute(c) {
00464       patchMap = PatchMap::Object();
00465       atomMap = AtomMap::Object();
00466       reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00467       
00468       SimParameters *params = Node::Object()->simParameters;
00469       accelMDdoDihe=false;
00470       if (params->accelMDOn) {
00471          if (params->accelMDdihe || params->accelMDdual) accelMDdoDihe=true;
00472       }
00473       if (params->pressureProfileOn) {
00474         pressureProfileSlabs = T::pressureProfileSlabs = 
00475           params->pressureProfileSlabs;
00476         int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
00477         pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00478           REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00479         int numAtomTypePairs = n*n;
00480         pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
00481       } else {
00482         pressureProfileReduction = NULL;
00483         pressureProfileData = NULL;
00484       }
00485       doLoadTuples = false;
00486       isBasePatch = 0;
00487 #ifdef USE_HOMETUPLES
00488       tuples = NULL;
00489 #endif
00490     }
00491 
00492     ComputeHomeTuples(ComputeID c, PatchIDList &pids) : Compute(c) {
00493       patchMap = PatchMap::Object();
00494       atomMap = AtomMap::Object();
00495       reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00496       SimParameters *params = Node::Object()->simParameters;
00497       accelMDdoDihe=false;
00498       if (params->accelMDOn) {
00499          if (params->accelMDdihe || params->accelMDdual) accelMDdoDihe=true;
00500       }
00501       if (params->pressureProfileOn) {
00502         pressureProfileSlabs = T::pressureProfileSlabs = 
00503           params->pressureProfileSlabs;
00504         int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
00505         pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00506           REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00507         int numAtomTypePairs = n*n;
00508         pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
00509       } else {
00510         pressureProfileReduction = NULL;
00511         pressureProfileData = NULL;
00512       }
00513       doLoadTuples = false;
00514       int nPatches = patchMap->numPatches();
00515       isBasePatch = new char[nPatches];
00516       int i;
00517       for (i=0; i<nPatches; ++i) { isBasePatch[i] = 0; }
00518       for (i=0; i<pids.size(); ++i) { isBasePatch[pids[i]] = 1; }
00519 #ifdef USE_HOMETUPLES
00520       tuples = NULL;
00521 #endif
00522     }
00523 
00524   public:
00525   
00526     virtual ~ComputeHomeTuples() {
00527       delete reduction;
00528       delete [] isBasePatch;
00529       delete pressureProfileReduction;
00530       delete pressureProfileData;
00531 #ifdef USE_HOMETUPLES
00532       if (tuples != NULL) delete tuples;
00533 #endif
00534     }
00535 
00536     //======================================================================
00537     // initialize() - Method is invoked only the first time
00538     // atom maps, patchmaps etc are ready and we are about to start computations
00539     //======================================================================
00540     virtual void initialize(void) {
00541 
00542 #ifdef NAMD_CUDA
00543       ProxyMgr *proxyMgr = ProxyMgr::Object();
00544 #endif
00545 
00546 #ifdef USE_HOMETUPLES
00547       tuples = new HomeTuples<T, S, P>();
00548 #endif
00549 
00550       // Start with empty list
00551       tuplePatchList.clear();
00552     
00553       int nPatches = patchMap->numPatches();
00554       int pid;
00555       for (pid=0; pid<nPatches; ++pid) {
00556         if ( isBasePatch[pid] ) {
00557 #ifdef NAMD_CUDA
00558           proxyMgr->createProxy(pid);
00559 #endif
00560           Patch *patch = patchMap->patch(pid);
00561           tuplePatchList.add(TuplePatchElem(patch, this));
00562         }
00563       }
00564     
00565       // Gather all proxy patches (neighbors, that is)
00566       PatchID neighbors[PatchMap::MaxOneOrTwoAway];
00567     
00568       for (pid=0; pid<nPatches; ++pid) if ( isBasePatch[pid] ) {
00569         int numNeighbors = patchMap->upstreamNeighbors(pid,neighbors);
00570         for ( int i = 0; i < numNeighbors; ++i ) {
00571           if ( ! tuplePatchList.find(TuplePatchElem(neighbors[i])) ) {
00572 #ifdef NAMD_CUDA
00573             proxyMgr->createProxy(neighbors[i]);
00574 #endif
00575             Patch *patch = patchMap->patch(neighbors[i]);
00576             tuplePatchList.add(TuplePatchElem(patch, this));
00577           }
00578         }
00579       }
00580       setNumPatches(tuplePatchList.size());
00581       doLoadTuples = true;
00582 
00583       basePriority = COMPUTE_PROXY_PRIORITY;  // no patch dependence
00584     }
00585 
00586     //======================================================================
00587     // atomUpdate() - Method is invoked after anytime that atoms have been
00588     // changed in patches used by this Compute object.
00589     //======================================================================
00590     void atomUpdate(void) {
00591       doLoadTuples = true;
00592     }
00593 
00594 //-------------------------------------------------------------------
00595 // Routine which is called by enqueued work msg.  It wraps
00596 // actualy Force computation with the apparatus needed
00597 // to get access to atom positions, return forces etc.
00598 //-------------------------------------------------------------------
00599     virtual void doWork(void) {
00600 
00601       LdbCoordinator::Object()->startWork(ldObjHandle);
00602 
00603       // Open Boxes - register that we are using Positions
00604       // and will be depositing Forces.
00605       UniqueSetIter<TuplePatchElem> ap(tuplePatchList);
00606       for (ap = ap.begin(); ap != ap.end(); ap++) {
00607         ap->x = ap->positionBox->open();
00608         ap->xExt = ap->p->getCompAtomExtInfo();
00609         if ( ap->p->flags.doMolly ) ap->x_avg = ap->avgPositionBox->open();
00610         ap->r = ap->forceBox->open();
00611         ap->f = ap->r->f[Results::normal];
00612         if (accelMDdoDihe) ap->af = ap->r->f[Results::amdf]; // for dihedral-only or dual-boost accelMD
00613       } 
00614     
00615       BigReal reductionData[T::reductionDataSize];
00616       int tupleCount = 0;
00617       int numAtomTypes = T::pressureProfileAtomTypes;
00618       int numAtomTypePairs = numAtomTypes*numAtomTypes;
00619     
00620       for ( int i = 0; i < T::reductionDataSize; ++i ) reductionData[i] = 0;
00621       if (pressureProfileData) {
00622         memset(pressureProfileData, 0, 3*pressureProfileSlabs*numAtomTypePairs*sizeof(BigReal));
00623         // Silly variable hiding of the previous iterator
00624         UniqueSetIter<TuplePatchElem> newap(tuplePatchList);
00625         newap = newap.begin();
00626         const Lattice &lattice = newap->p->lattice;
00627         T::pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
00628         T::pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
00629       }
00630 
00631       if ( ! Node::Object()->simParameters->commOnly ) {
00632       if ( doLoadTuples ) {
00633 #ifdef USE_HOMETUPLES
00634         tuples->loadTuples(tuplePatchList, isBasePatch, AtomMap::Object());
00635 #else
00636         loadTuples();
00637 #endif
00638         doLoadTuples = false;
00639       }
00640       // take triplet and pass with tuple info to force eval
00641 #ifdef USE_HOMETUPLES
00642       T *al = (T *)tuples->getTupleList();
00643       const int ntuple = tuples->getNumTuples();
00644 #else
00645       T *al = tupleList.begin();
00646       const int ntuple = tupleList.size();
00647 #endif
00648       if ( ntuple ) T::computeForce(al, ntuple, reductionData, pressureProfileData);
00649       tupleCount += ntuple;
00650       }
00651  
00652     LdbCoordinator::Object()->endWork(ldObjHandle);
00653 
00654       T::submitReductionData(reductionData,reduction);
00655       reduction->item(T::reductionChecksumLabel) += (BigReal)tupleCount;
00656       reduction->submit();
00657 
00658       if (pressureProfileReduction) {
00659         // For ease of calculation we stored interactions between types
00660         // i and j in (ni+j).  For efficiency now we coalesce the
00661         // cross interactions so that just i<=j are stored.
00662         const int arraysize = 3*pressureProfileSlabs;
00663         const BigReal *data = pressureProfileData;
00664         for (int i=0; i<numAtomTypes; i++) {
00665           for (int j=0; j<numAtomTypes; j++) {
00666             int ii=i;
00667             int jj=j;
00668             if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
00669             const int reductionOffset = 
00670               (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
00671             for (int k=0; k<arraysize; k++) {
00672               pressureProfileReduction->item(reductionOffset+k) += data[k];
00673             }
00674             data += arraysize;
00675           }
00676         }
00677         pressureProfileReduction->submit();
00678       }
00679     
00680       // Close boxes - i.e. signal we are done with Positions and
00681       // AtomProperties and that we are depositing Forces
00682       for (ap = ap.begin(); ap != ap.end(); ap++) {
00683         ap->positionBox->close(&(ap->x));
00684         if ( ap->p->flags.doMolly ) ap->avgPositionBox->close(&(ap->x_avg));
00685         ap->forceBox->close(&(ap->r));
00686       }
00687     }
00688 };
00689 
00690 
00691 #endif
00692 

Generated on Wed Jun 20 01:17:11 2018 for NAMD by  doxygen 1.4.7