00001
00007 #include "InfoStream.h"
00008 #include "Node.h"
00009 #include "PatchMap.h"
00010 #include "PatchMap.inl"
00011 #include "AtomMap.h"
00012 #include "ComputeFullDirect.h"
00013 #include "ComputeNonbondedUtil.h"
00014 #include "PatchMgr.h"
00015 #include "Molecule.h"
00016 #include "ReductionMgr.h"
00017 #include "Communicate.h"
00018 #include "Lattice.h"
00019
00020 #define MIN_DEBUG_LEVEL 3
00021 #include "Debug.h"
00022 #include "SimParameters.h"
00023
00024 ComputeFullDirect::ComputeFullDirect(ComputeID c) : ComputeHomePatches(c)
00025 {
00026 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00027 SimParameters *simParams = Node::Object()->simParameters;
00028 if ( simParams->accelMDOn ) {
00029 amd_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_AMD);
00030 } else {
00031 amd_reduction = NULL;
00032 }
00033 useAvgPositions = 1;
00034 }
00035
00036 ComputeFullDirect::~ComputeFullDirect()
00037 {
00038 delete reduction;
00039 delete amd_reduction;
00040 }
00041
00042 BigReal calc_fulldirect(BigReal *data1, BigReal *results1, int n1,
00043 BigReal *data2, BigReal *results2, int n2,
00044 int selfmode, Lattice *lattice, Tensor &virial)
00045 {
00046 if ( lattice->a_p() || lattice->b_p() || lattice->c_p() ) {
00047 #define FULLDIRECT_PERIODIC
00048 #include "ComputeFullDirectBase.h"
00049 } else {
00050 #undef FULLDIRECT_PERIODIC
00051 #include "ComputeFullDirectBase.h"
00052 }
00053 }
00054
00055 void ComputeFullDirect::doWork()
00056 {
00057 int numLocalAtoms;
00058 BigReal *localData;
00059 BigReal *localResults;
00060 BigReal *newLocalResults;
00061 register BigReal *local_ptr;
00062 Lattice *lattice;
00063
00064 int numWorkingPes = (PatchMap::Object())->numNodesWithPatches();
00065
00066 ResizeArrayIter<PatchElem> ap(patchList);
00067
00068
00069 if ( ! patchList[0].p->flags.doFullElectrostatics )
00070 {
00071 for (ap = ap.begin(); ap != ap.end(); ap++) {
00072 CompAtom *x = (*ap).positionBox->open();
00073 Results *r = (*ap).forceBox->open();
00074 (*ap).positionBox->close(&x);
00075 (*ap).forceBox->close(&r);
00076 }
00077 reduction->submit();
00078 if ( amd_reduction ) amd_reduction->submit();
00079 return;
00080 }
00081
00082
00083 numLocalAtoms = 0;
00084 for (ap = ap.begin(); ap != ap.end(); ap++) {
00085 numLocalAtoms += (*ap).p->getNumAtoms();
00086 }
00087
00088 localData = new BigReal[4*numLocalAtoms];
00089 localResults = new BigReal[3*numLocalAtoms];
00090 newLocalResults = new BigReal[3*numLocalAtoms];
00091
00092 lattice = &((*(ap.begin())).p->lattice);
00093
00094
00095 local_ptr = localData;
00096 for (ap = ap.begin(); ap != ap.end(); ap++) {
00097 CompAtom *x = (*ap).positionBox->open();
00098 if ( patchList[0].p->flags.doMolly ) {
00099 (*ap).positionBox->close(&x);
00100 x = (*ap).avgPositionBox->open();
00101 }
00102 int numAtoms = (*ap).p->getNumAtoms();
00103
00104 for(int i=0; i<numAtoms; ++i)
00105 {
00106 *(local_ptr++) = x[i].position.x;
00107 *(local_ptr++) = x[i].position.y;
00108 *(local_ptr++) = x[i].position.z;
00109 *(local_ptr++) = x[i].charge;
00110 }
00111
00112 if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
00113 else { (*ap).positionBox->close(&x); }
00114 }
00115
00116
00117 local_ptr = localResults;
00118 for(int j=0; j<numLocalAtoms; ++j)
00119 {
00120 *(local_ptr++) = 0.;
00121 *(local_ptr++) = 0.;
00122 *(local_ptr++) = 0.;
00123 }
00124
00125
00126 BigReal electEnergy = 0;
00127 Tensor virial;
00128
00129 #define PEMOD(N) (((N)+numWorkingPes)%numWorkingPes)
00130
00131 int numStages = numWorkingPes / 2 + 2;
00132 int lastStage = numStages - 2;
00133 int sendDataPE = PEMOD(CkMyPe()+1);
00134 int recvDataPE = PEMOD(CkMyPe()-1);
00135 int sendResultsPE = PEMOD(CkMyPe()-1);
00136 int recvResultsPE = PEMOD(CkMyPe()+1);
00137 int numRemoteAtoms = numLocalAtoms;
00138 int oldNumRemoteAtoms = 0;
00139 BigReal *remoteData = 0;
00140 BigReal *remoteResults = 0;
00141 register BigReal *remote_ptr;
00142 register BigReal *end_ptr;
00143
00144 MOStream *sendDataMsg=CkpvAccess(comm)->
00145 newOutputStream(sendDataPE, FULLTAG, BUFSIZE);
00146 MIStream *recvDataMsg=CkpvAccess(comm)->
00147 newInputStream(recvDataPE, FULLTAG);
00148
00149 for ( int stage = 0; stage < numStages; ++stage )
00150 {
00151
00152 if ( stage > 1 )
00153 {
00154 DebugM(4,"send remoteResults to sendResultsPE " << sendResultsPE << "\n");
00155 MOStream *msg=CkpvAccess(comm)->
00156 newOutputStream(sendResultsPE, FULLFORCETAG, BUFSIZE);
00157 msg->put(3*oldNumRemoteAtoms,remoteResults);
00158 delete [] remoteResults;
00159 msg->end();
00160 delete msg;
00161 sendResultsPE = PEMOD(sendResultsPE-1);
00162 }
00163
00164
00165 if ( stage < lastStage )
00166 {
00167 DebugM(4,"send remoteData to sendDataPE " << sendDataPE << "\n");
00168 sendDataMsg->put(numRemoteAtoms);
00169 sendDataMsg->put(4*numRemoteAtoms,(stage?remoteData:localData));
00170 sendDataMsg->end();
00171 }
00172
00173
00174 if ( stage > 0 && stage <= lastStage )
00175 {
00176 DebugM(4,"allocate new result storage\n");
00177 remoteResults = new BigReal[3*numRemoteAtoms];
00178 remote_ptr = remoteResults;
00179 end_ptr = remoteResults + 3*numRemoteAtoms;
00180 for ( ; remote_ptr != end_ptr; ++remote_ptr ) *remote_ptr = 0.;
00181 }
00182
00183
00184 if ( stage == 0 )
00185 {
00186 DebugM(4,"self interaction\n");
00187 electEnergy += calc_fulldirect(
00188 localData,localResults,numLocalAtoms,
00189 localData,localResults,numLocalAtoms,1,lattice,virial);
00190 }
00191 else if ( stage < lastStage ||
00192 ( stage == lastStage && ( numWorkingPes % 2 ) ) )
00193 {
00194 DebugM(4,"full other interaction\n");
00195 electEnergy += calc_fulldirect(
00196 localData,localResults,numLocalAtoms,
00197 remoteData,remoteResults,numRemoteAtoms,0,lattice,virial);
00198 }
00199 else if ( stage == lastStage )
00200 {
00201 DebugM(4,"half other interaction\n");
00202 if ( CkMyPe() < ( numWorkingPes / 2 ) )
00203 electEnergy += calc_fulldirect(
00204 localData,localResults,numLocalAtoms/2,
00205 remoteData,remoteResults,numRemoteAtoms,0,lattice,virial);
00206 else
00207 electEnergy += calc_fulldirect(
00208 localData,localResults,numLocalAtoms,
00209 remoteData + 4*(numRemoteAtoms/2),
00210 remoteResults + 3*(numRemoteAtoms/2),
00211 numRemoteAtoms - (numRemoteAtoms/2), 0,lattice,virial);
00212 }
00213
00214 delete [] remoteData; remoteData = 0;
00215 oldNumRemoteAtoms = numRemoteAtoms;
00216
00217
00218 if ( stage > 1 )
00219 {
00220 DebugM(4,"receive newLocalResults from recvResultsPE "
00221 << recvResultsPE << "\n");
00222 MIStream *msg=CkpvAccess(comm)->
00223 newInputStream(recvResultsPE, FULLFORCETAG);
00224 msg->get(3*numLocalAtoms,newLocalResults);
00225 delete msg;
00226 recvResultsPE = PEMOD(recvResultsPE+1);
00227 remote_ptr = newLocalResults;
00228 local_ptr = localResults;
00229 end_ptr = localResults + 3*numLocalAtoms;
00230 for ( ; local_ptr != end_ptr; ++local_ptr, ++remote_ptr )
00231 *local_ptr += *remote_ptr;
00232 }
00233
00234
00235 if ( stage < lastStage )
00236 {
00237 DebugM(4,"receive remoteData from recvDataPE "
00238 << recvDataPE << "\n");
00239 recvDataMsg->get(numRemoteAtoms);
00240 remoteData = new BigReal[4*numRemoteAtoms];
00241 recvDataMsg->get(4*numRemoteAtoms,remoteData);
00242 }
00243
00244 }
00245
00246 delete sendDataMsg;
00247 delete recvDataMsg;
00248
00249
00250 DebugM(4,"Full-electrostatics energy: " << electEnergy << "\n");
00251 reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += electEnergy;
00252 reduction->item(REDUCTION_VIRIAL_SLOW_XX) += virial.xx;
00253 reduction->item(REDUCTION_VIRIAL_SLOW_XY) += virial.xy;
00254 reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += virial.xz;
00255 reduction->item(REDUCTION_VIRIAL_SLOW_YX) += virial.yx;
00256 reduction->item(REDUCTION_VIRIAL_SLOW_YY) += virial.yy;
00257 reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += virial.yz;
00258 reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += virial.zx;
00259 reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += virial.zy;
00260 reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += virial.zz;
00261 reduction->submit();
00262
00263 if ( amd_reduction ) {
00264 amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += electEnergy;
00265 amd_reduction->item(REDUCTION_VIRIAL_SLOW_XX) += virial.xx;
00266 amd_reduction->item(REDUCTION_VIRIAL_SLOW_XY) += virial.xy;
00267 amd_reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += virial.xz;
00268 amd_reduction->item(REDUCTION_VIRIAL_SLOW_YX) += virial.yx;
00269 amd_reduction->item(REDUCTION_VIRIAL_SLOW_YY) += virial.yy;
00270 amd_reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += virial.yz;
00271 amd_reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += virial.zx;
00272 amd_reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += virial.zy;
00273 amd_reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += virial.zz;
00274 amd_reduction->submit();
00275 }
00276
00277
00278 local_ptr = localResults;
00279 for (ap = ap.begin(); ap != ap.end(); ap++) {
00280 Results *r = (*ap).forceBox->open();
00281 Force *f = r->f[Results::slow];
00282 int numAtoms = (*ap).p->getNumAtoms();
00283
00284 for(int i=0; i<numAtoms; ++i)
00285 {
00286 f[i].x += *(local_ptr++);
00287 f[i].y += *(local_ptr++);
00288 f[i].z += *(local_ptr++);
00289 }
00290
00291 (*ap).forceBox->close(&r);
00292 }
00293
00294
00295 delete [] localData;
00296 delete [] localResults;
00297 delete [] newLocalResults;
00298 }
00299