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