ComputeNonbondedSelf.C

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

Generated on Thu Nov 23 01:17:11 2017 for NAMD by  doxygen 1.4.7