OpenAtom  Version1.5a
AtomsCompute.C
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////=
2 ///////////////////////////////////////////////////////////////////////////////c
3 ///////////////////////////////////////////////////////////////////////////////=
4 /** \file AtomsCompute.C
5  * Processor group class Functions : Atoms
6  */
7 ///////////////////////////////////////////////////////////////////////////////=
8 #include "AtomsCache.h"
9 #include "AtomsCompute.h"
10 #include "energyGroup.h"
11 #include "eesCache.h"
12 #include "cp_state_ctrl/CP_State_GSpacePlane.h"
13 #include "fft_slab_ctrl/fftCacheSlab.h"
14 #include "utility/util.h"
15 #include "CPcharmParaInfoGrp.h"
16 #include "load_balance/IntMap.h"
17 #include "charm++.h"
18 #include "src_piny_physics_v1.0/include/class_defs/Interface_ctrl.h"
19 #include "PhysScratchCache.h"
20 
21 #include <cmath>
22 
23 //#include "CP_State_Plane.h"
24 
25 //----------------------------------------------------------------------------
26 #define CHARM_ON
27 #include "src_piny_physics_v1.0/include/class_defs/piny_constants.h"
28 #include "src_piny_physics_v1.0/include/class_defs/ATOM_OPERATIONS/class_atomintegrate.h"
29 #include "src_piny_physics_v1.0/include/class_defs/ATOM_OPERATIONS/class_atomoutput.h"
30 #include "src_piny_physics_v1.0/include/class_defs/CP_OPERATIONS/class_cprspaceion.h"
31 
32 //----------------------------------------------------------------------------
33 
34 extern CkVec <CProxy_PIBeadAtoms> UPIBeadAtomsProxy;
35 extern CkVec <IntMap2on2> GSImaptable;
36 extern CkVec <CProxy_EnergyGroup> UegroupProxy;
37 extern Config config;
38 extern CkVec <CProxy_CP_State_GSpacePlane> UgSpacePlaneProxy;
39 extern CkVec <CProxy_GSpaceDriver> UgSpaceDriverProxy;
40 extern CkVec <CProxy_AtomsCache> UatomsCacheProxy;
41 extern CkVec <CProxy_AtomsCompute> UatomsComputeProxy;
42 extern CkVec <CProxy_EnergyGroup> UegroupProxy;
43 extern CkVec <CProxy_StructFactCache> UsfCacheProxy;
44 extern CkVec <CProxy_eesCache> UeesCacheProxy;
45 extern CProxy_TemperController temperControllerProxy;
46 extern CProxy_InstanceController instControllerProxy;
47 extern CProxy_PhysScratchCache pScratchProxy;
48 
49 //----------------------------------------------------------------------------
50 
51 //#define _CP_DEBUG_PSI_OFF_
52 //#define _CP_ENERGY_GRP_VERBOSE_
53 //#define _CP_DEBUG_ATMS_
54 //#define _CP_DEBUG_ATMS_EXIT_
55 
56 ///////////////////////////////////////////////////////////////////////////////=
57 ///////////////////////////////////////////////////////////////////////////////c
58 ///////////////////////////////////////////////////////////////////////////////=
59 /**
60  *
61  * @addtogroup Atoms
62  * @{
63  *
64  */
65 ///////////////////////////////////////////////////////////////////////////////=
66 AtomsCompute::AtomsCompute(int n, int n_nl, int len_nhc_, int iextended_on_,int cp_min_opt_,
67  int cp_wave_opt_, int isokin_opt_,double kT_, Atom* a, AtomNHC *aNHC, int nChareAtoms_,
68  UberCollection _thisInstance) : natm(n), natm_nl(n_nl), len_nhc(len_nhc_), iextended_on(iextended_on_), cp_min_opt(cp_min_opt_), cp_wave_opt(cp_wave_opt_), isokin_opt(isokin_opt_), kT(kT_), nChareAtoms(nChareAtoms_), thisInstance(_thisInstance)
69 ///////////////////////////////////////////////////////////////////////////////=
70  {// begin routine
71 ///////////////////////////////////////////////////////////////////////////////=
72 /// parameters, options and energies
73  handleForcesCount=0;
74  ktemps = 0;
75  pot_ewd_rs = 0.0;
76  eKinetic = 0.0;
77  eKineticNhc = 0.0;
78  potNhc = 0.0;
79  potPIMDChain = 0.0;
80  countAtm = 0;
81  acceptCountfu = 0;
82  acceptCountX = 0;
83  acceptCountu = 0;
84  atomsCMrecv=atomsPIMDXrecv=false;
85  temperScreenFile = (UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch())->temperScreenFile;
86 
87 ///////////////////////////////////////////////////////////////////////////////=
88 /// Initial positions, forces, velocities
89  numPIMDBeads = config.UberImax;
90  PIBeadIndex = thisInstance.idxU.x;
91  TemperIndex = thisInstance.idxU.z;
92  atoms = new Atom[natm];
93  atomsNHC = new AtomNHC[natm];
94  AtomsCache *ag = UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
95  // we keep our own copy of this annoying thing to avoid data races
96  // with the force accumulation in the AtomsCache.
97  fastAtoms.q = new double[natm];
98  fastAtoms.x = new double[natm];
99  fastAtoms.y = new double[natm];
100  fastAtoms.z = new double[natm];
101  fastAtoms.fx = new double[natm];
102  fastAtoms.fy = new double[natm];
103  fastAtoms.fz = new double[natm];
104  fastAtoms.fxu = new double[natm];
105  fastAtoms.fyu = new double[natm];
106  fastAtoms.fzu = new double[natm];
107  for(int i=0;i<natm;i++){
108  fastAtoms.x[i] = a[i].x;
109  fastAtoms.y[i] = a[i].y;
110  fastAtoms.z[i] = a[i].z;
111  fastAtoms.fx[i] = a[i].fx;
112  fastAtoms.fy[i] = a[i].fy;
113  fastAtoms.fz[i] = a[i].fz;
114  }//endfor
115  iteration = &(ag->iteration);
116 
117  CmiMemcpy(atoms, a, natm * sizeof(Atom)); // atoms has no vectors
118  for(int i=0;i<natm;i++){atomsNHC[i].Init(&aNHC[i]);}
119  ftot = (double *)fftw_malloc((3*natm+2)*sizeof(double));
120 
121  zeroforces();
122  if(iextended_on==1 && cp_min_opt==0){
123  zeronhc();
124  }//endif
125 
126 ///////////////////////////////////////////////////////////////////////////////=
127 /// A copy of the atoms for fast routines : only need to copy charge once
128 
129  copyFastToSlow();
130  double *qq = fastAtoms.q;
131  for(int i=0;i<natm;i++){qq[i]=atoms[i].q;}
132 
133 ///////////////////////////////////////////////////////////////////////////////=
134 /// Number of messages to be received when atoms are moved : simple parallel
135 
136  int nproc = nChareAtoms;
137  int div = (natm / nproc);
138  int rem = (natm % nproc);
139 
140  nAtmMsgRecv = 0;
141  for(int myid=0;myid<nproc;myid++){
142  int natmNow = div;
143  if(myid< rem){natmNow++;}
144  if(natmNow>0){nAtmMsgRecv++;}
145  }//endfor
146 
147 ///////////////////////////////////////////////////////////////////////////////=
148 /// PIMD set up : Even if classical, this can run. It is harmless and
149 /// prevents careless bugs
150 
151  massPIMDScal = (double *)fftw_malloc(numPIMDBeads*sizeof(double));
152  initPIMD();
153 
154 //-----------------------------------------------------------------------------
155  }//end routine
156 ///////////////////////////////////////////////////////////////////////////////=
157 
158 
159 ///////////////////////////////////////////////////////////////////////////////=
160 ///////////////////////////////////////////////////////////////////////////////c
161 ///////////////////////////////////////////////////////////////////////////////=
162 /** Destructor
163  *
164  *
165  **/
166 ///////////////////////////////////////////////////////////////////////////////=
167 AtomsCompute::~AtomsCompute(){
168  delete [] atoms;
169  delete [] atomsNHC;
170  delete [] fastAtoms.q;
171  delete [] fastAtoms.x;
172  delete [] fastAtoms.y;
173  delete [] fastAtoms.z;
174  delete [] fastAtoms.fx;
175  delete [] fastAtoms.fy;
176  delete [] fastAtoms.fz;
177  delete [] fastAtoms.fxu;
178  delete [] fastAtoms.fyu;
179  delete [] fastAtoms.fzu;
180  fftw_free(PIMD_CM_Atoms.x);
181  fftw_free(PIMD_CM_Atoms.y);
182  fftw_free(PIMD_CM_Atoms.z);
183  fftw_free(ftot);
184 }
185 ///////////////////////////////////////////////////////////////////////////////=
186 
187 
188 ///////////////////////////////////////////////////////////////////////////=
189 /// Bead root is the 0th element. The AtomsCompute map will lookup the
190 /// GSP map to ensure that the 0th AtomCompute is co-mapped with GSP(0,0) in
191 /// each bead instance.
192 
193 ///////////////////////////////////////////////////////////////////////////=
194 ///////////////////////////////////////////////////////////////////////////c
195 ///////////////////////////////////////////////////////////////////////////=
196 
197 void AtomsCompute::init()
198 {
199 /// Bead root should be the 0th element. The AtomsCompute map will use
200 /// the GSP map to ensure that the 0th AtomCompute is co-mapped
201 /// with GSP(0,0) in each bead instance.
202  const int offset=thisInstance.getPO();
203  amBeadRoot = (CkMyPe()==GSImaptable[offset].get(0,0));
204  amZerothBead = (thisInstance.idxU.x==0);
205 
206  if(numPIMDBeads>1)
207  {
208  // make lists of Array IDs and their Commander PE
209  CkArrayID *atomsArrayids= new CkArrayID[numPIMDBeads];
210  CkArrayIndex **elems = new CkArrayIndex*[numPIMDBeads];
211  CkArrayIndex **elemsAll = new CkArrayIndex*[numPIMDBeads];
212  int *naelems = new int[numPIMDBeads];
213  int *naelemsAll = new int[numPIMDBeads];
214  for(int i=0;i<numPIMDBeads;i++){
215  elems[i]= new CkArrayIndex1D[1];
216  elemsAll[i]= new CkArrayIndex1D[nChareAtoms];
217  UberCollection instance=thisInstance;
218  naelems[i]=1;
219  naelemsAll[i]=nChareAtoms;
220  instance.idxU.x=i;
221  int offset=instance.calcPO();
222  elems[i][0]=CkArrayIndex1D(0);
223  for(int j=0;j<nChareAtoms;j++)
224  elemsAll[i][j]=CkArrayIndex1D(j);
225  atomsArrayids[i]=UatomsComputeProxy[offset].ckGetArrayID();
226  //CkPrintf("{%d}[%d] AtomsCompute::init elems[%d][0]=%d, atomsgrpids[%d]=%d amBeadRoot=%d\n",thisInstance.proxyOffset, CkMyPe(), i, elems[i][0], i, atomsgrpids[i], amBeadRoot);
227  }//endfor
228  // proxyHeadBeads=CProxySection_AtomsCompute(numPIMDBeads, atomsArrayids, elems, naelems);
229  proxyAllBeads=CProxySection_AtomsCompute(numPIMDBeads, atomsArrayids, elemsAll, naelemsAll);
230  delete [] naelems;
231  delete [] naelemsAll;
232  for(int i=0;i<numPIMDBeads;i++){delete [] elems[i]; delete [] elemsAll[i];}
233  delete [] elemsAll;
234  delete [] atomsArrayids;
235 
236  }//endif
237 
238 ///////////////////////////////////////////////////////////////////////////////=
239  }//end routine
240 ///////////////////////////////////////////////////////////////////////////////=
241 
242 
243 ///////////////////////////////////////////////////////////////////////////=
244 ///////////////////////////////////////////////////////////////////////////c
245 ///////////////////////////////////////////////////////////////////////////=
246 /** recvContribute
247  * Every Array member is sent all the forces at present.
248  */
249 ///////////////////////////////////////////////////////////////////////////=
250 void AtomsCompute::recvContribute(CkReductionMsg *msg) {
251 //////////////////////////////////////////////////////////////
252 #ifdef _CP_DEBUG_ATMS_
253  CkPrintf("{%d}[%d] AtomsCompute::recvContribute count %d\n ", thisInstance.proxyOffset, thisIndex, handleForcesCount);
254 #endif
255 ///////////////////////////////////////////////////////////////////////////=
256 /// Local pointers
257 
258  int i,j;
259  int myid = thisIndex;
260  contribMsg[handleForcesCount++]=msg;
261  EnergyGroup *eg = UegroupProxy[thisInstance.proxyOffset].ckLocalBranch();
262 
263 //////////////////////////////////////////////////////////////
264 /// Copy out the reduction of energy and forces and nuke msg
265  double *ftot = (double *) msg->getData();
266 
267  double pot_ewd_rs_loc = ftot[3*natm];
268  double potPerdCorrLoc = ftot[3*natm+1];
269 
270 //////////////////////////////////////////////////////////////
271 /// Tuck things that can be tucked.
272 
273  eg->estruct.eewald_real = pot_ewd_rs_loc;
274 
275 
276 #ifdef _CP_DEBUG_ATMS_EXIT_
277  if(myid==0){CkExit();}
278 #endif
279  if(handleForcesCount==2)
280  handleForces();
281 
282 ///////////////////////////////////////////////////////////////////////////=
283  }//end routine
284 ///////////////////////////////////////////////////////////////////////////=
285 
286 
287 ///////////////////////////////////////////////////////////////////////////=
288 ///////////////////////////////////////////////////////////////////////////c
289 ///////////////////////////////////////////////////////////////////////////=
290 /** recvContributeForces
291  * Every Array member has all the forces at present.
292  */
293 ///////////////////////////////////////////////////////////////////////////=
294 void AtomsCompute::recvContributeForces(CkReductionMsg *msg) {
295 //////////////////////////////////////////////////////////////
296 #ifdef _CP_DEBUG_ATMS_
297  CkPrintf("{%d}[%d] AtomsCompute::recvContributeForces count %d\n ", thisInstance.proxyOffset, thisIndex, handleForcesCount);
298 #endif
299 ///////////////////////////////////////////////////////////////////////////=
300 /// Local pointers
301  int myid = thisIndex;
302  contribMsg[handleForcesCount++]=msg;
303 
304 #ifdef _CP_DEBUG_ATMS_EXIT_
305  if(myid==0){CkExit();}
306 #endif
307 
308  if(handleForcesCount==2)
309  handleForces();
310 ///////////////////////////////////////////////////////////////////
311  }//end routine
312 ///////////////////////////////////////////////////////////////////
313 
314 ///////////////////////////////////////////////////////////////////
315 ///////////////////////////////////////////////////////////////////
316 ///////////////////////////////////////////////////////////////////
317 void AtomsCompute::handleForces() {
318 ///////////////////////////////////////////////////////////////////
319 #ifdef _CP_DEBUG_ATMS_
320  CkPrintf("{%d}[%d] AtomsCompute::handleForces \n ", thisInstance.proxyOffset, thisIndex);
321 #endif
322 ///////////////////////////////////////////////////////////////////
323  AtomsCache *ag = UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
324  int iteration = ag->iteration;
325 ///////////////////////////////////////////////////////////////////
326  copyFastToSlow();
327  zeroforces(); // we're now getting everyone's forces
328  double fmag=0.0;
329  handleForcesCount=0;
330  EnergyGroup *eg = UegroupProxy[thisInstance.proxyOffset].ckLocalBranch();
331 /// Model forces : This is fine for path integral checking too
332 
333  // sum the contributions
334  double *ftot0 = (double *) contribMsg[0]->getData();
335  double *ftot1 = (double *) contribMsg[1]->getData();
336  // no loop carried dependencies here, so a compiler worth its salt
337  // should easily vectorize these ops.
338  int i, j;
339  for(i=0,j=0;i<natm;i++,j+=3){
340  atoms[i].fx = ftot0[j] + ftot1[j];
341  atoms[i].fy = ftot0[j+1] + ftot1[j+1];
342  atoms[i].fz = ftot0[j+2] + ftot1[j+2];
343 #ifdef _CP_DEBUG_ATMS_
344  if(CkMyPe()==0)
345  {
346  CkPrintf("AtomCompute handleForces forces %d %.5g,%.5g,%.5g\n",i, atoms[i].fx, atoms[i].fy, atoms[i].fz);
347  }
348 #endif
349  fmag += atoms[i].fx * atoms[i].fx + atoms[i].fy * atoms[i].fy + atoms[i].fz * atoms[i].fz;
350  }// endfor
351 
352  fmag /= (double)(3*natm);
353  fmag = sqrt(fmag);
354  eg->estruct.fmag_atm = fmag;
355  delete contribMsg[0];
356  delete contribMsg[1];
357 
358 #ifdef _CP_DEBUG_PSI_OFF_
359  double omega = (0.0241888/15.0); // 15 fs^{-1}
360  double omega2 = omega*omega;
361  int npts = 200;
362  if(iteration==0){
363  px = new double *[natm];
364  for(i =0;i<natm;i++){
365  px[i] = new double[npts];
366  for(j=0;j<npts;j++){px[i][j]=0.0;}
367  }//endfor
368  for(i=0;i<natm;i++){
369  atoms[i].x = 0.0;
370  atoms[i].y = 0.0;
371  atoms[i].z = 0.0;
372  }//endfor
373  }//endif
374 
375  for(i=0;i<natm;i++){
376  atoms[i].fx = -omega2*atoms[i].x*atoms[i].m;
377  atoms[i].fy = -omega2*atoms[i].y*atoms[i].m;;
378  atoms[i].fz = -omega2*atoms[i].z*atoms[i].m;;
379  }//endfor
380 #endif
381 
382 ///////////////////////////////////////////////////////////////////
383 /// Check the forces
384 
385 #ifdef _GLENN_CHECK_FORCES
386  FILE *fp = fopen("forces.save","w");
387  for(i=0;i<natm;i++){
388  fprintf(fp,"%.12g %.12g %.12g\n",atoms[i].fx,atoms[i].fy,atoms[i].fz);
389  }
390  fclose(fp);
391  CkExit();
392 #endif
393 
394 
395 ///////////////////////////////////////////////////////////////////////////=
396 /// if classical go on to integration, otherwise Fx -> Fu
397 
398  if(numPIMDBeads>1){
399  if(iteration==0){
400  send_PIMD_Fx_and_x(); // atom integration must wait on both transforms
401  }else{
402  send_PIMD_Fx(); // atom integration must wait on Fx->Fu transformation
403  }//endif
404  }else{
405  integrateAtoms(); // Fx is all you need so go forth and integrate
406  }//endif
407 
408 ///////////////////////////////////////////////////////////////////////////=
409  }//end routine
410 ///////////////////////////////////////////////////////////////////////////=
411 
412 ///////////////////////////////////////////////////////////////////////////=
413 ///////////////////////////////////////////////////////////////////////////c
414 ///////////////////////////////////////////////////////////////////////////=
415 /**
416  * Integrate atoms. This is parallelized so that a subset of the atoms are
417  * computed on each processor and their results sent to AtomCompute->acceptAtoms().
418  */
419 ///////////////////////////////////////////////////////////////////////////=
420 void AtomsCompute::integrateAtoms(){
421 //////////////////////////////////////////////////////////////
422 #ifdef _CP_DEBUG_ATMS_
423  CkPrintf("GJM_DBG: Before atom integrate %d : %d\n",thisIndex,natm);
424 #endif
425 //////////////////////////////////////////////////////////////
426 /// Local pointers
427 
428  int mybead = thisInstance.idxU.x+1;
429  int myid = thisIndex;
430 
431  EnergyGroup *eg = UegroupProxy[thisInstance.proxyOffset].ckLocalBranch();
432 
433 ///////////////////////////////////////////////////////////////////
434 /// Atom parallelization: This decompostion is done dynamically because
435 /// its a pretty trivial choice and it allows us to support atom
436 /// migration into and out of the QM region as long that atom migration
437 /// mechanism updates natm and the atom array accordingly.
438 
439  int div = (natm / nChareAtoms);
440  int rem = (natm % nChareAtoms);
441  int natmStr = div*myid;
442  int natmNow = div;
443  if(myid>=rem){natmStr += rem;}
444  if(myid< rem){natmStr += myid;}
445  if(myid< rem){natmNow++;}
446  int natmEnd = natmNow+natmStr;
447 
448 //////////////////////////////////////////////////////////////
449 /// Zero some local copies energies and flags
450 
451  double eKinetic_loc = 0.0;
452  double eKineticNhc_loc = 0.0;
453  double potNhc_loc = 0.0;
454  double potPIMDChain_loc = 0.0;
455  int iwrite_atm = 0;
456  int myoutput_on = 0;
457 
458 //////////////////////////////////////////////////////////////
459 /// DEBUGGING : Compute the distribution function for model
460 
461 #ifdef _CP_DEBUG_PSI_OFF_
462  if(nChareAtoms>1 && numPIMDBeads>1){
463  CkPrintf("Harmonic oscillator debug test currently broken for PIMD\n");
464  CkPrintf("Need a reduction over beads of the chain energy\n");
465  CkExit();
466  }//endif
467  double omega = (0.0241888/15.0); // 15 fs^{-1} Ok for PIMD
468  double omega2 = omega*omega;
469  int npts = 200;
470  double sigma;
471  if(numPIMDBeads==1){
472  sigma = 1.0/sqrt(kT*atoms[0].m*omega2);
473  }else{
474  sigma = 1.0/sqrt(2.0*atoms[0].m*omega); // Assuming lots of beads
475  }//endif
476 
477  double dx = 6.0*sigma/(double)npts;
478  for(int i=natmStr;i<natmEnd;i++){
479  int i1 = (atoms[i].x+3.0*sigma)/dx;
480  int i2 = (atoms[i].y+3.0*sigma)/dx;
481  int i3 = (atoms[i].z+3.0*sigma)/dx;
482  i1 = (i1 >= 0 ? i1 : 0);
483  i2 = (i2 >= 0 ? i2 : 0);
484  i3 = (i3 >= 0 ? i3 : 0);
485  i1 = (i1 < npts ? i1 : npts-1);
486  i2 = (i2 < npts ? i2 : npts-1);
487  i3 = (i3 < npts ? i3 : npts-1);
488  px[i][i1] += 1.0;
489  px[i][i2] += 1.0;
490  px[i][i3] += 1.0;
491  }//endif
492 
493  if(((*iteration+1) % 1000)==0){
494  for(int i=natmStr;i<natmEnd;i++){
495  char fname[100];
496  sprintf(fname,"atmtest.dat.%d",i);
497  double anorm = sqrt(2.0/(sigma*sigma*acos(-1.0)));
498  FILE *fp = fopen(fname,"w");
499  for(int j =1;j<npts-1;j++){
500  double x = ((double)(2*j+1))*0.5*dx - 3.0*sigma;
501  double ans = anorm*exp(-0.5*x*x/(sigma*sigma));
502  double cnt = 3.0*((double)(*iteration+1));
503  fprintf(fp,"%g %g %g\n",x,(px[i][j]*anorm/cnt),ans);
504  }//endfor
505  fclose(fp);
506  }//endfor
507  }//endif
508 
509  double pot_harm = 0.0;
510  for(int i=natmStr;i<natmEnd;i++){
511  pot_harm += (atoms[i].m*omega2*(atoms[i].x*atoms[i].x+
512  atoms[i].y*atoms[i].y+
513  atoms[i].z*atoms[i].z));
514  }//endfor
515  pot_harm *= 0.5;
516 #endif
517 
518 //////////////////////////////////////////////////////////////
519 /// Path integral : Overwrite variables to keep integrator clean:
520 /// Add the chain force to transformed forces.
521 /// Scale the masses to fict bead masses
522 /// This is PIMD Respa implementation friendly
523 
524  if(numPIMDBeads>1 && cp_min_opt==0 && cp_wave_opt==0 && natmNow>0){
525  switchPIMDBeadForceMass(mybead,natmStr,natmEnd,&potPIMDChain_loc);
526  }//endif
527 
528 //////////////////////////////////////////////////////////////
529 /// Integrate the atoms : Path Integral Ready
530 
531 #ifndef _CP_DEBUG_SCALC_ONLY_
532  ATOMINTEGRATE::ctrl_atom_integrate(*iteration,natm,len_nhc,cp_min_opt,
533  cp_wave_opt,iextended_on,atoms,atomsNHC,myid,
534  &eKinetic_loc,&eKineticNhc_loc,&potNhc_loc,&iwrite_atm,
535  myoutput_on,natmNow,natmStr,natmEnd,mybead);
536 #endif
537 
538 //////////////////////////////////////////////////////////////
539 /// Path integral : Rescale to the physical masses
540 
541  if(numPIMDBeads>1 && cp_min_opt==0 && cp_wave_opt==0 && natmNow>0){
542  unswitchPIMDMass(mybead,natmStr,natmEnd);
543  }//endif
544 
545 //////////////////////////////////////////////////////////////
546 /// Debug output : Not correct for PIMD which needs chain PE
547 
548 #ifdef _CP_DEBUG_PSI_OFF_
549  double etot_atm;
550  if(isokin_opt==0){
551  etot_atm = eKinetic_loc+eKineticNhc_loc+potNhc_loc+pot_harm+potPIMDChain_loc;
552  CkPrintf("iteration %d : tot class energy %.12g on %d\n",*iteration,etot_atm,myid);
553  }else{
554  etot_atm = eKineticNhc_loc+potNhc_loc;
555  CkPrintf("iteration %d : tot class energy %.12g %.12g on %d\n",*iteration,etot_atm,
556  (eKinetic_loc+pot_harm),myid);
557  }//endif
558 #endif
559 
560 //////////////////////////////////////////////////////////////
561 /// Get ready for the next iteration :
562 /// zero forces, outputAtmEnergy, atomsDone
563 
564  zeroforces();
565 
566  if(cp_wave_opt==0 && cp_min_opt==0){
567  if(natmNow>0){
568  sendAtoms(eKinetic_loc,eKineticNhc_loc,potNhc_loc,potPIMDChain_loc,natmNow,natmStr,natmEnd);
569  }//endif
570  }else{
571  // in the minimization case the atoms don't move so we don't
572  // update the coordinates
573  eKinetic = 0.0;
574  eKineticNhc = 0.0;
575  potNhc = 0.0;
576  potPIMDChain = 0.0;
577  eg->estruct.eKinetic_atm = 0.0;
578  eg->estruct.eKineticNhc_atm = 0.0;
579  eg->estruct.potNhc_atm = 0.0;
580  eg->estruct.potPIMDChain = 0.0;
581 
582  copySlowToFast();
583  outputAtmEnergy();
584 
585  int i=0;
586  CkCallback cb(CkIndex_AtomsCompute::atomsDone(NULL), CkArrayIndex1D(0), UatomsComputeProxy[thisInstance.proxyOffset]);
587  contribute(sizeof(int),&i,CkReduction::sum_int,cb);
588  }//endif
589 
590 //-------------------------------------------------------------------------
591  }//end routine
592 ///////////////////////////////////////////////////////////////////////////=
593 
594 ///////////////////////////////////////////////////////////////////////////////=
595 ///////////////////////////////////////////////////////////////////////////////c
596 ///////////////////////////////////////////////////////////////////////////////=
597 /** trigger force computation
598  * based on real forces available in each processor's chare then contribute to
599  * global group reduction of all atom forces (within this bead) -> recvContribute
600  **/
601 ///////////////////////////////////////////////////////////////////////////////=
602 void AtomsCompute::startRealSpaceForces(){
603 ///////////////////////////////////////////////////////////////////////////=
604 #ifdef _CP_DEBUG_ATMS_
605  CkPrintf("{%d}[%d] AtomsCompute::startRealSpaceForces\n ", thisInstance.proxyOffset, thisIndex);
606 #endif
607 ///////////////////////////////////////////////////////////////////////////=
608 /// Atom parallelization : same as for integrate
609 
610  int myid = thisIndex;
611  int nproc = nChareAtoms;
612 
613 ///////////////////////////////////////////////////////////////////////////=
614 /// Get the real space atom forces plus these 4 quantities
615 
616  pot_ewd_rs = 0.0;
617  vself = 0.0;
618  vbgr = 0.0;
619  potPerdCorr = 0.0;
620 
621 #ifndef _CP_DEBUG_PSI_OFF_
622 #ifndef _CP_DEBUG_SCALC_ONLY_
623  CPRSPACEION::CP_getionforce(natm,&fastAtoms,myid,nproc,&pot_ewd_rs,&vself,&vbgr,&potPerdCorr, pScratchProxy.ckLocalBranch()->psscratch);
624 #endif
625 #endif
626 
627 ///////////////////////////////////////////////////////////////////////////=
628 /// Everybody contributes to the reduction (see contribute forces comment)
629 
630 #ifdef _CP_DEBUG_ATMS_
631  CkPrintf("GJM_DBG: calling contribute atm forces %d\n",myid);
632 #endif
633  // get the atomCache working to collect the other force contribution for us
634  if(thisIndex==0)
635  UatomsCacheProxy[thisInstance.proxyOffset].contributeforces();
636 
637  double *ftot = new double[(3*natm+2)];
638  for(int i=0,j=0; i<natm; i++,j+=3){
639  ftot[j] = fastAtoms.fx[i];
640  ftot[j+1] = fastAtoms.fy[i];
641  ftot[j+2] = fastAtoms.fz[i];
642  }//endfor
643  ftot[3*natm] =pot_ewd_rs;
644  ftot[3*natm+1]=potPerdCorr;
645  CkCallback cb(CkIndex_AtomsCompute::recvContribute(NULL), UatomsComputeProxy[thisInstance.proxyOffset]);
646  contribute((3*natm+2)*sizeof(double),ftot,CkReduction::sum_double,cb);
647  delete [] ftot;
648 ///////////////////////////////////////////////////////////////////////////=
649  }//end routine
650 ///////////////////////////////////////////////////////////////////////////=
651 
652 
653 ///////////////////////////////////////////////////////////////////////////=
654 /// Atom energy output
655 ///////////////////////////////////////////////////////////////////////////=
656 ///////////////////////////////////////////////////////////////////////////c
657 ///////////////////////////////////////////////////////////////////////////=
658 void AtomsCompute::outputAtmEnergy() {
659 ///////////////////////////////////////////////////////////////////////////=
660 
661  EnergyGroup *eg = UegroupProxy[thisInstance.proxyOffset].ckLocalBranch();
662  CPcharmParaInfo *sim = CPcharmParaInfo::get();
663  int myid = CkMyPe();
664  double eKinetic = eg->estruct.eKinetic_atm;
665  double eKineticNhc = eg->estruct.eKineticNhc_atm;
666  double fmag = eg->estruct.fmag_atm;
667  double pot_ewd_rs_now = eg->estruct.eewald_real;
668  double potNhc = eg->estruct.potNhc_atm;
669  double free_atm = 3*((double)natm);
670  int iperd = sim->iperd;
671 
672  if(myid==0){
673  if(iperd!=0){
674  fprintf(temperScreenFile,"Iter [%d] EWALD_REAL = %5.8lf\n", *iteration, pot_ewd_rs_now);
675  fprintf(temperScreenFile,"Iter [%d] EWALD_SELF = %5.8lf\n",*iteration, vself);
676  fprintf(temperScreenFile,"Iter [%d] EWALD_BGR = %5.8lf\n",*iteration, vbgr);
677  if(iperd!=3){
678  fprintf(temperScreenFile,"Iter [%d] EWALD_Perd = %5.8lf\n",*iteration, potPerdCorr);
679  }//endif
680  }else{
681  fprintf(temperScreenFile,"Iter [%d] ATM_COUL = %5.8lf\n",*iteration, pot_ewd_rs_now);
682  }//endif
683  if(cp_min_opt==0){
684  fprintf(temperScreenFile,"Iter [%d] atm eKin = %5.8lf\n",*iteration, eKinetic);
685  fprintf(temperScreenFile,"Iter [%d] atm Temp = %5.8lf\n",*iteration, (2.0*eKinetic*BOLTZ/free_atm));
686  fprintf(temperScreenFile,"Iter [%d] atm fmag = %5.8lf\n",*iteration, fmag);
687  if(iextended_on==1){
688  double free_Nhc;
689  if(isokin_opt==0){
690  free_Nhc = free_atm*((double)len_nhc);
691  }else{
692  free_Nhc = free_atm*((double)(len_nhc-1));
693  }//endif
694  fprintf(temperScreenFile,"Iter [%d] atm eKinNhc = %5.8lf\n",*iteration, eKineticNhc);
695  fprintf(temperScreenFile,"Iter [%d] atm TempNHC = %5.8lf\n",*iteration, (2.0*eKineticNhc*BOLTZ/free_Nhc));
696  fprintf(temperScreenFile,"Iter [%d] atm potNHC = %5.8lf\n",*iteration, potNhc);
697  }//endif
698  }else{
699  fprintf(temperScreenFile,"Iter [%d] atm fmag = %5.8lf\n",*iteration, fmag);
700  }//endif
701  }//endif
702 
703 //-------------------------------------------------------------------------
704  }//end routine
705 ///////////////////////////////////////////////////////////////////////////=
706 
707 
708 ///////////////////////////////////////////////////////////////////////////=
709 ///////////////////////////////////////////////////////////////////////////c
710 ///////////////////////////////////////////////////////////////////////////=
711 void AtomsCompute::send_PIMD_Fx(){
712 ///////////////////////////////////////////////////////////////////////////=
713 /// CkPrintf("{%d}[%d] AtomsCompute::send_PIMD_fx iteration %d\n ", thisInstance.proxyOffset, CkMyPe(), *iteration);
714 /// every BOC has all Fx might as well just bcast from beadroot
715 
716  if(amBeadRoot){
717  AtomXYZMsg *msg= new (natm, natm, natm, 8*sizeof(int)) AtomXYZMsg;
718  for(int atomI=0;atomI<natm;atomI++){
719  msg->x[atomI]=atoms[atomI].fx;
720  msg->y[atomI]=atoms[atomI].fy;
721  msg->z[atomI]=atoms[atomI].fz;
722  }//endfor
723  msg->index=PIBeadIndex;
724  UPIBeadAtomsProxy[thisInstance.proxyOffset].accept_PIMD_Fx(msg);
725  }else{
726  // everyone else should chill out for Fu
727  }//endif
728 //-------------------------------------------------------------------------
729  }//end routine
730 ///////////////////////////////////////////////////////////////////////////=
731 
732 
733 ///////////////////////////////////////////////////////////////////////////=
734 ///////////////////////////////////////////////////////////////////////////c
735 ///////////////////////////////////////////////////////////////////////////=
736 void AtomsCompute::send_PIMD_Fx_and_x(){
737 ///////////////////////////////////////////////////////////////////////////=
738 /// CkPrintf("{%d}[%d] AtomsCompute::send_PIMD_fx iteration %d\n ", thisInstance.proxyOffset, CkMyPe(), *iteration);
739 /// every BOC has all Fx might as well just bcast from beadroot
740 
741  if(amBeadRoot){
742  AtomXYZMsg *msg= new (2*natm, 2*natm, 2*natm, 8*sizeof(int)) AtomXYZMsg;
743  for(int atomI=0;atomI<natm;atomI++){
744  msg->x[atomI]=atoms[atomI].fx;
745  msg->y[atomI]=atoms[atomI].fy;
746  msg->z[atomI]=atoms[atomI].fz;
747  }//endfor
748  int ioff = natm;
749  for(int atomI=0;atomI<natm;atomI++){
750  msg->x[(atomI+ioff)]=atoms[atomI].x;
751  msg->y[(atomI+ioff)]=atoms[atomI].y;
752  msg->z[(atomI+ioff)]=atoms[atomI].z;
753  }//endfor
754  msg->index=PIBeadIndex;
755  UPIBeadAtomsProxy[thisInstance.proxyOffset].accept_PIMD_Fx_and_x(msg);
756  }else{
757  // everyone else should chill out for Fu
758  }//endif
759 //-------------------------------------------------------------------------
760  }//end routine
761 ///////////////////////////////////////////////////////////////////////////=
762 
763 
764 
765 ///////////////////////////////////////////////////////////////////////////=
766 ///////////////////////////////////////////////////////////////////////////c
767 ///////////////////////////////////////////////////////////////////////////=
768 /**
769  * This thing is sending its updated atoms to all elements.
770  * Effectively an all-to-all implementation of what is basically an
771  * all-gather. Cost of this isn't too bad for the array version if we
772  * use a section proxy for the multicast it constrains the number of
773  * destinations fairly managably. Formerly, this degraded to nAtoms
774  * global broadcasts in the old group implementation. If you've ever
775  * wondered why network fabrics quake in fear at the approach of an
776  * openatom run, this sort of thing should give you an inkling of the
777  * abuse we dish out. The only saving grace is that the number of
778  * atoms is relatively small.
779  */
780 void AtomsCompute::sendAtoms(double eKinetic_loc,double eKineticNhc_loc,double potNhc_loc,
781  double potPIMDChain_loc,int natmNow,int natmStr,int natmEnd){
782 ///////////////////////////////////////////////////////////////////////////=
783 /// Malloc the message
784 /// CkPrintf("{%d}[%d] AtomsCompute::sendAtoms.\n ", thisInstance.proxyOffset, CkMyPe());
785 
786  int nsize = 9*natmNow+4;
787  AtomMsg *msg = new (nsize,8*sizeof(int)) AtomMsg;
788  CkSetQueueing(msg, CK_QUEUEING_IFIFO);
789  *(int*)CkPriorityPtr(msg) = config.sfpriority-10;
790 
791  double *atmData = msg->data;
792  msg->nsize = nsize;
793  msg->natmStr = natmStr;
794  msg->natmEnd = natmEnd;
795 
796 ///////////////////////////////////////////////////////////////////////////=
797 /// pack atom positions : new for use : old for output
798 
799  for(int i=natmStr,j=0;i<natmEnd;i++,j+=9){
800  atmData[(j) ]=atoms[i].x; // for PIMD these are the XU
801  atmData[(j+1)]=atoms[i].y;
802  atmData[(j+2)]=atoms[i].z;
803  atmData[(j+3)]=atoms[i].xold; // for PIMD these are the X
804  atmData[(j+4)]=atoms[i].yold;
805  atmData[(j+5)]=atoms[i].zold;
806  atmData[(j+6)]=atoms[i].vxold; // for PIMD these are the Vxu
807  atmData[(j+7)]=atoms[i].vyold;
808  atmData[(j+8)]=atoms[i].vzold;
809  }//endfor
810 
811 ///////////////////////////////////////////////////////////////////////////=
812 /// pack the 3 energies
813 
814  atmData[(nsize-4)] = potPIMDChain_loc;
815  atmData[(nsize-3)] = eKinetic_loc;
816  atmData[(nsize-2)] = eKineticNhc_loc;
817  atmData[(nsize-1)] = potNhc_loc;
818 
819 ///////////////////////////////////////////////////////////////////////////=
820 /// Send the message to everyone in the group
821 
822 /// FIX THIS: should use delegated section to minimize stupid anytime
823 /// migration overhead. Also, this is really an all-gather.
824 
825  UatomsComputeProxy[thisInstance.proxyOffset].acceptAtoms(msg);
826 
827 //-------------------------------------------------------------------------
828  }//end routine
829 ///////////////////////////////////////////////////////////////////////////=
830 
831 
832 ///////////////////////////////////////////////////////////////////////////=
833 ///////////////////////////////////////////////////////////////////////////c
834 ///////////////////////////////////////////////////////////////////////////=
835 /**
836  * For dynamics we have to update all the caches with the new coordinates
837  */
838 void AtomsCompute::bcastAtomsToAtomCache()
839 {
840 ///////////////////////////////////////////////////////////////////////////=
841 /// Malloc the message
842 /// CkPrintf("{%d}[%d] AtomsCompute::bcastAtomsToAtomCache.\n ", thisInstance.proxyOffset, thisIndex);
843 
844  int nsize = 9*natm+4;
845  AtomMsg *msg = new (nsize,8*sizeof(int)) AtomMsg;
846  CkSetQueueing(msg, CK_QUEUEING_IFIFO);
847  *(int*)CkPriorityPtr(msg) = config.sfpriority-10;
848 
849  double *atmData = msg->data;
850  msg->nsize = nsize;
851  msg->natmStr = 0;
852  msg->natmEnd = natm;
853 
854 ///////////////////////////////////////////////////////////////////////////=
855 /// pack atom positions : new for use : old for output
856 
857  for(int i=0,j=0;i<natm;i++,j+=9){
858  atmData[(j) ]=atoms[i].x; // for PIMD these are the XU
859  atmData[(j+1)]=atoms[i].y;
860  atmData[(j+2)]=atoms[i].z;
861  atmData[(j+3)]=atoms[i].xold; // for PIMD these are the X
862  atmData[(j+4)]=atoms[i].yold;
863  atmData[(j+5)]=atoms[i].zold;
864  atmData[(j+6)]=atoms[i].vxold; // for PIMD these are the Vxu
865  atmData[(j+7)]=atoms[i].vyold;
866  atmData[(j+8)]=atoms[i].vzold;
867  }//endfor
868  EnergyGroup *eg = UegroupProxy[thisInstance.proxyOffset].ckLocalBranch();
869  atmData[nsize-4]=eg->estruct.eKinetic_atm;
870  atmData[nsize-3]=eg->estruct.eKineticNhc_atm;
871  atmData[nsize-2]=eg->estruct.potNhc_atm;
872  atmData[nsize-1]=eg->estruct.potPIMDChain;
873 
874 ///////////////////////////////////////////////////////////////////////////=
875 /// Send the message to everyone in the group
876 
877 /// FIX THIS: should use delegated section to minimize stupid anytime
878 /// migration overhead.
879 
880  UatomsCacheProxy[thisInstance.proxyOffset].acceptAtoms(msg);
881 
882 //-------------------------------------------------------------------------
883  }//end routine
884 ///////////////////////////////////////////////////////////////////////////=
885 
886 
887 ///////////////////////////////////////////////////////////////////////////=
888 ///////////////////////////////////////////////////////////////////////////c
889 ///////////////////////////////////////////////////////////////////////////=
890 /**
891  * Take packed message of a chunk of the atoms with updated positions. Update
892  * local copy of atoms. Update local energyGroup members. Print atom energies
893  * when we have all of them. Do file output of atoms if desired.
894  */
895 ///////////////////////////////////////////////////////////////////////////=
896  void AtomsCompute::acceptAtoms(AtomMsg *msg) {
897 ///////////////////////////////////////////////////////////////////////////=
898 /// CkPrintf("{%d}[%d] AtomsCompute::acceptAtoms.\n ", thisInstance.proxyOffset, CkMyPe());
899  double *atmData = msg->data;
900  int nsize = msg->nsize;
901  int natmStr = msg->natmStr;
902  int natmEnd = msg->natmEnd;
903 
904 ///////////////////////////////////////////////////////////////////////////=
905 /// unpack atom position and velocity : For PIMD x is reconstructed from xu
906 
907  if(numPIMDBeads>1){
908  for(int i=natmStr,j=0;i<natmEnd;i++,j+=9){
909  atoms[i].xu = atmData[(j) ];
910  atoms[i].yu = atmData[(j+1)];
911  atoms[i].zu = atmData[(j+2)];
912  atoms[i].xold = atmData[(j+3)];
913  atoms[i].yold = atmData[(j+4)];
914  atoms[i].zold = atmData[(j+5)];
915  atoms[i].vxold = atmData[(j+6)];
916  atoms[i].vyold = atmData[(j+7)];
917  atoms[i].vzold = atmData[(j+8)];
918  }//endfor
919  }else{
920  for(int i=natmStr,j=0;i<natmEnd;i++,j+=9){
921  atoms[i].x = atmData[(j) ];
922  atoms[i].y = atmData[(j+1)];
923  atoms[i].z = atmData[(j+2)];
924  atoms[i].xold = atmData[(j+3)];
925  atoms[i].yold = atmData[(j+4)];
926  atoms[i].zold = atmData[(j+5)];
927  atoms[i].vxold = atmData[(j+6)];
928  atoms[i].vyold = atmData[(j+7)];
929  atoms[i].vzold = atmData[(j+8)];
930  }//endfor
931  }//endif
932 #ifdef _CP_DEBUG_ATMS_
933  if(CkMyPe()==0)
934  {
935  for(int i=natmStr,j=0;i<natmEnd;i++,j+=9){
936 
937  CkPrintf("AtomCompute acceptAtoms[%d] updated to %.5g,%.5g,%.5g from %.5g,%.5g,%.5g\n",i, atoms[i].x, atoms[i].y, atoms[i].z, atoms[i].xold, atoms[i].yold, atoms[i].zold);
938  }
939  }
940 #endif
941 ///////////////////////////////////////////////////////////////////////////=
942 /// unpack energy
943 
944  countAtm++;
945  if(countAtm==1){
946  eKinetic = 0;
947  eKineticNhc = 0;
948  potNhc = 0;
949  potPIMDChain = 0;
950  }//endif
951 
952  potPIMDChain+= atmData[(nsize-4)];
953  eKinetic += atmData[(nsize-3)];
954  eKineticNhc += atmData[(nsize-2)];
955  potNhc += atmData[(nsize-1)];
956 
957 ///////////////////////////////////////////////////////////////////////////=
958 /// Delete the message
959 
960  delete msg;
961 
962 ///////////////////////////////////////////////////////////////////////////=
963 /// Copy to the fast vectors and phone home : If PIMD Transform U to X
964 
965  if(countAtm==nAtmMsgRecv){
966  countAtm = 0;
967 
968  //---------------------------------------------------------------------------------
969  // Do these quantities need to be updated on energy group chares that are not co-resident with an AtomsCompute?
970 
971  EnergyGroup *eg = UegroupProxy[thisInstance.proxyOffset].ckLocalBranch();
972  eg->estruct.eKinetic_atm = eKinetic;
973  eg->estruct.eKineticNhc_atm = eKineticNhc;
974  eg->estruct.potNhc_atm = potNhc;
975  eg->estruct.potPIMDChain = potPIMDChain;
976 
977  //---------------------------------------------------------------------------------
978  copySlowToFast(); // not complete for PIMD as you have xu not x but that's OK
979  outputAtmEnergy();
980 
981  //---------------------------------------------------------------------------------
982  // Output : Iteration is time of atoms[i].xold
983  // maxIter is 1 more than you need, slightly annoying but livable.
984 
985  int output_on = config.atmOutput;
986  if(output_on==1 && *iteration<=config.maxIter-1){
987  int pi_beads = 1;
988  int iwrite_atm = 0;
989  int myid = thisIndex;
990  ATOMOUTPUT::ctrl_piny_output(*iteration,natm,len_nhc,pi_beads,myid,atoms,atomsNHC,
991  &iwrite_atm,output_on,TemperIndex,PIBeadIndex);
992  if(myid==0 && iwrite_atm>0){
993  fprintf(temperScreenFile,"-----------------------------------\n");
994  fprintf(temperScreenFile, "Writing atoms to disk at time %d\n",*iteration);
995  fprintf(temperScreenFile,"-----------------------------------\n");
996  }//endif
997  }//endif
998 
999  //---------------------------------------------------------------------------------
1000  // Transform PIMD U to X
1001 
1002  if(numPIMDBeads>1){
1003  // CkPrintf("{%d}[%d] AtomComputep::acceptAtoms. numPIMDBeads >1 transform PIMD U to X iteration %d\n ",
1004  // thisInstance.proxyOffset, CkMyPe(), *iteration);
1005  if(amBeadRoot){send_PIMD_u();}
1006  if(amBeadRoot && amZerothBead){//For staging, this is the 1st bead not the CM but that's fine
1007  AtomXYZMsg *msg = new (natm,natm,natm) AtomXYZMsg;
1008  for(int atomI=0;atomI<natm;atomI++){
1009  msg->x[atomI]=atoms[atomI].xu;
1010  msg->y[atomI]=atoms[atomI].yu;
1011  msg->z[atomI]=atoms[atomI].zu;
1012  }//endfor
1013  // proxyHeadBeads.accept_PIMD_CM(msg);
1014  proxyAllBeads.accept_PIMD_CM(msg);
1015  }//endif : I am King of the Beads
1016  }else{
1017  //non PIMD case we're done
1018  int i=0;
1019  CkCallback cb(CkIndex_AtomsCompute::atomsDone(NULL), CkArrayIndex1D(0), UatomsComputeProxy[thisInstance.proxyOffset]);
1020  contribute(sizeof(int),&i,CkReduction::sum_int,cb);
1021  }//endif :
1022 
1023  }//endif : I have received all my messages
1024 
1025 //-------------------------------------------------------------------------
1026  }//end routine
1027 ///////////////////////////////////////////////////////////////////////////=
1028 
1029 void AtomsCompute::atomsDone(CkReductionMsg *msg)
1030 {
1031  // This means that all the atomCompute chares are done integrating
1032  // in minimization. Tell the atomCaches they can go ahead now.
1033  delete msg;
1034  if(cp_wave_opt==0 && cp_min_opt==0)
1035  bcastAtomsToAtomCache();
1036  else
1037  UatomsCacheProxy[thisInstance.proxyOffset].atomsDone();
1038 }
1039 ///////////////////////////////////////////////////////////////////////////////=
1040 
1041 
1042 ///////////////////////////////////////////////////////////////////////////////=
1043 ///////////////////////////////////////////////////////////////////////////////c
1044 ///////////////////////////////////////////////////////////////////////////////=
1045 void AtomsCompute::accept_PIMD_CM(AtomXYZMsg *msg){
1046 ///////////////////////////////////////////////////////////////////////////////=
1047 
1048  for(int atomnum=0; atomnum<natm; atomnum++){
1049  atoms[atomnum].xcm = msg->x[atomnum];
1050  atoms[atomnum].ycm = msg->y[atomnum];
1051  atoms[atomnum].zcm = msg->z[atomnum];
1052  }//endfor
1053 
1054  delete msg;
1055  atomsCMrecv=true;
1056 
1057  if(atomsPIMDXrecv){
1058 #ifdef _CP_DEBUG_ATMS_
1059  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_CM contributing to atomsDone atomsPIMDXrecv is %d iteration %d\n ", thisInstance.proxyOffset, thisIndex, atomsPIMDXrecv, *iteration);
1060 #endif
1061  int i=0;
1062  CkCallback cb(CkIndex_AtomsCompute::atomsDone(NULL), CkArrayIndex1D(0), UatomsComputeProxy[thisInstance.proxyOffset]);
1063  contribute(sizeof(int),&i,CkReduction::sum_int,cb);
1064  atomsCMrecv=atomsPIMDXrecv=false;
1065  }else{
1066 #ifdef _CP_DEBUG_ATMS_
1067  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_CM warning! atomsPIMDXrecv is %d iteration %d\n ", thisInstance.proxyOffset, thisIndex, atomsPIMDXrecv, *iteration);
1068 #endif
1069  }//endif
1070 
1071 ///////////////////////////////////////////////////////////////////////////////=
1072  }//end routine
1073 ///////////////////////////////////////////////////////////////////////////////=
1074 
1075 
1076 ///////////////////////////////////////////////////////////////////////////////=
1077 ///////////////////////////////////////////////////////////////////////////////c
1078 ///////////////////////////////////////////////////////////////////////////////=
1079 void AtomsCompute::send_PIMD_u(){
1080 ///////////////////////////////////////////////////////////////////////////////=
1081 /// CkPrintf("{%d}[%d] AtomsCompute::send_PIMD_u iteration %d\n ", thisInstance.proxyOffset, CkMyPe(), *iteration);
1082 
1083  for(int atomnum=0;atomnum<natm;atomnum++){
1084  UPIBeadAtomsProxy[thisInstance.proxyOffset][atomnum].accept_PIMD_u(
1085  atoms[atomnum].xu,atoms[atomnum].yu,atoms[atomnum].zu, PIBeadIndex);
1086  }//endfor
1087 
1088 }//end routine
1089 ///////////////////////////////////////////////////////////////////////////////=
1090 
1091 
1092 ///////////////////////////////////////////////////////////////////////////////=
1093 ///////////////////////////////////////////////////////////////////////////////c
1094 ///////////////////////////////////////////////////////////////////////////////=
1095 /// is broadcast to us
1096 void AtomsCompute::accept_PIMD_Fu(double _fxu, double _fyu, double _fzu, int atomI){
1097 ///////////////////////////////////////////////////////////////////////////////=
1098 
1099  atoms[atomI].fxu =_fxu; // FastAtoms not used for integration or output
1100  atoms[atomI].fyu =_fyu;
1101  atoms[atomI].fzu =_fzu;
1102 
1103  acceptCountfu++;
1104 #ifdef _CP_DEBUG_ATMS_
1105  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_fu (%d of %d) iteration %d %d %.5g %.5g %.5g\n",
1106  thisInstance.proxyOffset, thisIndex,acceptCountfu, natm, *iteration, atomI, _fxu, _fyu, _fzu);
1107 #endif
1108  if(acceptCountfu==natm){
1109  // CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_fu done calling integrator iteration %d\n ",
1110  // thisInstance.proxyOffset, CkMyPe(), *iteration);
1111  integrateAtoms();
1112  acceptCountfu=0;
1113  }//endif
1114 
1115 ///////////////////////////////////////////////////////////////////////////////=
1116  }//end routine
1117 ///////////////////////////////////////////////////////////////////////////////=
1118 
1119 ///////////////////////////////////////////////////////////////////////////////=
1120 ///////////////////////////////////////////////////////////////////////////////c
1121 ///////////////////////////////////////////////////////////////////////////////=
1122 /// is broadcast to us
1123 void AtomsCompute::accept_PIMD_Fu_and_u(double _fxu, double _fyu, double _fzu,
1124  double _xu, double _yu, double _zu, int atomI){
1125 ///////////////////////////////////////////////////////////////////////////////=
1126 
1127  atoms[atomI].fxu =_fxu; // FastAtoms not used for integration or output
1128  atoms[atomI].fyu =_fyu;
1129  atoms[atomI].fzu =_fzu;
1130 
1131  atoms[atomI].xu =_xu; // FastAtoms not used for integration or output
1132  atoms[atomI].yu =_yu;
1133  atoms[atomI].zu =_zu;
1134 
1135  acceptCountfu++;
1136 #ifdef _CP_DEBUG_ATMS_
1137  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_fu (%d of %d) iteration %d %d %.5g %.5g %.5g\n",
1138  thisInstance.proxyOffset, thisIndex,acceptCountfu, natm, *iteration, atomI, _fxu, _fyu, _fzu);
1139 #endif
1140  if(acceptCountfu==natm){
1141  // CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_fu done calling integrator iteration %d\n ",
1142  // thisInstance.proxyOffset, CkMyPe(), *iteration);
1143  integrateAtoms();
1144  acceptCountfu=0;
1145  }//endif
1146 
1147 ///////////////////////////////////////////////////////////////////////////////=
1148  }//end routine
1149 ///////////////////////////////////////////////////////////////////////////////=
1150 
1151 
1152 ///////////////////////////////////////////////////////////////////////////////=
1153 ///////////////////////////////////////////////////////////////////////////////c
1154 ///////////////////////////////////////////////////////////////////////////////=
1155 void AtomsCompute::accept_PIMD_x(double _x, double _y, double _z, int atomI){
1156 ///////////////////////////////////////////////////////////////////////////////=
1157 /// CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_x iteration %d\n ", thisInstance.proxyOffset, thisIndex, *iteration);
1158 ///////////////////////////////////////////////////////////////////////////////=
1159 
1160  atoms[atomI].x=_x;
1161  atoms[atomI].y=_y;
1162  atoms[atomI].z=_z;
1163 
1164  fastAtoms.x[atomI]=_x;
1165  fastAtoms.y[atomI]=_y;
1166  fastAtoms.z[atomI]=_z;
1167 
1168  acceptCountX++;
1169 
1170  if(acceptCountX==natm){
1171  acceptCountX=0;
1172  atomsPIMDXrecv=true;
1173  if(atomsCMrecv){
1174  int i=0;
1175 #ifdef _CP_DEBUG_ATMS_
1176  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_x contributing to atomsDone atomsCMrecv is %d iteration %d\n ", thisInstance.proxyOffset, thisIndex, atomsCMrecv, *iteration);
1177 #endif
1178  CkCallback cb(CkIndex_AtomsCompute::atomsDone(NULL), CkArrayIndex1D(0), UatomsComputeProxy[thisInstance.proxyOffset]);
1179 
1180  contribute(sizeof(int),&i,CkReduction::sum_int,cb);
1181  atomsCMrecv=atomsPIMDXrecv=false;
1182  }else{
1183 #ifdef _CP_DEBUG_ATMS_
1184  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_x warning! atomsCMrecv is %d iteration %d\n", thisInstance.proxyOffset, thisIndex, atomsCMrecv, *iteration);
1185 #endif
1186  }//endif
1187  }//endif : all data has arrived
1188 
1189 ///////////////////////////////////////////////////////////////////////////////=
1190  }//end routine
1191 ///////////////////////////////////////////////////////////////////////////////=
1192 
1193 
1194 ///////////////////////////////////////////////////////////////////////////////=
1195 ///////////////////////////////////////////////////////////////////////////////c
1196 ///////////////////////////////////////////////////////////////////////////////=
1197 /// done during initialization in 1st iteration
1198 ///////////////////////////////////////////////////////////////////////////////=
1199 void AtomsCompute::send_PIMD_x(){
1200 ///////////////////////////////////////////////////////////////////////////////=
1201 
1202  for(int atomnum=0;atomnum<natm;atomnum++){
1203  UPIBeadAtomsProxy[thisInstance.proxyOffset][atomnum].accept_PIMD_x(atoms[atomnum].x,atoms[atomnum].y,atoms[atomnum].z,
1204  PIBeadIndex);
1205  }//endfor
1206 
1207 ///////////////////////////////////////////////////////////////////////////////=
1208  }//endroutine
1209 ///////////////////////////////////////////////////////////////////////////////=
1210 
1211 void AtomsCompute::acceptNewTemperature(double temp)
1212 {
1213  // Hey GLENN do something with your new temperature here
1214  // when you're done
1215  int i=1;
1216  contribute(sizeof(int), &i, CkReduction::sum_int,
1217  CkCallback(CkIndex_InstanceController::atomsDoneNewTemp(NULL),CkArrayIndex1D(thisInstance.proxyOffset),instControllerProxy), thisInstance.proxyOffset);
1218 ///////////////////////////////////////////////////////////////////////////////=
1219  }//end routine
1220 ///////////////////////////////////////////////////////////////////////////////=
1221 
1222 ///////////////////////////////////////////////////////////////////////////////=
1223 ///////////////////////////////////////////////////////////////////////////////c
1224 ///////////////////////////////////////////////////////////////////////////////=
1225 /// done during initialization in 1st iteration
1226 ///////////////////////////////////////////////////////////////////////////////=
1227 void AtomsCompute::accept_PIMD_u(double _xu, double _yu, double _zu, int atomI){
1228 ///////////////////////////////////////////////////////////////////////////////=
1229 #ifdef _CP_DEBUG_ATMS_
1230  CkPrintf("{%d}[%d] AtomsCompute::accept_PIMD_u iteration %d\n", thisInstance.proxyOffset, thisIndex, *iteration);
1231 #endif
1232  atoms[atomI].xu =_xu;
1233  atoms[atomI].yu =_yu;
1234  atoms[atomI].zu =_zu;
1235 
1236  acceptCountu++;
1237  if(acceptCountu==natm){
1238  integrateAtoms();
1239  acceptCountu=0;
1240  }//endif
1241 
1242 ///////////////////////////////////////////////////////////////////////////////=
1243  }//end routine
1244 ///////////////////////////////////////////////////////////////////////////////=
1245 /*@}*/
1246 #include "Atoms.def.h"
1247 
AtomsCompute::outputAtmEnergy
void outputAtmEnergy()
= Atom energy output
Definition: AtomsCompute.C:658
UberCollection
holds the UberIndex and the offset for proxies
Definition: Uber.h:68
AtomsCompute::recvContributeForces
void recvContributeForces(CkReductionMsg *)
recvContributeForces Every Array member has all the forces at present.
Definition: AtomsCompute.C:294
AtomsCompute::accept_PIMD_Fu_and_u
void accept_PIMD_Fu_and_u(double _fxu, double _fyu, double _fzu, double _xu, double _yu, double _zu, int atomI)
= is broadcast to us
Definition: AtomsCompute.C:1123
AtomsCompute.h
Accepts reduction of forces from AtomsCache, integrates forces to produce new coordinates, distributes new coordinates to AtomsCache.
AtomsCompute::send_PIMD_u
void send_PIMD_u()
Definition: AtomsCompute.C:1079
instControllerProxy
CProxy_InstanceController instControllerProxy
Definition: cpaimd.ci:236
AtomsCompute::send_PIMD_Fx_and_x
void send_PIMD_Fx_and_x()
Definition: AtomsCompute.C:736
AtomsCompute::acceptAtoms
void acceptAtoms(AtomMsg *)
Take packed message of a chunk of the atoms with updated positions.
Definition: AtomsCompute.C:896
AtomsCompute::startRealSpaceForces
void startRealSpaceForces()
trigger force computation based on real forces available in each processor's chare then contribute to...
Definition: AtomsCompute.C:602
fftCacheSlab.h
Add type declarations for simulationConstants class (readonly vars) and once class for each type of o...
IntMap.h
Author: Eric J Bohm Date Created: June 4th, 2006.
AtomsCompute::accept_PIMD_u
void accept_PIMD_u(double _ux, double _uy, double _uz, int atomI)
= done during initialization in 1st iteration
Definition: AtomsCompute.C:1227
AtomsCompute::recvContribute
void recvContribute(CkReductionMsg *)
recvContribute Every Array member is sent all the forces at present.
Definition: AtomsCompute.C:250
AtomsCompute::~AtomsCompute
~AtomsCompute()
Destructor.
Definition: AtomsCompute.C:167
AtomsCompute::integrateAtoms
void integrateAtoms()
Integrate atoms.
Definition: AtomsCompute.C:420
AtomNHC
Definition: Atoms.h:129
Atom
Definition: Atoms.h:20
AtomsCompute::send_PIMD_x
void send_PIMD_x()
= done during initialization in 1st iteration
Definition: AtomsCompute.C:1199
AtomsCompute::sendAtoms
void sendAtoms(double, double, double, double, int, int, int)
This thing is sending its updated atoms to all elements.
Definition: AtomsCompute.C:780
AtomsCompute::accept_PIMD_Fu
void accept_PIMD_Fu(double _fxu, double _fyu, double _fzu, int atomI)
= is broadcast to us
Definition: AtomsCompute.C:1096
AtomsCompute::handleForces
void handleForces()
Definition: AtomsCompute.C:317
AtomsCompute::accept_PIMD_x
void accept_PIMD_x(double _x, double _y, double _z, int atomI)
Definition: AtomsCompute.C:1155
AtomsCompute::init
void init()
= Bead root is the 0th element.
Definition: AtomsCompute.C:197
AtomsCompute::send_PIMD_Fx
void send_PIMD_Fx()
Definition: AtomsCompute.C:711
EnergyGroup
EnergyGroup class.
Definition: energyGroup.h:14
AtomsCompute::bcastAtomsToAtomCache
void bcastAtomsToAtomCache()
For dynamics we have to update all the caches with the new coordinates.
Definition: AtomsCompute.C:838