NAMD
Controller.C
Go to the documentation of this file.
1 
7 /*****************************************************************************
8  * $Source: /home/cvs/namd/cvsroot/namd2/src/Controller.C,v $
9  * $Author: jim $
10  * $Date: 2017/03/29 21:29:03 $
11  * $Revision: 1.1326 $
12  *****************************************************************************/
13 
14 #include "InfoStream.h"
15 #include "converse.h"
16 #include "memusage.h"
17 #include "Node.h"
18 #include "Molecule.h"
19 #include "SimParameters.h"
20 #include "Controller.h"
21 #include "ReductionMgr.h"
22 #include "CollectionMaster.h"
23 #include "Output.h"
24 #include "strlib.h"
25 #include "BroadcastObject.h"
26 #include "NamdState.h"
27 #include "ScriptTcl.h"
28 #include "Broadcasts.h"
29 #include "LdbCoordinator.h"
30 #include "Thread.h"
31 #include <math.h>
32 #include <signal.h>
33 #include "NamdOneTools.h"
34 #include "PatchMap.h"
35 #include "PatchMap.inl"
36 #include "Random.h"
37 #include "imd.h"
38 #include "IMDOutput.h"
39 #include "BackEnd.h"
40 #include <fstream>
41 #include <iomanip>
42 #include <errno.h>
43 #include "qd.h"
44 
45 #if defined(USE_GETTIMEOFDAY)
46 #include <sys/time.h>
47 // Replace CmiWallTimer with our own based on gettimeofday
48 static double namdWallTimer() {
49  struct timeval t;
50  gettimeofday(&t, NULL);
51  double sec = double(t.tv_sec); // get absolute seconds
52  double usec = double (t.tv_usec); // plus microseconds
53  double totalSec = sec + 1e-6 * usec;
54  return totalSec;
55 }
56 #else
57 static double namdWallTimer() {
58  return CmiWallTimer();
59 }
60 #endif
61 
62 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
63 extern "C" void CApplicationDepositNode0Data(char *);
64 #endif
65 
66 #ifndef cbrt
67  // cbrt() not in math.h on goneril
68  #define cbrt(x) pow(x,(double)(1.0/3.0))
69 #endif
70 
71 //#define DEBUG_PRESSURE
72 #define MIN_DEBUG_LEVEL 3
73 //#define DEBUGM
74 #include "Debug.h"
75 
76 #define XXXBIGREAL 1.0e32
77 
78 //Modifications for alchemical fep
79 static char *FEPTITLE(int X)
80 {
81  static char tmp_string[21];
82  sprintf(tmp_string, "FepEnergy: %6d ",X);
83  return tmp_string;
84 }
85 
86 static char *FEPTITLE_BACK(int X)
87 {
88  static char tmp_string[21];
89  sprintf(tmp_string, "FepE_back: %6d ",X);
90  return tmp_string;
91 }
92 
93 static char *FEPTITLE2(int X)
94 {
95  static char tmp_string[21];
96  sprintf(tmp_string, "FEP: %7d",X);
97  return tmp_string;
98 }
99 
100 static char *TITITLE(int X)
101 {
102  static char tmp_string[21];
103  sprintf(tmp_string, "TI: %7d",X);
104  return tmp_string;
105 }
106 //fepe
107 
109 private:
110  RequireReduction *reduction;
111  const int nslabs;
112  const int freq;
113  int nelements;
114  int count;
115  char *name;
116 
117  BigReal *average;
118 
119 public:
120  PressureProfileReduction(int rtag, int numslabs, int numpartitions,
121  const char *myname, int outputfreq)
122  : nslabs(numslabs), freq(outputfreq) {
123  name = strdup(myname);
124  nelements = 3*nslabs * numpartitions;
125  reduction = ReductionMgr::Object()->willRequire(rtag,nelements);
126 
127  average = new BigReal[nelements];
128  count = 0;
129  }
131  delete [] average;
132  delete reduction;
133  free(name);
134  }
135  //
136  void getData(int firsttimestep, int step, const Lattice &lattice,
137  BigReal *total) {
138  reduction->require();
139 
140  int i;
141  double inv_volume = 1.0 / lattice.volume();
142  // accumulate the pressure profile computed for this step into the average.
143  int arraysize = 3*nslabs;
144  for (i=0; i<nelements; i++) {
145  BigReal val = reduction->item(i) * inv_volume;
146  average[i] += val;
147  total[i % arraysize] += val;
148  }
149  count++;
150  if (!(step % freq)) {
151  // convert NAMD internal virial to pressure in units of bar
152  BigReal scalefac = PRESSUREFACTOR * nslabs;
153 
154  iout << "PPROFILE" << name << ": " << step << " ";
155  if (step == firsttimestep) {
156  // output pressure profile for this step
157  for (i=0; i<nelements; i++)
158  iout << reduction->item(i)*scalefac*inv_volume << " ";
159  } else {
160  // output pressure profile averaged over the last Freq steps.
161  scalefac /= count;
162  for (i=0; i<nelements; i++)
163  iout << average[i]*scalefac << " ";
164  }
165  iout << "\n" << endi;
166  // Clear the average for the next block
167  memset(average, 0, nelements*sizeof(BigReal));
168  count = 0;
169  }
170  }
171 };
172 
174  CthThread thread;
175 };
176 
178  computeChecksum(0), marginViolations(0), pairlistWarnings(0),
179  simParams(Node::Object()->simParameters),
180  state(s),
181  collection(CollectionMaster::Object()),
182  startCTime(0),
183  firstCTime(CmiTimer()),
184  startWTime(0),
185  firstWTime(namdWallTimer()),
186  startBenchTime(0),
187  stepInFullRun(0),
188  ldbSteps(0),
189  fflush_count(3)
190 {
193  // for accelMD
194  if (simParams->accelMDOn) {
195  // REDUCTIONS_BASIC wil contain all potential energies and arrive first
197  // contents of amd_reduction be submitted to REDUCTIONS_AMD
199  // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
201  } else {
202  amd_reduction = NULL;
203  submit_reduction = NULL;
205  // do I need to set nodeReduction as well?
206  }
207  // pressure profile reductions
210  ppbonded = ppnonbonded = ppint = NULL;
212  int ntypes = simParams->pressureProfileAtomTypes;
214  // number of partitions for pairwise interactions
215  int npairs = (ntypes * (ntypes+1))/2;
216  pressureProfileAverage = new BigReal[3*nslabs];
217  memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
218  int freq = simParams->pressureProfileFreq;
221  nslabs, npairs, "BONDED", freq);
223  nslabs, npairs, "NONBONDED", freq);
224  // internal partitions by atom type, but not in a pairwise fashion
226  nslabs, ntypes, "INTERNAL", freq);
227  } else {
228  // just doing Ewald, so only need nonbonded
230  nslabs, npairs, "NONBONDED", freq);
231  }
232  }
234  random->split(0,PatchMap::Object()->numPatches()+1);
235 
236  heat = totalEnergy0 = 0.0;
239  stochRescale_count = 0;
240  if (simParams->stochRescaleOn) {
243  }
245  // strainRate tensor is symmetric to avoid rotation
248  if ( ! simParams->useFlexibleCell ) {
249  BigReal avg = trace(langevinPiston_strainRate) / 3.;
251  } else if ( simParams->useConstantRatio ) {
252 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
253  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
255 #undef AVGXY
256  }
258 
259  // Monte Carlo pressure control
265  for(int ax = 0; ax < MC_AXIS_TOTAL; ++ax) {
266  mc_trial[ax] = 0;
267  mc_accept[ax] = 0;
268  }
269  }
270 
271  if (simParams->multigratorOn) {
272  multigratorXi = 0.0;
275  Node *node = Node::Object();
276  Molecule *molecule = node->molecule;
277  BigReal Nf = molecule->num_deg_freedom();
279  multigratorNu.resize(n);
280  multigratorNuT.resize(n);
281  multigratorZeta.resize(n);
282  multigratorOmega.resize(n);
283  for (int i=0;i < n;i++) {
284  multigratorNu[i] = 0.0;
285  multigratorZeta[i] = 0.0;
286  if (i == 0) {
287  multigratorOmega[i] = Nf*kT0*tau*tau;
288  } else {
289  multigratorOmega[i] = kT0*tau*tau;
290  }
291  }
293  } else {
294  multigratorReduction = NULL;
295  }
296  origLattice = state->lattice;
298  temp_avg = 0;
299  pressure_avg = 0;
300  groupPressure_avg = 0;
301  avg_count = 0;
302  pressure_tavg = 0;
303  groupPressure_tavg = 0;
304  tavg_count = 0;
305  checkpoint_stored = 0;
306  drudeBondTemp = 0;
307  drudeBondTempAvg = 0;
308  cumAlchWork = 0;
309 
310  bondedEnergy_ti_1 = 0;
311  bondedEnergy_ti_2 = 0;
312  ljEnergy_ti_1 = 0;
313  ljEnergy_ti_2 = 0;
314  electEnergy_ti_1 = 0;
315  electEnergy_ti_2 = 0;
318  TIderivative = 0;
319 }
320 
322 {
323  delete broadcast;
324  delete reduction;
325  delete min_reduction;
326  delete amd_reduction;
327  delete submit_reduction;
328  delete ppbonded;
329  delete ppnonbonded;
330  delete ppint;
331  delete [] pressureProfileAverage;
332  delete random;
334 }
335 
336 void Controller::threadRun(Controller* arg)
337 {
338  arg->algorithm();
339 }
340 
341 void Controller::run(void)
342 {
343  // create a Thread and invoke it
344  DebugM(4, "Starting thread in controller on this=" << this << "\n");
345  threadWrapper = new CthThreadWrapper;
346  threadWrapper->thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
347  CthSetStrategyDefault(threadWrapper->thread);
348 #if CMK_BLUEGENE_CHARM
349  BgAttach(thread);
350 #endif
351  awaken();
352 }
353 
354 void Controller::awaken(void) {
355  CthAwaken(threadWrapper->thread);
356 }
357 
359  CthSuspend();
360 }
361 
363 {
364  int scriptTask;
365  int scriptSeq = 0;
366  BackEnd::awaken();
367  while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
368 #ifdef NODEGROUP_FORCE_REGISTER
369  CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
370  PatchData *patchData = cpdata.ckLocalBranch();
371 #endif
372  switch ( scriptTask ) {
373  case SCRIPT_OUTPUT:
376  break;
377  case SCRIPT_FORCEOUTPUT:
379  break;
380  case SCRIPT_MEASURE:
382  break;
383  case SCRIPT_REINITVELS:
384  iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
385  << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
386  break;
387  case SCRIPT_RESCALEVELS:
388  iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
389  << " BY " << simParams->scriptArg1 << "\n" << endi;
390  break;
392  // Parameter setting already reported in ScriptTcl
393  // Nothing to do!
394  break;
395  case SCRIPT_CHECKPOINT:
396  iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
397  << "\n" << endi;
398  checkpoint_stored = 1;
399  checkpoint_lattice = state->lattice;
401  break;
402  case SCRIPT_REVERT:
403  iout << "REVERTING AT STEP " << simParams->firstTimestep
404  << "\n" << endi;
405  if ( ! checkpoint_stored )
406  NAMD_die("Unable to revert, checkpoint was never called!");
407  state->lattice = checkpoint_lattice;
409  break;
411  iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
412  << " TO KEY " << simParams->scriptStringArg1;
413  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
414  iout << "\n" << endi;
415  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
416  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
418  }
420  cp.lattice = state->lattice;
421  cp.state = *(ControllerState*)this;
422  } else {
424  scriptTask, state->lattice, *(ControllerState*)this);
425  }
426  break;
428  iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
429  << " FROM KEY " << simParams->scriptStringArg1;
430  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
431  iout << "\n" << endi;
432  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
433  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
434  NAMD_die("Unable to load checkpoint, requested key was never stored.");
435  }
437  state->lattice = cp.lattice;
438  *(ControllerState*)this = cp.state;
439  } else {
441  scriptTask, state->lattice, *(ControllerState*)this);
442  }
443  break;
445  iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
446  << " FROM KEY " << simParams->scriptStringArg1;
447  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
448  iout << "\n" << endi;
449  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
450  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
451  NAMD_die("Unable to swap checkpoint, requested key was never stored.");
452  }
454  std::swap(state->lattice,cp.lattice);
455  std::swap(*(ControllerState*)this,cp.state);
456  } else {
458  scriptTask, state->lattice, *(ControllerState*)this);
459  }
460  break;
462  iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
463  << " FROM KEY " << simParams->scriptStringArg1;
464  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
465  iout << "\n" << endi;
466  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
467  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
468  NAMD_die("Unable to free checkpoint, requested key was never stored.");
469  }
472  } else {
474  scriptTask, state->lattice, *(ControllerState*)this);
475  }
476  break;
477  case SCRIPT_ATOMSENDRECV:
478  case SCRIPT_ATOMSEND:
479  case SCRIPT_ATOMRECV:
480  break;
481  case SCRIPT_MINIMIZE:
482  minimize();
483  break;
484  case SCRIPT_RUN:
485  case SCRIPT_CONTINUE:
486  integrate(scriptTask);
487  break;
488  default:
489  NAMD_bug("Unknown task in Controller::algorithm");
490  }
491  BackEnd::awaken();
492  }
495  terminate();
496 }
497 
498 
499 //
500 // XXX static and global variables are unsafe for shared memory builds.
501 // The use of global and static vars should be eliminated.
502 //
503 extern int eventEndOfTimeStep;
504 
505 // Handle SIGINT so that restart files get written completely.
506 static int gotsigint = 0;
507 static void my_sigint_handler(int sig) {
508  if (sig == SIGINT) gotsigint = 1;
509 }
510 extern "C" {
511  typedef void (*namd_sighandler_t)(int);
512 }
513 
514 void Controller::integrate(int scriptTask) {
515  char traceNote[24];
516 
517  int step = simParams->firstTimestep;
518 
519  const int numberOfSteps = simParams->N;
520  const int stepsPerCycle = simParams->stepsPerCycle;
521 
522  const int zeroMomentum = simParams->zeroMomentum;
523 
525  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
526  if (dofull)
528  else
530  if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
531 
533 
534  if ( scriptTask == SCRIPT_RUN ) {
535 
536  reassignVelocities(step); // only for full-step velecities
537  rescaleaccelMD(step);
538  adaptTempUpdate(step); // Init adaptive tempering;
539 
540  receivePressure(step);
541  if ( zeroMomentum && dofull && ! (step % slowFreq) )
542  correctMomentum(step);
543  printFepMessage(step);
544  printTiMessage(step);
545  printDynamicsEnergies(step);
546  outputFepEnergy(step);
547  outputTiEnergy(step);
548  if(traceIsOn()){
549  traceUserEvent(eventEndOfTimeStep);
550  sprintf(traceNote, "s:%d", step);
551  traceUserSuppliedNote(traceNote);
552  }
553  outputExtendedSystem(step);
554  rebalanceLoad(step);
555 
556  }
557 
558  // Handling SIGINT doesn't seem to be working on Lemieux, and it
559  // sometimes causes the net-xxx versions of NAMD to segfault on exit,
560  // so disable it for now.
561  // namd_sighandler_t oldhandler = signal(SIGINT,
562  // (namd_sighandler_t)my_sigint_handler);
563  for ( ++step ; step <= numberOfSteps; ++step )
564  {
565  adaptTempUpdate(step);
566  rescaleVelocities(step);
567  tcoupleVelocities(step);
569  berendsenPressure(step);
570  langevinPiston1(step);
571  rescaleaccelMD(step);
572  enqueueCollections(step); // after lattice scaling!
573  receivePressure(step);
574  if ( zeroMomentum && dofull && ! (step % slowFreq) )
575  correctMomentum(step);
576  langevinPiston2(step);
577  reassignVelocities(step);
578 
579  multigratorTemperature(step, 1);
580  multigratorPressure(step, 1);
581  multigratorPressure(step, 2);
582  multigratorTemperature(step, 2);
583 
584  printDynamicsEnergies(step);
585  outputFepEnergy(step);
586  outputTiEnergy(step);
587  if(traceIsOn()){
588  traceUserEvent(eventEndOfTimeStep);
589  sprintf(traceNote, "s:%d", step);
590  traceUserSuppliedNote(traceNote);
591  }
592  // if (gotsigint) {
593  // iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
594  // NAMD_quit();
595  // }
596  outputExtendedSystem(step);
597 #if CYCLE_BARRIER
598  cycleBarrier(!((step+1) % stepsPerCycle),step);
599 #elif PME_BARRIER
600  cycleBarrier(dofull && !(step%slowFreq),step);
601 #elif STEP_BARRIER
602  cycleBarrier(1, step);
603 #endif
604 
605  if(Node::Object()->specialTracing || simParams->statsOn){
606  int bstep = simParams->traceStartStep;
607  int estep = bstep + simParams->numTraceSteps;
608  if(step == bstep){
609  traceBarrier(1, step);
610  }
611  if(step == estep){
612  traceBarrier(0, step);
613  }
614  }
615 
616 #ifdef MEASURE_NAMD_WITH_PAPI
617  if(simParams->papiMeasure) {
618  int bstep = simParams->papiMeasureStartStep;
619  int estep = bstep + simParams->numPapiMeasureSteps;
620  if(step == bstep) {
621  papiMeasureBarrier(1, step);
622  }
623  if(step == estep) {
624  papiMeasureBarrier(0, step);
625  }
626  }
627 #endif
628 
629  rebalanceLoad(step);
630 
631 #if PME_BARRIER
632  cycleBarrier(dofull && !((step+1)%slowFreq),step); // step before PME
633 #endif
634  }
635  }else CthSuspend();
636  // signal(SIGINT, oldhandler);
637 }
638 
640  const int size = simParams->movingAverageWindowSize;
643  pressureAverage.reset(size);
645 }
646 
647 #ifdef NODEGROUP_FORCE_REGISTER
648 
649 void Controller::printStep(int step, NodeReduction *nr){
650  char traceNote[24];
651  this->nodeReduction = nr;
652  const int numberOfSteps = simParams->N;
653  const int stepsPerCycle = simParams->stepsPerCycle;
654  const int zeroMomentum = simParams->zeroMomentum;
656  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
657  if (dofull)
659  else
661  if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
662  // Ok, now I need a way to set reduction to the nodeReduction
663 
664 if(step == simParams->firstTimestep /* && scriptTask == SCRIPT_RUN */){
665  reassignVelocities(step); // only for full-step velecities
666  rescaleaccelMD(step);
667  adaptTempUpdate(step); // Init adaptive tempering;
668 
669  receivePressure(step);
670  if ( zeroMomentum && dofull && ! (step % slowFreq) )
671  correctMomentum(step);
672  printFepMessage(step);
673  printTiMessage(step);
674  printDynamicsEnergies(step);
675  outputFepEnergy(step);
676  outputTiEnergy(step);
677  if(traceIsOn()){
678  traceUserEvent(eventEndOfTimeStep);
679  sprintf(traceNote, "s:%d", step);
680  traceUserSuppliedNote(traceNote);
681  }
682  outputExtendedSystem(step);
683  }else{
684  adaptTempUpdate(step);
685  rescaleVelocities(step);
686  tcoupleVelocities(step);
687  // XXX for single GPU version, call directly from Sequencer
688  // stochRescaleVelocities(step);
689  // XXX for single GPU version, call directly from Sequencer
690  // berendsenPressure(step);
691  // XXX for single GPU version, call directly from Sequencer
692  // langevinPiston1(step);
693  // XXX for single GPU version, call directly from Sequencer
694  // monteCarloPressure_prepare(step);
695  // monteCarloPressure_accept(step);
696  rescaleaccelMD(step);
697  enqueueCollections(step); // after lattice scaling!
698  receivePressure(step);
699  if ( zeroMomentum && dofull && ! (step % slowFreq) )
700  correctMomentum(step);
701  langevinPiston2(step);
702  reassignVelocities(step);
703 
704  multigratorTemperature(step, 1);
705  multigratorPressure(step, 1);
706  multigratorPressure(step, 2);
707  multigratorTemperature(step, 2);
708 
709  printDynamicsEnergies(step);
710  outputFepEnergy(step);
711  outputTiEnergy(step);
712  if(traceIsOn()){
713  traceUserEvent(eventEndOfTimeStep);
714  sprintf(traceNote, "s:%d", step);
715  traceUserSuppliedNote(traceNote);
716  }
717  outputExtendedSystem(step);
718  }
719 }
720 
721 #endif
722 
723 
724 #define CALCULATE \
725  printMinimizeEnergies(step); \
726  outputExtendedSystem(step); \
727  rebalanceLoad(step); \
728  if ( step == numberOfSteps ) return; \
729  else ++step;
730 
731 #define MOVETO(X) \
732  if ( step == numberOfSteps ) { \
733  if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
734  if ( newDir || (mid.x-last.x) ) { \
735  broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
736  } else { \
737  broadcast->minimizeCoefficient.publish(minSeq++,0.); \
738  broadcast->minimizeCoefficient.publish(minSeq++,0.); \
739  min_reduction->require(); \
740  broadcast->minimizeCoefficient.publish(minSeq++,0.); \
741  } \
742  enqueueCollections(step); \
743  CALCULATE \
744  } else if ( (X)-last.x ) { \
745  broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
746  newDir = 0; \
747  last.x = (X); \
748  enqueueCollections(step); \
749  CALCULATE \
750  last.u = min_energy; \
751  last.dudx = -1. * min_f_dot_v; \
752  last.noGradient = min_huge_count; \
753  if ( minVerbose ) { \
754  iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
755  if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
756  iout << "\n" << endi; \
757  } \
758  }
759 
760 struct minpoint {
762  minpoint() : x(0), u(0), dudx(0), noGradient(1) { ; }
763 };
764 
766  // iout << "Controller::minimize() called.\n" << endi;
767 
768  const int minVerbose = simParams->minVerbose;
769  const int numberOfSteps = simParams->N;
770  int step = simParams->firstTimestep;
771  slowFreq = nbondFreq = 1;
772  BigReal linegoal = simParams->minLineGoal; // 1.0e-3
773  const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
774 
775  CALCULATE
776 
777  int minSeq = 0;
778 
779  // just move downhill until initial bad contacts go away or 100 steps
780  int old_min_huge_count = 2000000000;
781  int steps_at_same_min_huge_count = 0;
782  for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
783  if ( min_huge_count >= old_min_huge_count ) {
784  if ( ++steps_at_same_min_huge_count > 10 ) break;
785  } else {
786  old_min_huge_count = min_huge_count;
787  steps_at_same_min_huge_count = 0;
788  }
789  iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
790  broadcast->minimizeCoefficient.publish(minSeq++,1.);
791  enqueueCollections(step);
792  CALCULATE
793  }
794  if ( min_huge_count ) {
795  iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
796  }
797  iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
798 
799  int atStart = 2;
800  int errorFactor = 10;
801  BigReal old_f_dot_f = min_f_dot_f;
802  broadcast->minimizeCoefficient.publish(minSeq++,0.);
803  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
804  int newDir = 1;
807  while ( 1 ) {
808  // line minimization
809  // bracket minimum on line
810  minpoint lo,hi,mid,last;
811  BigReal x = 0;
812  lo.x = x;
813  lo.u = min_energy;
814  lo.dudx = -1. * min_f_dot_v;
816  mid = lo;
817  last = mid;
818  if ( minVerbose ) {
819  iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
820  if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
821  iout << "\n" << endi;
822  }
823  BigReal tol = fabs( linegoal * min_f_dot_v );
824  iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
825  fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
826  int start_with_huge = last.noGradient;
828  BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
829  x = maxstep; MOVETO(x);
830  // bracket minimum on line
831  while ( last.u < mid.u ) {
832  if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
833  x = 6 * maxstep; break;
834  }
835  lo = mid; mid = last;
836  x += maxstep;
837  if ( x > 5.5 * maxstep ) break;
838  MOVETO(x)
839  }
840  if ( x > 5.5 * maxstep ) {
841  iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
842  broadcast->minimizeCoefficient.publish(minSeq++,0.);
843  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
844  newDir = 1;
845  old_f_dot_f = min_f_dot_f;
848  continue;
849  }
850  hi = last;
851 #define PRINT_BRACKET \
852  iout << "LINE MINIMIZER BRACKET: DX " \
853  << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
854  " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
855  lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
857  // converge on minimum on line
858  int itcnt;
859  for ( itcnt = 10; itcnt > 0; --itcnt ) {
860  int progress = 1;
861  // select new position
862  if ( mid.noGradient ) {
863  if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) { // subdivide left side
864  x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
865  MOVETO(x)
866  if ( last.u <= mid.u ) {
867  hi = mid; mid = last;
868  } else {
869  lo = last;
870  }
871  } else { // subdivide right side
872  x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
873  MOVETO(x)
874  if ( last.u <= mid.u ) {
875  lo = mid; mid = last;
876  } else {
877  hi = last;
878  }
879  }
880  } else {
881  if ( mid.dudx > 0. ) { // subdivide left side
882  BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
883  BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
884  x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
885  BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
886  if ( xdiv ) x /= xdiv; else x = last.x;
887  if ( x > altxhi ) x = altxhi;
888  if ( x < altxlo ) x = altxlo;
889  if ( x-last.x == 0 ) break;
890  MOVETO(x)
891  if ( last.u <= mid.u ) {
892  hi = mid; mid = last;
893  } else {
894  if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
895  lo = last;
896  }
897  } else { // subdivide right side
898  BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
899  BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
900  x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
901  BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
902  if ( xdiv ) x /= xdiv; else x = last.x;
903  if ( x < altxlo ) x = altxlo;
904  if ( x > altxhi ) x = altxhi;
905  if ( x-last.x == 0 ) break;
906  MOVETO(x)
907  if ( last.u <= mid.u ) {
908  lo = mid; mid = last;
909  } else {
910  if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
911  hi = last;
912  }
913  }
914  }
916  if ( ! progress ) {
917  MOVETO(mid.x);
918  break;
919  }
920  if ( fabs(last.dudx) < tol ) break;
921  if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
922  if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
923  }
924  // new direction
925  broadcast->minimizeCoefficient.publish(minSeq++,0.);
926  BigReal c = min_f_dot_f / old_f_dot_f;
927  c = ( c > 1.5 ? 1.5 : c );
928  if ( atStart ) { c = 0; --atStart; }
929  if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
930  c = 0;
931  if ( errorFactor < 100 ) errorFactor += 10;
932  }
933  if ( c == 0 ) {
934  iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
935  }
936  broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
937  newDir = 1;
938  old_f_dot_f = min_f_dot_f;
941  }
942 
943 }
944 
945 #undef MOVETO
946 #undef CALCULATE
947 
948 // NOTE: Only isotropic case implemented here!
949 void Controller::multigratorPressure(int step, int callNumber) {
954  BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
957  {
958  // Compute new scaling factors and send them to Sequencer
959  BigReal V = state->lattice.volume();
960  BigReal Pinst = trace(controlPressure)/3.0;
961  BigReal PGsum = trace(momentumSqrSum);
962  //
963  multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
964  BigReal scale = exp(s*NGfac*multigratorXiT);
965  BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
966  // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
967  Tensor scaleTensor = Tensor::identity(scale);
968  Tensor volScaleTensor = Tensor::identity(scale);
969  Tensor velScaleTensor = Tensor::identity(velScale);
970  state->lattice.rescale(volScaleTensor);
971  if (callNumber == 1) {
972  broadcast->positionRescaleFactor.publish(step,scaleTensor);
973  broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
974  } else {
975  broadcast->positionRescaleFactor2.publish(step,scaleTensor);
976  broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);
977  }
978  }
979 
980  {
981  // Wait here for Sequencer to finish scaling and force computation
982  reduction->require();
983  Tensor virial_normal;
984  Tensor virial_nbond;
985  Tensor virial_slow;
986  Tensor intVirial_normal;
987  Tensor intVirial_nbond;
988  Tensor intVirial_slow;
989  Vector extForce_normal;
990  Vector extForce_nbond;
991  Vector extForce_slow;
992  GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
993  GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
994  GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
995  GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
996  GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
997  GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
998  GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
999  GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
1000  GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
1001  GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
1002  calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
1003  intVirial_normal, intVirial_nbond, intVirial_slow,
1004  extForce_normal, extForce_nbond, extForce_slow);
1005  if (callNumber == 2) {
1006  // Update temperature for the Temperature Cycle that is coming next
1009  }
1010  }
1011 
1012  {
1013  // Update pressure integrator
1014  BigReal V = state->lattice.volume();
1015  BigReal Pinst = trace(controlPressure)/3.0;
1016  BigReal PGsum = trace(momentumSqrSum);
1017  //
1018  multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
1019  }
1020 
1021  }
1022 }
1023 
1024 void Controller::multigratorTemperature(int step, int callNumber) {
1029  BigReal Nf = numDegFreedom;
1031  BigReal T1, T2, v;
1032  {
1033  T1 = temperature;
1034  BigReal kTinst = BOLTZMANN * temperature;
1035  for (int i=n-1;i >= 0;i--) {
1036  if (i == 0) {
1037  BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
1038  multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
1039  } else if (i == n-1) {
1040  multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
1041  } else {
1042  BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
1043  multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] +
1044  0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
1045  }
1046  }
1047  BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
1048  v = velScale;
1049  if (callNumber == 1)
1050  broadcast->velocityRescaleFactor.publish(step,velScale);
1051  else
1052  broadcast->velocityRescaleFactor2.publish(step,velScale);
1053  }
1054 
1055  {
1056  // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
1060  T2 = temperature;
1061  if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
1062  // If this is pressure cycle, receive new momentum product
1063  GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
1064  }
1065  }
1066 
1067  // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);
1068 
1069  {
1070  BigReal kTinst = BOLTZMANN * temperature;
1071  for (int i=0;i < n;i++) {
1072  if (i == 0) {
1073  BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
1074  multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
1075  } else if (i == n-1) {
1076  multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
1077  } else {
1078  BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
1079  multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] +
1080  0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
1081  }
1083  }
1084  }
1085 
1086  }
1087 }
1088 
1089 // Calculates Enthalpy for multigrator
1090 BigReal Controller::multigatorCalcEnthalpy(BigReal potentialEnergy, int step, int minimize) {
1091  Node *node = Node::Object();
1092  Molecule *molecule = node->molecule;
1093  SimParameters *simParameters = node->simParameters;
1094 
1095  BigReal V = state->lattice.volume();
1099  BigReal Nf = numDegFreedom;
1101  BigReal sumZeta = 0.0;
1102  for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
1103  sumZeta += multigratorZeta[i];
1104  }
1105  BigReal nuOmegaSum = 0.0;
1106  for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
1107  nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
1108  }
1109  BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
1110  BigReal eta = sqrt(kT0*W)*multigratorXi;
1111 
1112  BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);
1113 
1114 // if (!(step % 100))
1115  // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
1116  // kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));
1117 
1118  return enthalpy;
1119 }
1120 
1122 {
1124  if ( !(step % simParams->monteCarloPressureFreq) ) {
1125  Node *node = Node::Object();
1126  Molecule *molecule = node->molecule;
1127  // need to consider the fix atom case!
1128  int numAtom = molecule->numMolecules;
1129  BigReal volumeOld = mc_oldLattice.volume();
1130  BigReal volumeNew = state->lattice.volume();
1131  BigReal areaNew = (volumeNew / state->lattice.c().z);
1132  BigReal areaOld = (volumeOld / mc_oldLattice.c().z);
1134  BigReal pressureTarget = simParams->monteCarloPressureTarget;
1135  BigReal surfTensionTarget = simParams->surfaceTensionTarget;
1136 
1137  if(!simParams->CUDASOAintegrate) {
1138  // Wait here for Sequencer to finish scaling coordinates and energy/force computation
1139  reduction->require();
1140  }
1141 
1142  BigReal totalEnergyOld = mc_totalEnergyOld;
1143  BigReal totalEnergyNew = getTotalPotentialEnergy(step);
1144  BigReal weight = (totalEnergyNew - totalEnergyOld);
1145  weight += (pressureTarget * (volumeNew - volumeOld));
1146  weight -= (numAtom * kbT * log(volumeNew / volumeOld));
1147  weight -= (surfTensionTarget * (areaNew - areaOld));
1148  BigReal totalWeight = exp(-weight / kbT);
1149  int accepted = (random->uniform() < totalWeight);
1150  //accepted = 0;
1151  #if 0
1152  printf("MC-data (accept): step: %d, Pe: %d, numAtom: %d\n", step, CkMyPe(), numAtom);
1153  printf(" oldE: %f, newE: %f, deltaE: %f\n", totalEnergyOld, totalEnergyNew, totalEnergyNew - totalEnergyOld);
1154  printf(" oldV: %f, newV: %f, deltaV: %f\n", volumeOld, volumeNew, volumeNew - volumeOld);
1155  printf(" surfTen: %f, deltaA: %f \n", surfTensionTarget, areaNew - areaOld);
1156  printf(" weight: %f, accepted: %d\n\n", totalWeight, accepted);
1157  #endif
1158  if(accepted) {
1159  // what should I do if its accepted?
1160  ++mc_totalAccept;
1162  } else {
1163  // what should I do if its rejected?
1164  state->lattice = mc_oldLattice;
1165  }
1166  // update the maximum scale no matter if it's accepted or not
1167  {
1169  BigReal factor = 1.0;
1171  // safe check to make sure there is any accepted move
1172  if (accRate > 0.0) {
1173  factor = accRate / simParams->monteCarloAcceptanceRate;
1174  } else {
1175  factor = 0.5; // if no move was accepted, we scale by half
1176  }
1177  // reset the trial for picked axis
1179  BigReal maxVolume = state->lattice.volume() * 0.3;
1180  switch (mc_picked_axis) {
1181  case MC_X:
1182  monteCarloMaxVolume.x *= factor;
1183  if (monteCarloMaxVolume.x < 0.001) {
1184  iout << iWARN << "Small volume change in MC barostat!\n" << endi;
1185  iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
1186  monteCarloMaxVolume.x = 0.001;
1187  } else if (monteCarloMaxVolume.x > maxVolume ) {
1188  iout << iWARN << "Large volume change in MC barostat!\n" << endi;
1189  iout << iWARN << "Maximum volume change is set to " << maxVolume << " A^3.\n" << endi;
1190  monteCarloMaxVolume.x = maxVolume;
1191  }
1192  break;
1193 
1194  case MC_Y:
1195  monteCarloMaxVolume.y *= factor;
1196  if (monteCarloMaxVolume.y < 0.001) {
1197  iout << iWARN << "Small volume change in MC barostat!\n" << endi;
1198  iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
1199  monteCarloMaxVolume.y = 0.001;
1200  } else if (monteCarloMaxVolume.y > maxVolume ) {
1201  iout << iWARN << "Large volume change in MC barostat!\n" << endi;
1202  iout << iWARN << "Maximum volume change is set to " << maxVolume << " A^3.\n" << endi;
1203  monteCarloMaxVolume.y = maxVolume;
1204  }
1205  break;
1206 
1207  case MC_Z:
1208  monteCarloMaxVolume.z *= factor;
1209  if (monteCarloMaxVolume.z < 0.001) {
1210  iout << iWARN << "Small volume change in MC barostat!\n" << endi;
1211  iout << iWARN << "Maximum volume change is set to 0.001 A^3.\n" << endi;
1212  monteCarloMaxVolume.z = 0.001;
1213  } else if (monteCarloMaxVolume.z > maxVolume ) {
1214  iout << iWARN << "Large volume change in MC barostat!\n" << endi;
1215  iout << iWARN << "Maximum volume change is set to " << maxVolume << " A^3.\n" << endi;
1216  monteCarloMaxVolume.z = maxVolume;
1217  }
1218  break;
1219  }
1220  }
1221  }
1222 
1223  if ( !(step % simParams->outputEnergies) ) {
1224  BigReal acceptPerc = 100.0 * ((BigReal) mc_totalAccept /
1225  (BigReal) mc_totalTry);
1226  iout << "MC_BAROSTAT STEP: " << step << " TRIALS: " << mc_totalTry <<
1227  " ACCEPTED: " << mc_totalAccept << " %ACCEPTANCE: " << acceptPerc;
1228  if (simParams->useFlexibleCell) {
1229  if (simParams->useConstantArea) {
1230  iout << " MAX_VOLUME_Z: " << monteCarloMaxVolume.z << "\n\n";
1231  } else if (simParams->useConstantRatio) {
1232  iout << " MAX_VOLUME_XY: " << monteCarloMaxVolume.x <<
1233  " MAX_VOLUME_Z: " << monteCarloMaxVolume.z << "\n\n";
1234  } else {
1235  iout << " MAX_VOLUME_X: " << monteCarloMaxVolume.x
1236  << " MAX_VOLUME_Y: " << monteCarloMaxVolume.y
1237  << " MAX_VOLUME_Z: " << monteCarloMaxVolume.z
1238  << "\n\n";
1239  }
1240  } else {
1241  iout << " MAX_VOLUME_XYZ: " << monteCarloMaxVolume.x << "\n\n";
1242  }
1243  }
1244 
1246  #ifdef NODEGROUP_FORCE_REGISTER
1247  publishedMCAcceptance.insert({step, accepted});
1248  #endif
1249  }
1250  }
1251 }
1252 
1254 {
1256  if ( !(step % simParams->monteCarloPressureFreq) ) {
1257  if(!simParams->CUDASOAintegrate) {
1258  // Wait here for Sequencer to finish energy/force computation
1259  reduction->require();
1260  }
1261  mc_oldLattice = state->lattice;
1263  Tensor factor = Tensor::identity(1.0);
1264  // pick a random deltaV and calc scaling factor
1265  {
1266  BigReal oldVolume = mc_oldLattice.volume();
1267  bool flexible = simParams->useFlexibleCell;
1268  bool constArea = simParams->useConstantArea;
1269  bool constRatio = simParams->useConstantRatio;
1270  bool constZ, constY, constX;
1271  if (simParams->fixCellDims) {
1272  constX = simParams->fixCellDimX;
1273  constY = simParams->fixCellDimY;
1274  constZ = simParams->fixCellDimZ;
1275  } else {
1276  constZ = constY = constX = false;
1277  }
1278 
1279  if(flexible) {
1280  //Anisotropic fluctuation
1281  if(constArea) {
1282  // always pick z
1283  mc_picked_axis = MC_Z;
1284  BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
1285  factor.zz = (oldVolume + dV) / oldVolume;
1286  } else if (constRatio) {
1287  while(true) {
1288  // decide to scale x-y or scale z
1289  BigReal probAxis = 2.0 * random->uniform();
1290  if (probAxis < 1.0) {
1291  //picking x or y. We use MC_X for constant ratio
1292  mc_picked_axis = MC_X;
1293  BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
1294  BigReal scale = sqrt((oldVolume + dV)/oldVolume);
1295  factor.xx = factor.yy = scale;
1296  break;
1297  } else if (probAxis < 2.0 && !constZ) {
1298  //picking z
1299  mc_picked_axis = MC_Z;
1300  BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
1301  BigReal scale = (oldVolume + dV)/oldVolume;
1302  factor.zz = scale;
1303  break;
1304  }
1305  }
1306  } else {
1307  while(true) {
1308  // decide to scale x, y or z
1309  BigReal probAxis = 3.0 * random->uniform();
1310  if (probAxis < 1.0 && !constX) {
1311  // picking x
1312  mc_picked_axis = MC_X;
1313  BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
1314  factor.xx = (oldVolume + dV) / oldVolume;
1315  break;
1316  } else if (probAxis < 2.0 && !constY) {
1317  // picking y
1318  mc_picked_axis = MC_Y;
1319  BigReal dV = monteCarloMaxVolume.y * (2.0 * random->uniform() - 1.0);
1320  factor.yy = (oldVolume + dV) / oldVolume;
1321  break;
1322  } else if (probAxis < 3.0 && !constZ) {
1323  // picking z
1324  mc_picked_axis = MC_Z;
1325  BigReal dV = monteCarloMaxVolume.z * (2.0 * random->uniform() - 1.0);
1326  factor.zz = (oldVolume + dV) / oldVolume;
1327  break;
1328  }
1329  }
1330  }
1331  } else {
1332  //Isotropic fluctuation. We use MC_X for isotropic fluctuation
1333  mc_picked_axis = MC_X;
1334  BigReal dV = monteCarloMaxVolume.x * (2.0 * random->uniform() - 1.0);
1335  BigReal scale = pow((oldVolume + dV)/oldVolume, 1.0/3.0);
1336  factor.xx = factor.yy = factor.zz = scale;
1337  }
1338  // increament the MC trial
1339  ++mc_totalTry;
1341  }
1342 
1343  broadcast->positionRescaleFactor.publish(step,factor);
1344  state->lattice.rescale(factor);
1345  #ifdef NODEGROUP_FORCE_REGISTER
1346  publishedRescaleFactors.insert({step, factor});
1347  #endif
1348  }
1349  }
1350 }
1351 
1352 
1353 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
1354  if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
1355  if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
1356 
1358 {
1359  if ( simParams->berendsenPressureOn ) {
1362  const int freq = simParams->berendsenPressureFreq;
1363  if ( ! (berendsenPressure_count % freq) ) {
1367  // We only use on-diagonal terms (for now)
1368  factor = Tensor::diagonal(diagonal(factor));
1370  factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
1372  factor += Tensor::identity(1.0);
1373  int limited = 0;
1374  LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
1375  LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
1376  LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
1377  if ( limited ) {
1378  iout << iERROR << "Step " << step <<
1379  " cell rescaling factor limited.\n" << endi;
1380  }
1381  broadcast->positionRescaleFactor.publish(step,factor);
1382  state->lattice.rescale(factor);
1383 #ifdef NODEGROUP_FORCE_REGISTER
1384  publishedRescaleFactors.insert({step, factor});
1385 #endif
1386  }
1387  } else {
1390  }
1391 }
1392 
1393 #undef LIMIT_SCALING
1394 
1396 {
1397  if ( simParams->langevinPistonOn )
1398  {
1399  Tensor &strainRate = langevinPiston_strainRate;
1400  int cellDims = simParams->useFlexibleCell ? 1 : 3;
1401  BigReal dt = simParams->dt;
1402  BigReal dt_long = slowFreq * dt;
1405  BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
1406 
1407 #ifdef DEBUG_PRESSURE
1408  iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
1409 #endif
1410 
1411  if ( ! ( (step-1) % slowFreq ) )
1412  {
1413  BigReal gamma = 1 / simParams->langevinPistonDecay;
1414  BigReal f1 = exp( -0.5 * dt_long * gamma );
1415  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1416  strainRate *= f1;
1417  if ( simParams->useFlexibleCell ) {
1418  // We only use on-diagonal terms (for now)
1419  if ( simParams->useConstantRatio ) {
1420  BigReal r = f2 * random->gaussian();
1421  strainRate.xx += r;
1422  strainRate.yy += r;
1423  strainRate.zz += f2 * random->gaussian();
1424  } else {
1425  strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
1426  }
1427  } else {
1428  strainRate += f2 * Tensor::identity(random->gaussian());
1429  }
1430 
1431 #ifdef DEBUG_PRESSURE
1432  iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
1433 #endif
1434  }
1435 
1436  // Apply surface tension. If surfaceTensionTarget is zero, we get
1437  // the default (isotropic pressure) case.
1438 
1439  Tensor ptarget;
1440  ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
1441  if ( simParams->surfaceTensionTarget != 0. ) {
1442  ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
1443  simParams->surfaceTensionTarget / state->lattice.c().z;
1444  }
1445 
1446  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1447  ( controlPressure - ptarget );
1448 
1449  if (simParams->fixCellDims) {
1450  if (simParams->fixCellDimX) strainRate.xx = 0;
1451  if (simParams->fixCellDimY) strainRate.yy = 0;
1452  if (simParams->fixCellDimZ) strainRate.zz = 0;
1453  }
1454 
1455 
1456 #ifdef DEBUG_PRESSURE
1457  iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
1458 #endif
1459 
1461  if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
1462  {
1463  // We only use on-diagonal terms (for now)
1464  Tensor factor;
1465  if ( !simParams->useConstantArea ) {
1466  factor.xx = exp( dt_long * strainRate.xx );
1467  factor.yy = exp( dt_long * strainRate.yy );
1468  } else {
1469  factor.xx = factor.yy = 1;
1470  }
1471  factor.zz = exp( dt_long * strainRate.zz );
1472 #ifdef NODEGROUP_FORCE_REGISTER
1473  if(simParams->CUDASOAintegrate) // don't use broadcastMgr
1474  publishedRescaleFactors.insert({step, factor});
1475  else
1476  broadcast->positionRescaleFactor.publish(step,factor);
1477 #else
1478  broadcast->positionRescaleFactor.publish(step,factor);
1479 #endif
1480  state->lattice.rescale(factor);
1481 #ifdef DEBUG_PRESSURE
1482  iout << iINFO << "rescaling by: " << factor << "\n";
1483 #endif
1484  }
1485  } else { // langevinPistonBarrier
1486  if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
1487  {
1488  if ( ! positionRescaleFactor.xx ) { // first pass without MTS
1489  // We only use on-diagonal terms (for now)
1490  Tensor factor;
1491  if ( !simParams->useConstantArea ) {
1492  factor.xx = exp( dt_long * strainRate.xx );
1493  factor.yy = exp( dt_long * strainRate.yy );
1494  } else {
1495  factor.xx = factor.yy = 1;
1496  }
1497  factor.zz = exp( dt_long * strainRate.zz );
1498 #ifdef NODEGROUP_FORCE_REGISTER
1499  if(simParams->CUDASOAintegrate) // don't use broadcastMgr
1500  publishedRescaleFactors.insert({step, factor});
1501  else
1502  broadcast->positionRescaleFactor.publish(step,factor);
1503 #else
1504  broadcast->positionRescaleFactor.publish(step,factor);
1505 #endif
1506  positionRescaleFactor = factor;
1507  strainRate_old = strainRate;
1508  }
1510 #ifdef DEBUG_PRESSURE
1511  iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
1512 #endif
1513  }
1514  if ( ! ( (step-slowFreq/2) % slowFreq ) )
1515  {
1516  // We only use on-diagonal terms (for now)
1517  Tensor factor;
1518  if ( !simParams->useConstantArea ) {
1519  factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
1520  factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
1521  } else {
1522  factor.xx = factor.yy = 1;
1523  }
1524  factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
1525 #ifdef NODEGROUP_FORCE_REGISTER
1526  if(simParams->CUDASOAintegrate) // don't use broadcastMgr
1527  publishedRescaleFactors.insert({step+1, factor});
1528  else
1529  broadcast->positionRescaleFactor.publish(step+1,factor);
1530 #else
1531  broadcast->positionRescaleFactor.publish(step+1,factor);
1532 #endif
1533  positionRescaleFactor = factor;
1534  strainRate_old = strainRate;
1535  }
1536  }
1537 
1538  // corrections to integrator
1539  if ( ! ( step % nbondFreq ) )
1540  {
1541 #ifdef DEBUG_PRESSURE
1542  iout << iINFO << "correcting strain rate for nbond, ";
1543 #endif
1544  strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1545  ( (nbondFreq - 1) * controlPressure_nbond );
1546 #ifdef DEBUG_PRESSURE
1547  iout << "strain rate: " << strainRate << "\n";
1548 #endif
1549  }
1550  if ( ! ( step % slowFreq ) )
1551  {
1552 #ifdef DEBUG_PRESSURE
1553  iout << iINFO << "correcting strain rate for slow, ";
1554 #endif
1555  strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1556  ( (slowFreq - 1) * controlPressure_slow );
1557 #ifdef DEBUG_PRESSURE
1558  iout << "strain rate: " << strainRate << "\n";
1559 #endif
1560  }
1561  if (simParams->fixCellDims) {
1562  if (simParams->fixCellDimX) strainRate.xx = 0;
1563  if (simParams->fixCellDimY) strainRate.yy = 0;
1564  if (simParams->fixCellDimZ) strainRate.zz = 0;
1565  }
1566 
1567  }
1568 
1569 }
1570 
1572 {
1573  if ( simParams->langevinPistonOn )
1574  {
1575  Tensor &strainRate = langevinPiston_strainRate;
1576  int cellDims = simParams->useFlexibleCell ? 1 : 3;
1577  BigReal dt = simParams->dt;
1578  BigReal dt_long = slowFreq * dt;
1581  BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
1582 
1583  // corrections to integrator
1584  if ( ! ( step % nbondFreq ) )
1585  {
1586 #ifdef DEBUG_PRESSURE
1587  iout << iINFO << "correcting strain rate for nbond, ";
1588 #endif
1589  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1590  ( (nbondFreq - 1) * controlPressure_nbond );
1591 #ifdef DEBUG_PRESSURE
1592  iout << "strain rate: " << strainRate << "\n";
1593 #endif
1594  }
1595  if ( ! ( step % slowFreq ) )
1596  {
1597 #ifdef DEBUG_PRESSURE
1598  iout << iINFO << "correcting strain rate for slow, ";
1599 #endif
1600  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1601  ( (slowFreq - 1) * controlPressure_slow );
1602 #ifdef DEBUG_PRESSURE
1603  iout << "strain rate: " << strainRate << "\n";
1604 #endif
1605  }
1606 
1607  // Apply surface tension. If surfaceTensionTarget is zero, we get
1608  // the default (isotropic pressure) case.
1609 
1610  Tensor ptarget;
1611  ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
1612  if ( simParams->surfaceTensionTarget != 0. ) {
1613  ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
1614  simParams->surfaceTensionTarget / state->lattice.c().z;
1615  }
1616 
1617  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1618  ( controlPressure - ptarget );
1619 
1620 
1621 #ifdef DEBUG_PRESSURE
1622  iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
1623 #endif
1624 
1625  if ( ! ( step % slowFreq ) )
1626  {
1627  BigReal gamma = 1 / simParams->langevinPistonDecay;
1628  BigReal f1 = exp( -0.5 * dt_long * gamma );
1629  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1630  strainRate *= f1;
1631  if ( simParams->useFlexibleCell ) {
1632  // We only use on-diagonal terms (for now)
1633  if ( simParams->useConstantRatio ) {
1634  BigReal r = f2 * random->gaussian();
1635  strainRate.xx += r;
1636  strainRate.yy += r;
1637  strainRate.zz += f2 * random->gaussian();
1638  } else {
1639  strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
1640  }
1641  } else {
1642  strainRate += f2 * Tensor::identity(random->gaussian());
1643  }
1644 #ifdef DEBUG_PRESSURE
1645  iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
1646 #endif
1647  }
1648 
1649 #ifdef DEBUG_PRESSURE
1650  iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
1651 #endif
1652  if (simParams->fixCellDims) {
1653  if (simParams->fixCellDimX) strainRate.xx = 0;
1654  if (simParams->fixCellDimY) strainRate.yy = 0;
1655  if (simParams->fixCellDimZ) strainRate.zz = 0;
1656  }
1657  }
1658 
1659 
1660 }
1661 
1663 {
1664  const int rescaleFreq = simParams->rescaleFreq;
1665  if ( rescaleFreq > 0 ) {
1667  if ( rescaleVelocities_numTemps == rescaleFreq ) {
1669  BigReal rescaleTemp = simParams->rescaleTemp;
1671  (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
1672  rescaleTemp = adaptTempT;
1673  }
1674  BigReal factor = sqrt(rescaleTemp/avgTemp);
1675  broadcast->velocityRescaleFactor.publish(step,factor);
1676  //iout << "RESCALING VELOCITIES AT STEP " << step
1677  // << " FROM AVERAGE TEMPERATURE OF " << avgTemp
1678  // << " TO " << rescaleTemp << " KELVIN.\n" << endi;
1680  }
1681  }
1682 }
1683 
1685 
1686  Vector momentum;
1687  BigReal mass;
1688 #ifdef NODEGROUP_FORCE_REGISTER
1690  momentum.x = nodeReduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
1691  momentum.y = nodeReduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
1692  momentum.z = nodeReduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
1694  }else{
1695 #endif
1696  momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
1697  momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
1698  momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
1700 #ifdef NODEGROUP_FORCE_REGISTER
1701  }
1702 #endif
1703 
1704 
1705  if ( momentum.length2() == 0. )
1706  iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
1707  if ( mass == 0. )
1708  NAMD_die("Total mass is zero in Controller::correctMomentum");
1709 
1710  momentum *= (-1./mass);
1711  // What does this mean???
1712  broadcast->momentumCorrection.publish(step+slowFreq,momentum);
1713 }
1714 
1715 //Modifications for alchemical fep
1717 {
1719  && !simParams->alchLambdaFreq) {
1720  const BigReal alchLambda = simParams->alchLambda;
1721  const BigReal alchLambda2 = simParams->alchLambda2;
1722  const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
1723  const BigReal alchTemp = simParams->alchTemp;
1724  const int alchEquilSteps = simParams->alchEquilSteps;
1725  iout << "FEP: RESETTING FOR NEW FEP WINDOW "
1726  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
1727  if ( alchLambdaIDWS >= 0. ) {
1728  iout << " LAMBDA_IDWS " << alchLambdaIDWS;
1729  }
1730  iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
1731  << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
1732  << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
1733  << " K FOR FEP CALCULATION\n" << endi;
1734  }
1735 }
1737 {
1739  && !simParams->alchLambdaFreq) {
1740  const BigReal alchLambda = simParams->alchLambda;
1741  const int alchEquilSteps = simParams->alchEquilSteps;
1742  iout << "TI: RESETTING FOR NEW WINDOW "
1743  << "LAMBDA SET TO " << alchLambda
1744  << "\nTI: WINDOW TO HAVE " << alchEquilSteps
1745  << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
1746  }
1747 }
1748 //fepe
1749 
1751 {
1752  const int reassignFreq = simParams->reassignFreq;
1753  if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1754  BigReal newTemp = simParams->reassignTemp;
1755  newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
1756  if ( simParams->reassignIncr > 0.0 ) {
1757  if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
1758  newTemp = simParams->reassignHold;
1759  } else {
1760  if ( newTemp < simParams->reassignHold )
1761  newTemp = simParams->reassignHold;
1762  }
1763  iout << "REASSIGNING VELOCITIES AT STEP " << step
1764  << " TO " << newTemp << " KELVIN.\n" << endi;
1765  }
1766 }
1767 
1769 {
1770  if ( simParams->tCoupleOn )
1771  {
1772  const BigReal tCoupleTemp = simParams->tCoupleTemp;
1773  BigReal coefficient = 1.;
1774  if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
1775  broadcast->tcoupleCoefficient.publish(step,coefficient);
1776  }
1777 }
1778 
1791 {
1792  if ( simParams->stochRescaleOn ) {
1795  double coefficient = stochRescaleCoefficient();
1796  broadcast->stochRescaleCoefficient.publish(step,coefficient);
1797  stochRescale_count = 0;
1798  }
1799  }
1800 }
1801 
1807  const double stochRescaleTemp = simParams->stochRescaleTemp;
1808  double coefficient = 1;
1809  if ( temperature > 0 ) {
1810  double R1 = random->gaussian();
1811  // double gammaShape = 0.5*(numDegFreedom - 1);
1812  // R2sum is the sum of (numDegFreedom - 1) squared normal variables,
1813  // which is chi-squared distributed.
1814  // This is in turn a special case of the Gamma distribution,
1815  // which converges to a normal distribution in the limit of a
1816  // large shape parameter.
1817  // double R2sum = 2*(gammaShape+sqrt(gammaShape)*random->gaussian())+R1*R1;
1818  double R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
1819  double tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
1820 
1821  coefficient = sqrt(stochRescaleTimefactor +
1822  (1 - stochRescaleTimefactor)*tempfactor*R2sum +
1823  2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*
1825  }
1826  heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient-1);
1827  return coefficient;
1828 }
1829 
1830 static char *FORMAT(BigReal X, int decimal = 4)
1831 {
1832  static char tmp_string[50];
1833  static char format_string[50];
1834  const double maxnum = 99999999999.9999;
1835  if ( X > maxnum ) X = maxnum;
1836  if ( X < -maxnum ) X = -maxnum;
1837 
1838  int whole = (decimal <= 4 ? 14 : 10 + decimal);
1839  sprintf(format_string, " %%%d.%df", whole, decimal);
1840  sprintf(tmp_string, format_string, X);
1841  return tmp_string;
1842 }
1843 
1844 static char *FORMAT(const char *X, int decimal = 4)
1845 {
1846  static char tmp_string[50];
1847  static char format_string[50];
1848 
1849  int width = (decimal <= 4 ? 14 : 10 + decimal);
1850  sprintf(format_string, " %%%ds", width);
1851  sprintf(tmp_string, format_string, X);
1852  return tmp_string;
1853 }
1854 
1855 static char *ETITLE(int X)
1856 {
1857  static char tmp_string[21];
1858  sprintf(tmp_string,"ENERGY: %7d",X);
1859  return tmp_string;
1860 }
1861 
1862 void Controller::receivePressure(int step, int minimize)
1863 {
1864  Node *node = Node::Object();
1865  Molecule *molecule = node->molecule;
1866  Lattice &lattice = state->lattice;
1867  bool cudaIntegrator = simParams->CUDASOAintegrate;
1868 
1869  if(!cudaIntegrator) reduction->require();
1870 
1871  Tensor virial_normal;
1872  Tensor virial_nbond;
1873  Tensor virial_slow;
1874 #ifdef ALTVIRIAL
1875  Tensor altVirial_normal;
1876  Tensor altVirial_nbond;
1877  Tensor altVirial_slow;
1878 #endif
1879  Tensor intVirial_normal;
1880  Tensor intVirial_nbond;
1881  Tensor intVirial_slow;
1882  Vector extForce_normal;
1883  Vector extForce_nbond;
1884  Vector extForce_slow;
1885 
1886 #if 1
1887  numDegFreedom = molecule->num_deg_freedom();
1888  int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
1889  int numFixedGroups = molecule->num_fixed_groups();
1890  int numFixedAtoms = molecule->num_fixed_atoms();
1891 #endif
1892 #if 0
1893  int numAtoms = molecule->numAtoms;
1894  numDegFreedom = 3 * numAtoms;
1895  int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
1896  int numFixedAtoms =
1897  ( simParams->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
1898  int numLonepairs = molecule->numLonepairs;
1899  int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
1900  if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
1901  if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
1902  if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
1903  if ( ! ( numFixedAtoms || molecule->numConstraints
1904  || simParams->comMove || simParams->langevinOn ) ) {
1905  numDegFreedom -= 3;
1906  numGroupDegFreedom -= 3;
1907  }
1909  // this doesn't attempt to deal with fixed atoms or constraints
1910  numDegFreedom = 3 * molecule->numFepInitial;
1911  }
1912  int numRigidBonds = molecule->numRigidBonds;
1913  int numFixedRigidBonds =
1914  ( simParams->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
1915 
1916  // numLonepairs is subtracted here because all lonepairs have a rigid bond
1917  // to oxygen, but all of the LP degrees of freedom are dealt with above
1918  numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
1919 #endif
1920  BigReal groupKineticEnergyHalfstep, groupKineticEnergyCentered;
1921 #ifdef NODEGROUP_FORCE_REGISTER
1922  if(cudaIntegrator){
1925  groupKineticEnergyHalfstep = kineticEnergyHalfstep -
1927  groupKineticEnergyCentered = kineticEnergyCentered -
1929  }else{
1930 #endif
1933  groupKineticEnergyHalfstep = kineticEnergyHalfstep -
1935  groupKineticEnergyCentered = kineticEnergyCentered -
1937 #ifdef NODEGROUP_FORCE_REGISTER
1938  }
1939 #endif
1940 
1941 
1942 
1943  BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
1944  / ( numDegFreedom * BOLTZMANN );
1945  BigReal atomTempCentered = 2.0 * kineticEnergyCentered
1946  / ( numDegFreedom * BOLTZMANN );
1947  BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
1948  / ( numGroupDegFreedom * BOLTZMANN );
1949  BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
1950  / ( numGroupDegFreedom * BOLTZMANN );
1951 
1952  /* test code for comparing different temperatures
1953  iout << "TEMPTEST: " << step << " " <<
1954  atomTempHalfstep << " " <<
1955  atomTempCentered << " " <<
1956  groupTempHalfstep << " " <<
1957  groupTempCentered << "\n" << endi;
1958  iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
1959  numGroupDegFreedom << "\n" << endi;
1960  */
1961 #ifdef NODEGROUP_FORCE_REGISTER
1962  if(cudaIntegrator){
1963  GET_TENSOR(momentumSqrSum, nodeReduction,REDUCTION_MOMENTUM_SQUARED);
1964 
1965  GET_TENSOR(virial_normal,nodeReduction,REDUCTION_VIRIAL_NORMAL);
1966  GET_TENSOR(virial_nbond, nodeReduction,REDUCTION_VIRIAL_NBOND);
1967  GET_TENSOR(virial_slow, nodeReduction,REDUCTION_VIRIAL_SLOW);
1968  }else{
1969 #endif
1970  GET_TENSOR(momentumSqrSum, reduction,REDUCTION_MOMENTUM_SQUARED);
1971 
1972  GET_TENSOR(virial_normal, reduction,REDUCTION_VIRIAL_NORMAL);
1973  GET_TENSOR(virial_nbond, reduction,REDUCTION_VIRIAL_NBOND);
1974  GET_TENSOR(virial_slow, reduction,REDUCTION_VIRIAL_SLOW);
1975 #ifdef NODEGROUP_FORCE_REGISTER
1976  }
1977 #endif
1978 
1979 #ifdef ALTVIRIAL
1980 #ifdef NODEGROUP_FORCE_REGISTER
1981  if(cudaIntegrator){
1982  GET_TENSOR(altVirial_normal,nodeReduction,REDUCTION_ALT_VIRIAL_NORMAL);
1983  GET_TENSOR(altVirial_nbond, nodeReduction,REDUCTION_ALT_VIRIAL_NBOND);
1984  GET_TENSOR(altVirial_slow, nodeReduction,REDUCTION_ALT_VIRIAL_SLOW);
1985  }else{
1986 #endif
1987  GET_TENSOR(altVirial_normal, reduction,REDUCTION_ALT_VIRIAL_NORMAL);
1988  GET_TENSOR(altVirial_nbond, reduction,REDUCTION_ALT_VIRIAL_NBOND);
1989  GET_TENSOR(altVirial_slow, reduction,REDUCTION_ALT_VIRIAL_SLOW);
1990 #ifdef NODEGROUP_FORCE_REGISTER
1991  }
1992 #endif
1993 #endif
1994 
1995 #ifdef NODEGROUP_FORCE_REGISTER
1996  if(cudaIntegrator){
1997  GET_TENSOR(intVirial_normal,nodeReduction,REDUCTION_INT_VIRIAL_NORMAL);
1998  GET_TENSOR(intVirial_nbond, nodeReduction,REDUCTION_INT_VIRIAL_NBOND);
1999  GET_TENSOR(intVirial_slow, nodeReduction,REDUCTION_INT_VIRIAL_SLOW);
2000 
2001  GET_VECTOR(extForce_normal,nodeReduction,REDUCTION_EXT_FORCE_NORMAL);
2002  GET_VECTOR(extForce_nbond,nodeReduction,REDUCTION_EXT_FORCE_NBOND);
2003  GET_VECTOR(extForce_slow,nodeReduction,REDUCTION_EXT_FORCE_SLOW);
2004  }else{
2005 #endif
2006  GET_TENSOR(intVirial_normal, reduction,REDUCTION_INT_VIRIAL_NORMAL);
2007  GET_TENSOR(intVirial_nbond, reduction,REDUCTION_INT_VIRIAL_NBOND);
2008  GET_TENSOR(intVirial_slow, reduction,REDUCTION_INT_VIRIAL_SLOW);
2009 
2010  GET_VECTOR(extForce_normal, reduction,REDUCTION_EXT_FORCE_NORMAL);
2011  GET_VECTOR(extForce_nbond, reduction,REDUCTION_EXT_FORCE_NBOND);
2012  GET_VECTOR(extForce_slow, reduction,REDUCTION_EXT_FORCE_SLOW);
2013 #ifdef NODEGROUP_FORCE_REGISTER
2014  }
2015 #endif
2016 
2017  // APH NOTE: These four lines are now done in calcPressure()
2018  // Vector extPosition = lattice.origin();
2019  // virial_normal -= outer(extForce_normal,extPosition);
2020  // virial_nbond -= outer(extForce_nbond,extPosition);
2021  // virial_slow -= outer(extForce_slow,extPosition);
2022 
2025 
2026  if (simParams->drudeOn) {
2027  BigReal drudeComKE, drudeBondKE;
2028 #ifdef NODEGROUP_FORCE_REGISTER
2029  if(cudaIntegrator){
2030  drudeComKE
2032  drudeBondKE
2034  }else{
2035 #endif
2036  drudeComKE
2038  drudeBondKE
2040 #ifdef NODEGROUP_FORCE_REGISTER
2041  }
2042 #endif
2043  int g_bond = 3 * molecule->numDrudeAtoms;
2044  int g_com = numDegFreedom - g_bond;
2045  temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
2046  drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
2047  }
2048 
2049  // Calculate number of degrees of freedom (controlNumDegFreedom)
2050  if ( simParams->useGroupPressure )
2051  {
2052  controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
2053  if ( ! ( numFixedAtoms || molecule->numConstraints
2054  || simParams->comMove || simParams->langevinOn ) ) {
2055  controlNumDegFreedom -= 1;
2056  }
2057  }
2058  else
2059  {
2061  }
2062  if (simParams->fixCellDims) {
2066  }
2067 
2068  // Calculate pressure tensors using virials
2069  calcPressure(step, minimize,
2070  virial_normal, virial_nbond, virial_slow,
2071  intVirial_normal, intVirial_nbond, intVirial_slow,
2072  extForce_normal, extForce_nbond, extForce_slow);
2073 
2074 #ifdef DEBUG_PRESSURE
2075  iout << iINFO << "Control pressure = " << controlPressure <<
2076  " = " << controlPressure_normal << " + " <<
2077  controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
2078 #endif
2080  iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
2081  << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
2082  << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
2083  << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
2084  << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
2085  << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
2086  << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
2087  << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
2088  << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
2089  << endi;
2090  }
2091 }
2092 
2093 //
2094 // Calculates all pressure tensors using virials
2095 //
2096 // Sets variables:
2097 // pressure, pressure_normal, pressure_nbond, pressure_slow
2098 // groupPressure, groupPressure_normal, groupPressure_nbond, groupPressure_slow
2099 // controlPressure, controlPressure_normal, controlPressure_nbond, controlPressure_slow
2100 // pressure_amd
2101 //
2102 void Controller::calcPressure(int step, int minimize,
2103  const Tensor& virial_normal_in, const Tensor& virial_nbond_in, const Tensor& virial_slow_in,
2104  const Tensor& intVirial_normal, const Tensor& intVirial_nbond, const Tensor& intVirial_slow,
2105  const Vector& extForce_normal, const Vector& extForce_nbond, const Vector& extForce_slow) {
2106 
2107  Tensor virial_normal = virial_normal_in;
2108  Tensor virial_nbond = virial_nbond_in;
2109  Tensor virial_slow = virial_slow_in;
2110 
2111 #if 0
2112  Tensor tmp = virial_normal;
2113  fprintf(stderr, "VIRIAL_NORMAL %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2114  tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
2115 
2116  tmp = virial_nbond;
2117  fprintf(stderr, "VIRIAL_NBOND %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2118  tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
2119 
2120  tmp = virial_slow;
2121  fprintf(stderr, "VIRIAL_SLOW %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2122  tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
2123 
2124  tmp = intVirial_normal;
2125  fprintf(stderr, "INT_VIRIAL_NORMAL %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2126  tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
2127 
2128  tmp = intVirial_nbond;
2129  fprintf(stderr, "INT_VIRIAL_NBOND %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2130  tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
2131 
2132  tmp = intVirial_slow;
2133  fprintf(stderr, "INT_VIRIAL_SLOW %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
2134  tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
2135 
2136  Vector ext = extForce_normal;
2137  fprintf(stderr, "EXT_FORCE_NORMAL %1.2lf %1.2lf %1.2lf\n",
2138  ext.x, ext.y, ext.z);
2139 
2140  ext = extForce_nbond;
2141  fprintf(stderr, "EXT_FORCE_NBOND %1.2lf %1.2lf %1.2lf\n",
2142  ext.x, ext.y, ext.z);
2143 
2144  ext = extForce_slow;
2145  fprintf(stderr, "EXT_FORCE_SLOW %1.2lf %1.2lf %1.2lf\n",
2146  ext.x, ext.y, ext.z);
2147 
2148  //fprintf(stderr, "State Lattice for timestep %d\n", step);
2149  fprintf(stderr, "LATTICE %lf %lf %lf %lf %lf %lf %lf %lf %lf\n",
2150  state->lattice.a().x, state->lattice.a().y, state->lattice.a().z,
2151  state->lattice.b().x, state->lattice.b().y, state->lattice.b().z,
2152  state->lattice.c().x, state->lattice.c().y, state->lattice.c().z);
2153 
2154  fprintf(stderr, "VOLUME %3.10lf\n", state->lattice.volume());
2155 
2156 #endif
2157 
2158 
2159  Node *node = Node::Object();
2160  Molecule *molecule = node->molecule;
2161  Lattice &lattice = state->lattice;
2162  BigReal volume;
2163 
2164  Vector extPosition = lattice.origin();
2165  virial_normal -= outer(extForce_normal,extPosition);
2166  virial_nbond -= outer(extForce_nbond,extPosition);
2167  virial_slow -= outer(extForce_slow,extPosition);
2168 
2169  if ( (volume=lattice.volume()) != 0. )
2170  {
2171 
2172  if ((simParams->LJcorrection || simParams->LJcorrectionAlt) && volume) {
2173 #ifdef MEM_OPT_VERSION
2174  NAMD_bug("LJcorrection not supported in memory optimized build.");
2175 #else
2176  // Apply tail correction to pressure
2177  BigReal alchLambda = simParams->getCurrentLambda(step);
2178  virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
2179 #endif
2180  }
2181 
2182 
2183  // kinetic energy component included in virials
2184  if(simParams->CUDASOAintegrate && step != 0 && !simParams->isMultiTimeStepping()){
2185  pressure_normal = virial_normal / volume;
2186  // groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
2187  if (simParams->accelMDOn) {
2188  pressure_amd = virial_amd / volume;
2191  }
2192 
2193  if ( minimize || ! ( step % nbondFreq ) )
2194  {
2195  pressure_nbond = virial_nbond / volume;
2196  // groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
2197  }
2198 
2199  if ( minimize || ! ( step % slowFreq ) )
2200  {
2201  pressure_slow = virial_slow / volume;
2202  // groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
2203  }
2204 
2206  groupPressure = pressure - intVirial_normal / volume;
2207  // groupPressure = groupPressure_normal + groupPressure_nbond +
2208  // groupPressure_slow;
2209  }else{
2210  pressure_normal = virial_normal / volume;
2211 
2212  groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
2213  if (simParams->accelMDOn) {
2214  pressure_amd = virial_amd / volume;
2217  }
2218 
2219  if ( minimize || ! ( step % nbondFreq ) )
2220  {
2221  pressure_nbond = virial_nbond / volume;
2222  groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
2223  }
2224 
2225  if ( minimize || ! ( step % slowFreq ) )
2226  {
2227  pressure_slow = virial_slow / volume;
2228  groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
2229  }
2230 
2234  }
2235  }
2236  else
2237  {
2238  pressure = Tensor();
2239  groupPressure = Tensor();
2240  }
2241 
2242  if ( simParams->useGroupPressure )
2243  {
2248  }
2249  else
2250  {
2255  }
2256 
2257  if ( simParams->useFlexibleCell ) {
2258  // use symmetric pressure to control rotation
2259  // controlPressure_normal = symmetric(controlPressure_normal);
2260  // controlPressure_nbond = symmetric(controlPressure_nbond);
2261  // controlPressure_slow = symmetric(controlPressure_slow);
2262  // controlPressure = symmetric(controlPressure);
2263  // only use on-diagonal components for now
2268  if ( simParams->useConstantRatio ) {
2269 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
2270  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
2275 #undef AVGXY
2276  }
2277  } else {
2284  controlPressure =
2285  Tensor::identity(trace(controlPressure)/3.);
2286  }
2287 }
2288 
2290 (BigReal testV, int n,
2291  BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV){
2292  BigReal delta;
2293 
2294  if(testV > *Vmax){
2295  *Vmax = testV;
2296  }
2297  else if(testV < *Vmin){
2298  *Vmin = testV;
2299  }
2300 
2301  //mean and std calculated by Online Algorithm
2302  delta = testV - *Vavg;
2303  *Vavg += delta / (BigReal)n;
2304  *M2 += delta * (testV - (*Vavg));
2305 
2306  *sigmaV = sqrt(*M2/n);
2307 }
2308 
2310 (int iE, int step, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV,
2311  BigReal* k0, BigReal* k, BigReal* E, int *iEused, char *warn){
2312  switch(iE){
2313  case 2:
2314  *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
2315  // if k0 <=0 OR k0 > 1, switch to iE=1 mode
2316  if(*k0 > 1.)
2317  *warn |= 1;
2318  else if(*k0 <= 0.)
2319  *warn |= 2;
2320  //else stay in iE=2 mode
2321  else{
2322  *E = Vmin + (Vmax-Vmin)/(*k0);
2323  *iEused = 2;
2324  break;
2325  }
2326  case 1:
2327  *E = Vmax;
2328  *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
2329  if(*k0 > 1.)
2330  *k0 = 1.;
2331  *iEused = 1;
2332  break;
2333  }
2334 
2335  *k = *k0 / (Vmax - Vmin);
2336 }
2337 
2339 (BigReal k, BigReal E, BigReal testV, Tensor vir_orig,
2340  BigReal *dV, BigReal *factor, Tensor *vir){
2341  BigReal VE = testV - E;
2342  //if V < E, apply boost
2343  if(VE < 0){
2344  *factor = k * VE;
2345  *vir += vir_orig * (*factor);
2346  *dV += (*factor) * VE/2.;
2347  *factor += 1.;
2348  }
2349  else{
2350  *factor = 1.;
2351  }
2352 }
2353 
2355 (int step_n, char type, int V_n, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal M2, BigReal E, BigReal k, bool write_topic, bool lasttime){
2356  FILE *rest;
2357  char timestepstr[20];
2358  char *filename;
2359  int baselen;
2360  char *restart_name;
2361  const char *bsuffix;
2362 
2363  if(lasttime || simParams->restartFrequency == 0){
2364  filename = simParams->outputFilename;
2365  bsuffix = ".BAK";
2366  }
2367  else{
2368  filename = simParams->restartFilename;
2369  bsuffix = ".old";
2370  }
2371 
2372  baselen = strlen(filename);
2373  restart_name = new char[baselen+26];
2374 
2375  strcpy(restart_name, filename);
2376  if ( simParams->restartSave && !lasttime) {
2377  sprintf(timestepstr,".%d",step_n);
2378  strcat(restart_name, timestepstr);
2379  bsuffix = ".BAK";
2380  }
2381  strcat(restart_name, ".gamd");
2382 
2383  if(write_topic){
2384  NAMD_backup_file(restart_name,bsuffix);
2385 
2386  rest = fopen(restart_name, "w");
2387  if(rest){
2388  fprintf(rest, "# NAMD accelMDG restart file\n"
2389  "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
2390  "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
2391  type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
2392  fclose(rest);
2393  iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
2394  }
2395  else
2396  iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
2397  }
2398  else{
2399  rest = fopen(restart_name, "a");
2400  if(rest){
2401  fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
2402  type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
2403  fclose(rest);
2404  }
2405  else
2406  iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
2407  }
2408 
2409  delete[] restart_name;
2410 }
2411 
2412 
2413 void Controller::rescaleaccelMD(int step, int minimize)
2414 {
2415  if ( !simParams->accelMDOn ) return;
2416 
2418 
2419  // copy all to submit_reduction
2420  for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
2422  }
2424 
2425  if (step == simParams->firstTimestep) {
2426  accelMDdVAverage = 0;
2427  accelMDdV = 0;
2428  }
2429 // if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
2430  if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
2431 
2432  Node *node = Node::Object();
2433  Molecule *molecule = node->molecule;
2434  Lattice &lattice = state->lattice;
2435 
2436  const BigReal accelMDE = simParams->accelMDE;
2437  const BigReal accelMDalpha = simParams->accelMDalpha;
2438  const BigReal accelMDTE = simParams->accelMDTE;
2439  const BigReal accelMDTalpha = simParams->accelMDTalpha;
2440  const int accelMDOutFreq = simParams->accelMDOutFreq;
2441 
2442  //GaMD
2443  static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
2444  static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
2445  static BigReal k0P, kP, EP;
2446  static BigReal k0D, kD, ED;
2447  static int V_n = 1, iEusedD, iEusedP;
2448  static char warnD, warnP;
2449  static int restfreq;
2450 
2453  const int ntcmd = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
2454  const int ntebprep = ntcmd + simParams->accelMDGEquiPrepSteps;
2455  const int nteb = ntcmd + simParams->accelMDGEquiSteps;
2456  const int ntave = simParams->accelMDGStatWindow;
2457 
2458  BigReal volume;
2459  BigReal bondEnergy;
2460  BigReal angleEnergy;
2461  BigReal dihedralEnergy;
2462  BigReal improperEnergy;
2463  BigReal crosstermEnergy;
2464  BigReal boundaryEnergy;
2465  BigReal miscEnergy;
2466  BigReal amd_electEnergy;
2467  BigReal amd_ljEnergy;
2468  BigReal amd_ljEnergy_Corr = 0.;
2469  BigReal amd_electEnergySlow;
2470  BigReal amd_groLJEnergy;
2471  BigReal amd_groGaussEnergy;
2472  BigReal amd_goTotalEnergy;
2473  BigReal potentialEnergy;
2474  BigReal factor_dihe = 1;
2475  BigReal factor_tot = 1;
2476  BigReal testV;
2477  BigReal dV = 0.;
2478  Vector accelMDfactor;
2479  Tensor vir; //auto initialized to 0
2480  Tensor vir_dihe;
2481  Tensor vir_normal;
2482  Tensor vir_nbond;
2483  Tensor vir_slow;
2484 
2485  volume = lattice.volume();
2486 
2487  bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
2488  angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
2489  dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
2490  improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
2491  crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
2492  boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
2493  miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
2494 
2495  GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
2496  GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
2497  GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
2498  GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
2499 
2500  if ( !( step % nbondFreq ) ) {
2501  amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
2502  amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
2503  amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
2504  amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
2507  amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
2508  } else {
2509  amd_electEnergy = electEnergy;
2510  amd_ljEnergy = ljEnergy;
2511  amd_groLJEnergy = groLJEnergy;
2512  amd_groGaussEnergy = groGaussEnergy;
2513  amd_goTotalEnergy = goTotalEnergy;
2514  }
2515 
2516  if( (simParams->LJcorrection || simParams->LJcorrectionAlt) && volume ) {
2517  // Obtain tail correction (copied from printEnergies())
2518  // This value is only printed for sanity check
2519  // accelMD currently does not 'boost' tail correction
2520  amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
2521  }
2522 
2523  if ( !( step % slowFreq ) ) {
2524  amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
2525  } else {
2526  amd_electEnergySlow = electEnergySlow;
2527  }
2528 
2529  potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
2530  improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
2531  crosstermEnergy + boundaryEnergy + miscEnergy +
2532  amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
2533 
2534  //GaMD
2535  if(simParams->accelMDG){
2536  // if first time running accelMDG module
2537  if(step == firststep) {
2538  iEusedD = iEusedP = simParams->accelMDGiE;
2539  warnD = warnP = '\0';
2540 
2541  //restart file reading
2543  FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
2544  char line[256];
2545  int dihe_n=0, tot_n=0;
2546  if(!rest){
2547  sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
2548  NAMD_die(line);
2549  }
2550 
2551  while(fgets(line, 256, rest)){
2552  if(line[0] == 'D'){
2553  dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
2554  &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
2555  }
2556  else if(line[0] == 'T'){
2557  tot_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
2558  &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
2559  }
2560  }
2561 
2562  fclose(rest);
2563 
2564  V_n++;
2565 
2566  //check dihe settings
2568  if(dihe_n !=8)
2569  NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
2570  k0D = kD * (VmaxD - VminD);
2571  iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
2572  << " Vmax " << VmaxD << " Vmin " << VminD
2573  << " Vavg " << VavgD << " sigmaV " << sigmaVD
2574  << " E " << ED << " k " << kD
2575  << "\n" << endi;
2576  }
2577 
2578  //check tot settings
2580  if(tot_n !=8)
2581  NAMD_die("Invalid total potential energy format in accelMDG restart file");
2582  k0P = kP * (VmaxP - VminP);
2583  iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
2584  << " Vmax " << VmaxP << " Vmin " << VminP
2585  << " Vavg " << VavgP << " sigmaV " << sigmaVP
2586  << " E " << EP << " k " << kP
2587  << "\n" << endi;
2588  }
2589 
2590  iEusedD = (ED == VmaxD) ? 1 : 2;
2591  iEusedP = (EP == VmaxP) ? 1 : 2;
2592  }
2593  //local restfreq follows NAMD restartfreq (default: 500)
2595  }
2596 
2597  //for dihedral potential
2599  testV = dihedralEnergy + crosstermEnergy;
2600 
2601  //write restart file every restartfreq steps
2602  if(step > firststep && step % restfreq == 0)
2603  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2604  true, false);
2605  //write restart file at the end of the simulation
2606  if( simParams->accelMDLastStep > 0 ){
2607  if( step == simParams->accelMDLastStep )
2608  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2609  true, true);
2610  }
2611  else if(step == simParams->N)
2612  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2613  true, true);
2614 
2615  //conventional MD
2616  if(step < ntcmd){
2617  //very first step
2618  if(step == firststep && !simParams->accelMDGRestart){
2619  //initialize stat
2620  VmaxD = VminD = VavgD = testV;
2621  M2D = sigmaVD = 0.;
2622  }
2623  //first step after cmdprep
2624  else if(step == ntcmdprep && ntcmdprep != 0){
2625  //reset stat
2626  VmaxD = VminD = VavgD = testV;
2627  M2D = sigmaVD = 0.;
2628  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
2629  }
2630  //every ntave steps
2631  else if(ntave > 0 && step % ntave == 0){
2632  //update Vmax Vmin
2633  if(testV > VmaxD) VmaxD = testV;
2634  if(testV < VminD) VminD = testV;
2635  //reset avg and std
2636  VavgD = testV;
2637  M2D = sigmaVD = 0.;
2638  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
2639  }
2640  //normal steps
2641  else
2642  calc_accelMDG_mean_std(testV, V_n,
2643  &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2644 
2645  //last cmd step
2646  if(step == ntcmd - 1){
2648  VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2649  }
2650  }
2651  //equilibration
2652  else if(step < nteb){
2653  calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
2654  &dV, &factor_dihe, &vir);
2655 
2656  //first step after cmd
2657  if(step == ntcmd && simParams->accelMDGresetVaftercmd){
2658  //reset stat
2659  VmaxD = VminD = VavgD = testV;
2660  M2D = sigmaVD = 0.;
2661  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
2662  }
2663  //every ntave steps
2664  else if(ntave > 0 && step % ntave == 0){
2665  //update Vmax Vmin
2666  if(testV > VmaxD) VmaxD = testV;
2667  if(testV < VminD) VminD = testV;
2668  //reset avg and std
2669  VavgD = testV;
2670  M2D = sigmaVD = 0.;
2671  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
2672  }
2673  else
2674  calc_accelMDG_mean_std(testV, V_n,
2675  &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2676 
2677  //steps after ebprep
2678  if(step >= ntebprep)
2679  if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2681  VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2682  }
2683  //production
2684  else{
2685  calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
2686  &dV, &factor_dihe, &vir);
2687  }
2688 
2689  }
2690  //for total potential
2692  Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
2693  testV = potentialEnergy;
2694  if(simParams->accelMDdual){
2695  testV -= dihedralEnergy + crosstermEnergy;
2696  vir_tot -= vir_dihe;
2697  }
2698 
2699  //write restart file every restartfreq steps
2700  if(step > firststep && step % restfreq == 0)
2701  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2702  !simParams->accelMDdual, false);
2703  //write restart file at the end of simulation
2704  if( simParams->accelMDLastStep > 0 ){
2705  if( step == simParams->accelMDLastStep )
2706  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2707  !simParams->accelMDdual, true);
2708  }
2709  else if(step == simParams->N)
2710  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2711  !simParams->accelMDdual, true);
2712 
2713  //conventional MD
2714  if(step < ntcmd){
2715  //very first step
2716  if(step == firststep && !simParams->accelMDGRestart){
2717  //initialize stat
2718  VmaxP = VminP = VavgP = testV;
2719  M2P = sigmaVP = 0.;
2720  }
2721  //first step after cmdprep
2722  else if(step == ntcmdprep && ntcmdprep != 0){
2723  //reset stat
2724  VmaxP = VminP = VavgP = testV;
2725  M2P = sigmaVP = 0.;
2726  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
2727  }
2728  //every ntave steps
2729  else if(ntave > 0 && step % ntave == 0){
2730  //update Vmax Vmin
2731  if(testV > VmaxP) VmaxP = testV;
2732  if(testV < VminP) VminP = testV;
2733  //reset avg and std
2734  VavgP = testV;
2735  M2P = sigmaVP = 0.;
2736  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
2737  }
2738  //normal steps
2739  else
2740  calc_accelMDG_mean_std(testV, V_n,
2741  &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2742  //last cmd step
2743  if(step == ntcmd - 1){
2745  VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2746  }
2747  }
2748  //equilibration
2749  else if(step < nteb){
2750  calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
2751  &dV, &factor_tot, &vir);
2752 
2753  //first step after cmd
2754  if(step == ntcmd && simParams->accelMDGresetVaftercmd){
2755  //reset stat
2756  VmaxP = VminP = VavgP = testV;
2757  M2P = sigmaVP = 0.;
2758  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
2759  }
2760  //every ntave steps
2761  else if(ntave > 0 && step % ntave == 0){
2762  //update Vmax Vmin
2763  if(testV > VmaxP) VmaxP = testV;
2764  if(testV < VminP) VminP = testV;
2765  //reset avg and std
2766  VavgP = testV;
2767  M2P = sigmaVP = 0.;
2768  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
2769  }
2770  else
2771  calc_accelMDG_mean_std(testV, V_n,
2772  &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2773 
2774  //steps after ebprep
2775  if(step >= ntebprep)
2776  if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2778  VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2779  }
2780  //production
2781  else{
2782  calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
2783  &dV, &factor_tot, &vir);
2784  }
2785 
2786  }
2787  accelMDdVAverage += dV;
2788 
2789  //first step after ntcmdprep
2790  if((ntcmdprep > 0 && step == ntcmdprep) ||
2791  (step < nteb && ntave > 0 && step % ntave == 0) ||
2792  (simParams->accelMDGresetVaftercmd && step == ntcmd)){
2793  V_n = 1;
2794  }
2795 
2796  if(step < nteb)
2797  V_n++;
2798 
2799  }
2800  //aMD
2801  else{
2802  if (simParams->accelMDdihe) {
2803 
2804  testV = dihedralEnergy + crosstermEnergy;
2805  if ( testV < accelMDE ) {
2806  factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2807  factor_dihe *= factor_dihe;
2808  vir = vir_dihe * (factor_dihe - 1.0);
2809  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2810  accelMDdVAverage += dV;
2811  }
2812 
2813  } else if (simParams->accelMDdual) {
2814 
2815  testV = dihedralEnergy + crosstermEnergy;
2816  if ( testV < accelMDE ) {
2817  factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2818  factor_dihe *= factor_dihe;
2819  vir = vir_dihe * (factor_dihe - 1.0) ;
2820  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2821  }
2822 
2823  testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
2824  if ( testV < accelMDTE ) {
2825  factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
2826  factor_tot *= factor_tot;
2827  vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
2828  dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
2829  }
2830  accelMDdVAverage += dV;
2831 
2832  } else {
2833 
2834  testV = potentialEnergy;
2835  if ( testV < accelMDE ) {
2836  factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
2837  factor_tot *= factor_tot;
2838  vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
2839  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2840  accelMDdVAverage += dV;
2841  }
2842  }
2843  }
2844 
2845  accelMDfactor[0]=factor_dihe;
2846  accelMDfactor[1]=factor_tot;
2847  accelMDfactor[2]=1;
2848  broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
2849  virial_amd = vir;
2850  accelMDdV = dV;
2851 
2852  if ( factor_tot < 0.001 ) {
2853  iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
2854  << "\n" << endi;
2855  }
2856 
2857  if ( ! (step % accelMDOutFreq) ) {
2858  if ( !(step == simParams->firstTimestep) ) {
2859  accelMDdVAverage = accelMDdVAverage/accelMDOutFreq;
2860  }
2861  iout << "ACCELERATED MD: STEP " << step
2862  << " dV " << dV
2863  << " dVAVG " << accelMDdVAverage
2864  << " BOND " << bondEnergy
2865  << " ANGLE " << angleEnergy
2866  << " DIHED " << dihedralEnergy+crosstermEnergy
2867  << " IMPRP " << improperEnergy
2868  << " ELECT " << amd_electEnergy+amd_electEnergySlow
2869  << " VDW " << amd_ljEnergy
2870  << " POTENTIAL " << potentialEnergy;
2871  if(amd_ljEnergy_Corr)
2872  iout << " LJcorr " << amd_ljEnergy_Corr;
2873  iout << "\n" << endi;
2874  if(simParams->accelMDG){
2876  iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD
2877  << " Vmax " << VmaxD << " Vmin " << VminD
2878  << " Vavg " << VavgD << " sigmaV " << sigmaVD
2879  << " E " << ED << " k0 " << k0D << " k " << kD
2880  << "\n" << endi;
2881  if(warnD & 1)
2882  iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
2883  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2884  if(warnD & 2)
2885  iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
2886  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2888  iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP
2889  << " Vmax " << VmaxP << " Vmin " << VminP
2890  << " Vavg " << VavgP << " sigmaV " << sigmaVP
2891  << " E " << EP << " k0 " << k0P << " k " << kP
2892  << "\n" << endi;
2893  if(warnP & 1)
2894  iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
2895  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2896  if(warnP & 2)
2897  iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
2898  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2899  warnD = warnP = '\0';
2900  }
2901 
2902  accelMDdVAverage = 0;
2903 
2904  if (simParams->accelMDDebugOn) {
2905  Tensor p_normal;
2906  Tensor p_nbond;
2907  Tensor p_slow;
2908  Tensor p;
2909  if ( volume != 0. ) {
2910  p_normal = vir_normal/volume;
2911  p_nbond = vir_nbond/volume;
2912  p_slow = vir_slow/volume;
2913  p = vir/volume;
2914  }
2915  iout << " accelMD Scaling Factor: " << accelMDfactor << "\n"
2916  << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
2917  << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
2918  << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
2919  << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n"
2920  << endi;
2921  }
2922  }
2923 }
2924 
2926  if (!simParams->adaptTempOn) return;
2927  iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
2928  adaptTempDtMin = 0;
2929  adaptTempDtMax = 0;
2930  adaptTempAutoDt = false;
2931  if (simParams->adaptTempBins == 0) {
2932  iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
2933  std::ifstream adaptTempRead(simParams->adaptTempInFile);
2934  if (adaptTempRead) {
2935  int readInt;
2936  BigReal readReal;
2937  bool readBool;
2938  // step
2939  adaptTempRead >> readInt;
2940  // Start with min and max temperatures
2941  adaptTempRead >> adaptTempT; // KELVIN
2942  adaptTempRead >> adaptTempBetaMin; // KELVIN
2943  adaptTempRead >> adaptTempBetaMax; // KELVIN
2944  adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
2945  adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
2946  // In case file is manually edited
2948  readReal = adaptTempBetaMax;
2951  }
2952  adaptTempRead >> adaptTempBins;
2953  adaptTempRead >> adaptTempCg;
2954  adaptTempRead >> adaptTempDt;
2963  for(int j = 0; j < adaptTempBins; ++j) {
2964  adaptTempRead >> adaptTempPotEnergyAve[j];
2965  adaptTempRead >> adaptTempPotEnergyVar[j];
2966  adaptTempRead >> adaptTempPotEnergySamples[j];
2967  adaptTempRead >> adaptTempPotEnergyAveNum[j];
2968  adaptTempRead >> adaptTempPotEnergyVarNum[j];
2969  adaptTempRead >> adaptTempPotEnergyAveDen[j];
2971  }
2972  adaptTempRead.close();
2973  }
2974  else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
2975  }
2976  else {
2991  if (simParams->langevinOn)
2993  else if (simParams->rescaleFreq > 0)
2995  for(int j = 0; j < adaptTempBins; ++j){
2996  adaptTempPotEnergyAveNum[j] = 0.;
2997  adaptTempPotEnergyAveDen[j] = 0.;
2999  adaptTempPotEnergyVarNum[j] = 0.;
3000  adaptTempPotEnergyVar[j] = 0.;
3001  adaptTempPotEnergyAve[j] = 0.;
3003  }
3004  }
3005  if (simParams->adaptTempAutoDt > 0.0) {
3006  adaptTempAutoDt = true;
3009  }
3010  adaptTempDTave = 0;
3011  adaptTempDTavenum = 0;
3012  iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
3013  iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
3014  iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
3015  iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
3016  iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
3017  if (simParams->adaptTempRestartFreq > 0) {
3021  NAMD_die("Error opening restart file");
3022  }
3023 }
3024 
3027  adaptTempRestartFile.seekp(std::ios::beg);
3028  iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
3029  adaptTempRestartFile << step << " ";
3030  // Start with min and max temperatures
3031  adaptTempRestartFile << adaptTempT<< " "; // KELVIN
3032  adaptTempRestartFile << 1./adaptTempBetaMin << " "; // KELVIN
3033  adaptTempRestartFile << 1./adaptTempBetaMax << " "; // KELVIN
3037  adaptTempRestartFile << "\n" ;
3038  for(int j = 0; j < adaptTempBins; ++j) {
3045  adaptTempRestartFile << "\n";
3046  }
3048  }
3049 }
3050 
3051 void Controller::adaptTempUpdate(int step, int minimize)
3052 {
3053  //Beta = 1./T
3054  if ( !simParams->adaptTempOn ) return;
3055  int j = 0;
3056  if (step == simParams->firstTimestep) {
3057  adaptTempInit(step);
3058  return;
3059  }
3060  if ( minimize || (step < simParams->adaptTempFirstStep ) ||
3061  ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
3062  const int adaptTempOutFreq = simParams->adaptTempOutFreq;
3063  const bool adaptTempDebug = simParams->adaptTempDebug;
3064  //Calculate Current inverse temperature and bin
3065  BigReal adaptTempBeta = 1./adaptTempT;
3066  adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
3067 
3068  if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
3069  iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
3070  << " adaptTempBeta: " << adaptTempBeta
3071  << " adaptTempDBeta: " << adaptTempDBeta
3072  << " betaMin:" << adaptTempBetaMin
3073  << " betaMax: " << adaptTempBetaMax << "\n";
3076 
3077  BigReal potentialEnergy;
3078  BigReal potEnergyAverage;
3079  BigReal potEnergyVariance;
3080  potentialEnergy = totalEnergy-kineticEnergy;
3081 
3082  //calculate new bin average and variance using adaptive averaging
3085  adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
3086 
3089  potEnergyVariance -= potEnergyAverage*potEnergyAverage;
3090 
3091  adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
3092  adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
3093 
3094  // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
3095  // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
3096  // is constant for each bin. This is to estimate <E(beta)> where beta \in
3097  // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
3098  if ( ! ( step % simParams->adaptTempFreq ) ) {
3099  // If adaptTempBin not at the edge, calculate improved average:
3100  if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
3101  // Get Averaging Limits:
3102  BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
3103  BigReal betaPlus;
3104  BigReal betaMinus;
3105  int nMinus =0;
3106  int nPlus = 0;
3107  if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
3108  if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
3109  betaMinus = adaptTempBeta - deltaBeta;
3110  betaPlus = adaptTempBeta + deltaBeta;
3111  BigReal betaMinus2 = betaMinus*betaMinus;
3112  BigReal betaPlus2 = betaPlus*betaPlus;
3113  nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
3114  nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
3115  // Variables for <E(beta)> estimate:
3116  BigReal potEnergyAve0 = 0.0;
3117  BigReal potEnergyAve1 = 0.0;
3118  // Integral terms
3119  BigReal A0 = 0;
3120  BigReal A1 = 0;
3121  BigReal A2 = 0;
3122  //A0 phi_s integral for beta_minus < beta < beta_{i+1}
3123  BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
3124  BigReal DeltaE2Ave;
3125  BigReal DeltaE2AveSum = 0;
3126  BigReal betaj;
3127  for (j = nMinus+1; j <= adaptTempBin; ++j) {
3128  potEnergyAve0 += adaptTempPotEnergyAve[j];
3129  DeltaE2Ave = adaptTempPotEnergyVar[j];
3130  DeltaE2AveSum += DeltaE2Ave;
3131  A0 += j*DeltaE2Ave;
3132  }
3133  A0 *= adaptTempDBeta;
3134  A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
3135  A0 *= adaptTempDBeta;
3136  betaj = adaptTempBetaN[nMinus+1]-betaMinus;
3137  betaj *= betaj;
3138  A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
3139  A0 /= (betaNp1-betaMinus);
3140 
3141  //A1 phi_s integral for beta_{i+1} < beta < beta_plus
3142  DeltaE2AveSum = 0;
3143  for (j = adaptTempBin+1; j < nPlus; j++) {
3144  potEnergyAve1 += adaptTempPotEnergyAve[j];
3145  DeltaE2Ave = adaptTempPotEnergyVar[j];
3146  DeltaE2AveSum += DeltaE2Ave;
3147  A1 += j*DeltaE2Ave;
3148  }
3149  A1 *= adaptTempDBeta;
3150  A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
3151  A1 *= adaptTempDBeta;
3152  if ( nPlus < adaptTempBins ) {
3153  betaj = betaPlus-adaptTempBetaN[nPlus];
3154  betaj *= betaj;
3155  A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
3156  }
3157  A1 /= (betaPlus-betaNp1);
3158 
3159  //A2 phi_t integral for beta_i
3160  A2 = 0.5*adaptTempDBeta*potEnergyVariance;
3161 
3162  // Now calculate a+ and a-
3163  DeltaE2Ave = A0-A1;
3164  BigReal aplus = 0;
3165  BigReal aminus = 0;
3166  if (DeltaE2Ave != 0) {
3167  aplus = (A0-A2)/(A0-A1);
3168  if (aplus < 0) {
3169  aplus = 0;
3170  }
3171  if (aplus > 1) {
3172  aplus = 1;
3173  }
3174  aminus = 1-aplus;
3175  potEnergyAve0 *= adaptTempDBeta;
3176  potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
3177  potEnergyAve0 /= (betaNp1-betaMinus);
3178  potEnergyAve1 *= adaptTempDBeta;
3179  if ( nPlus < adaptTempBins ) {
3180  potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
3181  }
3182  potEnergyAve1 /= (betaPlus-betaNp1);
3183  potEnergyAverage = aminus*potEnergyAve0;
3184  potEnergyAverage += aplus*potEnergyAve1;
3185  }
3186  if (simParams->adaptTempDebug) {
3187  iout << "ADAPTEMP DEBUG:" << "\n"
3188  << " adaptTempBin: " << adaptTempBin << "\n"
3189  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
3190  << " adaptTempBeta: " << adaptTempBeta << "\n"
3191  << " adaptTemp: " << adaptTempT<< "\n"
3192  << " betaMin: " << adaptTempBetaMin << "\n"
3193  << " betaMax: " << adaptTempBetaMax << "\n"
3194  << " gammaAve: " << gammaAve << "\n"
3195  << " deltaBeta: " << deltaBeta << "\n"
3196  << " betaMinus: " << betaMinus << "\n"
3197  << " betaPlus: " << betaPlus << "\n"
3198  << " nMinus: " << nMinus << "\n"
3199  << " nPlus: " << nPlus << "\n"
3200  << " A0: " << A0 << "\n"
3201  << " A1: " << A1 << "\n"
3202  << " A2: " << A2 << "\n"
3203  << " a+: " << aplus << "\n"
3204  << " a-: " << aminus << "\n"
3205  << endi;
3206  }
3207  }
3208  else {
3209  if (simParams->adaptTempDebug) {
3210  iout << "ADAPTEMP DEBUG:" << "\n"
3211  << " adaptTempBin: " << adaptTempBin << "\n"
3212  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
3213  << " adaptTempBeta: " << adaptTempBeta << "\n"
3214  << " adaptTemp: " << adaptTempT<< "\n"
3215  << " betaMin: " << adaptTempBetaMin << "\n"
3216  << " betaMax: " << adaptTempBetaMax << "\n"
3217  << " gammaAve: " << gammaAve << "\n"
3218  << " deltaBeta: N/A\n"
3219  << " betaMinus: N/A\n"
3220  << " betaPlus: N/A\n"
3221  << " nMinus: N/A\n"
3222  << " nPlus: N/A\n"
3223  << " A0: N/A\n"
3224  << " A1: N/A\n"
3225  << " A2: N/A\n"
3226  << " a+: N/A\n"
3227  << " a-: N/A\n"
3228  << endi;
3229  }
3230  }
3231 
3232  //dT is new temperature
3233  BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
3234  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
3235  dT += adaptTempT;
3236 
3237  // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
3238  // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
3239  if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
3241  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
3242  dT += adaptTempT;
3243  // Check again, if not then keep original adaptTempTor assign random.
3244  if ( dT > 1./adaptTempBetaMin ) {
3245  if (!simParams->adaptTempRandom) {
3246  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
3247  // << " K higher than adaptTempTmax."
3248  // << " Keeping temperature at "
3249  // << adaptTempT<< "\n"<< endi;
3250  dT = adaptTempT;
3251  }
3252  else {
3253  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
3254  // << " K higher than adaptTempTmax."
3255  // << " Assigning random temperature in range\n" << endi;
3257  dT = 1./dT;
3258  }
3259  }
3260  else if ( dT < 1./adaptTempBetaMax ) {
3261  if (!simParams->adaptTempRandom) {
3262  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
3263  // << " K lower than adaptTempTmin."
3264  // << " Keeping temperature at " << adaptTempT<< "\n" << endi;
3265  dT = adaptTempT;
3266  }
3267  else {
3268  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
3269  // << " K lower than adaptTempTmin."
3270  // << " Assigning random temperature in range\n" << endi;
3272  dT = 1./dT;
3273  }
3274  }
3275  else if (adaptTempAutoDt) {
3276  //update temperature step size counter
3277  //FOR "TRUE" ADAPTIVE TEMPERING
3278  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
3279  if (adaptTempTdiff > 0) {
3280  adaptTempDTave += adaptTempTdiff;
3282 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
3283  }
3284  if(adaptTempDTavenum == 100){
3285  BigReal Frac;
3287  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
3288  Frac = adaptTempDTave/Frac;
3289  //if average temperature jump is > 3% of temperature range,
3290  //modify jump size to match 3%
3291  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
3292  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
3293  Frac = adaptTempDtMax/Frac;
3294  iout << "ADAPTEMP: Updating adaptTempDt to ";
3295  adaptTempDt *= Frac;
3296  iout << adaptTempDt << "\n" << endi;
3297  }
3298  adaptTempDTave = 0;
3299  adaptTempDTavenum = 0;
3300  }
3301  }
3302  }
3303  else if (adaptTempAutoDt) {
3304  //update temperature step size counter
3305  // FOR "TRUE" ADAPTIVE TEMPERING
3306  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
3307  if (adaptTempTdiff > 0) {
3308  adaptTempDTave += adaptTempTdiff;
3310 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
3311  }
3312  if(adaptTempDTavenum == 100){
3313  BigReal Frac;
3315  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
3316  Frac = adaptTempDTave/Frac;
3317  //if average temperature jump is > 3% of temperature range,
3318  //modify jump size to match 3%
3319  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
3320  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
3321  Frac = adaptTempDtMax/Frac;
3322  iout << "ADAPTEMP: Updating adaptTempDt to ";
3323  adaptTempDt *= Frac;
3324  iout << adaptTempDt << "\n" << endi;
3325  }
3326  adaptTempDTave = 0;
3327  adaptTempDTavenum = 0;
3328 
3329  }
3330 
3331  }
3332 
3333  adaptTempT = dT;
3335  }
3336  adaptTempWriteRestart(step);
3337  if ( ! (step % adaptTempOutFreq) ) {
3338  iout << "ADAPTEMP: STEP " << step
3339  << " TEMP " << adaptTempT
3340  << " ENERGY " << std::setprecision(10) << potentialEnergy
3341  << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
3342  << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
3343  iout << "\n" << endi;
3344  }
3345 
3346 }
3347 
3348 
3349 void Controller::compareChecksums(int step, int forgiving) {
3350  Node *node = Node::Object();
3351  Molecule *molecule = node->molecule;
3352  bool cudaIntegrator = simParams->CUDASOAintegrate;
3353 
3354  // Some basic correctness checking
3355  BigReal checksum, checksum_b;
3356  char errmsg[256];
3357 
3358 #ifdef NODEGROUP_FORCE_REGISTER
3359  if(cudaIntegrator)
3361  else{
3362 #endif
3363  checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
3364 #ifdef NODEGROUP_FORCE_REGISTER
3365  }
3366 #endif
3367  if ( ((int)checksum) != molecule->numAtoms && step != 0) {
3368  sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
3369  (int)checksum, molecule->numAtoms);
3370  NAMD_bug(errmsg);
3371  }
3372  // do we ever need this if we're on a GPU-based code?
3373 #ifdef NODEGROUP_FORCE_REGISTER
3374  if(cudaIntegrator) checksum = computeChecksum; //checksum = nodeReduction->item(REDUCTION_COMPUTE_CHECKSUM);
3375  else {
3376 #endif
3378 #ifdef NODEGROUP_FORCE_REGISTER
3379  }
3380 #endif
3381 
3382  if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
3383  if ( ! computeChecksum ) {
3384  computesPartitioned = 0;
3385  } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
3386  computesPartitioned = 1;
3387  } else {
3388  NAMD_bug("Bad global compute count!\n");
3389  }
3390  computeChecksum = ((int)checksum);
3391  }
3392 #ifdef NODEGROUP_FORCE_REGISTER
3393  if(cudaIntegrator)checksum_b = checksum = nodeReduction->item(REDUCTION_BOND_CHECKSUM);
3394  else{
3395 #endif
3396  checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
3397 #ifdef NODEGROUP_FORCE_REGISTER
3398  }
3399 #endif
3400  if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
3401  sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
3402  (int)checksum, molecule->numCalcBonds);
3403  if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
3404  iout << iWARN << errmsg << endi;
3405  else NAMD_bug(errmsg);
3406  }
3407 
3408 #ifdef NODEGROUP_FORCE_REGISTER
3409  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_ANGLE_CHECKSUM);
3410 
3411  else{
3412 #endif
3414 
3415 #ifdef NODEGROUP_FORCE_REGISTER
3416  }
3417 #endif
3418  if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
3419  sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
3420  (int)checksum, molecule->numCalcAngles);
3421  if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
3422  iout << iWARN << errmsg << endi;
3423  else NAMD_bug(errmsg);
3424  }
3425 #ifdef NODEGROUP_FORCE_REGISTER
3426  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_DIHEDRAL_CHECKSUM);
3427  else {
3428 #endif
3430 #ifdef NODEGROUP_FORCE_REGISTER
3431  }
3432 #endif
3433  if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
3434  sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
3435  (int)checksum, molecule->numCalcDihedrals);
3436  if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
3437  iout << iWARN << errmsg << endi;
3438  else NAMD_bug(errmsg);
3439  }
3440 
3441 #ifdef NODEGROUP_FORCE_REGISTER
3442  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_IMPROPER_CHECKSUM);
3443  else {
3444 #endif
3446 #ifdef NODEGROUP_FORCE_REGISTER
3447  }
3448 #endif
3449  if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
3450  sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
3451  (int)checksum, molecule->numCalcImpropers);
3452  if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
3453  iout << iWARN << errmsg << endi;
3454  else NAMD_bug(errmsg);
3455  }
3456 
3457 #ifdef NODEGROUP_FORCE_REGISTER
3458  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_THOLE_CHECKSUM);
3459  else{
3460 #endif
3462 #ifdef NODEGROUP_FORCE_REGISTER
3463  }
3464 #endif
3465  if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
3466  sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
3467  (int)checksum, molecule->numCalcTholes);
3468  if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
3469  iout << iWARN << errmsg << endi;
3470  else NAMD_bug(errmsg);
3471  }
3472 
3473 #ifdef NODEGROUP_FORCE_REGISTER
3474  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_ANISO_CHECKSUM);
3475  else {
3476 #endif
3478 #ifdef NODEGROUP_FORCE_REGISTER
3479  }
3480 #endif
3481  if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
3482  sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
3483  (int)checksum, molecule->numCalcAnisos);
3484  if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
3485  iout << iWARN << errmsg << endi;
3486  else NAMD_bug(errmsg);
3487  }
3488 
3489 #ifdef NODEGROUP_FORCE_REGISTER
3490  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_CROSSTERM_CHECKSUM);
3491  else{
3492 #endif
3494 #ifdef NODEGROUP_FORCE_REGISTER
3495  }
3496 #endif
3497  if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
3498  sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
3499  (int)checksum, molecule->numCalcCrossterms);
3500  if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
3501  iout << iWARN << errmsg << endi;
3502  else NAMD_bug(errmsg);
3503  }
3504 
3505 #ifdef NODEGROUP_FORCE_REGISTER
3506  if(cudaIntegrator )checksum = nodeReduction->item(REDUCTION_EXCLUSION_CHECKSUM);
3507  else{
3508 #endif
3510 #ifdef NODEGROUP_FORCE_REGISTER
3511  }
3512 #endif
3513  if ( ((int64)checksum) > molecule->numCalcExclusions &&
3514  ( ! simParams->mollyOn || step % slowFreq ) ) {
3515  if ( forgiving )
3516  iout << iWARN << "High global exclusion count ("
3517  << ((int64)checksum) << " vs "
3518  << molecule->numCalcExclusions << "), possible error!\n"
3519  << iWARN << "This warning is not unusual during minimization.\n"
3520  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
3521  else {
3522  char errmsg[256];
3523  sprintf(errmsg, "High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
3524  (long long)checksum, (long long)molecule->numCalcExclusions);
3525  NAMD_bug(errmsg);
3526  }
3527  }
3528 
3529  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
3531  if ( forgiving || ! simParams->fullElectFrequency ) {
3532  iout << iWARN << "Low global exclusion count! ("
3533  << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
3534  if ( forgiving ) {
3535  iout << "\n"
3536  << iWARN << "This warning is not unusual during minimization.\n"
3537  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
3538  } else {
3539  iout << " System unstable or pairlistdist or cutoff too small.\n";
3540  }
3541  iout << endi;
3542  } else {
3543  char errmsg[256];
3544  sprintf(errmsg, "Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
3545  (long long)checksum, (long long)molecule->numCalcExclusions);
3546  NAMD_bug(errmsg);
3547  }
3548  }
3549 
3550 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
3551 #ifdef NODEGROUP_FORCE_REGISTER
3552  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_EXCLUSION_CHECKSUM_CUDA);
3553  else{
3554 #endif
3556 #ifdef NODEGROUP_FORCE_REGISTER
3557  }
3558 #endif
3559  if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
3560  ( ! simParams->mollyOn || step % slowFreq ) ) {
3561  if ( forgiving )
3562  iout << iWARN << "High global CUDA exclusion count ("
3563  << ((int64)checksum) << " vs "
3564  << molecule->numCalcFullExclusions << "), possible error!\n"
3565  << iWARN << "This warning is not unusual during minimization.\n"
3566  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
3567  else {
3568  char errmsg[256];
3569  sprintf(errmsg, "High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
3570  (long long)checksum, (long long)molecule->numCalcFullExclusions);
3571  NAMD_bug(errmsg);
3572  }
3573  }
3574 
3575  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
3577  if ( forgiving || ! simParams->fullElectFrequency ) {
3578  iout << iWARN << "Low global CUDA exclusion count! ("
3579  << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ") on step " << step;
3580  if ( forgiving ) {
3581  iout << "\n"
3582  << iWARN << "This warning is not unusual during minimization.\n"
3583  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
3584  } else {
3585  iout << " System unstable or pairlistdist or cutoff too small.\n";
3586  }
3587  iout << endi;
3588  } else {
3589  char errmsg[256];
3590  sprintf(errmsg, "Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
3591  (long long)checksum, (long long)molecule->numCalcFullExclusions);
3592  NAMD_bug(errmsg);
3593  }
3594  }
3595 #endif
3596 
3597 #ifdef NODEGROUP_FORCE_REGISTER
3598  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_MARGIN_VIOLATIONS);
3599  else{
3600 #endif
3602 #ifdef NODEGROUP_FORCE_REGISTER
3603  }
3604 #endif
3605  if ( ((int)checksum) && ! marginViolations ) {
3606  iout << iERROR << "Margin is too small for " << ((int)checksum) <<
3607  " atoms during timestep " << step << ".\n" << iERROR <<
3608  "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
3609  }
3610  marginViolations += (int)checksum;
3611 
3612 #ifdef NODEGROUP_FORCE_REGISTER
3613  if(cudaIntegrator) checksum = nodeReduction->item(REDUCTION_PAIRLIST_WARNINGS);
3614  else {
3615 #endif
3617 #ifdef NODEGROUP_FORCE_REGISTER
3618  }
3619 #endif
3620  if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
3621  iout << iINFO <<
3622  "Pairlistdist is too small for " << ((int)checksum) <<
3623  " computes during timestep " << step << ".\n" << endi;
3624  }
3625  if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
3626 
3627 #ifdef NODEGROUP_FORCE_REGISTER
3628  if(cudaIntegrator)checksum = nodeReduction->item(REDUCTION_STRAY_CHARGE_ERRORS);
3629  else{
3630 #endif
3632 #ifdef NODEGROUP_FORCE_REGISTER
3633  }
3634 #endif
3635  if ( checksum ) {
3636  if ( forgiving )
3637  iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
3638  else NAMD_bug("Stray PME grid charges detected!\n");
3639  }
3640 }
3641 
3642 void Controller::printTiming(int step) {
3643 
3644  if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
3645  {
3646  const double endWTime = namdWallTimer() - firstWTime;
3647  const double endCTime = CmiTimer() - firstCTime;
3648 
3649  // fflush about once per minute
3650  if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
3651 
3652  const double elapsedW =
3653  (endWTime - startWTime) / simParams->outputTiming;
3654  const double elapsedC =
3655  (endCTime - startCTime) / simParams->outputTiming;
3656 
3657  const double remainingW = elapsedW * (simParams->N - step);
3658  const double remainingW_hours = remainingW / 3600;
3659 
3660  startWTime = endWTime;
3661  startCTime = endCTime;
3662 
3663 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
3664  if ( simParams->computeEnergies < 60 &&
3665  step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
3666  iout << iWARN << "Energy evaluation is expensive, increase computeEnergies to improve performance.\n" << endi;
3667  iout << iWARN << "If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n" << endi;
3668  }
3669 #endif
3670  if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
3671  if (simParams->nsPerDayOn) {
3672  BigReal ns = simParams->dt / 1000000.0;
3673  BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
3674  BigReal nsPerDay = ns / (elapsedW * days);
3675  CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
3676  ", %g ns/days"
3677  ", %g hours remaining, %f MB of memory in use.\n",
3678  step, endCTime, elapsedC, endWTime, elapsedW,
3679  nsPerDay,
3680  remainingW_hours, memusage_MB());
3681  }
3682  else {
3683  CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
3684  ", %g hours remaining, %f MB of memory in use.\n",
3685  step, endCTime, elapsedC, endWTime, elapsedW,
3686  remainingW_hours, memusage_MB());
3687  }
3689  // calculate running average and sample variance of
3690  // the measured time per step (elapsedW)
3691  perfstats.add_sample(elapsedW);
3692  double t_avg = perfstats.average();
3693  double t_std = perfstats.standard_deviation();
3694  // convert t_avg to average number of nanoseconds simulated per day
3695  double ns_per_day = (simParams->dt / t_avg) * (60 * 60 * 24 * 1e-6);
3696  CmiPrintf("PERFORMANCE: %d averaging %g ns/day, %g sec/step with"
3697  " standard deviation %g\n", step, ns_per_day, t_avg, t_std);
3698  }
3699  if ( fflush_count ) { --fflush_count; fflush(stdout); }
3700  double benchmarkTime = (double) simParams->benchmarkTime;
3701  if (benchmarkTime > 0 && benchmarkTime <= endWTime) {
3702  // terminate NAMD benchmark
3703  char s[100];
3704  sprintf(s, "Exceeded benchmark time limit %g seconds",
3705  benchmarkTime);
3706  NAMD_quit(s);
3707  }
3708  }
3709  }
3710 }
3711 
3713 
3714  rescaleaccelMD(step,1);
3715  receivePressure(step,1);
3716 
3717  Node *node = Node::Object();
3718  Molecule *molecule = node->molecule;
3719  compareChecksums(step,1);
3720 
3721  printEnergies(step,1);
3722 
3728 
3729 #ifdef DEBUG_MINIMIZE
3730  printf("%s, line %d\n", __FILE__, __LINE__);
3731  printf("Step %d:\n", step);
3732  printf(" min_energy = %f\n", min_energy);
3733  printf(" min_f_dot_f = %f\n", min_f_dot_f);
3734  printf(" min_f_dot_v = %f\n", min_f_dot_v);
3735  printf(" min_v_dot_v = %f\n", min_v_dot_v);
3736  printf(" min_huge_count = %d\n", min_huge_count);
3737  printf("\n");
3738 #endif
3739 }
3740 
3742  compareChecksums(step);
3743  printEnergies(step,0);
3744 }
3745 
3747 {
3748  compareChecksums(step, 0);
3749 
3750  Node *node = Node::Object();
3751  Molecule *molecule = node->molecule;
3752  Lattice &lattice = state->lattice;
3753  BigReal volume = lattice.volume();
3754  // store total energy
3755  BigReal totalPotentialEnergy = 0.0;
3756  BigReal bondEnergy, angleEnergy, dihedralEnergy;
3757  BigReal improperEnergy, crosstermEnergy;
3759  BigReal miscEnergy, bcEnergy;
3760  bool cudaIntegrator = simParams->CUDASOAintegrate;
3761 
3762 #ifdef NODEGROUP_FORCE_REGISTER
3763  if(cudaIntegrator){
3764  bondEnergy = nodeReduction->item(REDUCTION_BOND_ENERGY);
3765  angleEnergy = nodeReduction->item(REDUCTION_ANGLE_ENERGY);
3766  dihedralEnergy = nodeReduction->item(REDUCTION_DIHEDRAL_ENERGY);
3767  improperEnergy = nodeReduction->item(REDUCTION_IMPROPER_ENERGY);
3768  crosstermEnergy = nodeReduction->item(REDUCTION_CROSSTERM_ENERGY);
3769  bcEnergy = nodeReduction->item(REDUCTION_BC_ENERGY);
3770  miscEnergy = nodeReduction->item(REDUCTION_MISC_ENERGY);
3771  }else{
3772 #endif
3773  bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
3774  angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
3775  dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
3776  improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
3777  crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
3778  bcEnergy = reduction->item(REDUCTION_BC_ENERGY);
3779  miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
3780 #ifdef NODEGROUP_FORCE_REGISTER
3781  }
3782 #endif
3783 
3784  totalPotentialEnergy += bondEnergy + angleEnergy + dihedralEnergy +
3785  improperEnergy + crosstermEnergy + miscEnergy
3786  + bcEnergy;
3787 
3788  if (! ( step % nbondFreq ))
3789  {
3790 #ifdef NODEGROUP_FORCE_REGISTER
3791  if(cudaIntegrator){
3794 
3795  }else{
3796 #endif
3799 #ifdef NODEGROUP_FORCE_REGISTER
3800  }
3801 #endif
3802  totalPotentialEnergy += electEnergy + ljEnergy;
3803  }
3804 
3805  if (! ( step % slowFreq ))
3806  {
3807 #ifdef NODEGROUP_FORCE_REGISTER
3808  if(cudaIntegrator){
3810  }else{
3811 #endif
3813 #ifdef NODEGROUP_FORCE_REGISTER
3814  }
3815 #endif
3816  totalPotentialEnergy += electEnergySlow;
3817  }
3818 
3819  if (simParams->LJcorrection && volume) {
3820 #ifdef MEM_OPT_VERSION
3821  NAMD_bug("LJcorrection not supported in memory optimized build.");
3822 #else
3823  // Apply tail correction to energy.
3824  BigReal alchLambda = simParams->getCurrentLambda(step);
3825  totalPotentialEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
3826 #endif
3827  }
3828 
3829  #if 0
3830  printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, Bond: %f, Angle: %f, Dih: %f, IMP: %f, Cross: %f\n",
3831  step, CkMyPe(), bondEnergy, angleEnergy, dihedralEnergy, improperEnergy, crosstermEnergy);
3832  printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, LJ: %f, Real: %f, Recip: %f\n", step, CkMyPe(),
3834  printf("MC-data (TotalPotentialEnergy): step: %d, Pe: %d, MISC: %f, BC: %f, TOTAL: %f\n", step, CkMyPe(),
3835  miscEnergy, bcEnergy, totalPotentialEnergy);
3836  #endif
3837 
3838  return totalPotentialEnergy;
3839 }
3840 
3841 
3842 void Controller::printEnergies(int step, int minimize)
3843 {
3844  Node *node = Node::Object();
3845  Molecule *molecule = node->molecule;
3846  Lattice &lattice = state->lattice;
3847  bool cudaIntegrator = simParams->CUDASOAintegrate;
3848 
3849  // Drude model ANISO energy is added into BOND energy
3850  // and THOLE energy is added into ELECT energy
3851 
3852  BigReal bondEnergy;
3853  BigReal angleEnergy;
3854  BigReal dihedralEnergy;
3855  BigReal improperEnergy;
3856  BigReal crosstermEnergy;
3857  BigReal boundaryEnergy;
3858  BigReal miscEnergy;
3859  BigReal potentialEnergy;
3860  BigReal flatEnergy;
3861  BigReal smoothEnergy;
3862  BigReal work;
3863 
3864  Vector momentum;
3865  Vector angularMomentum;
3866  BigReal volume = lattice.volume();
3867 #ifdef NODEGROUP_FORCE_REGISTER
3868  if(cudaIntegrator){
3869  bondEnergy = nodeReduction->item(REDUCTION_BOND_ENERGY);
3870  angleEnergy = nodeReduction->item(REDUCTION_ANGLE_ENERGY);
3871  dihedralEnergy = nodeReduction->item(REDUCTION_DIHEDRAL_ENERGY);
3872  improperEnergy = nodeReduction->item(REDUCTION_IMPROPER_ENERGY);
3873  crosstermEnergy = nodeReduction->item(REDUCTION_CROSSTERM_ENERGY);
3874  boundaryEnergy = nodeReduction->item(REDUCTION_BC_ENERGY);
3875  miscEnergy = nodeReduction->item(REDUCTION_MISC_ENERGY);
3876  }else{
3877 #endif
3878  bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
3879  angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
3880  dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
3881  improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
3882  crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
3883  boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
3884  miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
3885 #ifdef NODEGROUP_FORCE_REGISTER
3886  }
3887 #endif
3888 
3889  if ( minimize || ! ( step % nbondFreq ) )
3890  {
3891 #ifdef NODEGROUP_FORCE_REGISTER
3892  if(cudaIntegrator){
3895 
3896  // JLai
3899 
3903 
3904  //fepb
3908 
3915  }else{
3916 #endif
3919 
3920  // JLai
3923 
3927 
3928  //fepb
3932 
3939 #ifdef NODEGROUP_FORCE_REGISTER
3940  }
3941 #endif
3942 //fepe
3943  }
3944 
3945  if ( minimize || ! ( step % slowFreq ) )
3946  {
3947 #ifdef NODEGROUP_FORCE_REGISTER
3948  if(cudaIntegrator){
3950  //fepb
3952 
3957  }else{
3958 #endif
3960  //fepb
3962 
3967 #ifdef NODEGROUP_FORCE_REGISTER
3968  }
3969 #endif
3970 //fepe
3971  }
3972 
3973  if ((simParams->LJcorrection || simParams->LJcorrectionAlt) && volume) {
3974 #ifdef MEM_OPT_VERSION
3975  NAMD_bug("LJcorrection not supported in memory optimized build.");
3976 #else
3977  // Apply tail correction to energy.
3978  BigReal alchLambda = simParams->getCurrentLambda(step);
3979  ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
3980 
3981  if (simParams->alchOn) {
3982  if (simParams->alchFepOn) {
3983  BigReal alchLambda2 = simParams->getCurrentLambda2(step);
3984  ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
3985  } else if (simParams->alchThermIntOn) {
3986  ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
3987  ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
3988  }
3989  }
3990 #endif
3991  }
3992 
3993 //fepb BKR - Compute alchemical work if using dynamic lambda. This is here so
3994 // that the cumulative work can be given during a callback.
3995  if (simParams->alchLambdaFreq > 0) {
3996  if (step <=
3998  cumAlchWork = 0.0;
3999  }
4000  alchWork = computeAlchWork(step);
4001  cumAlchWork += alchWork;
4002  }
4003 //fepe
4004 #ifdef NODEGROUP_FORCE_REGISTER
4005  if(cudaIntegrator){
4006  momentum.x = nodeReduction->item(REDUCTION_MOMENTUM_X);
4007  momentum.y = nodeReduction->item(REDUCTION_MOMENTUM_Y);
4008  momentum.z = nodeReduction->item(REDUCTION_MOMENTUM_Z);
4009  angularMomentum.x = nodeReduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
4010  angularMomentum.y = nodeReduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
4011  angularMomentum.z = nodeReduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
4012  }else{
4013 #endif
4014  momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
4015  momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
4016  momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
4017  angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
4018  angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
4019  angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
4020 #ifdef NODEGROUP_FORCE_REGISTER
4021  }
4022 #endif
4023 
4024  // Ported by JLai
4025  potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
4026  + improperEnergy + electEnergy + electEnergySlow + ljEnergy
4027  + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy
4029  // End of port
4030  totalEnergy = potentialEnergy + kineticEnergy;
4031  if ( ! cudaIntegrator) {
4032  flatEnergy = (totalEnergy
4034  if ( !(step%slowFreq) ) {
4035  // only adjust based on most accurate energies
4037  if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
4038  if ( step != simParams->firstTimestep ) {
4039  smooth2_avg *= 0.9375;
4040  smooth2_avg += 0.0625 * s;
4041  }
4042  }
4043  smoothEnergy = (flatEnergy + smooth2_avg
4045  }
4046 
4047  // Reset values for accumulated heat and work.
4048  if (step <= simParams->firstTimestep &&
4050  heat = 0.0;
4052  }
4053  if ( simParams->outputMomenta && ! minimize &&
4054  ! ( step % simParams->outputMomenta ) )
4055  {
4056  iout << "MOMENTUM: " << step
4057  << " P: " << momentum
4058  << " L: " << angularMomentum
4059  << "\n" << endi;
4060  }
4061 
4062  if ( simParams->outputPressure ) {
4063  if ( ! cudaIntegrator ) {
4066  }
4067  tavg_count += 1;
4068  if ( minimize || ! ( step % simParams->outputPressure ) ) {
4069  if (cudaIntegrator) {
4070  // tensors are valid for outputPressure step
4071  // update averages, then copy into tensor for output
4105  iout << "PRESSURE: " << step << " "
4106  << PRESSUREFACTOR * pressure << "\n"
4107  << "GPRESSURE: " << step << " "
4108  << PRESSUREFACTOR * groupPressure << "\n";
4109  // tavg_count makes output behaves the same as for standard case
4110  // but is no longer part of the average
4111  if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
4112  << (PRESSUREFACTOR) * pressure_tavg << "\n"
4113  << "GPRESSAVG: " << step << " "
4114  << (PRESSUREFACTOR) * groupPressure_tavg << "\n";
4115  iout << endi;
4116  tavg_count = 0;
4117  }
4118  else {
4119  iout << "PRESSURE: " << step << " "
4120  << PRESSUREFACTOR * pressure << "\n"
4121  << "GPRESSURE: " << step << " "
4122  << PRESSUREFACTOR * groupPressure << "\n";
4123  if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
4124  << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
4125  << "GPRESSAVG: " << step << " "
4127  iout << endi;
4128  pressure_tavg = 0;
4129  groupPressure_tavg = 0;
4130  tavg_count = 0;
4131  }
4132  }
4133  }
4134 
4135  // pressure profile reductions
4136  if (pressureProfileSlabs) {
4137  const int freq = simParams->pressureProfileFreq;
4138  const int arraysize = 3*pressureProfileSlabs;
4139 
4140  BigReal *total = new BigReal[arraysize];
4141  memset(total, 0, arraysize*sizeof(BigReal));
4142  const int first = simParams->firstTimestep;
4143 
4144  if (ppbonded) ppbonded->getData(first, step, lattice, total);
4145  if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
4146  if (ppint) ppint->getData(first, step, lattice, total);
4147  for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
4149 
4150  if (!(step % freq)) {
4151  // convert NAMD internal virial to pressure in units of bar
4153 
4154  iout << "PRESSUREPROFILE: " << step << " ";
4155  if (step == first) {
4156  // output pressure profile for this step
4157  for (int i=0; i<arraysize; i++) {
4158  iout << total[i] * scalefac << " ";
4159  }
4160  } else {
4161  // output pressure profile averaged over the last count steps.
4162  scalefac /= pressureProfileCount;
4163  for (int i=0; i<arraysize; i++)
4164  iout << pressureProfileAverage[i]*scalefac << " ";
4165  }
4166  iout << "\n" << endi;
4167 
4168  // Clear the average for the next block
4169  memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
4171  }
4172  delete [] total;
4173  }
4174 
4175  if ( step != simParams->firstTimestep || stepInFullRun == 0 ) { // skip repeated first step
4176  if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
4177  int benchPhase = stepInFullRun / simParams->firstLdbStep;
4178  if ( benchPhase > 0 && benchPhase < 10 ) {
4179 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
4180  if ( simParams->computeEnergies < 60 ) {
4181  iout << iWARN << "Energy evaluation is expensive, increase computeEnergies to improve performance.\n";
4182  iout << iWARN << "If computeEnergies is not defined, its value defaults to the same as outputEnergies, and increasing outputEnergies would be helpful to improve the performance.\n";
4183  }
4184 #endif
4185  iout << iINFO;
4186  if ( benchPhase < 4 ) iout << "Initial time: ";
4187  else iout << "Benchmark time: ";
4188  iout << CkNumPes() << " CPUs ";
4189 
4190  {
4191  BigReal wallPerStep =
4192  (namdWallTimer() - startBenchTime) / simParams->firstLdbStep;
4193  BigReal ns = simParams->dt / 1000000.0;
4194  BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
4195 
4196  if (simParams->nsPerDayOn) {
4197  BigReal nsPerDay = ns / (wallPerStep * days);
4198  iout << wallPerStep << " s/step " << nsPerDay << " ns/day ";
4199  }
4200  else {
4201  BigReal daysPerNano = wallPerStep * days / ns;
4202  iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
4203  }
4204  iout << memusage_MB() << " MB memory\n" << endi;
4205  }
4206 
4207  }
4208  startBenchTime = namdWallTimer();
4209  }
4210  ++stepInFullRun;
4211  }
4212 
4213  printTiming(step);
4214 
4215  Vector pairVDWForce, pairElectForce;
4216  if ( simParams->pairInteractionOn ){
4217 #ifdef NODEGROUP_FORCE_REGISTER
4218  if(cudaIntegrator){
4219  GET_VECTOR(pairVDWForce,nodeReduction,REDUCTION_PAIR_VDW_FORCE);
4220  GET_VECTOR(pairElectForce,nodeReduction,REDUCTION_PAIR_ELECT_FORCE);
4221  }else{
4222 #endif
4223  GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
4224  GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
4225 #ifdef NODEGROUP_FORCE_REGISTER
4226  }
4227 #endif
4228  }
4229 
4230  // Compute cumulative nonequilibrium work (including shadow work).
4232  work = totalEnergy - totalEnergy0 - heat;
4233  }
4234 
4235  // callback to Tcl with whatever we can
4236 #ifdef NAMD_TCL
4237 #define CALLBACKDATA(LABEL,VALUE) \
4238  labels << (LABEL) << " "; values << (VALUE) << " ";
4239 #define CALLBACKLIST(LABEL,VALUE) \
4240  labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
4241  if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
4242 
4243  std::ostringstream labels, values;
4244 #if CMK_BLUEGENEL
4245  // the normal version below gives a compiler error
4246  values.precision(16);
4247 #else
4248  values << std::setprecision(16);
4249 #endif
4250  CALLBACKDATA("TS",step);
4251  CALLBACKDATA("BOND",bondEnergy);
4252  CALLBACKDATA("ANGLE",angleEnergy);
4253  CALLBACKDATA("DIHED",dihedralEnergy);
4254  CALLBACKDATA("CROSS",crosstermEnergy);
4255  CALLBACKDATA("IMPRP",improperEnergy);
4257  CALLBACKDATA("VDW",ljEnergy);
4258  CALLBACKDATA("BOUNDARY",boundaryEnergy);
4259  CALLBACKDATA("MISC",miscEnergy);
4260  CALLBACKDATA("POTENTIAL",potentialEnergy);
4261  CALLBACKDATA("KINETIC",kineticEnergy);
4262  CALLBACKDATA("TOTAL",totalEnergy);
4263  CALLBACKDATA("TEMP",temperature);
4264  CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
4266  CALLBACKDATA("VOLUME",lattice.volume());
4267  CALLBACKLIST("CELL_A",lattice.a());
4268  CALLBACKLIST("CELL_B",lattice.b());
4269  CALLBACKLIST("CELL_C",lattice.c());
4270  CALLBACKLIST("CELL_O",lattice.origin());
4271  labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
4272  << lattice.b_p() << " " << lattice.c_p() << "} ";
4273  if (simParams->drudeOn) {
4274  CALLBACKDATA("DRUDEBOND",drudeBondTemp);
4275  }
4276  if ( simParams->pairInteractionOn ) {
4277  CALLBACKLIST("VDW_FORCE",pairVDWForce);
4278  CALLBACKLIST("ELECT_FORCE",pairElectForce);
4279  }
4281  CALLBACKLIST("HEAT",heat);
4282  CALLBACKLIST("WORK",work);
4283  }
4284  if (simParams->alchOn) {
4285  if (simParams->alchThermIntOn) {
4286  CALLBACKLIST("BOND1", bondedEnergy_ti_1);
4289  CALLBACKLIST("VDW1", ljEnergy_ti_1);
4290  CALLBACKLIST("BOND2", bondedEnergy_ti_2);
4293  CALLBACKLIST("VDW2", ljEnergy_ti_2);
4294  if (simParams->alchLambdaFreq > 0) {
4295  CALLBACKLIST("CUMALCHWORK", cumAlchWork);
4296  }
4297  } else if (simParams->alchFepOn) {
4298  CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
4299  improperEnergy + bondedEnergyDiff_f);
4301  CALLBACKLIST("VDW2", ljEnergy_f);
4302  }
4303  }
4304 
4305  labels << '\0'; values << '\0'; // insane but makes Linux work
4306  state->callback_labelstring = labels.str();
4307  state->callback_valuestring = values.str();
4308  // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
4309  }
4310 #undef CALLBACKDATA
4311 #endif
4312 
4313  if ( ! cudaIntegrator) {
4315  temp_avg += temperature;
4316  pressure_avg += trace(pressure)/3.;
4317  groupPressure_avg += trace(groupPressure)/3.;
4318  avg_count += 1;
4319  }
4320 
4322  ! (step % simParams->outputPairlists) ) {
4323  iout << iINFO << pairlistWarnings <<
4324  " pairlist warnings in past " << simParams->outputPairlists <<
4325  " steps.\n" << endi;
4326  pairlistWarnings = 0;
4327  }
4328 
4329  BigReal enthalpy;
4330  if (simParams->multigratorOn && ((step % simParams->computeEnergies) == 0)) {
4331  enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
4332  }
4333 
4334  // XXX
4335  // Important note: in CPU-only/GPU-offload modes
4336  // printEnergies is called EVERY STEP.
4337  // Reduction averages are summed above and then later are output
4338  // divided by count. Sums and counter are reset after printing.
4339  // XXX
4340  // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
4341  if ( ! minimize && step % simParams->outputEnergies ) return;
4342 
4343  if (cudaIntegrator) {
4346  pressureAverage.addSample((1./3)*trace(pressure));
4348  }
4349 
4350  // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
4351 
4352  if (simParams->IMDon && !(step % simParams->IMDfreq)) {
4353  IMDEnergies energies;
4354  energies.tstep = step;
4355  energies.T = temp_avg/avg_count;
4356  energies.Etot = totalEnergy;
4357  energies.Epot = potentialEnergy;
4358  energies.Evdw = ljEnergy;
4359  energies.Eelec = electEnergy + electEnergySlow;
4360  energies.Ebond = bondEnergy;
4361  energies.Eangle = angleEnergy;
4362  energies.Edihe = dihedralEnergy + crosstermEnergy;
4363  energies.Eimpr = improperEnergy;
4364  Node::Object()->imd->gather_energies(&energies);
4365  }
4366 
4367  if ( marginViolations ) {
4368  iout << iERROR << marginViolations <<
4369  " margin violations detected since previous energy output.\n" << endi;
4370  }
4371  marginViolations = 0;
4372 
4373  int precision = simParams->outputEnergiesPrecision;
4374  if ( (step % (10 * (minimize?1:simParams->outputEnergies) ) ) == 0 )
4375  {
4376  iout << "ETITLE: TS";
4377  iout << FORMAT("BOND", precision);
4378  iout << FORMAT("ANGLE", precision);
4379  iout << FORMAT("DIHED", precision);
4380  if ( ! simParams->mergeCrossterms ) iout << FORMAT("CROSS", precision);
4381  iout << FORMAT("IMPRP", precision);
4382  iout << " ";
4383  iout << FORMAT("ELECT", precision);
4384  iout << FORMAT("VDW", precision);
4385  iout << FORMAT("BOUNDARY", precision);
4386  iout << FORMAT("MISC", precision);
4387  iout << FORMAT("KINETIC", precision);
4388  iout << " ";
4389  iout << FORMAT("TOTAL", precision);
4390  iout << FORMAT("TEMP", precision);
4391  iout << FORMAT("POTENTIAL", precision);
4392  if (cudaIntegrator) {
4393  iout << FORMAT("TOTALAVG", precision);
4394  }
4395  else {
4396  // iout << FORMAT("TOTAL2", precision);
4397  iout << FORMAT("TOTAL3", precision);
4398  }
4399  iout << FORMAT("TEMPAVG", precision);
4400  if ( volume != 0. ) {
4401  iout << " ";
4402  iout << FORMAT("PRESSURE", precision);
4403  iout << FORMAT("GPRESSURE", precision);
4404  iout << FORMAT("VOLUME", precision);
4405  iout << FORMAT("PRESSAVG", precision);
4406  iout << FORMAT("GPRESSAVG", precision);
4407  }
4409  iout << " ";
4410  iout << FORMAT("HEAT", precision);
4411  iout << FORMAT("WORK", precision);
4412  }
4413  if (simParams->drudeOn) {
4414  iout << " ";
4415  iout << FORMAT("DRUDEBOND", precision);
4416  iout << FORMAT("DRBONDAVG", precision);
4417  }
4418  // Ported by JLai
4419  if (simParams->goGroPair) {
4420  iout << " ";
4421  iout << FORMAT("GRO_PAIR_LJ", precision);
4422  iout << FORMAT("GRO_PAIR_GAUSS", precision);
4423  }
4424 
4425  if (simParams->goForcesOn) {
4426  iout << " ";
4427  iout << FORMAT("NATIVE", precision);
4428  iout << FORMAT("NONNATIVE", precision);
4429  //iout << FORMAT("REL_NATIVE", precision);
4430  //iout << FORMAT("REL_NONNATIVE", precision);
4431  iout << FORMAT("GOTOTAL", precision);
4432  //iout << FORMAT("GOAVG", precision);
4433  }
4434  // End of port -- JLai
4435 
4436  if (simParams->alchOn) {
4437  if (simParams->alchThermIntOn) {
4438  iout << "\nTITITLE: TS";
4439  iout << FORMAT("BOND1", precision);
4440  iout << FORMAT("ELECT1", precision);
4441  iout << FORMAT("VDW1", precision);
4442  iout << FORMAT("BOND2", precision);
4443  iout << " ";
4444  iout << FORMAT("ELECT2", precision);
4445  iout << FORMAT("VDW2", precision);
4446  if (simParams->alchLambdaFreq > 0) {
4447  iout << FORMAT("LAMBDA", precision);
4448  iout << FORMAT("ALCHWORK", precision);
4449  iout << FORMAT("CUMALCHWORK", precision);
4450  }
4451  } else if (simParams->alchFepOn) {
4452  iout << "\nFEPTITLE: TS";
4453  iout << FORMAT("BOND2", precision);
4454  iout << FORMAT("ELECT2", precision);
4455  iout << FORMAT("VDW2", precision);
4456  if (simParams->alchLambdaFreq > 0) {
4457  iout << FORMAT("LAMBDA", precision);
4458  }
4459  }
4460  }
4461 
4462  iout << "\n\n" << endi;
4463 
4464  if (simParams->qmForcesOn) {
4465  iout << "QMETITLE: TS";
4466  iout << FORMAT("QMID", precision);
4467  iout << FORMAT("ENERGY", precision);
4468  if (simParams->PMEOn) iout << FORMAT("PMECORRENERGY", precision);
4469  iout << "\n\n" << endi;
4470  }
4471 
4472  }
4473 
4474  // N.B. HP's aCC compiler merges FORMAT calls in the same expression.
4475  // Need separate statements because data returned in static array.
4476  iout << ETITLE(step);
4477  iout << FORMAT(bondEnergy, precision);
4478  iout << FORMAT(angleEnergy, precision);
4479  if ( simParams->mergeCrossterms ) {
4480  iout << FORMAT(dihedralEnergy+crosstermEnergy, precision);
4481  } else {
4482  iout << FORMAT(dihedralEnergy, precision);
4483  iout << FORMAT(crosstermEnergy, precision);
4484  }
4485  iout << FORMAT(improperEnergy, precision);
4486  iout << " ";
4487  iout << FORMAT(electEnergy+electEnergySlow, precision);
4488  iout << FORMAT(ljEnergy, precision);
4489  iout << FORMAT(boundaryEnergy, precision);
4490  iout << FORMAT(miscEnergy, precision);
4491  iout << FORMAT(kineticEnergy, precision);
4492  iout << " ";
4493  iout << FORMAT(totalEnergy, precision);
4494  iout << FORMAT(temperature, precision);
4495  iout << FORMAT(potentialEnergy, precision);
4496  if (cudaIntegrator) {
4497  iout << FORMAT(totalEnergyAverage.average(), precision);
4498  iout << FORMAT(temperatureAverage.average(), precision);
4499  }
4500  else {
4501  // iout << FORMAT(flatEnergy, precision);
4502  iout << FORMAT(smoothEnergy, precision);
4503  iout << FORMAT(temp_avg/avg_count, precision);
4504  }
4505  if ( volume != 0. )
4506  {
4507  iout << " ";
4508  iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3., precision);
4509  iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3., precision);
4510  iout << FORMAT(volume, precision);
4511  if (cudaIntegrator) {
4514  }
4515  else {
4518  }
4519  }
4521  iout << " ";
4522  iout << FORMAT(heat, precision);
4523  iout << FORMAT(work, precision);
4524  }
4525  if (simParams->drudeOn) {
4526  iout << " ";
4527  iout << FORMAT(drudeBondTemp, precision);
4528  iout << FORMAT(drudeBondTempAvg/avg_count, precision);
4529  }
4530  // Ported by JLai
4531  if (simParams->goGroPair) {
4532  iout << " ";
4533  iout << FORMAT(groLJEnergy, precision);
4534  iout << FORMAT(groGaussEnergy, precision);
4535  }
4536 
4537  if (simParams->goForcesOn) {
4538  iout << " ";
4539  iout << FORMAT(goNativeEnergy, precision);
4540  iout << FORMAT(goNonnativeEnergy, precision);
4541  //iout << FORMAT(relgoNativeEnergy, precision);
4542  //iout << FORMAT(relgoNonnativeEnergy, precision);
4543  iout << FORMAT(goTotalEnergy, precision);
4544  //iout << FORMAT("not implemented", precision);
4545  } // End of port -- JLai
4546 
4547  if (simParams->alchOn) {
4548  if (simParams->alchThermIntOn) {
4549  iout << "\n";
4550  iout << TITITLE(step);
4551  iout << FORMAT(bondedEnergy_ti_1, precision);
4553  electEnergyPME_ti_1, precision);
4554  iout << FORMAT(ljEnergy_ti_1, precision);
4555  iout << FORMAT(bondedEnergy_ti_2, precision);
4556  iout << " ";
4558  electEnergyPME_ti_2, precision);
4559  iout << FORMAT(ljEnergy_ti_2, precision);
4560  if (simParams->alchLambdaFreq > 0) {
4561  iout << FORMAT(simParams->getCurrentLambda(step), precision);
4562  iout << FORMAT(alchWork, precision);
4563  iout << FORMAT(cumAlchWork, precision);
4564  }
4565  } else if (simParams->alchFepOn) {
4566  iout << "\n";
4567  iout << FEPTITLE2(step);
4568  iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy
4569  + improperEnergy + bondedEnergyDiff_f, precision);
4570  iout << FORMAT(electEnergy_f + electEnergySlow_f, precision);
4571  iout << FORMAT(ljEnergy_f, precision);
4572  if (simParams->alchLambdaFreq > 0) {
4573  iout << FORMAT(simParams->getCurrentLambda(step), precision);
4574  }
4575  }
4576  }