Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

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 
00013 #define MIN_DEBUG_LEVEL 4
00014 // #define DEBUGM
00015 #include "Debug.h"
00016 
00017 ComputeNonbondedPair::ComputeNonbondedPair(ComputeID c, PatchID pid[], int trans[],
00018                 ComputeNonbondedWorkArrays* _workArrays,
00019                 int minPartition, int maxPartition, int numPartitions)
00020   : ComputePatchPair(c,pid,trans), workArrays(_workArrays),
00021     minPart(minPartition), maxPart(maxPartition), numParts(numPartitions)
00022 {
00023   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00024   if (pressureProfileOn) {
00025     int n = pressureProfileAtomTypes; 
00026     pressureProfileData = new BigReal[3*n*n*pressureProfileSlabs];
00027     pressureProfileReduction = ReductionMgr::Object()->willSubmit(
00028         REDUCTIONS_PPROF_NONBONDED, 3*pressureProfileSlabs*((n*(n+1))/2));
00029   } else {
00030     pressureProfileReduction = NULL;
00031     pressureProfileData = NULL;
00032   }
00033   pairlistsValid = 0;
00034   pairlistTolerance = 0.;
00035 }
00036 
00037 void ComputeNonbondedPair::initialize() {
00038   ComputePatchPair::initialize();
00039   for (int i=0; i<2; i++) {
00040     avgPositionBox[i] = patch[i]->registerAvgPositionPickup(cid,trans[i]);
00041   }
00042 }
00043 
00044 ComputeNonbondedPair::~ComputeNonbondedPair()
00045 {
00046   delete reduction;
00047   delete pressureProfileReduction;
00048   delete [] pressureProfileData;
00049   for (int i=0; i<2; i++) {
00050     if (avgPositionBox[i] != NULL) {
00051       patch[i]->unregisterAvgPositionPickup(cid,&avgPositionBox[i]);
00052     }
00053   }
00054 }
00055 
00056 int ComputeNonbondedPair::noWork() {
00057 
00058   // return 0;  // for testing
00059   if ( numAtoms[0] && numAtoms[1] && patch[0]->flags.doNonbonded )
00060   {
00061     return 0;  // work to do, enqueue as usual
00062   } else
00063   {
00064     // Inform load balancer
00065     LdbCoordinator::Object()->startWork(cid,0); // Timestep not used
00066     // fake out patches and reduction system
00067 
00068     BigReal reductionData[reductionDataSize];
00069     int i;
00070     for ( i = 0; i < reductionDataSize; ++i ) reductionData[i] = 0;
00071     if (pressureProfileOn) {
00072       int n = pressureProfileAtomTypes;
00073       memset(pressureProfileData, 0, 3*n*n*pressureProfileSlabs*sizeof(BigReal));
00074     }
00075     CompAtom* p[2];
00076     CompAtom* p_avg[2];
00077     Results* r[2];
00078 
00079     // Open up positionBox, forceBox, and atomBox
00080     for (i=0; i<2; i++) {
00081       p[i] = positionBox[i]->open();
00082       r[i] = forceBox[i]->open();
00083       if ( patch[0]->flags.doMolly ) p_avg[i] = avgPositionBox[i]->open();
00084     }
00085 
00086     // Close up boxes
00087     for (i=0; i<2; i++) {
00088       positionBox[i]->close(&p[i]);
00089       forceBox[i]->close(&r[i]);
00090       if ( patch[0]->flags.doMolly ) avgPositionBox[i]->close(&p_avg[i]);
00091     }
00092 
00093     submitReductionData(reductionData,reduction);
00094     if (pressureProfileOn)
00095       submitPressureProfileData(pressureProfileData, pressureProfileReduction);
00096 
00097     // Inform load balancer
00098     LdbCoordinator::Object()->endWork(cid,0); // Timestep not used
00099 
00100     reduction->submit();
00101     if (pressureProfileOn) 
00102       pressureProfileReduction->submit();
00103 
00104     return 1;  // no work to do, do not enqueue
00105   }
00106 }
00107 
00108 #ifdef MEM_OPT_VERSION
00109 void ComputeNonbondedPair::doForce(CompAtom* p[2], CompAtomExt* pExt[2], Results* r[2])
00110 #else
00111 void ComputeNonbondedPair::doForce(CompAtom* p[2], Results* r[2])
00112 #endif
00113 {
00114   // Inform load balancer. 
00115   // I assume no threads will suspend until endWork is called
00116   LdbCoordinator::Object()->startWork(cid,0); // Timestep not used
00117 
00118 #ifdef TRACE_COMPUTE_OBJECTS
00119     double traceObjStartTime = CmiWallTimer();
00120 #endif
00121 
00122   DebugM(2,"doForce() called.\n");
00123   DebugM(2, numAtoms[0] << " patch #1 atoms and " <<
00124         numAtoms[1] << " patch #2 atoms\n");
00125 
00126   BigReal reductionData[reductionDataSize];
00127   for ( int i = 0; i < reductionDataSize; ++i ) reductionData[i] = 0;
00128   if (pressureProfileOn) {
00129     int n = pressureProfileAtomTypes; 
00130     memset(pressureProfileData, 0, 3*n*n*pressureProfileSlabs*sizeof(BigReal));
00131     // adjust lattice dimensions to allow constant pressure
00132     const Lattice &lattice = patch[0]->lattice;
00133     pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
00134     pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
00135   }
00136   if ( numAtoms[0] && numAtoms[1] )
00137   {
00138     nonbonded params;
00139     params.reduction = reductionData;
00140     params.pressureProfileReduction = pressureProfileData;
00141 
00142     params.minPart = minPart;
00143     params.maxPart = maxPart;
00144     params.numParts = numParts;
00145 
00146     params.workArrays = workArrays;
00147 
00148     params.pairlists = &pairlists;
00149     params.savePairlists = 0;
00150     params.usePairlists = 0;
00151     if ( patch[0]->flags.savePairlists ) {
00152       params.savePairlists = 1;
00153       params.usePairlists = 1;
00154     } else if ( patch[0]->flags.usePairlists && patch[1]->flags.usePairlists ) {
00155       if ( ! pairlistsValid ||
00156            ( patch[0]->flags.maxAtomMovement +
00157              patch[1]->flags.maxAtomMovement > pairlistTolerance ) ) {
00158         reductionData[pairlistWarningIndex] += 1;
00159       } else {
00160         params.usePairlists = 1;
00161       }
00162     }
00163     if ( ! params.usePairlists ) {
00164       pairlistsValid = 0;
00165     }
00166     params.plcutoff = cutoff;
00167     params.groupplcutoff = cutoff +
00168         patch[0]->flags.maxGroupRadius + patch[1]->flags.maxGroupRadius;
00169     if ( params.savePairlists ) {
00170       pairlistsValid = 1;
00171       pairlistTolerance = patch[0]->flags.pairlistTolerance +
00172                           patch[1]->flags.pairlistTolerance;
00173       params.plcutoff += pairlistTolerance;
00174       params.groupplcutoff += pairlistTolerance;
00175     }
00176 
00177     // swap to place more atoms in inner loop (second patch)
00178     int a = 0;  int b = 1;
00179     if ( numAtoms[0] > numAtoms[1] ) { a = 1; b = 0; }
00180 
00181     const Lattice &lattice = patch[0]->lattice;
00182     params.offset = lattice.offset(trans[a]) - lattice.offset(trans[b]);
00183 
00184     int doEnergy = patch[0]->flags.doEnergy;
00185       params.p[0] = p[a];
00186       params.p[1] = p[b];
00187     #ifdef MEM_OPT_VERSION      
00188       params.pExt[0] = pExt[a]; 
00189       params.pExt[1] = pExt[b];
00190     #endif
00191       params.ff[0] = r[a]->f[Results::nbond];
00192       params.ff[1] = r[b]->f[Results::nbond];
00193       params.numAtoms[0] = numAtoms[a];
00194       params.numAtoms[1] = numAtoms[b];
00195 
00196       // DMK - Atom Separation (water vs. non-water)
00197       #if NAMD_SeparateWaters != 0
00198         params.numWaterAtoms[0] = numWaterAtoms[a];
00199         params.numWaterAtoms[1] = numWaterAtoms[b];
00200       #endif
00201 
00202       if ( patch[0]->flags.doFullElectrostatics )
00203       {
00204         params.fullf[0] = r[a]->f[Results::slow];
00205         params.fullf[1] = r[b]->f[Results::slow];
00206         if ( patch[0]->flags.doMolly ) {
00207           if ( doEnergy ) calcPairEnergy(&params);
00208           else calcPair(&params);
00209           CompAtom *p_avg[2];
00210           p_avg[0] = avgPositionBox[0]->open();
00211           p_avg[1] = avgPositionBox[1]->open();
00212           params.p[0] = p_avg[a];
00213           params.p[1] = p_avg[b];
00214           if ( doEnergy ) calcSlowPairEnergy(&params);
00215           else calcSlowPair(&params);
00216           avgPositionBox[0]->close(&p_avg[0]);
00217           avgPositionBox[1]->close(&p_avg[1]);
00218         } else if ( patch[0]->flags.maxForceMerged == Results::slow ) {
00219           if ( doEnergy ) calcMergePairEnergy(&params);
00220           else calcMergePair(&params);
00221         } else {
00222           if ( doEnergy ) calcFullPairEnergy(&params);
00223           else calcFullPair(&params);
00224         }
00225       }
00226       else
00227         if ( doEnergy ) calcPairEnergy(&params);
00228         else calcPair(&params);
00229   }
00230 
00231   submitReductionData(reductionData,reduction);
00232   if (pressureProfileOn)
00233     submitPressureProfileData(pressureProfileData, pressureProfileReduction);
00234 
00235 #ifdef TRACE_COMPUTE_OBJECTS
00236     traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+cid, traceObjStartTime, CmiWallTimer());
00237 #endif
00238 
00239   // Inform load balancer
00240   LdbCoordinator::Object()->endWork(cid,0); // Timestep not used
00241 
00242   reduction->submit();
00243   if (pressureProfileOn)
00244     pressureProfileReduction->submit();
00245 }
00246 

Generated on Fri Aug 29 04:07:39 2008 for NAMD by  doxygen 1.3.9.1