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 #ifdef MEM_OPT_VERSION
00043 CompAtomExt *xExt;
00044 #endif
00045 CompAtom *x_avg;
00046 Results *r;
00047 Force *f;
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 #ifdef MEM_OPT_VERSION
00059 xExt = NULL;
00060 #endif
00061 x_avg = NULL;
00062 r = NULL;
00063 f = NULL;
00064 }
00065
00066 TuplePatchElem(Patch *p_param, ComputeID cid) {
00067 patchID = p_param->getPatchID();
00068 p = p_param;
00069 positionBox = p_param->registerPositionPickup(cid);
00070 avgPositionBox = p_param->registerAvgPositionPickup(cid);
00071 forceBox = p_param->registerForceDeposit(cid);
00072 x = NULL;
00073 #ifdef MEM_OPT_VERSION
00074 xExt = NULL;
00075 #endif
00076 x_avg = NULL;
00077 r = NULL;
00078 f = NULL;
00079 }
00080
00081 ~TuplePatchElem() {};
00082
00083 int operator==(const TuplePatchElem &elem) const {
00084 return (elem.patchID == patchID);
00085 }
00086
00087 int operator<(const TuplePatchElem &elem) const {
00088 return (patchID < elem.patchID);
00089 }
00090 };
00091
00092 typedef UniqueSet<TuplePatchElem> TuplePatchList;
00093 typedef UniqueSetIter<TuplePatchElem> TuplePatchListIter;
00094
00095 class AtomMap;
00096 class ReductionMgr;
00097
00098 template <class T, class S, class P> class ComputeHomeTuples : public Compute {
00099
00100 protected:
00101
00102 virtual void loadTuples(void) {
00103 int numTuples;
00104
00105 #ifdef MEM_OPT_VERSION
00106 AtomSignature *allSigs;
00107 #else
00108 int32 **tuplesByAtom;
00109 S *tupleStructs;
00110 #endif
00111
00112 const P *tupleValues;
00113 Node *node = Node::Object();
00114
00115 #ifdef MEM_OPT_VERSION
00116 allSigs = node->molecule->atomSigPool;
00117 #else
00118 T::getMoleculePointers(node->molecule,
00119 &numTuples, &tuplesByAtom, &tupleStructs);
00120 #endif
00121
00122 T::getParameterPointers(node->parameters, &tupleValues);
00123
00124 tupleList.resize(0);
00125
00126 LocalID aid[T::size];
00127
00128 const int lesOn = node->simParameters->lesOn;
00129 Real invLesFactor = lesOn ?
00130 1.0/node->simParameters->lesFactor :
00131 1.0;
00132
00133
00134 TuplePatchListIter ai(tuplePatchList);
00135
00136 for ( ai = ai.begin(); ai != ai.end(); ai++ )
00137 {
00138 CompAtom *atom = (*ai).x;
00139 Patch *patch = (*ai).p;
00140 int numAtoms = patch->getNumAtoms();
00141 #ifdef MEM_OPT_VERSION
00142 CompAtomExt *atomExt = (*ai).xExt;
00143 #endif
00144
00145
00146 for (int j=0; j < numAtoms; j++)
00147 {
00148
00149 #ifdef MEM_OPT_VERSION
00150 AtomSignature *thisAtomSig = &allSigs[atomExt[j].sigId];
00151 TupleSignature *allTuples;
00152 T::getTupleInfo(thisAtomSig, &numTuples, &allTuples);
00153 for(int k=0; k<numTuples; k++) {
00154 T t(atom[j].id, &allTuples[k], tupleValues);
00155 #else
00156
00157 int32 *curTuple = tuplesByAtom[atom[j].id];
00158 for( ; *curTuple != -1; ++curTuple) {
00159 T t(&tupleStructs[*curTuple],tupleValues);
00160 #endif
00161 register int i;
00162 aid[0] = atomMap->localID(t.atomID[0]);
00163 int homepatch = aid[0].pid;
00164 int samepatch = 1;
00165 int has_les = lesOn && node->molecule->get_fep_type(t.atomID[0]);
00166 for (i=1; i < T::size; i++) {
00167 aid[i] = atomMap->localID(t.atomID[i]);
00168 samepatch = samepatch && ( homepatch == aid[i].pid );
00169 has_les |= lesOn && node->molecule->get_fep_type(t.atomID[i]);
00170 }
00171 if ( samepatch ) continue;
00172 t.scale = has_les ? invLesFactor : 1;
00173 for (i=1; i < T::size; i++) {
00174 homepatch = patchMap->downstream(homepatch,aid[i].pid);
00175 }
00176 if ( homepatch != notUsed && isBasePatch[homepatch] ) {
00177 for (i=0; i < T::size; i++) {
00178 TuplePatchElem *p;
00179 t.p[i] = p = tuplePatchList.find(TuplePatchElem(aid[i].pid));
00180 if ( ! p ) {
00181 #ifdef MEM_OPT_VERSION
00182 iout << iWARN << "Tuple with atoms ";
00183 #else
00184 iout << iWARN << "Tuple " << *curTuple << " with atoms ";
00185 #endif
00186 int erri;
00187 for( erri = 0; erri < T::size; erri++ ) {
00188 iout << t.atomID[erri] << "(" << aid[erri].pid << ") ";
00189 }
00190 iout << "missing patch " << aid[i].pid << "\n" << endi;
00191
00192 NAMD_die("Patch needed for tuple is missing.\n");
00193 }
00194 t.localIndex[i] = aid[i].index;
00195 }
00196 tupleList.add(t);
00197 }
00198 }
00199 }
00200 }
00201 }
00202
00203 int doLoadTuples;
00204
00205 protected:
00206
00207 ResizeArray<T> tupleList;
00208 TuplePatchList tuplePatchList;
00209
00210 PatchMap *patchMap;
00211 AtomMap *atomMap;
00212 SubmitReduction *reduction;
00213 SubmitReduction *pressureProfileReduction;
00214 BigReal *pressureProfileData;
00215 int pressureProfileSlabs;
00216 char *isBasePatch;
00217
00218 ComputeHomeTuples(ComputeID c) : Compute(c) {
00219 patchMap = PatchMap::Object();
00220 atomMap = AtomMap::Object();
00221 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00222
00223 SimParameters *params = Node::Object()->simParameters;
00224 if (params->pressureProfileOn) {
00225 pressureProfileSlabs = T::pressureProfileSlabs =
00226 params->pressureProfileSlabs;
00227 int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
00228 pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00229 REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00230 int numAtomTypePairs = n*n;
00231 pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
00232 } else {
00233 pressureProfileReduction = NULL;
00234 pressureProfileData = NULL;
00235 }
00236 doLoadTuples = false;
00237 isBasePatch = 0;
00238 }
00239
00240 ComputeHomeTuples(ComputeID c, PatchIDList pids) : Compute(c) {
00241 patchMap = PatchMap::Object();
00242 atomMap = AtomMap::Object();
00243 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00244 SimParameters *params = Node::Object()->simParameters;
00245 if (params->pressureProfileOn) {
00246 pressureProfileSlabs = T::pressureProfileSlabs =
00247 params->pressureProfileSlabs;
00248 int n = T::pressureProfileAtomTypes = params->pressureProfileAtomTypes;
00249 pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00250 REDUCTIONS_PPROF_BONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00251 int numAtomTypePairs = n*n;
00252 pressureProfileData = new BigReal[3*pressureProfileSlabs*numAtomTypePairs];
00253 } else {
00254 pressureProfileReduction = NULL;
00255 pressureProfileData = NULL;
00256 }
00257 doLoadTuples = false;
00258 int nPatches = patchMap->numPatches();
00259 isBasePatch = new char[nPatches];
00260 int i;
00261 for (i=0; i<nPatches; ++i) { isBasePatch[i] = 0; }
00262 for (i=0; i<pids.size(); ++i) { isBasePatch[pids[i]] = 1; }
00263 }
00264
00265 public:
00266
00267 virtual ~ComputeHomeTuples() {
00268 delete reduction;
00269 delete [] isBasePatch;
00270 delete pressureProfileReduction;
00271 delete pressureProfileData;
00272 }
00273
00274
00275
00276
00277
00278 virtual void initialize(void) {
00279
00280
00281 tuplePatchList.clear();
00282
00283 int nPatches = patchMap->numPatches();
00284 int pid;
00285 for (pid=0; pid<nPatches; ++pid) {
00286 if ( isBasePatch[pid] ) {
00287 Patch *patch = patchMap->patch(pid);
00288 tuplePatchList.add(TuplePatchElem(patch, cid));
00289 }
00290 }
00291
00292
00293 PatchID neighbors[PatchMap::MaxOneOrTwoAway];
00294
00295 for (pid=0; pid<nPatches; ++pid) if ( isBasePatch[pid] ) {
00296 int numNeighbors = patchMap->upstreamNeighbors(pid,neighbors);
00297 for ( int i = 0; i < numNeighbors; ++i ) {
00298 if ( ! tuplePatchList.find(TuplePatchElem(neighbors[i])) ) {
00299 Patch *patch = patchMap->patch(neighbors[i]);
00300 tuplePatchList.add(TuplePatchElem(patch, cid));
00301 }
00302 }
00303 }
00304 setNumPatches(tuplePatchList.size());
00305 doLoadTuples = true;
00306
00307 basePriority = COMPUTE_PROXY_PRIORITY;
00308 }
00309
00310
00311
00312
00313
00314 void atomUpdate(void) {
00315 doLoadTuples = true;
00316 }
00317
00318
00319
00320
00321
00322
00323 virtual void doWork(void) {
00324
00325
00326
00327 UniqueSetIter<TuplePatchElem> ap(tuplePatchList);
00328 for (ap = ap.begin(); ap != ap.end(); ap++) {
00329 ap->x = ap->positionBox->open();
00330 #ifdef MEM_OPT_VERSION
00331 ap->xExt = ap->p->getCompAtomExtInfo();
00332 #endif
00333 if ( ap->p->flags.doMolly ) ap->x_avg = ap->avgPositionBox->open();
00334 ap->r = ap->forceBox->open();
00335 ap->f = ap->r->f[Results::normal];
00336 }
00337
00338 BigReal reductionData[T::reductionDataSize];
00339 int tupleCount = 0;
00340 int numAtomTypes = T::pressureProfileAtomTypes;
00341 int numAtomTypePairs = numAtomTypes*numAtomTypes;
00342
00343 if ( ! Node::Object()->simParameters->commOnly ) {
00344 if ( doLoadTuples ) {
00345 loadTuples();
00346 doLoadTuples = false;
00347 }
00348
00349 for ( int i = 0; i < T::reductionDataSize; ++i ) reductionData[i] = 0;
00350 if (pressureProfileData) {
00351 memset(pressureProfileData, 0, 3*pressureProfileSlabs*numAtomTypePairs*sizeof(BigReal));
00352
00353 UniqueSetIter<TuplePatchElem> newap(tuplePatchList);
00354 newap = newap.begin();
00355 const Lattice &lattice = newap->p->lattice;
00356 T::pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
00357 T::pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
00358 }
00359
00360 T *al = tupleList.begin();
00361 const int ntuple = tupleList.size();
00362 for (int i=0; i<ntuple; ++i) {
00363 al[i].computeForce(reductionData, pressureProfileData);
00364 }
00365 tupleCount += ntuple;
00366 }
00367
00368 T::submitReductionData(reductionData,reduction);
00369 reduction->item(T::reductionChecksumLabel) += (BigReal)tupleCount;
00370 reduction->submit();
00371
00372 if (pressureProfileReduction) {
00373
00374
00375
00376 const int arraysize = 3*pressureProfileSlabs;
00377 const BigReal *data = pressureProfileData;
00378 for (int i=0; i<numAtomTypes; i++) {
00379 for (int j=0; j<numAtomTypes; j++) {
00380 int ii=i;
00381 int jj=j;
00382 if (ii > jj) { int tmp=ii; ii=jj; jj=tmp; }
00383 const int reductionOffset =
00384 (ii*numAtomTypes - (ii*(ii+1))/2 + jj)*arraysize;
00385 for (int k=0; k<arraysize; k++) {
00386 pressureProfileReduction->item(reductionOffset+k) += data[k];
00387 }
00388 data += arraysize;
00389 }
00390 }
00391 pressureProfileReduction->submit();
00392 }
00393
00394
00395
00396 for (ap = ap.begin(); ap != ap.end(); ap++) {
00397 ap->positionBox->close(&(ap->x));
00398 if ( ap->p->flags.doMolly ) ap->avgPositionBox->close(&(ap->x_avg));
00399 ap->forceBox->close(&(ap->r));
00400 }
00401 }
00402 };
00403
00404
00405 #endif
00406