NAMD
ComputeFullDirect.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "Node.h"
9 #include "PatchMap.h"
10 #include "PatchMap.inl"
11 #include "AtomMap.h"
12 #include "ComputeFullDirect.h"
13 #include "ComputeNonbondedUtil.h"
14 #include "PatchMgr.h"
15 #include "Molecule.h"
16 #include "ReductionMgr.h"
17 #include "Communicate.h"
18 #include "Lattice.h"
19 //#define DEBUGM
20 #define MIN_DEBUG_LEVEL 3
21 #include "Debug.h"
22 #include "SimParameters.h"
23 
25 {
28  useAvgPositions = 1;
29 }
30 
32 {
33  delete reduction;
34 }
35 
36 BigReal calc_fulldirect(BigReal *data1, BigReal *results1, int n1,
37  BigReal *data2, BigReal *results2, int n2,
38  int selfmode, Lattice *lattice, Tensor &virial)
39 {
40  if ( lattice->a_p() || lattice->b_p() || lattice->c_p() ) {
41  #define FULLDIRECT_PERIODIC
42  #include "ComputeFullDirectBase.h"
43  } else {
44  #undef FULLDIRECT_PERIODIC
45  #include "ComputeFullDirectBase.h"
46  }
47 }
48 
50 {
51  int numLocalAtoms;
52  BigReal *localData;
53  BigReal *localResults;
54  BigReal *newLocalResults;
55  register BigReal *local_ptr;
56  Lattice *lattice;
57 
58  int numWorkingPes = (PatchMap::Object())->numNodesWithPatches();
59 
60  if ( numWorkingPes > 1 ) NAMD_die("FullDirect not supported for parallel runs.");
61 
63 
64  // Skip computations if nothing to do.
65  if ( ! patchList[0].p->flags.doFullElectrostatics )
66  {
67  for (ap = ap.begin(); ap != ap.end(); ap++) {
68  CompAtom *x = (*ap).positionBox->open();
69  Results *r = (*ap).forceBox->open();
70  (*ap).positionBox->close(&x);
71  (*ap).forceBox->close(&r);
72  }
73  reduction->submit();
74  return;
75  }
76 
77  // allocate storage
78  numLocalAtoms = 0;
79  for (ap = ap.begin(); ap != ap.end(); ap++) {
80  numLocalAtoms += (*ap).p->getNumAtoms();
81  }
82 
83  localData = new BigReal[4*numLocalAtoms]; // freed at end of this method
84  localResults = new BigReal[3*numLocalAtoms]; // freed at end of this method
85  newLocalResults = new BigReal[3*numLocalAtoms]; // freed at end of this method
86 
87  lattice = &((*(ap.begin())).p->lattice);
88 
89  // get positions and charges
90  local_ptr = localData;
91  for (ap = ap.begin(); ap != ap.end(); ap++) {
92  CompAtom *x = (*ap).positionBox->open();
93  if ( patchList[0].p->flags.doMolly ) {
94  (*ap).positionBox->close(&x);
95  x = (*ap).avgPositionBox->open();
96  }
97  int numAtoms = (*ap).p->getNumAtoms();
98 
99  for(int i=0; i<numAtoms; ++i)
100  {
101  *(local_ptr++) = x[i].position.x;
102  *(local_ptr++) = x[i].position.y;
103  *(local_ptr++) = x[i].position.z;
104  *(local_ptr++) = x[i].charge;
105  }
106 
107  if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
108  else { (*ap).positionBox->close(&x); }
109  }
110 
111  // zero out forces
112  local_ptr = localResults;
113  for(int j=0; j<numLocalAtoms; ++j)
114  {
115  *(local_ptr++) = 0.;
116  *(local_ptr++) = 0.;
117  *(local_ptr++) = 0.;
118  }
119 
120  // perform calculations
121  BigReal electEnergy = 0;
122  Tensor virial;
123 
124 #define PEMOD(N) (((N)+numWorkingPes)%numWorkingPes)
125 
126  int numStages = numWorkingPes / 2 + 2;
127  int lastStage = numStages - 2;
128  int sendDataPE = PEMOD(CkMyPe()+1);
129  int recvDataPE = PEMOD(CkMyPe()-1);
130  int sendResultsPE = PEMOD(CkMyPe()-1);
131  int recvResultsPE = PEMOD(CkMyPe()+1);
132  int numRemoteAtoms = numLocalAtoms;
133  int oldNumRemoteAtoms = 0;
134  BigReal *remoteData = 0;
135  BigReal *remoteResults = 0;
136  register BigReal *remote_ptr;
137  register BigReal *end_ptr;
138 
139  MOStream *sendDataMsg=CkpvAccess(comm)->
140  newOutputStream(sendDataPE, FULLTAG, BUFSIZE);
141  MIStream *recvDataMsg=CkpvAccess(comm)->
142  newInputStream(recvDataPE, FULLTAG);
143 
144  for ( int stage = 0; stage < numStages; ++stage )
145  {
146  // send remoteResults to sendResultsPE
147  if ( stage > 1 )
148  {
149  DebugM(4,"send remoteResults to sendResultsPE " << sendResultsPE << "\n");
150  MOStream *msg=CkpvAccess(comm)->
151  newOutputStream(sendResultsPE, FULLFORCETAG, BUFSIZE);
152  msg->put(3*oldNumRemoteAtoms,remoteResults);
153  delete [] remoteResults;
154  msg->end();
155  delete msg;
156  sendResultsPE = PEMOD(sendResultsPE-1);
157  }
158 
159  // send remoteData to sendDataPE
160  if ( stage < lastStage )
161  {
162  DebugM(4,"send remoteData to sendDataPE " << sendDataPE << "\n");
163  sendDataMsg->put(numRemoteAtoms);
164  sendDataMsg->put(4*numRemoteAtoms,(stage?remoteData:localData));
165  sendDataMsg->end();
166  }
167 
168  // allocate new result storage
169  if ( stage > 0 && stage <= lastStage )
170  {
171  DebugM(4,"allocate new result storage\n");
172  remoteResults = new BigReal[3*numRemoteAtoms];
173  remote_ptr = remoteResults;
174  end_ptr = remoteResults + 3*numRemoteAtoms;
175  for ( ; remote_ptr != end_ptr; ++remote_ptr ) *remote_ptr = 0.;
176  }
177 
178  // do calculation
179  if ( stage == 0 )
180  { // self interaction
181  DebugM(4,"self interaction\n");
182  electEnergy += calc_fulldirect(
183  localData,localResults,numLocalAtoms,
184  localData,localResults,numLocalAtoms,1,lattice,virial);
185  }
186  else if ( stage < lastStage ||
187  ( stage == lastStage && ( numWorkingPes % 2 ) ) )
188  { // full other interaction
189  DebugM(4,"full other interaction\n");
190  electEnergy += calc_fulldirect(
191  localData,localResults,numLocalAtoms,
192  remoteData,remoteResults,numRemoteAtoms,0,lattice,virial);
193  }
194  else if ( stage == lastStage )
195  { // half other interaction
196  DebugM(4,"half other interaction\n");
197  if ( CkMyPe() < ( numWorkingPes / 2 ) )
198  electEnergy += calc_fulldirect(
199  localData,localResults,numLocalAtoms/2,
200  remoteData,remoteResults,numRemoteAtoms,0,lattice,virial);
201  else
202  electEnergy += calc_fulldirect(
203  localData,localResults,numLocalAtoms,
204  remoteData + 4*(numRemoteAtoms/2),
205  remoteResults + 3*(numRemoteAtoms/2),
206  numRemoteAtoms - (numRemoteAtoms/2), 0,lattice,virial);
207  }
208 
209  delete [] remoteData; remoteData = 0;
210  oldNumRemoteAtoms = numRemoteAtoms;
211 
212  // receive newLocalResults from recvResultsPE
213  if ( stage > 1 )
214  {
215  DebugM(4,"receive newLocalResults from recvResultsPE "
216  << recvResultsPE << "\n");
217  MIStream *msg=CkpvAccess(comm)->
218  newInputStream(recvResultsPE, FULLFORCETAG);
219  msg->get(3*numLocalAtoms,newLocalResults);
220  delete msg;
221  recvResultsPE = PEMOD(recvResultsPE+1);
222  remote_ptr = newLocalResults;
223  local_ptr = localResults;
224  end_ptr = localResults + 3*numLocalAtoms;
225  for ( ; local_ptr != end_ptr; ++local_ptr, ++remote_ptr )
226  *local_ptr += *remote_ptr;
227  }
228 
229  // receive remoteData from recvDataPE
230  if ( stage < lastStage )
231  {
232  DebugM(4,"receive remoteData from recvDataPE "
233  << recvDataPE << "\n");
234  recvDataMsg->get(numRemoteAtoms);
235  remoteData = new BigReal[4*numRemoteAtoms];
236  recvDataMsg->get(4*numRemoteAtoms,remoteData);
237  }
238 
239  }
240 
241  delete sendDataMsg;
242  delete recvDataMsg;
243 
244  // send out reductions
245  DebugM(4,"Full-electrostatics energy: " << electEnergy << "\n");
247  reduction->item(REDUCTION_VIRIAL_SLOW_XX) += virial.xx;
248  reduction->item(REDUCTION_VIRIAL_SLOW_XY) += virial.xy;
249  reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += virial.xz;
250  reduction->item(REDUCTION_VIRIAL_SLOW_YX) += virial.yx;
251  reduction->item(REDUCTION_VIRIAL_SLOW_YY) += virial.yy;
252  reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += virial.yz;
253  reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += virial.zx;
254  reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += virial.zy;
255  reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += virial.zz;
256  reduction->submit();
257 
258  // add in forces
259  local_ptr = localResults;
260  for (ap = ap.begin(); ap != ap.end(); ap++) {
261  Results *r = (*ap).forceBox->open();
262  Force *f = r->f[Results::slow];
263  int numAtoms = (*ap).p->getNumAtoms();
264 
265  for(int i=0; i<numAtoms; ++i)
266  {
267  f[i].x += *(local_ptr++);
268  f[i].y += *(local_ptr++);
269  f[i].z += *(local_ptr++);
270  }
271 
272  (*ap).forceBox->close(&r);
273  }
274 
275  // free storage
276  delete [] localData; // allocated at beginning of this method
277  delete [] localResults; // allocated at beginning of this method
278  delete [] newLocalResults; // allocated at beginning of this method
279 }
280 
static Node * Object()
Definition: Node.h:86
ComputeFullDirect(ComputeID c)
BigReal zy
Definition: Tensor.h:19
void end(void)
Definition: MStream.C:176
#define PEMOD(N)
BigReal xz
Definition: Tensor.h:17
int ComputeID
Definition: NamdTypes.h:183
static PatchMap * Object()
Definition: PatchMap.h:27
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
register BigReal electEnergy
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
BigReal yz
Definition: Tensor.h:18
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
MIStream * get(char &data)
Definition: MStream.h:29
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
Charge charge
Definition: NamdTypes.h:54
ResizeArrayIter< T > end(void) const
BigReal yx
Definition: Tensor.h:18
virtual ~ComputeFullDirect()
Force * f[maxNumForces]
Definition: PatchTypes.h:67
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:83
#define BUFSIZE
Definition: Communicate.h:15
BigReal xx
Definition: Tensor.h:17
BigReal zz
Definition: Tensor.h:19
#define simParams
Definition: Output.C:127
Definition: Tensor.h:15
BigReal xy
Definition: Tensor.h:17
BigReal y
Definition: Vector.h:66
#define FULLFORCETAG
Definition: common.h:154
BigReal yy
Definition: Tensor.h:18
BigReal calc_fulldirect(BigReal *data1, BigReal *results1, int n1, BigReal *data2, BigReal *results2, int n2, int selfmode, Lattice *lattice, Tensor &virial)
MOStream * put(char data)
Definition: MStream.h:112
void submit(void)
Definition: ReductionMgr.h:323
int b_p() const
Definition: Lattice.h:274
gridSize x
int a_p() const
Definition: Lattice.h:273
BigReal zx
Definition: Tensor.h:19
#define FULLTAG
Definition: common.h:153
ResizeArrayIter< T > begin(void) const
double BigReal
Definition: common.h:112
int c_p() const
Definition: Lattice.h:275