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
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 #ifdef NAMD_CUDA
00043 register_cuda_compute_pair(cid, patchID, trans);
00044 #endif
00045 }
00046
00047 ComputeNonbondedPair::~ComputeNonbondedPair()
00048 {
00049 delete reduction;
00050 delete pressureProfileReduction;
00051 delete [] pressureProfileData;
00052 for (int i=0; i<2; i++) {
00053 if (avgPositionBox[i] != NULL) {
00054 patch[i]->unregisterAvgPositionPickup(cid,&avgPositionBox[i]);
00055 }
00056 }
00057 }
00058
00059 int ComputeNonbondedPair::noWork() {
00060
00061
00062 if ( numAtoms[0] && numAtoms[1] && patch[0]->flags.doNonbonded
00063 #ifdef NAMD_CUDA
00064 && patch[0]->flags.doEnergy
00065 #endif
00066 )
00067 {
00068 return 0;
00069 } else
00070 {
00071
00072 #ifndef NAMD_CUDA
00073 LdbCoordinator::Object()->startWork(cid,0);
00074 #endif
00075
00076
00077 BigReal reductionData[reductionDataSize];
00078 int i;
00079 for ( i = 0; i < reductionDataSize; ++i ) reductionData[i] = 0;
00080 if (pressureProfileOn) {
00081 int n = pressureProfileAtomTypes;
00082 memset(pressureProfileData, 0, 3*n*n*pressureProfileSlabs*sizeof(BigReal));
00083 }
00084 CompAtom* p[2];
00085 CompAtom* p_avg[2];
00086 Results* r[2];
00087
00088
00089 for (i=0; i<2; i++) {
00090 p[i] = positionBox[i]->open();
00091 r[i] = forceBox[i]->open();
00092 if ( patch[0]->flags.doMolly ) p_avg[i] = avgPositionBox[i]->open();
00093 }
00094
00095
00096 for (i=0; i<2; i++) {
00097 positionBox[i]->close(&p[i]);
00098 forceBox[i]->close(&r[i]);
00099 if ( patch[0]->flags.doMolly ) avgPositionBox[i]->close(&p_avg[i]);
00100 }
00101
00102 submitReductionData(reductionData,reduction);
00103 if (pressureProfileOn)
00104 submitPressureProfileData(pressureProfileData, pressureProfileReduction);
00105
00106
00107 #ifndef NAMD_CUDA
00108 LdbCoordinator::Object()->endWork(cid,0);
00109 #endif
00110
00111 reduction->submit();
00112 if (pressureProfileOn)
00113 pressureProfileReduction->submit();
00114
00115 return 1;
00116 }
00117 }
00118
00119 void ComputeNonbondedPair::doForce(CompAtom* p[2], CompAtomExt* pExt[2], Results* r[2])
00120 {
00121
00122
00123 #ifndef NAMD_CUDA
00124 LdbCoordinator::Object()->startWork(cid,0);
00125 #endif
00126
00127 #ifdef TRACE_COMPUTE_OBJECTS
00128 double traceObjStartTime = CmiWallTimer();
00129 #endif
00130
00131 DebugM(2,"doForce() called.\n");
00132 DebugM(2, numAtoms[0] << " patch #1 atoms and " <<
00133 numAtoms[1] << " patch #2 atoms\n");
00134
00135 BigReal reductionData[reductionDataSize];
00136 for ( int i = 0; i < reductionDataSize; ++i ) reductionData[i] = 0;
00137 if (pressureProfileOn) {
00138 int n = pressureProfileAtomTypes;
00139 memset(pressureProfileData, 0, 3*n*n*pressureProfileSlabs*sizeof(BigReal));
00140
00141 const Lattice &lattice = patch[0]->lattice;
00142 pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
00143 pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
00144 }
00145 if ( numAtoms[0] && numAtoms[1] )
00146 {
00147 nonbonded params;
00148 params.reduction = reductionData;
00149 params.pressureProfileReduction = pressureProfileData;
00150
00151 params.minPart = minPart;
00152 params.maxPart = maxPart;
00153 params.numParts = numParts;
00154
00155 params.workArrays = workArrays;
00156
00157 params.pairlists = &pairlists;
00158 params.savePairlists = 0;
00159 params.usePairlists = 0;
00160 if ( patch[0]->flags.savePairlists ) {
00161 params.savePairlists = 1;
00162 params.usePairlists = 1;
00163 } else if ( patch[0]->flags.usePairlists && patch[1]->flags.usePairlists ) {
00164 if ( ! pairlistsValid ||
00165 ( patch[0]->flags.maxAtomMovement +
00166 patch[1]->flags.maxAtomMovement > pairlistTolerance ) ) {
00167 reductionData[pairlistWarningIndex] += 1;
00168 } else {
00169 params.usePairlists = 1;
00170 }
00171 }
00172 if ( ! params.usePairlists ) {
00173 pairlistsValid = 0;
00174 }
00175 params.plcutoff = cutoff;
00176 params.groupplcutoff = cutoff +
00177 patch[0]->flags.maxGroupRadius + patch[1]->flags.maxGroupRadius;
00178 if ( params.savePairlists ) {
00179 pairlistsValid = 1;
00180 pairlistTolerance = patch[0]->flags.pairlistTolerance +
00181 patch[1]->flags.pairlistTolerance;
00182 params.plcutoff += pairlistTolerance;
00183 params.groupplcutoff += pairlistTolerance;
00184 }
00185
00186
00187 int a = 0; int b = 1;
00188 if ( numAtoms[0] > numAtoms[1] ) { a = 1; b = 0; }
00189
00190 const Lattice &lattice = patch[0]->lattice;
00191 params.offset = lattice.offset(trans[a]) - lattice.offset(trans[b]);
00192
00193
00194
00195
00196 #if NAMD_ComputeNonbonded_SortAtoms != 0
00197
00198
00199 PatchMap* patchMap = PatchMap::Object();
00200 ScaledPosition p_a_center, p_b_center;
00201 p_a_center.x = patchMap->min_a(patchID[a]);
00202 p_a_center.y = patchMap->min_b(patchID[a]);
00203 p_a_center.z = patchMap->min_c(patchID[a]);
00204 p_b_center.x = patchMap->min_a(patchID[b]);
00205 p_b_center.y = patchMap->min_b(patchID[b]);
00206 p_b_center.z = patchMap->min_c(patchID[b]);
00207 p_a_center.x += patchMap->max_a(patchID[a]);
00208 p_a_center.y += patchMap->max_b(patchID[a]);
00209 p_a_center.z += patchMap->max_c(patchID[a]);
00210 p_b_center.x += patchMap->max_a(patchID[b]);
00211 p_b_center.y += patchMap->max_b(patchID[b]);
00212 p_b_center.z += patchMap->max_c(patchID[b]);
00213 p_a_center *= (BigReal)0.5;
00214 p_b_center *= (BigReal)0.5;
00215 p_a_center = lattice.unscale(p_a_center);
00216 p_b_center = lattice.unscale(p_b_center);
00217
00218
00219 p_a_center.x += params.offset.x;
00220 p_a_center.y += params.offset.y;
00221 p_a_center.z += params.offset.z;
00222
00223
00224 params.projLineVec = p_b_center - p_a_center;
00225 params.projLineVec /= params.projLineVec.length();
00226
00227 #endif
00228
00229 int doEnergy = patch[0]->flags.doEnergy;
00230 params.p[0] = p[a];
00231 params.p[1] = p[b];
00232 params.pExt[0] = pExt[a];
00233 params.pExt[1] = pExt[b];
00234 params.ff[0] = r[a]->f[Results::nbond];
00235 params.ff[1] = r[b]->f[Results::nbond];
00236 params.numAtoms[0] = numAtoms[a];
00237 params.numAtoms[1] = numAtoms[b];
00238
00239
00240 #if NAMD_SeparateWaters != 0
00241 params.numWaterAtoms[0] = numWaterAtoms[a];
00242 params.numWaterAtoms[1] = numWaterAtoms[b];
00243 #endif
00244
00245 if ( patch[0]->flags.doFullElectrostatics )
00246 {
00247 params.fullf[0] = r[a]->f[Results::slow];
00248 params.fullf[1] = r[b]->f[Results::slow];
00249 if ( patch[0]->flags.doMolly ) {
00250 if ( doEnergy ) calcPairEnergy(¶ms);
00251 else calcPair(¶ms);
00252 CompAtom *p_avg[2];
00253 p_avg[0] = avgPositionBox[0]->open();
00254 p_avg[1] = avgPositionBox[1]->open();
00255 params.p[0] = p_avg[a];
00256 params.p[1] = p_avg[b];
00257 if ( doEnergy ) calcSlowPairEnergy(¶ms);
00258 else calcSlowPair(¶ms);
00259 avgPositionBox[0]->close(&p_avg[0]);
00260 avgPositionBox[1]->close(&p_avg[1]);
00261 } else if ( patch[0]->flags.maxForceMerged == Results::slow ) {
00262 if ( doEnergy ) calcMergePairEnergy(¶ms);
00263 else calcMergePair(¶ms);
00264 } else {
00265 if ( doEnergy ) calcFullPairEnergy(¶ms);
00266 else calcFullPair(¶ms);
00267 }
00268 }
00269 else
00270 if ( doEnergy ) calcPairEnergy(¶ms);
00271 else calcPair(¶ms);
00272 }
00273
00274 submitReductionData(reductionData,reduction);
00275 if (pressureProfileOn)
00276 submitPressureProfileData(pressureProfileData, pressureProfileReduction);
00277
00278 #ifdef TRACE_COMPUTE_OBJECTS
00279 traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+cid, traceObjStartTime, CmiWallTimer());
00280 #endif
00281
00282
00283 #ifndef NAMD_CUDA
00284 LdbCoordinator::Object()->endWork(cid,0);
00285 #endif
00286
00287 reduction->submit();
00288 if (pressureProfileOn)
00289 pressureProfileReduction->submit();
00290 }
00291