ComputeNonbondedPair.C

Go to the documentation of this file.
00001 
00007 #include "ComputeNonbondedPair.h"
00008 #include "ReductionMgr.h"
00009 #include "Patch.h"
00010 #include "LdbCoordinator.h"
00011 #include "PatchMap.h"
00012 #include "ComputeMgr.h"
00013 #include "Molecule.h"
00014 
00015 #include "Node.h"
00016 #include "SimParameters.h"
00017 
00018 #define MIN_DEBUG_LEVEL 4
00019 // #define DEBUGM
00020 #include "Debug.h"
00021 
00022 ComputeNonbondedPair::ComputeNonbondedPair(ComputeID c, PatchID pid[], int trans[],
00023                 ComputeNonbondedWorkArrays* _workArrays,
00024                 int minPartition, int maxPartition, int numPartitions)
00025   : ComputePatchPair(c,pid,trans), workArrays(_workArrays),
00026     minPart(minPartition), maxPart(maxPartition), numParts(numPartitions)
00027 {
00028   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00029   if (pressureProfileOn) {
00030     int n = pressureProfileAtomTypes; 
00031     pressureProfileData = new BigReal[3*n*n*pressureProfileSlabs];
00032     pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00033         REDUCTIONS_PPROF_NONBONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00034   } else {
00035     pressureProfileReduction = NULL;
00036     pressureProfileData = NULL;
00037   }
00038   pairlistsValid = 0;
00039   pairlistTolerance = 0.;
00040   params.simParameters = Node::Object()->simParameters;
00041   params.parameters = Node::Object()->parameters;
00042   params.random = Node::Object()->rand;
00043 }
00044 
00045 void ComputeNonbondedPair::initialize() {
00046   ComputePatchPair::initialize();
00047   for (int i=0; i<2; i++) {
00048     avgPositionBox[i] = patch[i]->registerAvgPositionPickup(this);
00049     // BEGIN LA
00050     velocityBox[i] = patch[i]->registerVelocityPickup(this);
00051     // END LA
00052 
00053     psiSumBox[i] = patch[i]->registerPsiSumDeposit(this);
00054     intRadBox[i] = patch[i]->registerIntRadPickup(this);
00055     bornRadBox[i] = patch[i]->registerBornRadPickup(this);
00056     dEdaSumBox[i] = patch[i]->registerDEdaSumDeposit(this);
00057     dHdrPrefixBox[i] = patch[i]->registerDHdrPrefixPickup(this);
00058   }
00059 #ifdef NAMD_CUDA
00060   register_cuda_compute_pair(cid, patchID, trans);
00061 #endif
00062 }
00063 
00064 ComputeNonbondedPair::~ComputeNonbondedPair()
00065 {
00066   delete reduction;
00067   delete pressureProfileReduction;
00068   delete [] pressureProfileData;
00069   for (int i=0; i<2; i++) {
00070     if (avgPositionBox[i] != NULL) {
00071       patch[i]->unregisterAvgPositionPickup(this,&avgPositionBox[i]);
00072     }
00073     // BEGIN LA
00074     if (velocityBox[i] != NULL) {
00075       patch[i]->unregisterVelocityPickup(this,&velocityBox[i]);
00076     }
00077     // END LA
00078 
00079     if (psiSumBox[i] != NULL) {
00080       patch[i]->unregisterPsiSumDeposit(this,&psiSumBox[i]);
00081     }
00082     if (intRadBox[i] != NULL) {
00083       patch[i]->unregisterIntRadPickup(this,&intRadBox[i]);
00084     }
00085     if (bornRadBox[i] != NULL) {
00086       patch[i]->unregisterBornRadPickup(this,&bornRadBox[i]);
00087     }
00088     if (dEdaSumBox[i] != NULL) {
00089       patch[i]->unregisterDEdaSumDeposit(this,&dEdaSumBox[i]);
00090     }
00091     if (dHdrPrefixBox[i] != NULL) {
00092       patch[i]->unregisterDHdrPrefixPickup(this,&dHdrPrefixBox[i]);
00093     }
00094   }
00095 }
00096 
00097 int ComputeNonbondedPair::noWork() {
00098 
00099   if (patch[0]->flags.doGBIS) {
00100     gbisPhase = 1 + (gbisPhase % 3);//1->2->3->1...
00101   }
00102 
00103 #ifndef NAMD_CUDA
00104   if ( patch[0]->flags.doNonbonded && (numAtoms[0] && numAtoms[1]) ) {
00105     return 0;  // work to do, enqueue as usual
00106   } else {
00107 #else
00108   {
00109 #endif
00110 
00111     if (patch[0]->flags.doGBIS) {
00112      if (gbisPhase == 1) {
00113       for (int i=0; i<2; i++) {
00114         psiSumBox[i]->skip();
00115         intRadBox[i]->skip();
00116       }
00117       if (patch[0]->flags.doNonbonded) return 1;
00118       else gbisPhase = 2;
00119      }
00120      if (gbisPhase == 2) {
00121       for (int i=0; i<2; i++) {
00122         bornRadBox[i]->skip();
00123         dEdaSumBox[i]->skip();
00124       }
00125       if (patch[0]->flags.doNonbonded) return 1;
00126       else gbisPhase = 3;
00127      }
00128      if (gbisPhase == 3) {
00129       for (int i=0; i<2; i++) {
00130         dHdrPrefixBox[i]->skip();
00131       }
00132      }
00133     }
00134 
00135     // skip all boxes
00136     for (int i=0; i<2; i++) {
00137       positionBox[i]->skip();
00138       forceBox[i]->skip();
00139       if ( patch[0]->flags.doMolly ) avgPositionBox[i]->skip();
00140       // BEGIN LA
00141       if (patch[0]->flags.doLoweAndersen) velocityBox[i]->skip();
00142       // END LA
00143     }
00144 
00145     reduction->item(REDUCTION_COMPUTE_CHECKSUM) += 1.;
00146     reduction->submit();
00147     if (pressureProfileOn) 
00148       pressureProfileReduction->submit();
00149 
00150 #ifndef NAMD_CUDA
00151     // Inform load balancer
00152     LdbCoordinator::Object()->skipWork(ldObjHandle);
00153 #endif
00154 
00155     return 1;  // no work to do, do not enqueue
00156   }
00157 }
00158 
00159 void ComputeNonbondedPair::doForce(CompAtom* p[2], CompAtomExt* pExt[2], Results* r[2])
00160 {
00161   // Inform load balancer. 
00162   // I assume no threads will suspend until endWork is called
00163 
00164   //single phase declarations
00165   int doEnergy = patch[0]->flags.doEnergy;
00166   int a = 0;  int b = 1;
00167   // swap to place more atoms in inner loop (second patch)
00168   if ( numAtoms[0] > numAtoms[1] ) { a = 1; b = 0; }
00169   CompAtom* v[2];
00170 
00171 
00172 /*******************************************************************************
00173  * Prepare Parameters
00174 *******************************************************************************/
00175   if (!patch[0]->flags.doGBIS || gbisPhase == 1) {
00176 
00177 #ifdef TRACE_COMPUTE_OBJECTS
00178     double traceObjStartTime = CmiWallTimer();
00179 #endif
00180 
00181   DebugM(2,"doForce() called.\n");
00182   DebugM(2, numAtoms[0] << " patch #1 atoms and " <<
00183         numAtoms[1] << " patch #2 atoms\n");
00184 
00185 
00186   for ( int i = 0; i < reductionDataSize; ++i )
00187     reductionData[i] = 0;
00188   if (pressureProfileOn) {
00189     int n = pressureProfileAtomTypes; 
00190     memset(pressureProfileData, 0, 3*n*n*pressureProfileSlabs*sizeof(BigReal));
00191     // adjust lattice dimensions to allow constant pressure
00192     const Lattice &lattice = patch[0]->lattice;
00193     pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
00194     pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
00195   }
00196 
00197     params.reduction = reductionData;
00198     params.pressureProfileReduction = pressureProfileData;
00199 
00200     params.minPart = minPart;
00201     params.maxPart = maxPart;
00202     params.numParts = numParts;
00203 
00204     params.workArrays = workArrays;
00205 
00206     params.pairlists = &pairlists;
00207     params.savePairlists = 0;
00208     params.usePairlists = 0;
00209     if ( patch[0]->flags.savePairlists ) {
00210       params.savePairlists = 1;
00211       params.usePairlists = 1;
00212     } else if ( patch[0]->flags.usePairlists && patch[1]->flags.usePairlists ) {
00213       if ( ! pairlistsValid ||
00214            ( patch[0]->flags.maxAtomMovement +
00215              patch[1]->flags.maxAtomMovement > pairlistTolerance ) ) {
00216         reductionData[pairlistWarningIndex] += 1;
00217       } else {
00218         params.usePairlists = 1;
00219       }
00220     }
00221     if ( ! params.usePairlists ) {
00222       pairlistsValid = 0;
00223     }
00224     params.plcutoff = cutoff;
00225     params.groupplcutoff = cutoff +
00226         patch[0]->flags.maxGroupRadius + patch[1]->flags.maxGroupRadius;
00227     if ( params.savePairlists ) {
00228       pairlistsValid = 1;
00229       pairlistTolerance = patch[0]->flags.pairlistTolerance +
00230                           patch[1]->flags.pairlistTolerance;
00231       params.plcutoff += pairlistTolerance;
00232       params.groupplcutoff += pairlistTolerance;
00233     }
00234 
00235 
00236     const Lattice &lattice = patch[0]->lattice;
00237     params.offset = lattice.offset(trans[a]) - lattice.offset(trans[b]);
00238 
00239     PatchMap* patchMap = PatchMap::Object();
00240     params.offset_f = params.offset + lattice.unscale(patchMap->center(patchID[a]))
00241                                     - lattice.unscale(patchMap->center(patchID[b]));
00242 
00243     // Atom Sorting : If we are sorting the atoms along the line connecting
00244     //   the patch centers, then calculate a normalized vector pointing from
00245     //   patch a to patch b (i.e. outer loop patch to inner loop patch).
00246     #if NAMD_ComputeNonbonded_SortAtoms != 0
00247 
00248       params.projLineVec = params.offset_f * ( -1. / params.offset_f.length() );
00249 
00250     #endif
00251 
00252       params.p[0] = p[a];
00253       params.p[1] = p[b];
00254       params.pExt[0] = pExt[a]; 
00255       params.pExt[1] = pExt[b];
00256 #ifdef NAMD_KNL
00257       params.pFlt[0] = patch[a]->getCompAtomFlt();
00258       params.pFlt[1] = patch[b]->getCompAtomFlt();
00259 #endif
00260       // BEGIN LA
00261       params.doLoweAndersen = patch[0]->flags.doLoweAndersen;
00262       if (params.doLoweAndersen) {
00263           DebugM(4, "opening velocity boxes\n");
00264           v[0] = velocityBox[0]->open();
00265           v[1] = velocityBox[1]->open();
00266           params.v[0] = v[a];
00267           params.v[1] = v[b];
00268       }
00269       // END LA
00270 #ifndef NAMD_CUDA
00271       params.ff[0] = r[a]->f[Results::nbond_virial];
00272       params.ff[1] = r[b]->f[Results::nbond_virial];
00273 #endif
00274       params.numAtoms[0] = numAtoms[a];
00275       params.numAtoms[1] = numAtoms[b];
00276       params.step = patch[0]->flags.step;
00277 
00278       // DMK - Atom Separation (water vs. non-water)
00279       #if NAMD_SeparateWaters != 0
00280         params.numWaterAtoms[0] = numWaterAtoms[a];
00281         params.numWaterAtoms[1] = numWaterAtoms[b];
00282       #endif
00283 
00284 
00285 /*******************************************************************************
00286  * Call Nonbonded Functions
00287 *******************************************************************************/
00288       if (numAtoms[0] && numAtoms[1]) {//only do if has atoms since gbis noWork doesn't account for no atoms
00289 
00290       //force calculation calls
00291       if ( patch[0]->flags.doFullElectrostatics )
00292       {
00293 #ifndef NAMD_CUDA
00294         params.fullf[0] = r[a]->f[Results::slow_virial];
00295         params.fullf[1] = r[b]->f[Results::slow_virial];
00296 #endif
00297         if ( patch[0]->flags.doMolly ) {
00298           if ( doEnergy )
00299             calcPairEnergy(&params);
00300           else calcPair(&params);
00301           CompAtom *p_avg[2];
00302           p_avg[0] = avgPositionBox[0]->open();
00303           p_avg[1] = avgPositionBox[1]->open();
00304           params.p[0] = p_avg[a];
00305           params.p[1] = p_avg[b];
00306           if ( doEnergy ) calcSlowPairEnergy(&params);
00307           else calcSlowPair(&params);
00308           avgPositionBox[0]->close(&p_avg[0]);
00309           avgPositionBox[1]->close(&p_avg[1]);
00310         } else if ( patch[0]->flags.maxForceMerged == Results::slow ) {
00311           if ( doEnergy ) calcMergePairEnergy(&params);
00312     else calcMergePair(&params);
00313   } else {
00314     if ( doEnergy ) calcFullPairEnergy(&params);
00315     else calcFullPair(&params);
00316   }
00317       }
00318       else
00319         if ( doEnergy ) calcPairEnergy(&params);
00320         else calcPair(&params);
00321 
00322       }//end if has atoms
00323       
00324       // BEGIN LA
00325       if (params.doLoweAndersen) {
00326           DebugM(4, "closing velocity boxes\n");
00327           velocityBox[0]->close(&v[0]);
00328           velocityBox[1]->close(&v[1]);
00329       }
00330       // END LA
00331     }// end not gbis
00332 
00333 /*******************************************************************************
00334  * gbis Loop
00335 *******************************************************************************/
00336 if (patch[0]->flags.doGBIS) {
00337   SimParameters *simParams = Node::Object()->simParameters;
00338   gbisParams.sequence = sequence();
00339   gbisParams.doGBIS = patch[0]->flags.doGBIS;
00340   gbisParams.numPatches = 2;//pair
00341   gbisParams.gbisPhase = gbisPhase;
00342   gbisParams.doFullElectrostatics = patch[0]->flags.doFullElectrostatics;
00343   gbisParams.epsilon_s = simParams->solvent_dielectric;
00344   gbisParams.epsilon_p = simParams->dielectric;
00345   gbisParams.rho_0 = simParams->coulomb_radius_offset;
00346   gbisParams.kappa = simParams->kappa;
00347   gbisParams.cutoff = simParams->cutoff;
00348   gbisParams.doSmoothing = simParams->switchingActive;
00349   gbisParams.a_cut = simParams->alpha_cutoff;
00350   gbisParams.delta = simParams->gbis_delta;
00351   gbisParams.beta = simParams->gbis_beta;
00352   gbisParams.gamma = simParams->gbis_gamma;
00353   gbisParams.alpha_max = simParams->alpha_max;
00354   gbisParams.cid = cid;
00355   gbisParams.patchID[0] = patch[a]->getPatchID();
00356   gbisParams.patchID[1] = patch[b]->getPatchID();
00357   gbisParams.maxGroupRadius = patch[0]->flags.maxGroupRadius;
00358   if (patch[1]->flags.maxGroupRadius > gbisParams.maxGroupRadius)
00359     gbisParams.maxGroupRadius = patch[1]->flags.maxGroupRadius;
00360   gbisParams.doEnergy = doEnergy;
00361   gbisParams.fsMax = simParams->fsMax;
00362   for (int i = 0; i < numGBISPairlists; i++)
00363     gbisParams.gbisStepPairlists[i] = &gbisStepPairlists[i];
00364 
00365   //open boxes
00366   if (gbisPhase == 1) {
00367       gbisParams.intRad[0] = intRadBox[a]->open();
00368       gbisParams.intRad[1] = intRadBox[b]->open();
00369       gbisParams.psiSum[0] = psiSumBox[a]->open();
00370       gbisParams.psiSum[1] = psiSumBox[b]->open();
00371       gbisParams.gbInterEnergy=0;
00372       gbisParams.gbSelfEnergy=0;
00373 
00374   } else if (gbisPhase == 2) {
00375       gbisParams.bornRad[0] = bornRadBox[a]->open();
00376       gbisParams.bornRad[1] = bornRadBox[b]->open();
00377       gbisParams.dEdaSum[0] = dEdaSumBox[a]->open();
00378       gbisParams.dEdaSum[1] = dEdaSumBox[b]->open();
00379   } else if (gbisPhase == 3) {
00380       gbisParams.dHdrPrefix[0] = dHdrPrefixBox[a]->open();
00381       gbisParams.dHdrPrefix[1] = dHdrPrefixBox[b]->open();
00382   }
00383 
00384   //make call to calculate GBIS
00385   if ( !ComputeNonbondedUtil::commOnly ) {
00386     calcGBIS(&params,&gbisParams);
00387   }
00388 
00389   //close boxes
00390   if (gbisPhase == 1) {
00391       psiSumBox[0]->close(&(gbisParams.psiSum[a]));
00392       psiSumBox[1]->close(&(gbisParams.psiSum[b]));
00393   } else if (gbisPhase == 2) {
00394       dEdaSumBox[0]->close(&(gbisParams.dEdaSum[a]));
00395       dEdaSumBox[1]->close(&(gbisParams.dEdaSum[b]));
00396 
00397 
00398   } else if (gbisPhase == 3) {
00399       bornRadBox[0]->close(&(gbisParams.bornRad[a]));
00400       bornRadBox[1]->close(&(gbisParams.bornRad[b]));
00401       reduction->item(REDUCTION_ELECT_ENERGY) += gbisParams.gbInterEnergy;
00402       reduction->item(REDUCTION_ELECT_ENERGY) += gbisParams.gbSelfEnergy;
00403       intRadBox[0]->close(&(gbisParams.intRad[a]));
00404       intRadBox[1]->close(&(gbisParams.intRad[b]));
00405       dHdrPrefixBox[0]->close(&(gbisParams.dHdrPrefix[a]));
00406       dHdrPrefixBox[1]->close(&(gbisParams.dHdrPrefix[b]));
00407   }
00408 
00409 }//end if doGBIS
00410 
00411 
00412   if (!patch[0]->flags.doGBIS || gbisPhase == 3) {
00413   submitReductionData(reductionData,reduction);
00414   if (pressureProfileOn)
00415     submitPressureProfileData(pressureProfileData, pressureProfileReduction);
00416 
00417 #ifdef TRACE_COMPUTE_OBJECTS
00418     traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+cid, traceObjStartTime, CmiWallTimer());
00419 #endif
00420 
00421 
00422   reduction->submit();
00423   if (pressureProfileOn)
00424     pressureProfileReduction->submit();
00425   }//end gbis end phase
00426 
00427 }//end do Force
00428 

Generated on Fri Sep 22 01:17:11 2017 for NAMD by  doxygen 1.4.7