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 "memusage.h"
16 #include "Node.h"
17 #include "Molecule.h"
18 #include "SimParameters.h"
19 #include "Controller.h"
20 #include "ReductionMgr.h"
21 #include "CollectionMaster.h"
22 #include "Output.h"
23 #include "strlib.h"
24 #include "BroadcastObject.h"
25 #include "NamdState.h"
26 #include "ScriptTcl.h"
27 #include "Broadcasts.h"
28 #include "LdbCoordinator.h"
29 #include "Thread.h"
30 #include <math.h>
31 #include <signal.h>
32 #include "NamdOneTools.h"
33 #include "PatchMap.h"
34 #include "PatchMap.inl"
35 #include "Random.h"
36 #include "imd.h"
37 #include "IMDOutput.h"
38 #include "BackEnd.h"
39 #include <fstream>
40 #include <iomanip>
41 #include <errno.h>
42 #include "qd.h"
43 
45 
46 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
47 extern "C" void CApplicationDepositNode0Data(char *);
48 #endif
49 
50 #ifndef cbrt
51  // cbrt() not in math.h on goneril
52  #define cbrt(x) pow(x,(double)(1.0/3.0))
53 #endif
54 
55 //#define DEBUG_PRESSURE
56 #define MIN_DEBUG_LEVEL 3
57 //#define DEBUGM
58 #include "Debug.h"
59 
60 #define XXXBIGREAL 1.0e32
61 
63 private:
64  RequireReduction *reduction;
65  const int nslabs;
66  const int freq;
67  int nelements;
68  int count;
69  char *name;
70 
71  BigReal *average;
72 
73 public:
74  PressureProfileReduction(int rtag, int numslabs, int numpartitions,
75  const char *myname, int outputfreq)
76  : nslabs(numslabs), freq(outputfreq) {
77  name = strdup(myname);
78  nelements = 3*nslabs * numpartitions;
79  reduction = ReductionMgr::Object()->willRequire(rtag,nelements);
80 
81  average = new BigReal[nelements];
82  count = 0;
83  }
85  delete [] average;
86  delete reduction;
87  free(name);
88  }
89  //
90  void getData(int firsttimestep, int step, const Lattice &lattice,
91  BigReal *total) {
92  reduction->require();
93 
94  int i;
95  double inv_volume = 1.0 / lattice.volume();
96  // accumulate the pressure profile computed for this step into the average.
97  int arraysize = 3*nslabs;
98  for (i=0; i<nelements; i++) {
99  BigReal val = reduction->item(i) * inv_volume;
100  average[i] += val;
101  total[i % arraysize] += val;
102  }
103  count++;
104  if (!(step % freq)) {
105  // convert NAMD internal virial to pressure in units of bar
106  BigReal scalefac = PRESSUREFACTOR * nslabs;
107 
108  iout << "PPROFILE" << name << ": " << step << " ";
109  if (step == firsttimestep) {
110  // output pressure profile for this step
111  for (i=0; i<nelements; i++)
112  iout << reduction->item(i)*scalefac*inv_volume << " ";
113  } else {
114  // output pressure profile averaged over the last Freq steps.
115  scalefac /= count;
116  for (i=0; i<nelements; i++)
117  iout << average[i]*scalefac << " ";
118  }
119  iout << "\n" << endi;
120  // Clear the average for the next block
121  memset(average, 0, nelements*sizeof(BigReal));
122  count = 0;
123  }
124  }
125 };
126 
127 
129  computeChecksum(0), marginViolations(0), pairlistWarnings(0),
130  simParams(Node::Object()->simParameters),
131  state(s),
132  collection(CollectionMaster::Object()),
133  startCTime(0),
134  firstCTime(CmiTimer()),
135  startWTime(0),
136  firstWTime(CmiWallTimer()),
137  startBenchTime(0),
138  stepInFullRun(0),
139  ldbSteps(0),
140  fflush_count(3)
141 {
144  // for accelMD
145  if (simParams->accelMDOn) {
146  // REDUCTIONS_BASIC wil contain all potential energies and arrive first
148  // contents of amd_reduction be submitted to REDUCTIONS_AMD
150  // REDUCTIONS_AMD will contain Sequencer contributions and arrive second
152  } else {
153  amd_reduction = NULL;
154  submit_reduction = NULL;
156  }
157  // pressure profile reductions
160  ppbonded = ppnonbonded = ppint = NULL;
162  int ntypes = simParams->pressureProfileAtomTypes;
164  // number of partitions for pairwise interactions
165  int npairs = (ntypes * (ntypes+1))/2;
166  pressureProfileAverage = new BigReal[3*nslabs];
167  memset(pressureProfileAverage, 0, 3*nslabs*sizeof(BigReal));
168  int freq = simParams->pressureProfileFreq;
171  nslabs, npairs, "BONDED", freq);
173  nslabs, npairs, "NONBONDED", freq);
174  // internal partitions by atom type, but not in a pairwise fashion
176  nslabs, ntypes, "INTERNAL", freq);
177  } else {
178  // just doing Ewald, so only need nonbonded
180  nslabs, npairs, "NONBONDED", freq);
181  }
182  }
185 
186  heat = totalEnergy0 = 0.0;
189  stochRescale_count = 0;
190  if (simParams->stochRescaleOn) {
193  }
195  // strainRate tensor is symmetric to avoid rotation
198  if ( ! simParams->useFlexibleCell ) {
199  BigReal avg = trace(langevinPiston_strainRate) / 3.;
201  } else if ( simParams->useConstantRatio ) {
202 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
203  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
205 #undef AVGXY
206  }
208  if (simParams->multigratorOn) {
209  multigratorXi = 0.0;
212  Node *node = Node::Object();
213  Molecule *molecule = node->molecule;
214  BigReal Nf = molecule->num_deg_freedom();
216  multigratorNu.resize(n);
217  multigratorNuT.resize(n);
218  multigratorZeta.resize(n);
219  multigratorOmega.resize(n);
220  for (int i=0;i < n;i++) {
221  multigratorNu[i] = 0.0;
222  multigratorZeta[i] = 0.0;
223  if (i == 0) {
224  multigratorOmega[i] = Nf*kT0*tau*tau;
225  } else {
226  multigratorOmega[i] = kT0*tau*tau;
227  }
228  }
230  } else {
231  multigratorReduction = NULL;
232  }
233  origLattice = state->lattice;
235  temp_avg = 0;
236  pressure_avg = 0;
237  groupPressure_avg = 0;
238  avg_count = 0;
239  pressure_tavg = 0;
240  groupPressure_tavg = 0;
241  tavg_count = 0;
242  checkpoint_stored = 0;
243  drudeBondTemp = 0;
244  drudeBondTempAvg = 0;
245  cumAlchWork = 0;
246 }
247 
249 {
250  delete broadcast;
251  delete reduction;
252  delete min_reduction;
253  delete amd_reduction;
254  delete submit_reduction;
255  delete ppbonded;
256  delete ppnonbonded;
257  delete ppint;
258  delete [] pressureProfileAverage;
259  delete random;
261 }
262 
263 void Controller::threadRun(Controller* arg)
264 {
265  arg->algorithm();
266 }
267 
268 void Controller::run(void)
269 {
270  // create a Thread and invoke it
271  DebugM(4, "Starting thread in controller on this=" << this << "\n");
272  thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),CTRL_STK_SZ);
273  CthSetStrategyDefault(thread);
274 #if CMK_BLUEGENE_CHARM
275  BgAttach(thread);
276 #endif
277  awaken();
278 }
279 
280 
282 {
283  int scriptTask;
284  int scriptSeq = 0;
285  BackEnd::awaken();
286  while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END) {
287  switch ( scriptTask ) {
288  case SCRIPT_OUTPUT:
291  break;
292  case SCRIPT_FORCEOUTPUT:
294  break;
295  case SCRIPT_MEASURE:
297  break;
298  case SCRIPT_REINITVELS:
299  iout << "REINITIALIZING VELOCITIES AT STEP " << simParams->firstTimestep
300  << " TO " << simParams->initialTemp << " KELVIN.\n" << endi;
301  break;
302  case SCRIPT_RESCALEVELS:
303  iout << "RESCALING VELOCITIES AT STEP " << simParams->firstTimestep
304  << " BY " << simParams->scriptArg1 << "\n" << endi;
305  break;
307  // Parameter setting already reported in ScriptTcl
308  // Nothing to do!
309  break;
310  case SCRIPT_CHECKPOINT:
311  iout << "CHECKPOINTING AT STEP " << simParams->firstTimestep
312  << "\n" << endi;
313  checkpoint_stored = 1;
314  checkpoint_lattice = state->lattice;
316  break;
317  case SCRIPT_REVERT:
318  iout << "REVERTING AT STEP " << simParams->firstTimestep
319  << "\n" << endi;
320  if ( ! checkpoint_stored )
321  NAMD_die("Unable to revert, checkpoint was never called!");
322  state->lattice = checkpoint_lattice;
324  break;
326  iout << "STORING CHECKPOINT AT STEP " << simParams->firstTimestep
327  << " TO KEY " << simParams->scriptStringArg1;
328  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
329  iout << "\n" << endi;
330  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
331  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
333  }
335  cp.lattice = state->lattice;
336  cp.state = *(ControllerState*)this;
337  } else {
339  scriptTask, state->lattice, *(ControllerState*)this);
340  }
341  break;
343  iout << "LOADING CHECKPOINT AT STEP " << simParams->firstTimestep
344  << " FROM KEY " << simParams->scriptStringArg1;
345  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
346  iout << "\n" << endi;
347  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
348  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
349  NAMD_die("Unable to load checkpoint, requested key was never stored.");
350  }
352  state->lattice = cp.lattice;
353  *(ControllerState*)this = cp.state;
354  } else {
356  scriptTask, state->lattice, *(ControllerState*)this);
357  }
358  break;
360  iout << "SWAPPING CHECKPOINT AT STEP " << simParams->firstTimestep
361  << " FROM KEY " << simParams->scriptStringArg1;
362  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
363  iout << "\n" << endi;
364  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
365  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
366  NAMD_die("Unable to swap checkpoint, requested key was never stored.");
367  }
369  std::swap(state->lattice,cp.lattice);
370  std::swap(*(ControllerState*)this,cp.state);
371  } else {
373  scriptTask, state->lattice, *(ControllerState*)this);
374  }
375  break;
377  iout << "FREEING CHECKPOINT AT STEP " << simParams->firstTimestep
378  << " FROM KEY " << simParams->scriptStringArg1;
379  if ( CmiNumPartitions() > 1 ) iout << " ON REPLICA " << simParams->scriptIntArg1;
380  iout << "\n" << endi;
381  if ( simParams->scriptIntArg1 == CmiMyPartition() ) {
382  if ( ! checkpoints.count(simParams->scriptStringArg1) ) {
383  NAMD_die("Unable to free checkpoint, requested key was never stored.");
384  }
387  } else {
389  scriptTask, state->lattice, *(ControllerState*)this);
390  }
391  break;
392  case SCRIPT_ATOMSENDRECV:
393  case SCRIPT_ATOMSEND:
394  case SCRIPT_ATOMRECV:
395  break;
396  case SCRIPT_MINIMIZE:
397  minimize();
398  break;
399  case SCRIPT_RUN:
400  case SCRIPT_CONTINUE:
401  integrate(scriptTask);
402  break;
403  default:
404  NAMD_bug("Unknown task in Controller::algorithm");
405  }
406  BackEnd::awaken();
407  }
410  terminate();
411 }
412 
413 
414 //
415 // XXX static and global variables are unsafe for shared memory builds.
416 // The use of global and static vars should be eliminated.
417 //
418 extern int eventEndOfTimeStep;
419 
420 // Handle SIGINT so that restart files get written completely.
421 static int gotsigint = 0;
422 static void my_sigint_handler(int sig) {
423  if (sig == SIGINT) gotsigint = 1;
424 }
425 extern "C" {
426  typedef void (*namd_sighandler_t)(int);
427 }
428 
429 void Controller::integrate(int scriptTask) {
430  char traceNote[24];
431 
432  int step = simParams->firstTimestep;
433 
434  const int numberOfSteps = simParams->N;
435  const int stepsPerCycle = simParams->stepsPerCycle;
436 
437  const int zeroMomentum = simParams->zeroMomentum;
438 
440  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
441  if (dofull)
443  else
445  if ( step >= numberOfSteps ) slowFreq = nbondFreq = 1;
446 
447  if ( scriptTask == SCRIPT_RUN ) {
448 
449  reassignVelocities(step); // only for full-step velecities
450  rescaleaccelMD(step);
451  adaptTempUpdate(step); // Init adaptive tempering;
452 
453  receivePressure(step);
454  if ( zeroMomentum && dofull && ! (step % slowFreq) )
455  correctMomentum(step);
456  printFepMessage(step);
457  printTiMessage(step);
458  printDynamicsEnergies(step);
459  outputFepEnergy(step);
460  outputTiEnergy(step);
461  if(traceIsOn()){
462  traceUserEvent(eventEndOfTimeStep);
463  sprintf(traceNote, "s:%d", step);
464  traceUserSuppliedNote(traceNote);
465  }
466  outputExtendedSystem(step);
467  rebalanceLoad(step);
468 
469  }
470 
471  // Handling SIGINT doesn't seem to be working on Lemieux, and it
472  // sometimes causes the net-xxx versions of NAMD to segfault on exit,
473  // so disable it for now.
474  // namd_sighandler_t oldhandler = signal(SIGINT,
475  // (namd_sighandler_t)my_sigint_handler);
476  for ( ++step ; step <= numberOfSteps; ++step )
477  {
478 
479  adaptTempUpdate(step);
480  rescaleVelocities(step);
481  tcoupleVelocities(step);
483  berendsenPressure(step);
484  langevinPiston1(step);
485  rescaleaccelMD(step);
486  enqueueCollections(step); // after lattice scaling!
487  receivePressure(step);
488  if ( zeroMomentum && dofull && ! (step % slowFreq) )
489  correctMomentum(step);
490  langevinPiston2(step);
491  reassignVelocities(step);
492 
493  multigratorTemperature(step, 1);
494  multigratorPressure(step, 1);
495  multigratorPressure(step, 2);
496  multigratorTemperature(step, 2);
497 
498  printDynamicsEnergies(step);
499  outputFepEnergy(step);
500  outputTiEnergy(step);
501  if(traceIsOn()){
502  traceUserEvent(eventEndOfTimeStep);
503  sprintf(traceNote, "s:%d", step);
504  traceUserSuppliedNote(traceNote);
505  }
506  // if (gotsigint) {
507  // iout << iINFO << "Received SIGINT; shutting down.\n" << endi;
508  // NAMD_quit();
509  // }
510  outputExtendedSystem(step);
511 #if CYCLE_BARRIER
512  cycleBarrier(!((step+1) % stepsPerCycle),step);
513 #elif PME_BARRIER
514  cycleBarrier(dofull && !(step%slowFreq),step);
515 #elif STEP_BARRIER
516  cycleBarrier(1, step);
517 #endif
518 
519  if(Node::Object()->specialTracing || simParams->statsOn){
520  int bstep = simParams->traceStartStep;
521  int estep = bstep + simParams->numTraceSteps;
522  if(step == bstep){
523  traceBarrier(1, step);
524  }
525  if(step == estep){
526  traceBarrier(0, step);
527  }
528  }
529 
530 #ifdef MEASURE_NAMD_WITH_PAPI
531  if(simParams->papiMeasure) {
532  int bstep = simParams->papiMeasureStartStep;
533  int estep = bstep + simParams->numPapiMeasureSteps;
534  if(step == bstep) {
535  papiMeasureBarrier(1, step);
536  }
537  if(step == estep) {
538  papiMeasureBarrier(0, step);
539  }
540  }
541 #endif
542 
543  rebalanceLoad(step);
544 
545 #if PME_BARRIER
546  cycleBarrier(dofull && !((step+1)%slowFreq),step); // step before PME
547 #endif
548  }
549  // signal(SIGINT, oldhandler);
550 }
551 
552 
553 #define CALCULATE \
554  printMinimizeEnergies(step); \
555  outputExtendedSystem(step); \
556  rebalanceLoad(step); \
557  if ( step == numberOfSteps ) return; \
558  else ++step;
559 
560 #define MOVETO(X) \
561  if ( step == numberOfSteps ) { \
562  if ( minVerbose ) { iout << "LINE MINIMIZER: RETURNING TO " << mid.x << " FROM " << last.x << "\n" << endi; } \
563  if ( newDir || (mid.x-last.x) ) { \
564  broadcast->minimizeCoefficient.publish(minSeq++,mid.x-last.x); \
565  } else { \
566  broadcast->minimizeCoefficient.publish(minSeq++,0.); \
567  broadcast->minimizeCoefficient.publish(minSeq++,0.); \
568  min_reduction->require(); \
569  broadcast->minimizeCoefficient.publish(minSeq++,0.); \
570  } \
571  enqueueCollections(step); \
572  CALCULATE \
573  } else if ( (X)-last.x ) { \
574  broadcast->minimizeCoefficient.publish(minSeq++,(X)-last.x); \
575  newDir = 0; \
576  last.x = (X); \
577  enqueueCollections(step); \
578  CALCULATE \
579  last.u = min_energy; \
580  last.dudx = -1. * min_f_dot_v; \
581  last.noGradient = min_huge_count; \
582  if ( minVerbose ) { \
583  iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx; \
584  if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient; \
585  iout << "\n" << endi; \
586  } \
587  }
588 
589 struct minpoint {
591  minpoint() : x(0), u(0), dudx(0), noGradient(1) { ; }
592 };
593 
595  // iout << "Controller::minimize() called.\n" << endi;
596 
597  const int minVerbose = simParams->minVerbose;
598  const int numberOfSteps = simParams->N;
599  int step = simParams->firstTimestep;
600  slowFreq = nbondFreq = 1;
601  BigReal linegoal = simParams->minLineGoal; // 1.0e-3
602  const BigReal goldenRatio = 0.5 * ( sqrt(5.0) - 1.0 );
603 
604  CALCULATE
605 
606  int minSeq = 0;
607 
608  // just move downhill until initial bad contacts go away or 100 steps
609  int old_min_huge_count = 2000000000;
610  int steps_at_same_min_huge_count = 0;
611  for ( int i_slow = 0; min_huge_count && i_slow < 100; ++i_slow ) {
612  if ( min_huge_count >= old_min_huge_count ) {
613  if ( ++steps_at_same_min_huge_count > 10 ) break;
614  } else {
615  old_min_huge_count = min_huge_count;
616  steps_at_same_min_huge_count = 0;
617  }
618  iout << "MINIMIZER SLOWLY MOVING " << min_huge_count << " ATOMS WITH BAD CONTACTS DOWNHILL\n" << endi;
619  broadcast->minimizeCoefficient.publish(minSeq++,1.);
620  enqueueCollections(step);
621  CALCULATE
622  }
623  if ( min_huge_count ) {
624  iout << "MINIMIZER GIVING UP ON " << min_huge_count << " ATOMS WITH BAD CONTACTS\n" << endi;
625  }
626  iout << "MINIMIZER STARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
627 
628  int atStart = 2;
629  int errorFactor = 10;
630  BigReal old_f_dot_f = min_f_dot_f;
631  broadcast->minimizeCoefficient.publish(minSeq++,0.);
632  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
633  int newDir = 1;
636  while ( 1 ) {
637  // line minimization
638  // bracket minimum on line
639  minpoint lo,hi,mid,last;
640  BigReal x = 0;
641  lo.x = x;
642  lo.u = min_energy;
643  lo.dudx = -1. * min_f_dot_v;
645  mid = lo;
646  last = mid;
647  if ( minVerbose ) {
648  iout << "LINE MINIMIZER: POSITION " << last.x << " ENERGY " << last.u << " GRADIENT " << last.dudx;
649  if ( last.noGradient ) iout << " HUGECOUNT " << last.noGradient;
650  iout << "\n" << endi;
651  }
652  BigReal tol = fabs( linegoal * min_f_dot_v );
653  iout << "LINE MINIMIZER REDUCING GRADIENT FROM " <<
654  fabs(min_f_dot_v) << " TO " << tol << "\n" << endi;
655  int start_with_huge = last.noGradient;
657  BigReal maxstep = 0.1 / sqrt(min_reduction->item(0));
658  x = maxstep; MOVETO(x);
659  // bracket minimum on line
660  while ( last.u < mid.u ) {
661  if ( last.dudx < mid.dudx * (5.5 - x/maxstep)/5. ) {
662  x = 6 * maxstep; break;
663  }
664  lo = mid; mid = last;
665  x += maxstep;
666  if ( x > 5.5 * maxstep ) break;
667  MOVETO(x)
668  }
669  if ( x > 5.5 * maxstep ) {
670  iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM DUE TO POOR PROGRESS\n" << endi;
671  broadcast->minimizeCoefficient.publish(minSeq++,0.);
672  broadcast->minimizeCoefficient.publish(minSeq++,0.); // v = f
673  newDir = 1;
674  old_f_dot_f = min_f_dot_f;
675  min_f_dot_v = min_f_dot_f;
677  continue;
678  }
679  hi = last;
680 #define PRINT_BRACKET \
681  iout << "LINE MINIMIZER BRACKET: DX " \
682  << (mid.x-lo.x) << " " << (hi.x-mid.x) << \
683  " DU " << (mid.u-lo.u) << " " << (hi.u-mid.u) << " DUDX " << \
684  lo.dudx << " " << mid.dudx << " " << hi.dudx << " \n" << endi;
686  // converge on minimum on line
687  int itcnt;
688  for ( itcnt = 10; itcnt > 0; --itcnt ) {
689  int progress = 1;
690  // select new position
691  if ( mid.noGradient ) {
692  if ( ( mid.x - lo.x ) > ( hi.x - mid.x ) ) { // subdivide left side
693  x = (1.0 - goldenRatio) * lo.x + goldenRatio * mid.x;
694  MOVETO(x)
695  if ( last.u <= mid.u ) {
696  hi = mid; mid = last;
697  } else {
698  lo = last;
699  }
700  } else { // subdivide right side
701  x = (1.0 - goldenRatio) * hi.x + goldenRatio * mid.x;
702  MOVETO(x)
703  if ( last.u <= mid.u ) {
704  lo = mid; mid = last;
705  } else {
706  hi = last;
707  }
708  }
709  } else {
710  if ( mid.dudx > 0. ) { // subdivide left side
711  BigReal altxhi = 0.1 * lo.x + 0.9 * mid.x;
712  BigReal altxlo = 0.9 * lo.x + 0.1 * mid.x;
713  x = mid.dudx*(mid.x*mid.x-lo.x*lo.x) + 2*mid.x*(lo.u-mid.u);
714  BigReal xdiv = 2*(mid.dudx*(mid.x-lo.x)+(lo.u-mid.u));
715  if ( xdiv ) x /= xdiv; else x = last.x;
716  if ( x > altxhi ) x = altxhi;
717  if ( x < altxlo ) x = altxlo;
718  if ( x-last.x == 0 ) break;
719  MOVETO(x)
720  if ( last.u <= mid.u ) {
721  hi = mid; mid = last;
722  } else {
723  if ( lo.dudx < 0. && last.dudx > 0. ) progress = 0;
724  lo = last;
725  }
726  } else { // subdivide right side
727  BigReal altxlo = 0.1 * hi.x + 0.9 * mid.x;
728  BigReal altxhi = 0.9 * hi.x + 0.1 * mid.x;
729  x = mid.dudx*(mid.x*mid.x-hi.x*hi.x) + 2*mid.x*(hi.u-mid.u);
730  BigReal xdiv = 2*(mid.dudx*(mid.x-hi.x)+(hi.u-mid.u));
731  if ( xdiv ) x /= xdiv; else x = last.x;
732  if ( x < altxlo ) x = altxlo;
733  if ( x > altxhi ) x = altxhi;
734  if ( x-last.x == 0 ) break;
735  MOVETO(x)
736  if ( last.u <= mid.u ) {
737  lo = mid; mid = last;
738  } else {
739  if ( hi.dudx > 0. && last.dudx < 0. ) progress = 0;
740  hi = last;
741  }
742  }
743  }
745  if ( ! progress ) {
746  MOVETO(mid.x);
747  break;
748  }
749  if ( fabs(last.dudx) < tol ) break;
750  if ( lo.x != mid.x && (lo.u-mid.u) < tol * (mid.x-lo.x) ) break;
751  if ( mid.x != hi.x && (hi.u-mid.u) < tol * (hi.x-mid.x) ) break;
752  }
753  // new direction
754  broadcast->minimizeCoefficient.publish(minSeq++,0.);
755  BigReal c = min_f_dot_f / old_f_dot_f;
756  c = ( c > 1.5 ? 1.5 : c );
757  if ( atStart ) { c = 0; --atStart; }
758  if ( c*c*min_v_dot_v > errorFactor*min_f_dot_f ) {
759  c = 0;
760  if ( errorFactor < 100 ) errorFactor += 10;
761  }
762  if ( c == 0 ) {
763  iout << "MINIMIZER RESTARTING CONJUGATE GRADIENT ALGORITHM\n" << endi;
764  }
765  broadcast->minimizeCoefficient.publish(minSeq++,c); // v = c*v+f
766  newDir = 1;
767  old_f_dot_f = min_f_dot_f;
768  min_f_dot_v = c * min_f_dot_v + min_f_dot_f;
769  min_v_dot_v = c*c*min_v_dot_v + 2*c*min_f_dot_v + min_f_dot_f;
770  }
771 
772 }
773 
774 #undef MOVETO
775 #undef CALCULATE
776 
777 // NOTE: Only isotropic case implemented here!
778 void Controller::multigratorPressure(int step, int callNumber) {
783  BigReal NGfac = 1.0/sqrt(3.0*NG + 1.0);
786  {
787  // Compute new scaling factors and send them to Sequencer
788  BigReal V = state->lattice.volume();
789  BigReal Pinst = trace(controlPressure)/3.0;
790  BigReal PGsum = trace(momentumSqrSum);
791  //
792  multigratorXiT = multigratorXi + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
793  BigReal scale = exp(s*NGfac*multigratorXiT);
794  BigReal velScale = exp(-s*NGfac*(1.0 + 1.0/NG)*multigratorXiT);
795  // fprintf(stderr, "%d | T %lf P %lf V %1.3lf\n", step, temperature, Pinst*PRESSUREFACTOR, V);
796  Tensor scaleTensor = Tensor::identity(scale);
797  Tensor volScaleTensor = Tensor::identity(scale);
798  Tensor velScaleTensor = Tensor::identity(velScale);
799  state->lattice.rescale(volScaleTensor);
800  if (callNumber == 1) {
801  broadcast->positionRescaleFactor.publish(step,scaleTensor);
802  broadcast->velocityRescaleTensor.publish(step,velScaleTensor);
803  } else {
804  broadcast->positionRescaleFactor2.publish(step,scaleTensor);
805  broadcast->velocityRescaleTensor2.publish(step,velScaleTensor);
806  }
807  }
808 
809  {
810  // Wait here for Sequencer to finish scaling and force computation
811  reduction->require();
812  Tensor virial_normal;
813  Tensor virial_nbond;
814  Tensor virial_slow;
815  Tensor intVirial_normal;
816  Tensor intVirial_nbond;
817  Tensor intVirial_slow;
818  Vector extForce_normal;
819  Vector extForce_nbond;
820  Vector extForce_slow;
821  GET_TENSOR(momentumSqrSum, reduction, REDUCTION_MOMENTUM_SQUARED);
822  GET_TENSOR(virial_normal, reduction, REDUCTION_VIRIAL_NORMAL);
823  GET_TENSOR(virial_nbond, reduction, REDUCTION_VIRIAL_NBOND);
824  GET_TENSOR(virial_slow, reduction, REDUCTION_VIRIAL_SLOW);
825  GET_TENSOR(intVirial_normal, reduction, REDUCTION_INT_VIRIAL_NORMAL);
826  GET_TENSOR(intVirial_nbond, reduction, REDUCTION_INT_VIRIAL_NBOND);
827  GET_TENSOR(intVirial_slow, reduction, REDUCTION_INT_VIRIAL_SLOW);
828  GET_VECTOR(extForce_normal, reduction, REDUCTION_EXT_FORCE_NORMAL);
829  GET_VECTOR(extForce_nbond, reduction, REDUCTION_EXT_FORCE_NBOND);
830  GET_VECTOR(extForce_slow, reduction, REDUCTION_EXT_FORCE_SLOW);
831  calcPressure(step, 0, virial_normal, virial_nbond, virial_slow,
832  intVirial_normal, intVirial_nbond, intVirial_slow,
833  extForce_normal, extForce_nbond, extForce_slow);
834  if (callNumber == 2) {
835  // Update temperature for the Temperature Cycle that is coming next
837  temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
838  }
839  }
840 
841  {
842  // Update pressure integrator
843  BigReal V = state->lattice.volume();
844  BigReal Pinst = trace(controlPressure)/3.0;
845  BigReal PGsum = trace(momentumSqrSum);
846  //
847  multigratorXi = multigratorXiT + 0.5*s*NGfac/kT0*( 3.0*V*(Pinst - P0) + PGsum/NG );
848  }
849 
850  }
851 }
852 
853 void Controller::multigratorTemperature(int step, int callNumber) {
858  BigReal Nf = numDegFreedom;
860  BigReal T1, T2, v;
861  {
862  T1 = temperature;
863  BigReal kTinst = BOLTZMANN * temperature;
864  for (int i=n-1;i >= 0;i--) {
865  if (i == 0) {
866  BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
867  multigratorNuT[0] = exp(-0.5*t*NuOmega)*multigratorNu[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
868  } else if (i == n-1) {
869  multigratorNuT[n-1] = multigratorNu[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
870  } else {
871  BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
872  multigratorNuT[i] = exp(-0.5*t*NuOmega)*multigratorNu[i] +
873  0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
874  }
875  }
876  BigReal velScale = exp(-t*multigratorNuT[0]/multigratorOmega[0]);
877  v = velScale;
878  if (callNumber == 1)
879  broadcast->velocityRescaleFactor.publish(step,velScale);
880  else
881  broadcast->velocityRescaleFactor2.publish(step,velScale);
882  }
883 
884  {
885  // Wait here for Sequencer to finish scaling and re-calculating kinetic energy
888  temperature = 2.0 * kineticEnergy / ( numDegFreedom * BOLTZMANN );
889  T2 = temperature;
890  if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
891  // If this is pressure cycle, receive new momentum product
892  GET_TENSOR(momentumSqrSum, multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED);
893  }
894  }
895 
896  // fprintf(stderr, "%d | T %lf scale %lf T' %lf\n", step, T1, v, T2);
897 
898  {
899  BigReal kTinst = BOLTZMANN * temperature;
900  for (int i=0;i < n;i++) {
901  if (i == 0) {
902  BigReal NuOmega = (n > 1) ? multigratorNuT[1]/multigratorOmega[1] : 0.0;
903  multigratorNu[0] = exp(-0.5*t*NuOmega)*multigratorNuT[0] + 0.5*Nf*t*exp(-0.25*t*NuOmega)*(kTinst - kT0);
904  } else if (i == n-1) {
905  multigratorNu[n-1] = multigratorNuT[n-1] + 0.5*t*(pow(multigratorNu[n-2],2.0)/multigratorOmega[n-2] - kT0);
906  } else {
907  BigReal NuOmega = multigratorNuT[i+1]/multigratorOmega[i+1];
908  multigratorNu[i] = exp(-0.5*t*NuOmega)*multigratorNuT[i] +
909  0.5*t*exp(-0.25*t*NuOmega)*(pow(multigratorNu[i-1],2.0)/multigratorOmega[i-1] - kT0);
910  }
912  }
913  }
914 
915  }
916 }
917 
918 // Calculates Enthalpy for multigrator
919 BigReal Controller::multigatorCalcEnthalpy(BigReal potentialEnergy, int step, int minimize) {
920  Node *node = Node::Object();
921  Molecule *molecule = node->molecule;
922  SimParameters *simParameters = node->simParameters;
923 
924  BigReal V = state->lattice.volume();
928  BigReal Nf = numDegFreedom;
930  BigReal sumZeta = 0.0;
931  for (int i=1;i < simParams->multigratorNoseHooverChainLength;i++) {
932  sumZeta += multigratorZeta[i];
933  }
934  BigReal nuOmegaSum = 0.0;
935  for (int i=0;i < simParams->multigratorNoseHooverChainLength;i++) {
936  nuOmegaSum += multigratorNu[i]*multigratorNu[i]/(2.0*multigratorOmega[i]);
937  }
938  BigReal W = (3.0*NG + 1.0)*kT0*tauV*tauV;
939  BigReal eta = sqrt(kT0*W)*multigratorXi;
940 
941  BigReal enthalpy = kineticEnergy + potentialEnergy + eta*eta/(2.0*W) + P0*V + nuOmegaSum + kT0*(Nf*multigratorZeta[0] + sumZeta);
942 
943 // if (!(step % 100))
944  // fprintf(stderr, "enthalpy %lf %lf %lf %lf %lf %lf %lf\n", enthalpy,
945  // kineticEnergy, potentialEnergy, eta*eta/(2.0*W), P0*V, nuOmegaSum, kT0*(Nf*multigratorZeta[0] + sumZeta));
946 
947  return enthalpy;
948 }
949 
951 {
955  const int freq = simParams->berendsenPressureFreq;
956  if ( ! (berendsenPressure_count % freq) ) {
960  // We only use on-diagonal terms (for now)
961  factor = Tensor::diagonal(diagonal(factor));
963  factor *= ( ( simParams->berendsenPressureCompressibility / 3.0 ) *
965  factor += Tensor::identity(1.0);
966 #define LIMIT_SCALING(VAR,MIN,MAX,FLAG) {\
967  if ( VAR < (MIN) ) { VAR = (MIN); FLAG = 1; } \
968  if ( VAR > (MAX) ) { VAR = (MAX); FLAG = 1; } }
969  int limited = 0;
970  LIMIT_SCALING(factor.xx,1./1.03,1.03,limited)
971  LIMIT_SCALING(factor.yy,1./1.03,1.03,limited)
972  LIMIT_SCALING(factor.zz,1./1.03,1.03,limited)
973 #undef LIMIT_SCALING
974  if ( limited ) {
975  iout << iERROR << "Step " << step <<
976  " cell rescaling factor limited.\n" << endi;
977  }
979  state->lattice.rescale(factor);
980  }
981  } else {
984  }
985 }
986 
988 {
990  {
991  Tensor &strainRate = langevinPiston_strainRate;
992  int cellDims = simParams->useFlexibleCell ? 1 : 3;
993  BigReal dt = simParams->dt;
994  BigReal dt_long = slowFreq * dt;
997  BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
998 
999 #ifdef DEBUG_PRESSURE
1000  iout << iINFO << "entering langevinPiston1, strain rate: " << strainRate << "\n";
1001 #endif
1002 
1003  if ( ! ( (step-1) % slowFreq ) )
1004  {
1005  BigReal gamma = 1 / simParams->langevinPistonDecay;
1006  BigReal f1 = exp( -0.5 * dt_long * gamma );
1007  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1008  strainRate *= f1;
1009  if ( simParams->useFlexibleCell ) {
1010  // We only use on-diagonal terms (for now)
1011  if ( simParams->useConstantRatio ) {
1012  BigReal r = f2 * random->gaussian();
1013  strainRate.xx += r;
1014  strainRate.yy += r;
1015  strainRate.zz += f2 * random->gaussian();
1016  } else {
1017  strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
1018  }
1019  } else {
1020  strainRate += f2 * Tensor::identity(random->gaussian());
1021  }
1022 
1023 #ifdef DEBUG_PRESSURE
1024  iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
1025 #endif
1026  }
1027 
1028  // Apply surface tension. If surfaceTensionTarget is zero, we get
1029  // the default (isotropic pressure) case.
1030 
1031  Tensor ptarget;
1032  ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
1033  if ( simParams->surfaceTensionTarget != 0. ) {
1034  ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
1035  simParams->surfaceTensionTarget / state->lattice.c().z;
1036  }
1037 
1038  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1039  ( controlPressure - ptarget );
1040 
1041  if (simParams->fixCellDims) {
1042  if (simParams->fixCellDimX) strainRate.xx = 0;
1043  if (simParams->fixCellDimY) strainRate.yy = 0;
1044  if (simParams->fixCellDimZ) strainRate.zz = 0;
1045  }
1046 
1047 
1048 #ifdef DEBUG_PRESSURE
1049  iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
1050 #endif
1051 
1053  if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
1054  {
1055  // We only use on-diagonal terms (for now)
1056  Tensor factor;
1057  if ( !simParams->useConstantArea ) {
1058  factor.xx = exp( dt_long * strainRate.xx );
1059  factor.yy = exp( dt_long * strainRate.yy );
1060  } else {
1061  factor.xx = factor.yy = 1;
1062  }
1063  factor.zz = exp( dt_long * strainRate.zz );
1064  broadcast->positionRescaleFactor.publish(step,factor);
1065  state->lattice.rescale(factor);
1066 #ifdef DEBUG_PRESSURE
1067  iout << iINFO << "rescaling by: " << factor << "\n";
1068 #endif
1069  }
1070  } else { // langevinPistonBarrier
1071  if ( ! ( (step-1-slowFreq/2) % slowFreq ) )
1072  {
1073  if ( ! positionRescaleFactor.xx ) { // first pass without MTS
1074  // We only use on-diagonal terms (for now)
1075  Tensor factor;
1076  if ( !simParams->useConstantArea ) {
1077  factor.xx = exp( dt_long * strainRate.xx );
1078  factor.yy = exp( dt_long * strainRate.yy );
1079  } else {
1080  factor.xx = factor.yy = 1;
1081  }
1082  factor.zz = exp( dt_long * strainRate.zz );
1083  broadcast->positionRescaleFactor.publish(step,factor);
1084  positionRescaleFactor = factor;
1085  strainRate_old = strainRate;
1086  }
1088 #ifdef DEBUG_PRESSURE
1089  iout << iINFO << "rescaling by: " << positionRescaleFactor << "\n";
1090 #endif
1091  }
1092  if ( ! ( (step-slowFreq/2) % slowFreq ) )
1093  {
1094  // We only use on-diagonal terms (for now)
1095  Tensor factor;
1096  if ( !simParams->useConstantArea ) {
1097  factor.xx = exp( (dt+dt_long) * strainRate.xx - dt * strainRate_old.xx );
1098  factor.yy = exp( (dt+dt_long) * strainRate.yy - dt * strainRate_old.yy );
1099  } else {
1100  factor.xx = factor.yy = 1;
1101  }
1102  factor.zz = exp( (dt+dt_long) * strainRate.zz - dt * strainRate_old.zz );
1103  broadcast->positionRescaleFactor.publish(step+1,factor);
1104  positionRescaleFactor = factor;
1105  strainRate_old = strainRate;
1106  }
1107  }
1108 
1109  // corrections to integrator
1110  if ( ! ( step % nbondFreq ) )
1111  {
1112 #ifdef DEBUG_PRESSURE
1113  iout << iINFO << "correcting strain rate for nbond, ";
1114 #endif
1115  strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1116  ( (nbondFreq - 1) * controlPressure_nbond );
1117 #ifdef DEBUG_PRESSURE
1118  iout << "strain rate: " << strainRate << "\n";
1119 #endif
1120  }
1121  if ( ! ( step % slowFreq ) )
1122  {
1123 #ifdef DEBUG_PRESSURE
1124  iout << iINFO << "correcting strain rate for slow, ";
1125 #endif
1126  strainRate -= ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1127  ( (slowFreq - 1) * controlPressure_slow );
1128 #ifdef DEBUG_PRESSURE
1129  iout << "strain rate: " << strainRate << "\n";
1130 #endif
1131  }
1132  if (simParams->fixCellDims) {
1133  if (simParams->fixCellDimX) strainRate.xx = 0;
1134  if (simParams->fixCellDimY) strainRate.yy = 0;
1135  if (simParams->fixCellDimZ) strainRate.zz = 0;
1136  }
1137 
1138  }
1139 
1140 }
1141 
1143 {
1144  if ( simParams->langevinPistonOn )
1145  {
1146  Tensor &strainRate = langevinPiston_strainRate;
1147  int cellDims = simParams->useFlexibleCell ? 1 : 3;
1148  BigReal dt = simParams->dt;
1149  BigReal dt_long = slowFreq * dt;
1152  BigReal mass = controlNumDegFreedom * kT * tau * tau * cellDims;
1153 
1154  // corrections to integrator
1155  if ( ! ( step % nbondFreq ) )
1156  {
1157 #ifdef DEBUG_PRESSURE
1158  iout << iINFO << "correcting strain rate for nbond, ";
1159 #endif
1160  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1161  ( (nbondFreq - 1) * controlPressure_nbond );
1162 #ifdef DEBUG_PRESSURE
1163  iout << "strain rate: " << strainRate << "\n";
1164 #endif
1165  }
1166  if ( ! ( step % slowFreq ) )
1167  {
1168 #ifdef DEBUG_PRESSURE
1169  iout << iINFO << "correcting strain rate for slow, ";
1170 #endif
1171  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1172  ( (slowFreq - 1) * controlPressure_slow );
1173 #ifdef DEBUG_PRESSURE
1174  iout << "strain rate: " << strainRate << "\n";
1175 #endif
1176  }
1177 
1178  // Apply surface tension. If surfaceTensionTarget is zero, we get
1179  // the default (isotropic pressure) case.
1180 
1181  Tensor ptarget;
1182  ptarget.xx = ptarget.yy = ptarget.zz = simParams->langevinPistonTarget;
1183  if ( simParams->surfaceTensionTarget != 0. ) {
1184  ptarget.xx = ptarget.yy = simParams->langevinPistonTarget -
1185  simParams->surfaceTensionTarget / state->lattice.c().z;
1186  }
1187 
1188  strainRate += ( 0.5 * dt * cellDims * state->lattice.volume() / mass ) *
1189  ( controlPressure - ptarget );
1190 
1191 
1192 #ifdef DEBUG_PRESSURE
1193  iout << iINFO << "integrating half step, strain rate: " << strainRate << "\n";
1194 #endif
1195 
1196  if ( ! ( step % slowFreq ) )
1197  {
1198  BigReal gamma = 1 / simParams->langevinPistonDecay;
1199  BigReal f1 = exp( -0.5 * dt_long * gamma );
1200  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kT / mass );
1201  strainRate *= f1;
1202  if ( simParams->useFlexibleCell ) {
1203  // We only use on-diagonal terms (for now)
1204  if ( simParams->useConstantRatio ) {
1205  BigReal r = f2 * random->gaussian();
1206  strainRate.xx += r;
1207  strainRate.yy += r;
1208  strainRate.zz += f2 * random->gaussian();
1209  } else {
1210  strainRate += f2 * Tensor::diagonal(random->gaussian_vector());
1211  }
1212  } else {
1213  strainRate += f2 * Tensor::identity(random->gaussian());
1214  }
1215 #ifdef DEBUG_PRESSURE
1216  iout << iINFO << "applying langevin, strain rate: " << strainRate << "\n";
1217 #endif
1218  }
1219 
1220 #ifdef DEBUG_PRESSURE
1221  iout << iINFO << "exiting langevinPiston2, strain rate: " << strainRate << "\n";
1222 #endif
1223  if (simParams->fixCellDims) {
1224  if (simParams->fixCellDimX) strainRate.xx = 0;
1225  if (simParams->fixCellDimY) strainRate.yy = 0;
1226  if (simParams->fixCellDimZ) strainRate.zz = 0;
1227  }
1228  }
1229 
1230 
1231 }
1232 
1234 {
1235  const int rescaleFreq = simParams->rescaleFreq;
1236  if ( rescaleFreq > 0 ) {
1238  if ( rescaleVelocities_numTemps == rescaleFreq ) {
1240  BigReal rescaleTemp = simParams->rescaleTemp;
1242  (!(simParams->adaptTempLastStep > 0) || step < simParams->adaptTempLastStep )) {
1243  rescaleTemp = adaptTempT;
1244  }
1245  BigReal factor = sqrt(rescaleTemp/avgTemp);
1246  broadcast->velocityRescaleFactor.publish(step,factor);
1247  //iout << "RESCALING VELOCITIES AT STEP " << step
1248  // << " FROM AVERAGE TEMPERATURE OF " << avgTemp
1249  // << " TO " << rescaleTemp << " KELVIN.\n" << endi;
1251  }
1252  }
1253 }
1254 
1256 
1257  Vector momentum;
1258  momentum.x = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_X);
1259  momentum.y = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Y);
1260  momentum.z = reduction->item(REDUCTION_HALFSTEP_MOMENTUM_Z);
1262 
1263  if ( momentum.length2() == 0. )
1264  iout << iERROR << "Momentum is exactly zero; probably a bug.\n" << endi;
1265  if ( mass == 0. )
1266  NAMD_die("Total mass is zero in Controller::correctMomentum");
1267 
1268  momentum *= (-1./mass);
1269 
1270  broadcast->momentumCorrection.publish(step+slowFreq,momentum);
1271 }
1272 
1273 //Modifications for alchemical fep
1275 {
1277  && !simParams->alchLambdaFreq) {
1278  const BigReal alchLambda = simParams->alchLambda;
1279  const BigReal alchLambda2 = simParams->alchLambda2;
1280  const BigReal alchLambdaIDWS = simParams->alchLambdaIDWS;
1281  const BigReal alchTemp = simParams->alchTemp;
1282  const int alchEquilSteps = simParams->alchEquilSteps;
1283  iout << "FEP: RESETTING FOR NEW FEP WINDOW "
1284  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 " << alchLambda2;
1285  if ( alchLambdaIDWS >= 0. ) {
1286  iout << " LAMBDA_IDWS " << alchLambdaIDWS;
1287  }
1288  iout << "\nFEP: WINDOW TO HAVE " << alchEquilSteps
1289  << " STEPS OF EQUILIBRATION PRIOR TO FEP DATA COLLECTION.\n"
1290  << "FEP: USING CONSTANT TEMPERATURE OF " << alchTemp
1291  << " K FOR FEP CALCULATION\n" << endi;
1292  }
1293 }
1295 {
1297  && !simParams->alchLambdaFreq) {
1298  const BigReal alchLambda = simParams->alchLambda;
1299  const int alchEquilSteps = simParams->alchEquilSteps;
1300  iout << "TI: RESETTING FOR NEW WINDOW "
1301  << "LAMBDA SET TO " << alchLambda
1302  << "\nTI: WINDOW TO HAVE " << alchEquilSteps
1303  << " STEPS OF EQUILIBRATION PRIOR TO TI DATA COLLECTION.\n" << endi;
1304  }
1305 }
1306 //fepe
1307 
1309 {
1310  const int reassignFreq = simParams->reassignFreq;
1311  if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1312  BigReal newTemp = simParams->reassignTemp;
1313  newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
1314  if ( simParams->reassignIncr > 0.0 ) {
1315  if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
1316  newTemp = simParams->reassignHold;
1317  } else {
1318  if ( newTemp < simParams->reassignHold )
1319  newTemp = simParams->reassignHold;
1320  }
1321  iout << "REASSIGNING VELOCITIES AT STEP " << step
1322  << " TO " << newTemp << " KELVIN.\n" << endi;
1323  }
1324 }
1325 
1327 {
1328  if ( simParams->tCoupleOn )
1329  {
1330  const BigReal tCoupleTemp = simParams->tCoupleTemp;
1331  BigReal coefficient = 1.;
1332  if ( temperature > 0. ) coefficient = tCoupleTemp/temperature - 1.;
1333  broadcast->tcoupleCoefficient.publish(step,coefficient);
1334  }
1335 }
1336 
1349 {
1350  if ( simParams->stochRescaleOn ) {
1353  double coefficient = stochRescaleCoefficient();
1354  broadcast->stochRescaleCoefficient.publish(step,coefficient);
1355  stochRescale_count = 0;
1356  }
1357  }
1358 }
1359 
1365  const double stochRescaleTemp = simParams->stochRescaleTemp;
1366  double coefficient = 1;
1367  if ( temperature > 0 ) {
1368  double R1 = random->gaussian();
1369  // double gammaShape = 0.5*(numDegFreedom - 1);
1370  // R2sum is the sum of (numDegFreedom - 1) squared normal variables,
1371  // which is chi-squared distributed.
1372  // This is in turn a special case of the Gamma distribution,
1373  // which converges to a normal distribution in the limit of a
1374  // large shape parameter.
1375  // double R2sum = 2*(gammaShape+sqrt(gammaShape)*random->gaussian())+R1*R1;
1376  double R2sum = random->sum_of_squared_gaussians(numDegFreedom-1);
1377  double tempfactor = stochRescaleTemp/(temperature*numDegFreedom);
1378 
1379  coefficient = sqrt(stochRescaleTimefactor +
1380  (1 - stochRescaleTimefactor)*tempfactor*R2sum +
1381  2*R1*sqrt(tempfactor*(1 - stochRescaleTimefactor)*
1383  }
1384  heat += 0.5*numDegFreedom*BOLTZMANN*temperature*(coefficient*coefficient-1);
1385  return coefficient;
1386 }
1387 
1388 static char *FORMAT(BigReal X, int decimal = 4)
1389 {
1390  static char tmp_string[50];
1391  static char format_string[50];
1392  const double maxnum = 99999999999.9999;
1393  if ( X > maxnum ) X = maxnum;
1394  if ( X < -maxnum ) X = -maxnum;
1395 
1396  int whole = (decimal <= 4 ? 14 : 10 + decimal);
1397  sprintf(format_string, " %%%d.%df", whole, decimal);
1398  sprintf(tmp_string, format_string, X);
1399  return tmp_string;
1400 }
1401 
1402 static char *FORMAT(const char *X, int decimal = 4)
1403 {
1404  static char tmp_string[50];
1405  static char format_string[50];
1406 
1407  int width = (decimal <= 4 ? 14 : 10 + decimal);
1408  sprintf(format_string, " %%%ds", width);
1409  sprintf(tmp_string, format_string, X);
1410  return tmp_string;
1411 }
1412 
1413 static char *ETITLE(int X)
1414 {
1415  static char tmp_string[21];
1416  sprintf(tmp_string,"ENERGY: %7d",X);
1417  return tmp_string;
1418 }
1419 
1420 void Controller::receivePressure(int step, int minimize)
1421 {
1422  Node *node = Node::Object();
1423  Molecule *molecule = node->molecule;
1424  SimParameters *simParameters = node->simParameters;
1425  Lattice &lattice = state->lattice;
1426 
1427  reduction->require();
1428 
1429  Tensor virial_normal;
1430  Tensor virial_nbond;
1431  Tensor virial_slow;
1432 #ifdef ALTVIRIAL
1433  Tensor altVirial_normal;
1434  Tensor altVirial_nbond;
1435  Tensor altVirial_slow;
1436 #endif
1437  Tensor intVirial_normal;
1438  Tensor intVirial_nbond;
1439  Tensor intVirial_slow;
1440  Vector extForce_normal;
1441  Vector extForce_nbond;
1442  Vector extForce_slow;
1443 
1444 #if 1
1445  numDegFreedom = molecule->num_deg_freedom();
1446  int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
1447  int numFixedGroups = molecule->num_fixed_groups();
1448  int numFixedAtoms = molecule->num_fixed_atoms();
1449 #endif
1450 #if 0
1451  int numAtoms = molecule->numAtoms;
1452  numDegFreedom = 3 * numAtoms;
1453  int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
1454  int numFixedAtoms =
1455  ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
1456  int numLonepairs = molecule->numLonepairs;
1457  int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
1458  if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
1459  if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
1460  if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
1461  if ( ! ( numFixedAtoms || molecule->numConstraints
1462  || simParameters->comMove || simParameters->langevinOn ) ) {
1463  numDegFreedom -= 3;
1464  numGroupDegFreedom -= 3;
1465  }
1466  if (simParameters->pairInteractionOn) {
1467  // this doesn't attempt to deal with fixed atoms or constraints
1468  numDegFreedom = 3 * molecule->numFepInitial;
1469  }
1470  int numRigidBonds = molecule->numRigidBonds;
1471  int numFixedRigidBonds =
1472  ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
1473 
1474  // numLonepairs is subtracted here because all lonepairs have a rigid bond
1475  // to oxygen, but all of the LP degrees of freedom are dealt with above
1476  numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
1477 #endif
1478 
1481 
1482  BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
1484  BigReal groupKineticEnergyCentered = kineticEnergyCentered -
1486 
1487  BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
1488  / ( numDegFreedom * BOLTZMANN );
1489  BigReal atomTempCentered = 2.0 * kineticEnergyCentered
1490  / ( numDegFreedom * BOLTZMANN );
1491  BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
1492  / ( numGroupDegFreedom * BOLTZMANN );
1493  BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
1494  / ( numGroupDegFreedom * BOLTZMANN );
1495 
1496  /* test code for comparing different temperatures
1497  iout << "TEMPTEST: " << step << " " <<
1498  atomTempHalfstep << " " <<
1499  atomTempCentered << " " <<
1500  groupTempHalfstep << " " <<
1501  groupTempCentered << "\n" << endi;
1502  iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
1503  numGroupDegFreedom << "\n" << endi;
1504  */
1505 
1506  GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
1507 
1508  GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
1509  GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
1510  GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
1511 
1512 #ifdef ALTVIRIAL
1513  GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
1514  GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
1515  GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
1516 #endif
1517 
1518  GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
1519  GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
1520  GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
1521 
1522  GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
1523  GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
1524  GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
1525  // APH NOTE: These four lines are now done in calcPressure()
1526  // Vector extPosition = lattice.origin();
1527  // virial_normal -= outer(extForce_normal,extPosition);
1528  // virial_nbond -= outer(extForce_nbond,extPosition);
1529  // virial_slow -= outer(extForce_slow,extPosition);
1530 
1533 
1534  if (simParameters->drudeOn) {
1535  BigReal drudeComKE
1537  BigReal drudeBondKE
1539  int g_bond = 3 * molecule->numDrudeAtoms;
1540  int g_com = numDegFreedom - g_bond;
1541  temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
1542  drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
1543  }
1544 
1545  // Calculate number of degrees of freedom (controlNumDegFreedom)
1546  if ( simParameters->useGroupPressure )
1547  {
1548  controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
1549  if ( ! ( numFixedAtoms || molecule->numConstraints
1550  || simParameters->comMove || simParameters->langevinOn ) ) {
1551  controlNumDegFreedom -= 1;
1552  }
1553  }
1554  else
1555  {
1557  }
1558  if (simParameters->fixCellDims) {
1559  if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
1560  if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
1561  if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
1562  }
1563 
1564  // Calculate pressure tensors using virials
1565  calcPressure(step, minimize,
1566  virial_normal, virial_nbond, virial_slow,
1567  intVirial_normal, intVirial_nbond, intVirial_slow,
1568  extForce_normal, extForce_nbond, extForce_slow);
1569 
1570 #ifdef DEBUG_PRESSURE
1571  iout << iINFO << "Control pressure = " << controlPressure <<
1572  " = " << controlPressure_normal << " + " <<
1573  controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
1574 #endif
1575  if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
1576  iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
1577  << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
1578  << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
1579  << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
1580  << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
1581  << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
1582  << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
1583  << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
1584  << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
1585  << endi;
1586  }
1587 }
1588 
1589 //
1590 // Calculates all pressure tensors using virials
1591 //
1592 // Sets variables:
1593 // pressure, pressure_normal, pressure_nbond, pressure_slow
1594 // groupPressure, groupPressure_normal, groupPressure_nbond, groupPressure_slow
1595 // controlPressure, controlPressure_normal, controlPressure_nbond, controlPressure_slow
1596 // pressure_amd
1597 //
1598 void Controller::calcPressure(int step, int minimize,
1599  const Tensor& virial_normal_in, const Tensor& virial_nbond_in, const Tensor& virial_slow_in,
1600  const Tensor& intVirial_normal, const Tensor& intVirial_nbond, const Tensor& intVirial_slow,
1601  const Vector& extForce_normal, const Vector& extForce_nbond, const Vector& extForce_slow) {
1602 
1603  Tensor virial_normal = virial_normal_in;
1604  Tensor virial_nbond = virial_nbond_in;
1605  Tensor virial_slow = virial_slow_in;
1606 
1607  // Tensor tmp = virial_normal;
1608  // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
1609  // tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
1610 
1611  Node *node = Node::Object();
1612  Molecule *molecule = node->molecule;
1613  SimParameters *simParameters = node->simParameters;
1614  Lattice &lattice = state->lattice;
1615 
1616  BigReal volume;
1617 
1618  Vector extPosition = lattice.origin();
1619  virial_normal -= outer(extForce_normal,extPosition);
1620  virial_nbond -= outer(extForce_nbond,extPosition);
1621  virial_slow -= outer(extForce_slow,extPosition);
1622 
1623  if ( (volume=lattice.volume()) != 0. )
1624  {
1625 
1626  if (simParameters->LJcorrection && volume) {
1627 #ifdef MEM_OPT_VERSION
1628  NAMD_bug("LJcorrection not supported in memory optimized build.");
1629 #else
1630  // Apply tail correction to pressure
1631  BigReal alchLambda = simParameters->getCurrentLambda(step);
1632  virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
1633 #endif
1634  }
1635 
1636  // kinetic energy component included in virials
1637  pressure_normal = virial_normal / volume;
1638  groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
1639 
1640  if (simParameters->accelMDOn) {
1641  pressure_amd = virial_amd / volume;
1644  }
1645 
1646  if ( minimize || ! ( step % nbondFreq ) )
1647  {
1648  pressure_nbond = virial_nbond / volume;
1649  groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
1650  }
1651 
1652  if ( minimize || ! ( step % slowFreq ) )
1653  {
1654  pressure_slow = virial_slow / volume;
1655  groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
1656  }
1657 
1661  }
1662  else
1663  {
1664  pressure = Tensor();
1665  groupPressure = Tensor();
1666  }
1667 
1668  if ( simParameters->useGroupPressure )
1669  {
1674  }
1675  else
1676  {
1681  }
1682 
1683  if ( simParameters->useFlexibleCell ) {
1684  // use symmetric pressure to control rotation
1685  // controlPressure_normal = symmetric(controlPressure_normal);
1686  // controlPressure_nbond = symmetric(controlPressure_nbond);
1687  // controlPressure_slow = symmetric(controlPressure_slow);
1688  // controlPressure = symmetric(controlPressure);
1689  // only use on-diagonal components for now
1694  if ( simParameters->useConstantRatio ) {
1695 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
1696  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
1701 #undef AVGXY
1702  }
1703  } else {
1710  controlPressure =
1711  Tensor::identity(trace(controlPressure)/3.);
1712  }
1713 }
1714 
1716 (BigReal testV, int n,
1717  BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV){
1718  BigReal delta;
1719 
1720  if(testV > *Vmax){
1721  *Vmax = testV;
1722  }
1723  else if(testV < *Vmin){
1724  *Vmin = testV;
1725  }
1726 
1727  //mean and std calculated by Online Algorithm
1728  delta = testV - *Vavg;
1729  *Vavg += delta / (BigReal)n;
1730  *M2 += delta * (testV - (*Vavg));
1731 
1732  *sigmaV = sqrt(*M2/n);
1733 }
1734 
1736 (int iE, int step, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV,
1737  BigReal* k0, BigReal* k, BigReal* E, int *iEused, char *warn){
1738  switch(iE){
1739  case 2:
1740  *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
1741  // if k0 <=0 OR k0 > 1, switch to iE=1 mode
1742  if(*k0 > 1.)
1743  *warn |= 1;
1744  else if(*k0 <= 0.)
1745  *warn |= 2;
1746  //else stay in iE=2 mode
1747  else{
1748  *E = Vmin + (Vmax-Vmin)/(*k0);
1749  *iEused = 2;
1750  break;
1751  }
1752  case 1:
1753  *E = Vmax;
1754  *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
1755  if(*k0 > 1.)
1756  *k0 = 1.;
1757  *iEused = 1;
1758  break;
1759  }
1760 
1761  *k = *k0 / (Vmax - Vmin);
1762 }
1763 
1765 (BigReal k, BigReal E, BigReal testV, Tensor vir_orig,
1766  BigReal *dV, BigReal *factor, Tensor *vir){
1767  BigReal VE = testV - E;
1768  //if V < E, apply boost
1769  if(VE < 0){
1770  *factor = k * VE;
1771  *vir += vir_orig * (*factor);
1772  *dV += (*factor) * VE/2.;
1773  *factor += 1.;
1774  }
1775  else{
1776  *factor = 1.;
1777  }
1778 }
1779 
1781 (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){
1782  FILE *rest;
1783  char timestepstr[20];
1784  char *filename;
1785  int baselen;
1786  char *restart_name;
1787  const char *bsuffix;
1788 
1789  if(lasttime || simParams->restartFrequency == 0){
1790  filename = simParams->outputFilename;
1791  bsuffix = ".BAK";
1792  }
1793  else{
1794  filename = simParams->restartFilename;
1795  bsuffix = ".old";
1796  }
1797 
1798  baselen = strlen(filename);
1799  restart_name = new char[baselen+26];
1800 
1801  strcpy(restart_name, filename);
1802  if ( simParams->restartSave && !lasttime) {
1803  sprintf(timestepstr,".%d",step_n);
1804  strcat(restart_name, timestepstr);
1805  bsuffix = ".BAK";
1806  }
1807  strcat(restart_name, ".gamd");
1808 
1809  if(write_topic){
1810  NAMD_backup_file(restart_name,bsuffix);
1811 
1812  rest = fopen(restart_name, "w");
1813  if(rest){
1814  fprintf(rest, "# NAMD accelMDG restart file\n"
1815  "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
1816  "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
1817  type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
1818  fclose(rest);
1819  iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
1820  }
1821  else
1822  iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
1823  }
1824  else{
1825  rest = fopen(restart_name, "a");
1826  if(rest){
1827  fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
1828  type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
1829  fclose(rest);
1830  }
1831  else
1832  iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
1833  }
1834 
1835  delete[] restart_name;
1836 }
1837 
1838 
1839 void Controller::rescaleaccelMD(int step, int minimize)
1840 {
1841  if ( !simParams->accelMDOn ) return;
1842 
1844 
1845  // copy all to submit_reduction
1846  for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
1848  }
1850 
1851  if (step == simParams->firstTimestep) {
1852  accelMDdVAverage = 0;
1853  accelMDdV = 0;
1854  }
1855 // if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
1856  if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
1857 
1858  Node *node = Node::Object();
1859  Molecule *molecule = node->molecule;
1860  Lattice &lattice = state->lattice;
1861 
1862  const BigReal accelMDE = simParams->accelMDE;
1863  const BigReal accelMDalpha = simParams->accelMDalpha;
1864  const BigReal accelMDTE = simParams->accelMDTE;
1865  const BigReal accelMDTalpha = simParams->accelMDTalpha;
1866  const int accelMDOutFreq = simParams->accelMDOutFreq;
1867 
1868  //GaMD
1869  static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
1870  static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
1871  static BigReal k0P, kP, EP;
1872  static BigReal k0D, kD, ED;
1873  static int V_n = 1, iEusedD, iEusedP;
1874  static char warnD, warnP;
1875  static int restfreq;
1876 
1879  const int ntcmd = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
1880  const int ntebprep = ntcmd + simParams->accelMDGEquiPrepSteps;
1881  const int nteb = ntcmd + simParams->accelMDGEquiSteps;
1882  const int ntave = simParams->accelMDGStatWindow;
1883 
1884  BigReal volume;
1885  BigReal bondEnergy;
1886  BigReal angleEnergy;
1887  BigReal dihedralEnergy;
1888  BigReal improperEnergy;
1889  BigReal crosstermEnergy;
1890  BigReal boundaryEnergy;
1891  BigReal miscEnergy;
1892  BigReal amd_electEnergy;
1893  BigReal amd_ljEnergy;
1894  BigReal amd_ljEnergy_Corr = 0.;
1895  BigReal amd_electEnergySlow;
1896  BigReal amd_groLJEnergy;
1897  BigReal amd_groGaussEnergy;
1898  BigReal amd_goTotalEnergy;
1899  BigReal potentialEnergy;
1900  BigReal factor_dihe = 1;
1901  BigReal factor_tot = 1;
1902  BigReal testV;
1903  BigReal dV = 0.;
1904  Vector accelMDfactor;
1905  Tensor vir; //auto initialized to 0
1906  Tensor vir_dihe;
1907  Tensor vir_normal;
1908  Tensor vir_nbond;
1909  Tensor vir_slow;
1910 
1911  volume = lattice.volume();
1912 
1913  bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
1914  angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
1915  dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
1916  improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
1917  crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
1918  boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
1919  miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
1920 
1921  GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
1922  GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
1923  GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
1924  GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
1925 
1926  if ( !( step % nbondFreq ) ) {
1927  amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
1928  amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
1929  amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
1930  amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
1933  amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
1934  } else {
1935  amd_electEnergy = electEnergy;
1936  amd_ljEnergy = ljEnergy;
1937  amd_groLJEnergy = groLJEnergy;
1938  amd_groGaussEnergy = groGaussEnergy;
1939  amd_goTotalEnergy = goTotalEnergy;
1940  }
1941 
1942  if( simParams->LJcorrection && volume ) {
1943  // Obtain tail correction (copied from printEnergies())
1944  // This value is only printed for sanity check
1945  // accelMD currently does not 'boost' tail correction
1946  amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
1947  }
1948 
1949  if ( !( step % slowFreq ) ) {
1950  amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
1951  } else {
1952  amd_electEnergySlow = electEnergySlow;
1953  }
1954 
1955  potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
1956  improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
1957  crosstermEnergy + boundaryEnergy + miscEnergy +
1958  amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
1959 
1960  //GaMD
1961  if(simParams->accelMDG){
1962  // if first time running accelMDG module
1963  if(step == firststep) {
1964  iEusedD = iEusedP = simParams->accelMDGiE;
1965  warnD = warnP = '\0';
1966 
1967  //restart file reading
1969  FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
1970  char line[256];
1971  int dihe_n=0, tot_n=0;
1972  if(!rest){
1973  sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
1974  NAMD_die(line);
1975  }
1976 
1977  while(fgets(line, 256, rest)){
1978  if(line[0] == 'D'){
1979  dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
1980  &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
1981  }
1982  else if(line[0] == 'T'){
1983  tot_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
1984  &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
1985  }
1986  }
1987 
1988  fclose(rest);
1989 
1990  V_n++;
1991 
1992  //check dihe settings
1994  if(dihe_n !=8)
1995  NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
1996  k0D = kD * (VmaxD - VminD);
1997  iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
1998  << " Vmax " << VmaxD << " Vmin " << VminD
1999  << " Vavg " << VavgD << " sigmaV " << sigmaVD
2000  << " E " << ED << " k " << kD
2001  << "\n" << endi;
2002  }
2003 
2004  //check tot settings
2006  if(tot_n !=8)
2007  NAMD_die("Invalid total potential energy format in accelMDG restart file");
2008  k0P = kP * (VmaxP - VminP);
2009  iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
2010  << " Vmax " << VmaxP << " Vmin " << VminP
2011  << " Vavg " << VavgP << " sigmaV " << sigmaVP
2012  << " E " << EP << " k " << kP
2013  << "\n" << endi;
2014  }
2015 
2016  iEusedD = (ED == VmaxD) ? 1 : 2;
2017  iEusedP = (EP == VmaxP) ? 1 : 2;
2018  }
2019  //local restfreq follows NAMD restartfreq (default: 500)
2021  }
2022 
2023  //for dihedral potential
2025  testV = dihedralEnergy + crosstermEnergy;
2026 
2027  //write restart file every restartfreq steps
2028  if(step > firststep && step % restfreq == 0)
2029  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2030  true, false);
2031  //write restart file at the end of the simulation
2032  if( simParams->accelMDLastStep > 0 ){
2033  if( step == simParams->accelMDLastStep )
2034  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2035  true, true);
2036  }
2037  else if(step == simParams->N)
2038  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2039  true, true);
2040 
2041  //conventional MD
2042  if(step < ntcmd){
2043  //very first step
2044  if(step == firststep && !simParams->accelMDGRestart){
2045  //initialize stat
2046  VmaxD = VminD = VavgD = testV;
2047  M2D = sigmaVD = 0.;
2048  }
2049  //first step after cmdprep
2050  else if(step == ntcmdprep && ntcmdprep != 0){
2051  //reset stat
2052  VmaxD = VminD = VavgD = testV;
2053  M2D = sigmaVD = 0.;
2054  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
2055  }
2056  //every ntave steps
2057  else if(ntave > 0 && step % ntave == 0){
2058  //update Vmax Vmin
2059  if(testV > VmaxD) VmaxD = testV;
2060  if(testV < VminD) VminD = testV;
2061  //reset avg and std
2062  VavgD = testV;
2063  M2D = sigmaVD = 0.;
2064  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
2065  }
2066  //normal steps
2067  else
2068  calc_accelMDG_mean_std(testV, V_n,
2069  &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2070 
2071  //last cmd step
2072  if(step == ntcmd - 1){
2074  VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2075  }
2076  }
2077  //equilibration
2078  else if(step < nteb){
2079  calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
2080  &dV, &factor_dihe, &vir);
2081 
2082  //first step after cmd
2083  if(step == ntcmd && simParams->accelMDGresetVaftercmd){
2084  //reset stat
2085  VmaxD = VminD = VavgD = testV;
2086  M2D = sigmaVD = 0.;
2087  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
2088  }
2089  //every ntave steps
2090  else if(ntave > 0 && step % ntave == 0){
2091  //update Vmax Vmin
2092  if(testV > VmaxD) VmaxD = testV;
2093  if(testV < VminD) VminD = testV;
2094  //reset avg and std
2095  VavgD = testV;
2096  M2D = sigmaVD = 0.;
2097  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
2098  }
2099  else
2100  calc_accelMDG_mean_std(testV, V_n,
2101  &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2102 
2103  //steps after ebprep
2104  if(step >= ntebprep)
2105  if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2107  VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2108  }
2109  //production
2110  else{
2111  calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
2112  &dV, &factor_dihe, &vir);
2113  }
2114 
2115  }
2116  //for total potential
2118  Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
2119  testV = potentialEnergy;
2120  if(simParams->accelMDdual){
2121  testV -= dihedralEnergy + crosstermEnergy;
2122  vir_tot -= vir_dihe;
2123  }
2124 
2125  //write restart file every restartfreq steps
2126  if(step > firststep && step % restfreq == 0)
2127  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2128  !simParams->accelMDdual, false);
2129  //write restart file at the end of simulation
2130  if( simParams->accelMDLastStep > 0 ){
2131  if( step == simParams->accelMDLastStep )
2132  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2133  !simParams->accelMDdual, true);
2134  }
2135  else if(step == simParams->N)
2136  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2137  !simParams->accelMDdual, true);
2138 
2139  //conventional MD
2140  if(step < ntcmd){
2141  //very first step
2142  if(step == firststep && !simParams->accelMDGRestart){
2143  //initialize stat
2144  VmaxP = VminP = VavgP = testV;
2145  M2P = sigmaVP = 0.;
2146  }
2147  //first step after cmdprep
2148  else if(step == ntcmdprep && ntcmdprep != 0){
2149  //reset stat
2150  VmaxP = VminP = VavgP = testV;
2151  M2P = sigmaVP = 0.;
2152  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
2153  }
2154  //every ntave steps
2155  else if(ntave > 0 && step % ntave == 0){
2156  //update Vmax Vmin
2157  if(testV > VmaxP) VmaxP = testV;
2158  if(testV < VminP) VminP = testV;
2159  //reset avg and std
2160  VavgP = testV;
2161  M2P = sigmaVP = 0.;
2162  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
2163  }
2164  //normal steps
2165  else
2166  calc_accelMDG_mean_std(testV, V_n,
2167  &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2168  //last cmd step
2169  if(step == ntcmd - 1){
2171  VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2172  }
2173  }
2174  //equilibration
2175  else if(step < nteb){
2176  calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
2177  &dV, &factor_tot, &vir);
2178 
2179  //first step after cmd
2180  if(step == ntcmd && simParams->accelMDGresetVaftercmd){
2181  //reset stat
2182  VmaxP = VminP = VavgP = testV;
2183  M2P = sigmaVP = 0.;
2184  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
2185  }
2186  //every ntave steps
2187  else if(ntave > 0 && step % ntave == 0){
2188  //update Vmax Vmin
2189  if(testV > VmaxP) VmaxP = testV;
2190  if(testV < VminP) VminP = testV;
2191  //reset avg and std
2192  VavgP = testV;
2193  M2P = sigmaVP = 0.;
2194  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
2195  }
2196  else
2197  calc_accelMDG_mean_std(testV, V_n,
2198  &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2199 
2200  //steps after ebprep
2201  if(step >= ntebprep)
2202  if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2204  VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2205  }
2206  //production
2207  else{
2208  calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
2209  &dV, &factor_tot, &vir);
2210  }
2211 
2212  }
2213  accelMDdVAverage += dV;
2214 
2215  //first step after ntcmdprep
2216  if((ntcmdprep > 0 && step == ntcmdprep) ||
2217  (step < nteb && ntave > 0 && step % ntave == 0) ||
2218  (simParams->accelMDGresetVaftercmd && step == ntcmd)){
2219  V_n = 1;
2220  }
2221 
2222  if(step < nteb)
2223  V_n++;
2224 
2225  }
2226  //aMD
2227  else{
2228  if (simParams->accelMDdihe) {
2229 
2230  testV = dihedralEnergy + crosstermEnergy;
2231  if ( testV < accelMDE ) {
2232  factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2233  factor_dihe *= factor_dihe;
2234  vir = vir_dihe * (factor_dihe - 1.0);
2235  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2236  accelMDdVAverage += dV;
2237  }
2238 
2239  } else if (simParams->accelMDdual) {
2240 
2241  testV = dihedralEnergy + crosstermEnergy;
2242  if ( testV < accelMDE ) {
2243  factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2244  factor_dihe *= factor_dihe;
2245  vir = vir_dihe * (factor_dihe - 1.0) ;
2246  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2247  }
2248 
2249  testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
2250  if ( testV < accelMDTE ) {
2251  factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
2252  factor_tot *= factor_tot;
2253  vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
2254  dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
2255  }
2256  accelMDdVAverage += dV;
2257 
2258  } else {
2259 
2260  testV = potentialEnergy;
2261  if ( testV < accelMDE ) {
2262  factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
2263  factor_tot *= factor_tot;
2264  vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
2265  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2266  accelMDdVAverage += dV;
2267  }
2268  }
2269  }
2270 
2271  accelMDfactor[0]=factor_dihe;
2272  accelMDfactor[1]=factor_tot;
2273  accelMDfactor[2]=1;
2274  broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
2275  virial_amd = vir;
2276  accelMDdV = dV;
2277 
2278  if ( factor_tot < 0.001 ) {
2279  iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
2280  << "\n" << endi;
2281  }
2282 
2283  if ( ! (step % accelMDOutFreq) ) {
2284  if ( !(step == simParams->firstTimestep) ) {
2285  accelMDdVAverage = accelMDdVAverage/accelMDOutFreq;
2286  }
2287  iout << "ACCELERATED MD: STEP " << step
2288  << " dV " << dV
2289  << " dVAVG " << accelMDdVAverage
2290  << " BOND " << bondEnergy
2291  << " ANGLE " << angleEnergy
2292  << " DIHED " << dihedralEnergy+crosstermEnergy
2293  << " IMPRP " << improperEnergy
2294  << " ELECT " << amd_electEnergy+amd_electEnergySlow
2295  << " VDW " << amd_ljEnergy
2296  << " POTENTIAL " << potentialEnergy;
2297  if(amd_ljEnergy_Corr)
2298  iout << " LJcorr " << amd_ljEnergy_Corr;
2299  iout << "\n" << endi;
2300  if(simParams->accelMDG){
2302  iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD
2303  << " Vmax " << VmaxD << " Vmin " << VminD
2304  << " Vavg " << VavgD << " sigmaV " << sigmaVD
2305  << " E " << ED << " k0 " << k0D << " k " << kD
2306  << "\n" << endi;
2307  if(warnD & 1)
2308  iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
2309  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2310  if(warnD & 2)
2311  iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
2312  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2314  iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP
2315  << " Vmax " << VmaxP << " Vmin " << VminP
2316  << " Vavg " << VavgP << " sigmaV " << sigmaVP
2317  << " E " << EP << " k0 " << k0P << " k " << kP
2318  << "\n" << endi;
2319  if(warnP & 1)
2320  iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
2321  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2322  if(warnP & 2)
2323  iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
2324  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2325  warnD = warnP = '\0';
2326  }
2327 
2328  accelMDdVAverage = 0;
2329 
2330  if (simParams->accelMDDebugOn) {
2331  Tensor p_normal;
2332  Tensor p_nbond;
2333  Tensor p_slow;
2334  Tensor p;
2335  if ( volume != 0. ) {
2336  p_normal = vir_normal/volume;
2337  p_nbond = vir_nbond/volume;
2338  p_slow = vir_slow/volume;
2339  p = vir/volume;
2340  }
2341  iout << " accelMD Scaling Factor: " << accelMDfactor << "\n"
2342  << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
2343  << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
2344  << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
2345  << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n"
2346  << endi;
2347  }
2348  }
2349 }
2350 
2352  if (!simParams->adaptTempOn) return;
2353  iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
2354  adaptTempDtMin = 0;
2355  adaptTempDtMax = 0;
2356  adaptTempAutoDt = false;
2357  if (simParams->adaptTempBins == 0) {
2358  iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
2359  std::ifstream adaptTempRead(simParams->adaptTempInFile);
2360  if (adaptTempRead) {
2361  int readInt;
2362  BigReal readReal;
2363  bool readBool;
2364  // step
2365  adaptTempRead >> readInt;
2366  // Start with min and max temperatures
2367  adaptTempRead >> adaptTempT; // KELVIN
2368  adaptTempRead >> adaptTempBetaMin; // KELVIN
2369  adaptTempRead >> adaptTempBetaMax; // KELVIN
2370  adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
2371  adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
2372  // In case file is manually edited
2373  if (adaptTempBetaMin > adaptTempBetaMax){
2374  readReal = adaptTempBetaMax;
2375  adaptTempBetaMax = adaptTempBetaMin;
2376  adaptTempBetaMin = adaptTempBetaMax;
2377  }
2378  adaptTempRead >> adaptTempBins;
2379  adaptTempRead >> adaptTempCg;
2380  adaptTempRead >> adaptTempDt;
2388  adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
2389  for(int j = 0; j < adaptTempBins; ++j) {
2390  adaptTempRead >> adaptTempPotEnergyAve[j];
2391  adaptTempRead >> adaptTempPotEnergyVar[j];
2392  adaptTempRead >> adaptTempPotEnergySamples[j];
2393  adaptTempRead >> adaptTempPotEnergyAveNum[j];
2394  adaptTempRead >> adaptTempPotEnergyVarNum[j];
2395  adaptTempRead >> adaptTempPotEnergyAveDen[j];
2396  adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
2397  }
2398  adaptTempRead.close();
2399  }
2400  else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
2401  }
2402  else {
2417  if (simParams->langevinOn)
2419  else if (simParams->rescaleFreq > 0)
2421  for(int j = 0; j < adaptTempBins; ++j){
2422  adaptTempPotEnergyAveNum[j] = 0.;
2423  adaptTempPotEnergyAveDen[j] = 0.;
2425  adaptTempPotEnergyVarNum[j] = 0.;
2426  adaptTempPotEnergyVar[j] = 0.;
2427  adaptTempPotEnergyAve[j] = 0.;
2429  }
2430  }
2431  if (simParams->adaptTempAutoDt > 0.0) {
2432  adaptTempAutoDt = true;
2435  }
2436  adaptTempDTave = 0;
2437  adaptTempDTavenum = 0;
2438  iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
2439  iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
2440  iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
2441  iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
2442  iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
2443  if (simParams->adaptTempRestartFreq > 0) {
2447  NAMD_die("Error opening restart file");
2448  }
2449 }
2450 
2453  adaptTempRestartFile.seekp(std::ios::beg);
2454  iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
2455  adaptTempRestartFile << step << " ";
2456  // Start with min and max temperatures
2457  adaptTempRestartFile << adaptTempT<< " "; // KELVIN
2458  adaptTempRestartFile << 1./adaptTempBetaMin << " "; // KELVIN
2459  adaptTempRestartFile << 1./adaptTempBetaMax << " "; // KELVIN
2463  adaptTempRestartFile << "\n" ;
2464  for(int j = 0; j < adaptTempBins; ++j) {
2471  adaptTempRestartFile << "\n";
2472  }
2474  }
2475 }
2476 
2477 void Controller::adaptTempUpdate(int step, int minimize)
2478 {
2479  //Beta = 1./T
2480  if ( !simParams->adaptTempOn ) return;
2481  int j = 0;
2482  if (step == simParams->firstTimestep) {
2483  adaptTempInit(step);
2484  return;
2485  }
2486  if ( minimize || (step < simParams->adaptTempFirstStep ) ||
2487  ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
2488  const int adaptTempOutFreq = simParams->adaptTempOutFreq;
2489  const bool adaptTempDebug = simParams->adaptTempDebug;
2490  //Calculate Current inverse temperature and bin
2491  BigReal adaptTempBeta = 1./adaptTempT;
2492  adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
2493 
2494  if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
2495  iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
2496  << " adaptTempBeta: " << adaptTempBeta
2497  << " adaptTempDBeta: " << adaptTempDBeta
2498  << " betaMin:" << adaptTempBetaMin
2499  << " betaMax: " << adaptTempBetaMax << "\n";
2502 
2503  BigReal potentialEnergy;
2504  BigReal potEnergyAverage;
2505  BigReal potEnergyVariance;
2506  potentialEnergy = totalEnergy-kineticEnergy;
2507 
2508  //calculate new bin average and variance using adaptive averaging
2511  adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
2512 
2515  potEnergyVariance -= potEnergyAverage*potEnergyAverage;
2516 
2517  adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
2518  adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
2519 
2520  // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
2521  // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
2522  // is constant for each bin. This is to estimate <E(beta)> where beta \in
2523  // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
2524  if ( ! ( step % simParams->adaptTempFreq ) ) {
2525  // If adaptTempBin not at the edge, calculate improved average:
2526  if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
2527  // Get Averaging Limits:
2528  BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
2529  BigReal betaPlus;
2530  BigReal betaMinus;
2531  int nMinus =0;
2532  int nPlus = 0;
2533  if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
2534  if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
2535  betaMinus = adaptTempBeta - deltaBeta;
2536  betaPlus = adaptTempBeta + deltaBeta;
2537  BigReal betaMinus2 = betaMinus*betaMinus;
2538  BigReal betaPlus2 = betaPlus*betaPlus;
2539  nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
2540  nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
2541  // Variables for <E(beta)> estimate:
2542  BigReal potEnergyAve0 = 0.0;
2543  BigReal potEnergyAve1 = 0.0;
2544  // Integral terms
2545  BigReal A0 = 0;
2546  BigReal A1 = 0;
2547  BigReal A2 = 0;
2548  //A0 phi_s integral for beta_minus < beta < beta_{i+1}
2549  BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
2550  BigReal DeltaE2Ave;
2551  BigReal DeltaE2AveSum = 0;
2552  BigReal betaj;
2553  for (j = nMinus+1; j <= adaptTempBin; ++j) {
2554  potEnergyAve0 += adaptTempPotEnergyAve[j];
2555  DeltaE2Ave = adaptTempPotEnergyVar[j];
2556  DeltaE2AveSum += DeltaE2Ave;
2557  A0 += j*DeltaE2Ave;
2558  }
2559  A0 *= adaptTempDBeta;
2560  A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
2561  A0 *= adaptTempDBeta;
2562  betaj = adaptTempBetaN[nMinus+1]-betaMinus;
2563  betaj *= betaj;
2564  A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
2565  A0 /= (betaNp1-betaMinus);
2566 
2567  //A1 phi_s integral for beta_{i+1} < beta < beta_plus
2568  DeltaE2AveSum = 0;
2569  for (j = adaptTempBin+1; j < nPlus; j++) {
2570  potEnergyAve1 += adaptTempPotEnergyAve[j];
2571  DeltaE2Ave = adaptTempPotEnergyVar[j];
2572  DeltaE2AveSum += DeltaE2Ave;
2573  A1 += j*DeltaE2Ave;
2574  }
2575  A1 *= adaptTempDBeta;
2576  A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
2577  A1 *= adaptTempDBeta;
2578  if ( nPlus < adaptTempBins ) {
2579  betaj = betaPlus-adaptTempBetaN[nPlus];
2580  betaj *= betaj;
2581  A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
2582  }
2583  A1 /= (betaPlus-betaNp1);
2584 
2585  //A2 phi_t integral for beta_i
2586  A2 = 0.5*adaptTempDBeta*potEnergyVariance;
2587 
2588  // Now calculate a+ and a-
2589  DeltaE2Ave = A0-A1;
2590  BigReal aplus = 0;
2591  BigReal aminus = 0;
2592  if (DeltaE2Ave != 0) {
2593  aplus = (A0-A2)/(A0-A1);
2594  if (aplus < 0) {
2595  aplus = 0;
2596  }
2597  if (aplus > 1) {
2598  aplus = 1;
2599  }
2600  aminus = 1-aplus;
2601  potEnergyAve0 *= adaptTempDBeta;
2602  potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
2603  potEnergyAve0 /= (betaNp1-betaMinus);
2604  potEnergyAve1 *= adaptTempDBeta;
2605  if ( nPlus < adaptTempBins ) {
2606  potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
2607  }
2608  potEnergyAve1 /= (betaPlus-betaNp1);
2609  potEnergyAverage = aminus*potEnergyAve0;
2610  potEnergyAverage += aplus*potEnergyAve1;
2611  }
2612  if (simParams->adaptTempDebug) {
2613  iout << "ADAPTEMP DEBUG:" << "\n"
2614  << " adaptTempBin: " << adaptTempBin << "\n"
2615  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
2616  << " adaptTempBeta: " << adaptTempBeta << "\n"
2617  << " adaptTemp: " << adaptTempT<< "\n"
2618  << " betaMin: " << adaptTempBetaMin << "\n"
2619  << " betaMax: " << adaptTempBetaMax << "\n"
2620  << " gammaAve: " << gammaAve << "\n"
2621  << " deltaBeta: " << deltaBeta << "\n"
2622  << " betaMinus: " << betaMinus << "\n"
2623  << " betaPlus: " << betaPlus << "\n"
2624  << " nMinus: " << nMinus << "\n"
2625  << " nPlus: " << nPlus << "\n"
2626  << " A0: " << A0 << "\n"
2627  << " A1: " << A1 << "\n"
2628  << " A2: " << A2 << "\n"
2629  << " a+: " << aplus << "\n"
2630  << " a-: " << aminus << "\n"
2631  << endi;
2632  }
2633  }
2634  else {
2635  if (simParams->adaptTempDebug) {
2636  iout << "ADAPTEMP DEBUG:" << "\n"
2637  << " adaptTempBin: " << adaptTempBin << "\n"
2638  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
2639  << " adaptTempBeta: " << adaptTempBeta << "\n"
2640  << " adaptTemp: " << adaptTempT<< "\n"
2641  << " betaMin: " << adaptTempBetaMin << "\n"
2642  << " betaMax: " << adaptTempBetaMax << "\n"
2643  << " gammaAve: " << gammaAve << "\n"
2644  << " deltaBeta: N/A\n"
2645  << " betaMinus: N/A\n"
2646  << " betaPlus: N/A\n"
2647  << " nMinus: N/A\n"
2648  << " nPlus: N/A\n"
2649  << " A0: N/A\n"
2650  << " A1: N/A\n"
2651  << " A2: N/A\n"
2652  << " a+: N/A\n"
2653  << " a-: N/A\n"
2654  << endi;
2655  }
2656  }
2657 
2658  //dT is new temperature
2659  BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
2660  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
2661  dT += adaptTempT;
2662 
2663  // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
2664  // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
2665  if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
2667  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
2668  dT += adaptTempT;
2669  // Check again, if not then keep original adaptTempTor assign random.
2670  if ( dT > 1./adaptTempBetaMin ) {
2671  if (!simParams->adaptTempRandom) {
2672  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
2673  // << " K higher than adaptTempTmax."
2674  // << " Keeping temperature at "
2675  // << adaptTempT<< "\n"<< endi;
2676  dT = adaptTempT;
2677  }
2678  else {
2679  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
2680  // << " K higher than adaptTempTmax."
2681  // << " Assigning random temperature in range\n" << endi;
2682  dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
2683  dT = 1./dT;
2684  }
2685  }
2686  else if ( dT < 1./adaptTempBetaMax ) {
2687  if (!simParams->adaptTempRandom) {
2688  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
2689  // << " K lower than adaptTempTmin."
2690  // << " Keeping temperature at " << adaptTempT<< "\n" << endi;
2691  dT = adaptTempT;
2692  }
2693  else {
2694  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
2695  // << " K lower than adaptTempTmin."
2696  // << " Assigning random temperature in range\n" << endi;
2697  dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
2698  dT = 1./dT;
2699  }
2700  }
2701  else if (adaptTempAutoDt) {
2702  //update temperature step size counter
2703  //FOR "TRUE" ADAPTIVE TEMPERING
2704  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
2705  if (adaptTempTdiff > 0) {
2706  adaptTempDTave += adaptTempTdiff;
2708 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
2709  }
2710  if(adaptTempDTavenum == 100){
2711  BigReal Frac;
2713  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
2714  Frac = adaptTempDTave/Frac;
2715  //if average temperature jump is > 3% of temperature range,
2716  //modify jump size to match 3%
2717  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
2718  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
2719  Frac = adaptTempDtMax/Frac;
2720  iout << "ADAPTEMP: Updating adaptTempDt to ";
2721  adaptTempDt *= Frac;
2722  iout << adaptTempDt << "\n" << endi;
2723  }
2724  adaptTempDTave = 0;
2725  adaptTempDTavenum = 0;
2726  }
2727  }
2728  }
2729  else if (adaptTempAutoDt) {
2730  //update temperature step size counter
2731  // FOR "TRUE" ADAPTIVE TEMPERING
2732  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
2733  if (adaptTempTdiff > 0) {
2734  adaptTempDTave += adaptTempTdiff;
2736 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
2737  }
2738  if(adaptTempDTavenum == 100){
2739  BigReal Frac;
2741  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
2742  Frac = adaptTempDTave/Frac;
2743  //if average temperature jump is > 3% of temperature range,
2744  //modify jump size to match 3%
2745  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
2746  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
2747  Frac = adaptTempDtMax/Frac;
2748  iout << "ADAPTEMP: Updating adaptTempDt to ";
2749  adaptTempDt *= Frac;
2750  iout << adaptTempDt << "\n" << endi;
2751  }
2752  adaptTempDTave = 0;
2753  adaptTempDTavenum = 0;
2754 
2755  }
2756 
2757  }
2758 
2759  adaptTempT = dT;
2761  }
2762  adaptTempWriteRestart(step);
2763  if ( ! (step % adaptTempOutFreq) ) {
2764  iout << "ADAPTEMP: STEP " << step
2765  << " TEMP " << adaptTempT
2766  << " ENERGY " << std::setprecision(10) << potentialEnergy
2767  << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
2768  << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
2769  iout << "\n" << endi;
2770  }
2771 
2772 }
2773 
2774 
2775 void Controller::compareChecksums(int step, int forgiving) {
2776  Node *node = Node::Object();
2777  Molecule *molecule = node->molecule;
2778 
2779  // Some basic correctness checking
2780  BigReal checksum, checksum_b;
2781  char errmsg[256];
2782 
2783  checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
2784  if ( ((int)checksum) != molecule->numAtoms ) {
2785  sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
2786  (int)checksum, molecule->numAtoms);
2787  NAMD_bug(errmsg);
2788  }
2789 
2791  if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
2792  if ( ! computeChecksum ) {
2793  computesPartitioned = 0;
2794  } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
2795  computesPartitioned = 1;
2796  } else {
2797  NAMD_bug("Bad global compute count!\n");
2798  }
2799  computeChecksum = ((int)checksum);
2800  }
2801 
2802  checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
2803  if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
2804  sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
2805  (int)checksum, molecule->numCalcBonds);
2806  if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
2807  iout << iWARN << errmsg << endi;
2808  else NAMD_bug(errmsg);
2809  }
2810 
2812  if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
2813  sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
2814  (int)checksum, molecule->numCalcAngles);
2815  if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
2816  iout << iWARN << errmsg << endi;
2817  else NAMD_bug(errmsg);
2818  }
2819 
2821  if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
2822  sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
2823  (int)checksum, molecule->numCalcDihedrals);
2824  if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
2825  iout << iWARN << errmsg << endi;
2826  else NAMD_bug(errmsg);
2827  }
2828 
2830  if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
2831  sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
2832  (int)checksum, molecule->numCalcImpropers);
2833  if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
2834  iout << iWARN << errmsg << endi;
2835  else NAMD_bug(errmsg);
2836  }
2837 
2839  if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
2840  sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
2841  (int)checksum, molecule->numCalcTholes);
2842  if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
2843  iout << iWARN << errmsg << endi;
2844  else NAMD_bug(errmsg);
2845  }
2846 
2848  if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
2849  sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
2850  (int)checksum, molecule->numCalcAnisos);
2851  if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
2852  iout << iWARN << errmsg << endi;
2853  else NAMD_bug(errmsg);
2854  }
2855 
2857  if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
2858  sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
2859  (int)checksum, molecule->numCalcCrossterms);
2860  if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
2861  iout << iWARN << errmsg << endi;
2862  else NAMD_bug(errmsg);
2863  }
2864 
2866  if ( ((int64)checksum) > molecule->numCalcExclusions &&
2867  ( ! simParams->mollyOn || step % slowFreq ) ) {
2868  if ( forgiving )
2869  iout << iWARN << "High global exclusion count ("
2870  << ((int64)checksum) << " vs "
2871  << molecule->numCalcExclusions << "), possible error!\n"
2872  << iWARN << "This warning is not unusual during minimization.\n"
2873  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
2874  else {
2875  char errmsg[256];
2876  sprintf(errmsg, "High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2877  (long long)checksum, (long long)molecule->numCalcExclusions);
2878  NAMD_bug(errmsg);
2879  }
2880  }
2881 
2882  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
2884  if ( forgiving || ! simParams->fullElectFrequency ) {
2885  iout << iWARN << "Low global exclusion count! ("
2886  << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
2887  if ( forgiving ) {
2888  iout << "\n"
2889  << iWARN << "This warning is not unusual during minimization.\n"
2890  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
2891  } else {
2892  iout << " System unstable or pairlistdist or cutoff too small.\n";
2893  }
2894  iout << endi;
2895  } else {
2896  char errmsg[256];
2897  sprintf(errmsg, "Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2898  (long long)checksum, (long long)molecule->numCalcExclusions);
2899  NAMD_bug(errmsg);
2900  }
2901  }
2902 
2903 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
2905  if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
2906  ( ! simParams->mollyOn || step % slowFreq ) ) {
2907  if ( forgiving )
2908  iout << iWARN << "High global CUDA exclusion count ("
2909  << ((int64)checksum) << " vs "
2910  << molecule->numCalcFullExclusions << "), possible error!\n"
2911  << iWARN << "This warning is not unusual during minimization.\n"
2912  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
2913  else {
2914  char errmsg[256];
2915  sprintf(errmsg, "High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2916  (long long)checksum, (long long)molecule->numCalcFullExclusions);
2917  NAMD_bug(errmsg);
2918  }
2919  }
2920 
2921  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
2923  if ( forgiving || ! simParams->fullElectFrequency ) {
2924  iout << iWARN << "Low global CUDA exclusion count! ("
2925  << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
2926  if ( forgiving ) {
2927  iout << "\n"
2928  << iWARN << "This warning is not unusual during minimization.\n"
2929  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
2930  } else {
2931  iout << " System unstable or pairlistdist or cutoff too small.\n";
2932  }
2933  iout << endi;
2934  } else {
2935  char errmsg[256];
2936  sprintf(errmsg, "Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2937  (long long)checksum, (long long)molecule->numCalcFullExclusions);
2938  NAMD_bug(errmsg);
2939  }
2940  }
2941 #endif
2942 
2944  if ( ((int)checksum) && ! marginViolations ) {
2945  iout << iERROR << "Margin is too small for " << ((int)checksum) <<
2946  " atoms during timestep " << step << ".\n" << iERROR <<
2947  "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
2948  }
2949  marginViolations += (int)checksum;
2950 
2952  if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
2953  iout << iINFO <<
2954  "Pairlistdist is too small for " << ((int)checksum) <<
2955  " computes during timestep " << step << ".\n" << endi;
2956  }
2957  if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
2958 
2960  if ( checksum ) {
2961  if ( forgiving )
2962  iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
2963  else NAMD_bug("Stray PME grid charges detected!\n");
2964  }
2965 }
2966 
2967 void Controller::printTiming(int step) {
2968 
2969  if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
2970  {
2971  const double endWTime = CmiWallTimer() - firstWTime;
2972  const double endCTime = CmiTimer() - firstCTime;
2973 
2974  // fflush about once per minute
2975  if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
2976 
2977  const double elapsedW =
2978  (endWTime - startWTime) / simParams->outputTiming;
2979  const double elapsedC =
2980  (endCTime - startCTime) / simParams->outputTiming;
2981 
2982  const double remainingW = elapsedW * (simParams->N - step);
2983  const double remainingW_hours = remainingW / 3600;
2984 
2985  startWTime = endWTime;
2986  startCTime = endCTime;
2987 
2988 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
2989  if ( simParams->outputEnergies < 60 &&
2990  step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
2991  iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
2992  }
2993 #endif
2994  if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
2995  CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
2996  ", %g hours remaining, %f MB of memory in use.\n",
2997  step, endCTime, elapsedC, endWTime, elapsedW,
2998  remainingW_hours, memusage_MB());
2999  if ( fflush_count ) { --fflush_count; fflush(stdout); }
3000  }
3001  }
3002 }
3003 
3005 
3006  rescaleaccelMD(step,1);
3007  receivePressure(step,1);
3008 
3009  Node *node = Node::Object();
3010  Molecule *molecule = node->molecule;
3011  compareChecksums(step,1);
3012 
3013  printEnergies(step,1);
3014 
3020 
3021 }
3022 
3024 
3025  Node *node = Node::Object();
3026  Molecule *molecule = node->molecule;
3027  SimParameters *simParameters = node->simParameters;
3028  Lattice &lattice = state->lattice;
3029 
3030  compareChecksums(step);
3031 
3032  printEnergies(step,0);
3033 }
3034 
3035 void Controller::printEnergies(int step, int minimize)
3036 {
3037  Node *node = Node::Object();
3038  Molecule *molecule = node->molecule;
3039  SimParameters *simParameters = node->simParameters;
3040  Lattice &lattice = state->lattice;
3041 
3042  // Drude model ANISO energy is added into BOND energy
3043  // and THOLE energy is added into ELECT energy
3044 
3045  BigReal bondEnergy;
3046  BigReal angleEnergy;
3047  BigReal dihedralEnergy;
3048  BigReal improperEnergy;
3049  BigReal crosstermEnergy;
3050  BigReal boundaryEnergy;
3051  BigReal miscEnergy;
3052  BigReal potentialEnergy;
3053  BigReal flatEnergy;
3054  BigReal smoothEnergy;
3055  BigReal work;
3056 
3057  Vector momentum;
3058  Vector angularMomentum;
3059  BigReal volume = lattice.volume();
3060 
3061  bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
3062  angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
3063  dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
3064  improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
3065  crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
3066  boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
3067  miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
3068 
3069  if ( minimize || ! ( step % nbondFreq ) )
3070  {
3073 
3074  // JLai
3077 
3081 
3082 //fepb
3086 
3093 //fepe
3094  }
3095 
3096  if ( minimize || ! ( step % slowFreq ) )
3097  {
3099 //fepb
3101 
3106 //fepe
3107  }
3108 
3109  if (simParameters->LJcorrection && volume) {
3110 #ifdef MEM_OPT_VERSION
3111  NAMD_bug("LJcorrection not supported in memory optimized build.");
3112 #else
3113  // Apply tail correction to energy.
3114  BigReal alchLambda = simParameters->getCurrentLambda(step);
3115  ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
3116 
3117  if (simParameters->alchOn) {
3118  if (simParameters->alchFepOn) {
3119  BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
3120  ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
3121  } else if (simParameters->alchThermIntOn) {
3122  ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
3123  ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
3124  }
3125  }
3126 #endif
3127  }
3128 
3129 //fepb BKR - Compute alchemical work if using dynamic lambda. This is here so
3130 // that the cumulative work can be given during a callback.
3131  if (simParameters->alchLambdaFreq > 0) {
3132  if (step <=
3133  simParameters->firstTimestep + simParameters->alchEquilSteps) {
3134  cumAlchWork = 0.0;
3135  }
3136  alchWork = computeAlchWork(step);
3137  cumAlchWork += alchWork;
3138  }
3139 //fepe
3140 
3141  momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
3142  momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
3143  momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
3144  angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
3145  angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
3146  angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
3147 
3148  // Ported by JLai
3149  potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
3150  + improperEnergy + electEnergy + electEnergySlow + ljEnergy
3151  + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy
3153  // End of port
3154  totalEnergy = potentialEnergy + kineticEnergy;
3155  flatEnergy = (totalEnergy
3157  if ( !(step%slowFreq) ) {
3158  // only adjust based on most accurate energies
3160  if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
3161  if ( step != simParams->firstTimestep ) {
3162  smooth2_avg *= 0.9375;
3163  smooth2_avg += 0.0625 * s;
3164  }
3165  }
3166  smoothEnergy = (flatEnergy + smooth2_avg
3168 
3169  // Reset values for accumulated heat and work.
3170  if (step <= simParameters->firstTimestep &&
3171  (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
3172  heat = 0.0;
3174  }
3175  if ( simParameters->outputMomenta && ! minimize &&
3176  ! ( step % simParameters->outputMomenta ) )
3177  {
3178  iout << "MOMENTUM: " << step
3179  << " P: " << momentum
3180  << " L: " << angularMomentum
3181  << "\n" << endi;
3182  }
3183 
3184  if ( simParameters->outputPressure ) {
3187  tavg_count += 1;
3188  if ( minimize || ! ( step % simParameters->outputPressure ) ) {
3189  iout << "PRESSURE: " << step << " "
3190  << PRESSUREFACTOR * pressure << "\n"
3191  << "GPRESSURE: " << step << " "
3192  << PRESSUREFACTOR * groupPressure << "\n";
3193  if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
3194  << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
3195  << "GPRESSAVG: " << step << " "
3197  iout << endi;
3198  pressure_tavg = 0;
3199  groupPressure_tavg = 0;
3200  tavg_count = 0;
3201  }
3202  }
3203 
3204  // pressure profile reductions
3205  if (pressureProfileSlabs) {
3206  const int freq = simParams->pressureProfileFreq;
3207  const int arraysize = 3*pressureProfileSlabs;
3208 
3209  BigReal *total = new BigReal[arraysize];
3210  memset(total, 0, arraysize*sizeof(BigReal));
3211  const int first = simParams->firstTimestep;
3212 
3213  if (ppbonded) ppbonded->getData(first, step, lattice, total);
3214  if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
3215  if (ppint) ppint->getData(first, step, lattice, total);
3216  for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
3218 
3219  if (!(step % freq)) {
3220  // convert NAMD internal virial to pressure in units of bar
3222 
3223  iout << "PRESSUREPROFILE: " << step << " ";
3224  if (step == first) {
3225  // output pressure profile for this step
3226  for (int i=0; i<arraysize; i++) {
3227  iout << total[i] * scalefac << " ";
3228  }
3229  } else {
3230  // output pressure profile averaged over the last count steps.
3231  scalefac /= pressureProfileCount;
3232  for (int i=0; i<arraysize; i++)
3233  iout << pressureProfileAverage[i]*scalefac << " ";
3234  }
3235  iout << "\n" << endi;
3236 
3237  // Clear the average for the next block
3238  memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
3240  }
3241  delete [] total;
3242  }
3243 
3244  if ( step != simParams->firstTimestep || stepInFullRun == 0 ) { // skip repeated first step
3245  if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
3246  int benchPhase = stepInFullRun / simParams->firstLdbStep;
3247  if ( benchPhase > 0 && benchPhase < 10 ) {
3248 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
3249  if ( simParams->outputEnergies < 60 ) {
3250  iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
3251  }
3252 #endif
3253  iout << iINFO;
3254  if ( benchPhase < 4 ) iout << "Initial time: ";
3255  else iout << "Benchmark time: ";
3256  iout << CkNumPes() << " CPUs ";
3257  {
3258  BigReal wallPerStep =
3259  (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
3260  BigReal ns = simParams->dt / 1000000.0;
3261  BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
3262  BigReal daysPerNano = wallPerStep * days / ns;
3263  iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
3264  iout << memusage_MB() << " MB memory\n" << endi;
3265  }
3266  }
3267  startBenchTime = CmiWallTimer();
3268  }
3269  ++stepInFullRun;
3270  }
3271 
3272  printTiming(step);
3273 
3274  Vector pairVDWForce, pairElectForce;
3275  if ( simParameters->pairInteractionOn ) {
3276  GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
3277  GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
3278  }
3279 
3280  // Compute cumulative nonequilibrium work (including shadow work).
3281  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3282  work = totalEnergy - totalEnergy0 - heat;
3283  }
3284 
3285  // callback to Tcl with whatever we can
3286 #ifdef NAMD_TCL
3287 #define CALLBACKDATA(LABEL,VALUE) \
3288  labels << (LABEL) << " "; values << (VALUE) << " ";
3289 #define CALLBACKLIST(LABEL,VALUE) \
3290  labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
3291  if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
3292  std::ostringstream labels, values;
3293 #if CMK_BLUEGENEL
3294  // the normal version below gives a compiler error
3295  values.precision(16);
3296 #else
3297  values << std::setprecision(16);
3298 #endif
3299  CALLBACKDATA("TS",step);
3300  CALLBACKDATA("BOND",bondEnergy);
3301  CALLBACKDATA("ANGLE",angleEnergy);
3302  CALLBACKDATA("DIHED",dihedralEnergy);
3303  CALLBACKDATA("CROSS",crosstermEnergy);
3304  CALLBACKDATA("IMPRP",improperEnergy);
3306  CALLBACKDATA("VDW",ljEnergy);
3307  CALLBACKDATA("BOUNDARY",boundaryEnergy);
3308  CALLBACKDATA("MISC",miscEnergy);
3309  CALLBACKDATA("POTENTIAL",potentialEnergy);
3310  CALLBACKDATA("KINETIC",kineticEnergy);
3311  CALLBACKDATA("TOTAL",totalEnergy);
3312  CALLBACKDATA("TEMP",temperature);
3313  CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
3314  CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
3315  CALLBACKDATA("VOLUME",lattice.volume());
3316  CALLBACKLIST("CELL_A",lattice.a());
3317  CALLBACKLIST("CELL_B",lattice.b());
3318  CALLBACKLIST("CELL_C",lattice.c());
3319  CALLBACKLIST("CELL_O",lattice.origin());
3320  labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
3321  << lattice.b_p() << " " << lattice.c_p() << "} ";
3322  if (simParameters->drudeOn) {
3323  CALLBACKDATA("DRUDEBOND",drudeBondTemp);
3324  }
3325  if ( simParameters->pairInteractionOn ) {
3326  CALLBACKLIST("VDW_FORCE",pairVDWForce);
3327  CALLBACKLIST("ELECT_FORCE",pairElectForce);
3328  }
3329  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3330  CALLBACKLIST("HEAT",heat);
3331  CALLBACKLIST("WORK",work);
3332  }
3333  if (simParameters->alchOn) {
3334  if (simParameters->alchThermIntOn) {
3335  CALLBACKLIST("BOND1", bondedEnergy_ti_1);
3338  CALLBACKLIST("VDW1", ljEnergy_ti_1);
3339  CALLBACKLIST("BOND2", bondedEnergy_ti_2);
3342  CALLBACKLIST("VDW2", ljEnergy_ti_2);
3343  if (simParameters->alchLambdaFreq > 0) {
3344  CALLBACKLIST("CUMALCHWORK", cumAlchWork);
3345  }
3346  } else if (simParameters->alchFepOn) {
3347  CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
3348  improperEnergy + bondedEnergyDiff_f);
3350  CALLBACKLIST("VDW2", ljEnergy_f);
3351  }
3352  }
3353 
3354  labels << '\0'; values << '\0'; // insane but makes Linux work
3355  state->callback_labelstring = labels.str();
3356  state->callback_valuestring = values.str();
3357  // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
3358  }
3359 #undef CALLBACKDATA
3360 #endif
3361 
3363 
3364  temp_avg += temperature;
3365  pressure_avg += trace(pressure)/3.;
3366  groupPressure_avg += trace(groupPressure)/3.;
3367  avg_count += 1;
3368 
3370  ! (step % simParams->outputPairlists) ) {
3371  iout << iINFO << pairlistWarnings <<
3372  " pairlist warnings in past " << simParams->outputPairlists <<
3373  " steps.\n" << endi;
3374  pairlistWarnings = 0;
3375  }
3376 
3377  BigReal enthalpy;
3378  if (simParameters->multigratorOn && ((step % simParameters->outputEnergies) == 0)) {
3379  enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
3380  }
3381 
3382  // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
3383  if ( ! minimize && step % simParameters->outputEnergies ) return;
3384  // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
3385 
3386  if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
3387  IMDEnergies energies;
3388  energies.tstep = step;
3389  energies.T = temp_avg/avg_count;
3390  energies.Etot = totalEnergy;
3391  energies.Epot = potentialEnergy;
3392  energies.Evdw = ljEnergy;
3393  energies.Eelec = electEnergy + electEnergySlow;
3394  energies.Ebond = bondEnergy;
3395  energies.Eangle = angleEnergy;
3396  energies.Edihe = dihedralEnergy + crosstermEnergy;
3397  energies.Eimpr = improperEnergy;
3398  Node::Object()->imd->gather_energies(&energies);
3399  }
3400 
3401  if ( marginViolations ) {
3402  iout << iERROR << marginViolations <<
3403  " margin violations detected since previous energy output.\n" << endi;
3404  }
3405  marginViolations = 0;
3406 
3407  int precision = simParameters->outputEnergiesPrecision;
3408  if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
3409  {
3410  iout << "ETITLE: TS";
3411  iout << FORMAT("BOND", precision);
3412  iout << FORMAT("ANGLE", precision);
3413  iout << FORMAT("DIHED", precision);
3414  if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS", precision);
3415  iout << FORMAT("IMPRP", precision);
3416  iout << " ";
3417  iout << FORMAT("ELECT", precision);
3418  iout << FORMAT("VDW", precision);
3419  iout << FORMAT("BOUNDARY", precision);
3420  iout << FORMAT("MISC", precision);
3421  iout << FORMAT("KINETIC", precision);
3422  iout << " ";
3423  iout << FORMAT("TOTAL", precision);
3424  iout << FORMAT("TEMP", precision);
3425  iout << FORMAT("POTENTIAL", precision);
3426  // iout << FORMAT("TOTAL2", precision);
3427  iout << FORMAT("TOTAL3", precision);
3428  iout << FORMAT("TEMPAVG", precision);
3429  if ( volume != 0. ) {
3430  iout << " ";
3431  iout << FORMAT("PRESSURE", precision);
3432  iout << FORMAT("GPRESSURE", precision);
3433  iout << FORMAT("VOLUME", precision);
3434  iout << FORMAT("PRESSAVG", precision);
3435  iout << FORMAT("GPRESSAVG", precision);
3436  }
3437  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3438  iout << " ";
3439  iout << FORMAT("HEAT", precision);
3440  iout << FORMAT("WORK", precision);
3441  }
3442  if (simParameters->drudeOn) {
3443  iout << " ";
3444  iout << FORMAT("DRUDEBOND", precision);
3445  iout << FORMAT("DRBONDAVG", precision);
3446  }
3447  // Ported by JLai
3448  if (simParameters->goGroPair) {
3449  iout << " ";
3450  iout << FORMAT("GRO_PAIR_LJ", precision);
3451  iout << FORMAT("GRO_PAIR_GAUSS", precision);
3452  }
3453 
3454  if (simParameters->goForcesOn) {
3455  iout << " ";
3456  iout << FORMAT("NATIVE", precision);
3457  iout << FORMAT("NONNATIVE", precision);
3458  //iout << FORMAT("REL_NATIVE", precision);
3459  //iout << FORMAT("REL_NONNATIVE", precision);
3460  iout << FORMAT("GOTOTAL", precision);
3461  //iout << FORMAT("GOAVG", precision);
3462  }
3463  // End of port -- JLai
3464 
3465  if (simParameters->alchOn) {
3466  if (simParameters->alchThermIntOn) {
3467  iout << "\nTITITLE: TS";
3468  iout << FORMAT("BOND1", precision);
3469  iout << FORMAT("ELECT1", precision);
3470  iout << FORMAT("VDW1", precision);
3471  iout << FORMAT("BOND2", precision);
3472  iout << " ";
3473  iout << FORMAT("ELECT2", precision);
3474  iout << FORMAT("VDW2", precision);
3475  if (simParameters->alchLambdaFreq > 0) {
3476  iout << FORMAT("LAMBDA", precision);
3477  iout << FORMAT("ALCHWORK", precision);
3478  iout << FORMAT("CUMALCHWORK", precision);
3479  }
3480  } else if (simParameters->alchFepOn) {
3481  iout << "\nFEPTITLE: TS";
3482  iout << FORMAT("BOND2", precision);
3483  iout << FORMAT("ELECT2", precision);
3484  iout << FORMAT("VDW2", precision);
3485  if (simParameters->alchLambdaFreq > 0) {
3486  iout << FORMAT("LAMBDA", precision);
3487  }
3488  }
3489  }
3490 
3491  iout << "\n\n" << endi;
3492 
3493  if (simParameters->qmForcesOn) {
3494  iout << "QMETITLE: TS";
3495  iout << FORMAT("QMID", precision);
3496  iout << FORMAT("ENERGY", precision);
3497  if (simParameters->PMEOn) iout << FORMAT("PMECORRENERGY", precision);
3498  iout << "\n\n" << endi;
3499  }
3500 
3501  }
3502 
3503  // N.B. HP's aCC compiler merges FORMAT calls in the same expression.
3504  // Need separate statements because data returned in static array.
3505  iout << ETITLE(step);
3506  iout << FORMAT(bondEnergy, precision);
3507  iout << FORMAT(angleEnergy, precision);
3508  if ( simParameters->mergeCrossterms ) {
3509  iout << FORMAT(dihedralEnergy+crosstermEnergy, precision);
3510  } else {
3511  iout << FORMAT(dihedralEnergy, precision);
3512  iout << FORMAT(crosstermEnergy, precision);
3513  }
3514  iout << FORMAT(improperEnergy, precision);
3515  iout << " ";
3516  iout << FORMAT(electEnergy+electEnergySlow, precision);
3517  iout << FORMAT(ljEnergy, precision);
3518  iout << FORMAT(boundaryEnergy, precision);
3519  iout << FORMAT(miscEnergy, precision);
3520  iout << FORMAT(kineticEnergy, precision);
3521  iout << " ";
3522  iout << FORMAT(totalEnergy, precision);
3523  iout << FORMAT(temperature, precision);
3524  iout << FORMAT(potentialEnergy, precision);
3525  // iout << FORMAT(flatEnergy, precision);
3526  iout << FORMAT(smoothEnergy, precision);
3527  iout << FORMAT(temp_avg/avg_count, precision);
3528  if ( volume != 0. )
3529  {
3530  iout << " ";
3531  iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3., precision);
3532  iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3., precision);
3533  iout << FORMAT(volume, precision);
3536  }
3537  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3538  iout << " ";
3539  iout << FORMAT(heat, precision);
3540  iout << FORMAT(work, precision);
3541  }
3542  if (simParameters->drudeOn) {
3543  iout << " ";
3544  iout << FORMAT(drudeBondTemp, precision);
3545  iout << FORMAT(drudeBondTempAvg/avg_count, precision);
3546  }
3547  // Ported by JLai
3548  if (simParameters->goGroPair) {
3549  iout << " ";
3550  iout << FORMAT(groLJEnergy, precision);
3551  iout << FORMAT(groGaussEnergy, precision);
3552  }
3553 
3554  if (simParameters->goForcesOn) {
3555  iout << " ";
3556  iout << FORMAT(goNativeEnergy, precision);
3557  iout << FORMAT(goNonnativeEnergy, precision);
3558  //iout << FORMAT(relgoNativeEnergy, precision);
3559  //iout << FORMAT(relgoNonnativeEnergy, precision);
3560  iout << FORMAT(goTotalEnergy, precision);
3561  //iout << FORMAT("not implemented", precision);
3562  } // End of port -- JLai
3563 
3564  if (simParameters->alchOn) {
3565  if (simParameters->alchThermIntOn) {
3566  iout << "\n";
3567  iout << TITITLE(step);
3568  iout << FORMAT(bondedEnergy_ti_1, precision);
3570  electEnergyPME_ti_1, precision);
3571  iout << FORMAT(ljEnergy_ti_1, precision);
3572  iout << FORMAT(bondedEnergy_ti_2, precision);
3573  iout << " ";
3575  electEnergyPME_ti_2, precision);
3576  iout << FORMAT(ljEnergy_ti_2, precision);
3577  if (simParameters->alchLambdaFreq > 0) {
3578  iout << FORMAT(simParameters->getCurrentLambda(step), precision);
3579  iout << FORMAT(alchWork, precision);
3580  iout << FORMAT(cumAlchWork, precision);
3581  }
3582  } else if (simParameters->alchFepOn) {
3583  iout << "\n";
3584  iout << FEPTITLE2(step);
3585  iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy
3586  + improperEnergy + bondedEnergyDiff_f, precision);
3587  iout << FORMAT(electEnergy_f + electEnergySlow_f, precision);
3588  iout << FORMAT(ljEnergy_f, precision);
3589  if (simParameters->alchLambdaFreq > 0) {
3590  iout << FORMAT(simParameters->getCurrentLambda(step), precision);
3591  }
3592  }
3593  }
3594 
3595  iout << "\n\n" << endi;
3596 
3597 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
3598  char webout[80];
3599  sprintf(webout,"%d %d %d %d",(int)totalEnergy,
3600  (int)(potentialEnergy),
3601  (int)kineticEnergy,(int)temperature);
3602  CApplicationDepositNode0Data(webout);
3603 #endif
3604 
3605  if (simParameters->pairInteractionOn) {
3606  iout << "PAIR INTERACTION:";
3607  iout << " STEP: " << step;
3608  iout << " VDW_FORCE: ";
3609  iout << FORMAT(pairVDWForce.x, precision);
3610  iout << FORMAT(pairVDWForce.y, precision);
3611  iout << FORMAT(pairVDWForce.z, precision);
3612  iout << " ELECT_FORCE: ";
3613  iout << FORMAT(pairElectForce.x, precision);
3614  iout << FORMAT(pairElectForce.y, precision);
3615  iout << FORMAT(pairElectForce.z, precision);
3616  iout << "\n" << endi;
3617  }
3618  drudeBondTempAvg = 0;
3619  temp_avg = 0;
3620  pressure_avg = 0;
3621  groupPressure_avg = 0;
3622  avg_count = 0;
3623 
3624  if ( fflush_count ) {
3625  --fflush_count;
3626  fflush(stdout);
3627  }
3628 }
3629 
3631  Lattice &lattice = state->lattice;
3632  file << "#$LABELS step";
3633  if ( lattice.a_p() ) file << " a_x a_y a_z";
3634  if ( lattice.b_p() ) file << " b_x b_y b_z";
3635  if ( lattice.c_p() ) file << " c_x c_y c_z";
3636  file << " o_x o_y o_z";
3637  if ( simParams->langevinPistonOn ) {
3638  file << " s_x s_y s_z s_u s_v s_w";
3639  }
3640  file << std::endl;
3641 }
3642 
3644  Lattice &lattice = state->lattice;
3645  file.precision(12);
3646  file << step;
3647  if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
3648  if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
3649  if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
3650  file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
3651  if ( simParams->langevinPistonOn ) {
3652  Vector strainRate = diagonal(langevinPiston_strainRate);
3653  Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
3654  file << " " << strainRate.x;
3655  file << " " << strainRate.y;
3656  file << " " << strainRate.z;
3657  file << " " << strainRate2.x;
3658  file << " " << strainRate2.y;
3659  file << " " << strainRate2.z;
3660  }
3661  file << std::endl;
3662 }
3663 
3665 {
3666  if ( Output::coordinateNeeded(timestep) ){
3667  collection->enqueuePositions(timestep,state->lattice);
3668  }
3669  if ( Output::velocityNeeded(timestep) )
3670  collection->enqueueVelocities(timestep);
3671  if ( Output::forceNeeded(timestep) )
3672  collection->enqueueForces(timestep);
3673 }
3674 
3675 //Modifications for alchemical fep
3677  if (simParams->alchOn && simParams->alchFepOn) {
3678  const int stepInRun = step - simParams->firstTimestep;
3679  const int alchEquilSteps = simParams->alchEquilSteps;
3680  const BigReal alchLambda = simParams->alchLambda;
3681  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
3682 
3683  if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
3684  FepNo = 0;
3685  exp_dE_ByRT = 0.0;
3686  net_dE = 0.0;
3687  }
3691 
3692  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
3693  (step / simParams->alchIDWSFreq) % 2 == 1 )) {
3694  // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
3695  FepNo++;
3696  exp_dE_ByRT += exp(-dE/RT);
3697  net_dE += dE;
3698  }
3699 
3700  if (simParams->alchOutFreq) {
3701  if (stepInRun == 0) {
3702  if (!fepFile.is_open()) {
3705  iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
3706  if(alchEnsembleAvg){
3707  fepSum = 0.0;
3708  fepFile << "# STEP Elec "
3709  << "vdW dE dE_avg Temp dG\n"
3710  << "# l l+dl "
3711  << " l l+dl E(l+dl)-E(l)" << std::endl;
3712  }
3713  else{
3714  fepFile << "# STEP Elec "
3715  << "vdW dE Temp\n"
3716  << "# l l+dl "
3717  << " l l+dl E(l+dl)-E(l)" << std::endl;
3718  }
3719  }
3720  if(!step){
3721  fepFile << "#NEW FEP WINDOW: "
3722  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
3723  << simParams->alchLambda2;
3724  if ( simParams->alchLambdaIDWS >= 0. ) {
3725  fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
3726  }
3727  fepFile << std::endl;
3728  }
3729  }
3730  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3731  fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3732  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3733  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3734  }
3735  if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
3736  writeFepEnergyData(step, fepFile);
3737  fepFile.flush();
3738  }
3739  if (alchEnsembleAvg && (step == simParams->N)) {
3740  fepSum = fepSum + dG;
3741  fepFile << "#Free energy change for lambda window [ " << alchLambda
3742  << " " << simParams->alchLambda2 << " ] is " << dG
3743  << " ; net change until now is " << fepSum << std::endl;
3744  fepFile.flush();
3745  }
3746  }
3747  }
3748 }
3749 
3752  const int stepInRun = step - simParams->firstTimestep;
3753  const int alchEquilSteps = simParams->alchEquilSteps;
3754  const int alchLambdaFreq = simParams->alchLambdaFreq;
3755 
3756  if (stepInRun == 0 || stepInRun == alchEquilSteps) {
3757  TiNo = 0;
3758  net_dEdl_bond_1 = 0;
3759  net_dEdl_bond_2 = 0;
3760  net_dEdl_elec_1 = 0;
3761  net_dEdl_elec_2 = 0;
3762  net_dEdl_lj_1 = 0;
3763  net_dEdl_lj_2 = 0;
3764  }
3765  TiNo++;
3774 
3775  // Don't accumulate block averages (you'll get a divide by zero!) or even
3776  // open the TI output if alchOutFreq is zero.
3777  if (simParams->alchOutFreq) {
3778  if (stepInRun == 0 || stepInRun == alchEquilSteps
3779  || (! ((step - 1) % simParams->alchOutFreq))) {
3780  // output of instantaneous dU/dl now replaced with running average
3781  // over last alchOutFreq steps (except for step 0)
3782  recent_TiNo = 0;
3783  recent_dEdl_bond_1 = 0;
3784  recent_dEdl_bond_2 = 0;
3785  recent_dEdl_elec_1 = 0;
3786  recent_dEdl_elec_2 = 0;
3787  recent_dEdl_lj_1 = 0;
3788  recent_dEdl_lj_2 = 0;
3789  recent_alchWork = 0;
3790  }
3791  recent_TiNo++;
3795  + electEnergyPME_ti_1);
3797  + electEnergyPME_ti_2);
3801 
3802  if (stepInRun == 0) {
3803  if (!tiFile.is_open()) {
3806  /* BKR - This has been rather drastically updated to better match
3807  stdout. This was necessary for several reasons:
3808  1) PME global scaling is obsolete (now removed)
3809  2) scaling of bonded terms was added
3810  3) alchemical work is now accumulated when switching is active
3811  */
3812  iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
3813  tiFile << "#TITITLE: TS";
3814  tiFile << FORMAT("BOND1");
3815  tiFile << FORMAT("AVGBOND1");
3816  tiFile << FORMAT("ELECT1");
3817  tiFile << FORMAT("AVGELECT1");
3818  tiFile << " ";
3819  tiFile << FORMAT("VDW1");
3820  tiFile << FORMAT("AVGVDW1");
3821  tiFile << FORMAT("BOND2");
3822  tiFile << FORMAT("AVGBOND2");
3823  tiFile << FORMAT("ELECT2");
3824  tiFile << " ";
3825  tiFile << FORMAT("AVGELECT2");
3826  tiFile << FORMAT("VDW2");
3827  tiFile << FORMAT("AVGVDW2");
3828  if (alchLambdaFreq > 0) {
3829  tiFile << FORMAT("ALCHWORK");
3830  tiFile << FORMAT("CUMALCHWORK");
3831  }
3832  tiFile << std::endl;
3833  }
3834 
3835  if (alchLambdaFreq > 0) {
3836  tiFile << "#ALCHEMICAL SWITCHING ACTIVE "
3837  << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
3838  << "\n#LAMBDA SCHEDULE: "
3839  << "dL: " << simParams->getLambdaDelta()
3840  << " Freq: " << alchLambdaFreq;
3841  }
3842  else {
3843  const BigReal alchLambda = simParams->alchLambda;
3844  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
3845  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
3846  const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
3847  const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
3848  const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
3849  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
3850  tiFile << "#NEW TI WINDOW: "
3851  << "LAMBDA " << alchLambda
3852  << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
3853  << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
3854  << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
3855  << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
3856  }
3857  tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
3858  << std::endl;
3859  }
3860 
3861  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3862  tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3863  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3864  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3865  }
3866  if ((step%simParams->alchOutFreq) == 0) {
3867  writeTiEnergyData(step, tiFile);
3868  tiFile.flush();
3869  }
3870  }
3871  }
3872 }
3873 
3874 /*
3875  Work is accumulated whenever alchLambda changes. In the current scheme we
3876  always increment lambda _first_, then integrate in time. Therefore the work
3877  is wrt the "old" lambda before the increment.
3878 */
3880  // alchemical scaling factors for groups 1/2 at the previous lambda
3881  const BigReal oldLambda = simParams->getCurrentLambda(step-1);
3882  const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
3883  const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
3884  const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
3885  const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
3886  const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
3887  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
3888  // alchemical scaling factors for groups 1/2 at the new lambda
3889  const BigReal alchLambda = simParams->getCurrentLambda(step);
3890  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
3891  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
3892  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
3893  const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
3894  const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
3895  const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda);
3896 
3897  return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
3898  (elec_lambda_12 - elec_lambda_1)*
3900  (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
3901  (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
3902  (elec_lambda_22 - elec_lambda_2)*
3904  (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
3905  );
3906 }
3907 
3912 
3913  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
3914  const int stepInRun = step - simParams->firstTimestep;
3915  const BigReal alchLambda = simParams->alchLambda;
3916 
3917  if(stepInRun){
3918  if ( simParams->alchLambdaIDWS >= 0. &&
3919  (step / simParams->alchIDWSFreq) % 2 == 0 ) {
3920  // IDWS is active and we are on a "backward" timestep
3921  fepFile << FEPTITLE_BACK(step);
3922  } else {
3923  // "forward" timestep
3924  fepFile << FEPTITLE(step);
3925  }
3926  fepFile << FORMAT(eeng);
3927  fepFile << FORMAT(eeng_f);
3928  fepFile << FORMAT(ljEnergy);
3930  fepFile << FORMAT(dE);
3931  if(alchEnsembleAvg){
3932  BigReal dE_avg = net_dE / FepNo;
3933  fepFile << FORMAT(dE_avg);
3934  }
3936  if(alchEnsembleAvg){
3937  dG = -(RT * log(exp_dE_ByRT / FepNo));
3938  fepFile << FORMAT(dG);
3939  }
3940  fepFile << std::endl;
3941  }
3942 }
3943 
3945  tiFile << TITITLE(step);
3950  tiFile << " ";
3956  tiFile << " ";
3960  if (simParams->alchLambdaFreq > 0) {
3963  }
3964  tiFile << std::endl;
3965 }
3966 //fepe
3967 
3969 {
3970 
3971  if ( step >= 0 ) {
3972 
3973  // Write out eXtended System Trajectory (XST) file
3974  if ( simParams->xstFrequency &&
3975  ((step % simParams->xstFrequency) == 0) )
3976  {
3977  if ( ! xstFile.is_open() )
3978  {
3979  iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
3982  while (!xstFile) {
3983  if ( errno == EINTR ) {
3984  CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
3985  xstFile.clear();
3987  continue;
3988  }
3989  char err_msg[257];
3990  sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
3991  NAMD_err(err_msg);
3992  }
3993  xstFile << "# NAMD extended system trajectory file" << std::endl;
3995  }
3997  xstFile.flush();
3998  }
3999 
4000  // Write out eXtended System Configuration (XSC) files
4001  // Output a restart file
4002  if ( simParams->restartFrequency &&
4003  ((step % simParams->restartFrequency) == 0) &&
4004  (step != simParams->firstTimestep) )
4005  {
4006  iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
4007  << step << "\n" << endi;
4008  char fname[140];
4009  const char *bsuffix = ".old";
4010  strcpy(fname, simParams->restartFilename);
4011  if ( simParams->restartSave ) {
4012  char timestepstr[20];
4013  sprintf(timestepstr,".%d",step);
4014  strcat(fname, timestepstr);
4015  bsuffix = ".BAK";
4016  }
4017  strcat(fname, ".xsc");
4018  NAMD_backup_file(fname,bsuffix);
4019  ofstream_namd xscFile(fname);
4020  while (!xscFile) {
4021  if ( errno == EINTR ) {
4022  CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
4023  xscFile.clear();
4024  xscFile.open(fname);
4025  continue;
4026  }
4027  char err_msg[257];
4028  sprintf(err_msg, "Error opening XSC restart file %s",fname);
4029  NAMD_err(err_msg);
4030  }
4031  xscFile << "# NAMD extended system configuration restart file" << std::endl;
4032  writeExtendedSystemLabels(xscFile);
4033  writeExtendedSystemData(step,xscFile);
4034  if (!xscFile) {
4035  char err_msg[257];
4036  sprintf(err_msg, "Error writing XSC restart file %s",fname);
4037  NAMD_err(err_msg);
4038  }
4039  }
4040 
4041  }
4042 
4043  // Output final coordinates
4044  if (step == FILE_OUTPUT || step == END_OF_RUN)
4045  {
4046  int realstep = ( step == FILE_OUTPUT ?
4048  iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
4049  << realstep << "\n" << endi;
4050  static char fname[140];
4051  strcpy(fname, simParams->outputFilename);
4052  strcat(fname, ".xsc");
4053  NAMD_backup_file(fname);
4054  ofstream_namd xscFile(fname);
4055  while (!xscFile) {
4056  if ( errno == EINTR ) {
4057  CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
4058  xscFile.clear();
4059  xscFile.open(fname);
4060  continue;
4061  }
4062  char err_msg[257];
4063  sprintf(err_msg, "Error opening XSC output file %s",fname);
4064  NAMD_err(err_msg);
4065  }
4066  xscFile << "# NAMD extended system configuration output file" << std::endl;
4067  writeExtendedSystemLabels(xscFile);
4068  writeExtendedSystemData(realstep,xscFile);
4069  if (!xscFile) {
4070  char err_msg[257];
4071  sprintf(err_msg, "Error writing XSC output file %s",fname);
4072  NAMD_err(err_msg);
4073  }
4074  }
4075 
4076  // Close trajectory file
4077  if (step == END_OF_RUN) {
4078  if ( xstFile.is_open() ) {
4079  xstFile.close();
4080  iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
4081  }
4082  }
4083 
4084 }
4085 
4086 
4087 void Controller::recvCheckpointReq(const char *key, int task, checkpoint &cp) { // responding replica
4088  switch ( task ) {
4090  if ( ! checkpoints.count(key) ) {
4091  checkpoints[key] = new checkpoint;
4092  }
4093  *checkpoints[key] = cp;
4094  break;
4096  if ( ! checkpoints.count(key) ) {
4097  NAMD_die("Unable to load checkpoint, requested key was never stored.");
4098  }
4099  cp = *checkpoints[key];
4100  break;
4102  if ( ! checkpoints.count(key) ) {
4103  NAMD_die("Unable to swap checkpoint, requested key was never stored.");
4104  }
4105  std::swap(cp,*checkpoints[key]);
4106  break;
4108  if ( ! checkpoints.count(key) ) {
4109  NAMD_die("Unable to free checkpoint, requested key was never stored.");
4110  }
4111  delete checkpoints[key];
4112  checkpoints.erase(key);
4113  break;
4114  }
4115 }
4116 
4117 void Controller::recvCheckpointAck(checkpoint &cp) { // initiating replica
4119  state->lattice = cp.lattice;
4120  *(ControllerState*)this = cp.state;
4121  }
4122  CkpvAccess(_qd)->process();
4123 }
4124 
4125 
4127 {
4128  if ( ! ldbSteps ) {
4130  }
4131  if ( ! --ldbSteps ) {
4132  startBenchTime -= CmiWallTimer();
4133  Node::Object()->outputPatchComputeMaps("before_ldb", step);
4135  startBenchTime += CmiWallTimer();
4136  fflush_count = 3;
4137  }
4138 }
4139 
4140 void Controller::cycleBarrier(int doBarrier, int step) {
4141 #if USE_BARRIER
4142  if (doBarrier) {
4143  broadcast->cycleBarrier.publish(step,1);
4144  CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
4145  CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4146  }
4147 #endif
4148 }
4149 
4150 void Controller::traceBarrier(int turnOnTrace, int step) {
4151  CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4152  CProxy_Node nd(CkpvAccess(BOCclass_group).node);
4153  nd.traceBarrier(turnOnTrace, step);
4154  CthSuspend();
4155 }
4156 
4158  broadcast->traceBarrier.publish(step,1);
4159  CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4160  awaken();
4161 }
4162 
4163 #ifdef MEASURE_NAMD_WITH_PAPI
4164 void Controller::papiMeasureBarrier(int turnOnMeasure, int step){
4165  CkPrintf("Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4166  CProxy_Node nd(CkpvAccess(BOCclass_group).node);
4167  nd.papiMeasureBarrier(turnOnMeasure, step);
4168  CthSuspend();
4169 }
4170 
4171 void Controller::resumeAfterPapiMeasureBarrier(int step){
4172  broadcast->papiMeasureBarrier.publish(step,1);
4173  CkPrintf("Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4174  awaken();
4175 }
4176 #endif
4177 
4179  BackEnd::awaken();
4180  CthFree(thread);
4181  CthSuspend();
4182 }
4183 
ofstream_namd tiFile
Definition: Controller.h:262
static Node * Object()
Definition: Node.h:86
Vector gaussian_vector(void)
Definition: Random.h:208
BigReal adaptTempBetaMin
Definition: Controller.h:317
BigReal adaptTempTmax
int checkpoint_stored
Definition: Controller.h:268
int adaptTempRestartFreq
int numFixedGroups
Definition: Molecule.h:604
Bool accelMDGresetVaftercmd
BigReal berendsenPressureCompressibility
BigReal net_dE
Definition: Controller.h:120
void enqueueCollections(int)
Definition: Controller.C:3664
#define AVGXY(T)
Bool berendsenPressureOn
Tensor controlPressure_slow
Definition: Controller.h:78
BigReal * adaptTempBetaN
Definition: Controller.h:313
char scriptStringArg1[128]
BigReal berendsenPressureRelaxationTime
int adaptTempBins
Definition: Controller.h:320
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:81
#define XXXBIGREAL
Definition: Controller.C:60
int pressureProfileSlabs
Definition: Controller.h:245
BigReal electEnergy_ti_1
Definition: Controller.h:128
#define CALLBACKDATA(LABEL, VALUE)
void recvCheckpointReq(const char *key, int task, checkpoint &cp)
Definition: Controller.C:4087
void cycleBarrier(int, int)
Definition: Controller.C:4140
BigReal net_dEdl_lj_2
Definition: Controller.h:139
void rescaleVelocities(int)
Definition: Controller.C:1233
BigReal smooth2_avg
Definition: Controller.h:36
void rescaleaccelMD(int step, int minimize=0)
Definition: Controller.C:1839
BigReal dG
Definition: Controller.h:121
int nbondFreq
Definition: Controller.h:79
BigReal accelMDE
BigReal min_v_dot_v
Definition: Controller.h:97
static void awaken(void)
Definition: BackEnd.C:316
static int velocityNeeded(int)
Definition: Output.C:365
void write_accelMDG_rest_file(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)
Definition: Controller.C:1781
BigReal getLambdaDelta(void)
int numCalcBonds
Definition: Molecule.h:622
ControllerState state
Definition: Controller.h:274
BigReal adaptTempCgamma
int numCalcAnisos
Definition: Molecule.h:632
void compareChecksums(int, int=0)
Definition: Controller.C:2775
BigReal noGradient
Definition: Controller.C:590
BigReal uniform(void)
Definition: Random.h:109
Tensor groupPressure_nbond
Definition: Controller.h:74
void NAMD_err(const char *err_msg)
Definition: common.C:106
BigReal goNativeEnergy
Definition: Controller.h:109
std::map< std::string, checkpoint * > checkpoints
Definition: Controller.h:276
BigReal langevinPistonTemp
SimpleBroadcastObject< int > traceBarrier
Definition: Broadcasts.h:84
BigReal accelMDLastStep
void calcPressure(int step, int minimize, const Tensor &virial_normal_in, const Tensor &virial_nbond_in, const Tensor &virial_slow_in, const Tensor &intVirial_normal, const Tensor &intVirial_nbond, const Tensor &intVirial_slow, const Vector &extForce_normal, const Vector &extForce_nbond, const Vector &extForce_slow)
Definition: Controller.C:1598
int eventEndOfTimeStep
Definition: Node.C:286
BigReal fepSum
Definition: Controller.h:124
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
BigReal adaptTempTmin
BigReal net_dEdl_lj_1
Definition: Controller.h:138
#define BOLTZMANN
Definition: common.h:47
static int coordinateNeeded(int)
Definition: Output.C:198
Definition: Node.h:78
void printTiMessage(int)
Definition: Controller.C:1294
BigReal temp_avg
Definition: Controller.h:81
#define FILE_OUTPUT
Definition: Output.h:25
IMDOutput * imd
Definition: Node.h:183
int fflush_count
Definition: Controller.h:222
BigReal gaussian(void)
Definition: Random.h:116
BigReal adaptTempT
Definition: Controller.h:314
void adaptTempWriteRestart(int step)
Definition: Controller.C:2451
int numHydrogenGroups
Definition: Molecule.h:600
static __global__ void(const patch_pair *patch_pairs, const atom *atoms, const atom_param *atom_params, const int *vdw_types, unsigned int *plist, float4 *tmpforces, float4 *slow_tmpforces, float4 *forces, float4 *slow_forces, float *tmpvirials, float *slow_tmpvirials, float *virials, float *slow_virials, unsigned int *global_counters, int *force_ready_queue, const unsigned int *overflow_exclusions, const int npatches, const int block_begin, const int total_block_count, int *block_order, exclmask *exclmasks, const int lj_table_size, const float3 lata, const float3 latb, const float3 latc, const float cutoff2, const float plcutoff2, const int doSlow)
SimpleBroadcastObject< Vector > momentumCorrection
Definition: Broadcasts.h:79
ScriptTcl * getScript(void)
Definition: Node.h:192
virtual void algorithm(void)
Definition: Controller.C:281
BigReal adaptTempDt
Definition: Controller.h:323
static PatchMap * Object()
Definition: PatchMap.h:27
BigReal ljEnergy
Definition: Controller.h:106
void calc_accelMDG_mean_std(BigReal testV, int step_n, BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV)
Definition: Controller.C:1716
Tensor controlPressure
Definition: Controller.h:170
void adaptTempInit(int step)
Definition: Controller.C:2351
void minimize()
Definition: Controller.C:594
BigReal totalEnergy0
Definition: Controller.h:164
#define EVAL_MEASURE
Definition: Output.h:27
char adaptTempInFile[NAMD_FILENAME_BUFFER_SIZE]
PressureProfileReduction(int rtag, int numslabs, int numpartitions, const char *myname, int outputfreq)
Definition: Controller.C:74
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
void printMinimizeEnergies(int)
Definition: Controller.C:3004
void integrate(int)
Definition: Controller.C:429
BigReal multigratorPressureTarget
#define PRESSUREFACTOR
Definition: common.h:49
int min_huge_count
Definition: Controller.h:98
int accelMDGEquiPrepSteps
double stochRescaleCoefficient()
Definition: Controller.C:1364
BigReal net_dEdl_bond_2
Definition: Controller.h:135
BigReal surfaceTensionTarget
void writeExtendedSystemLabels(ofstream_namd &file)
Definition: Controller.C:3630
BigReal reassignTemp
BigReal pressure_avg
Definition: Controller.h:82
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
void enqueueVelocities(int seq)
std::vector< BigReal > multigratorOmega
Definition: Controller.h:215
std::ostream & endi(std::ostream &s)
Definition: InfoStream.C:54
BigReal z
Definition: Vector.h:66
BigReal alchLambda2
BigReal accelMDTE
BigReal sum_of_squared_gaussians(int64_t num_gaussians)
Return the sum of i.i.d. squared Gaussians.
Definition: Random.h:178
float Eelec
Definition: imd.h:32
BigReal electEnergySlow_f
Definition: Controller.h:115
#define X
Definition: msm_defn.h:29
BigReal minLineGoal
std::vector< BigReal > multigratorNuT
Definition: Controller.h:214
int berendsenPressureFreq
std::ostream & iWARN(std::ostream &s)
Definition: InfoStream.C:82
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
ofstream_namd & flush()
Definition: fstream_namd.C:17
BigReal recent_alchWork
Definition: Controller.h:150
if(ComputeNonbondedUtil::goMethod==2)
BigReal u
Definition: Controller.C:590
SimpleBroadcastObject< BigReal > adaptTemperature
Definition: Broadcasts.h:86
Tensor groupPressure_tavg
Definition: Controller.h:86
Bool langevinPistonBarrier
BigReal * adaptTempPotEnergyAveNum
Definition: Controller.h:307
static int gotsigint
Definition: Controller.C:421
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:51
Tensor pressure_amd
Definition: Controller.h:71
BigReal accelMDalpha
int pressureProfileSlabs
int stochRescale_count
Definition: Controller.h:192
#define GET_TENSOR(O, R, A)
Definition: ReductionMgr.h:59
int num_fixed_atoms() const
Definition: Molecule.h:498
BigReal tail_corr_ener
Definition: Molecule.h:471
#define CTRL_STK_SZ
Definition: Thread.h:12
BigReal adaptTempDBeta
Definition: Controller.h:321
BigReal recent_dEdl_bond_2
Definition: Controller.h:145
void enqueuePositions(int seq, Lattice &lattice)
BigReal electEnergySlow
Definition: Controller.h:105
Vector origin() const
Definition: Lattice.h:262
SimpleBroadcastObject< BigReal > tcoupleCoefficient
Definition: Broadcasts.h:76
void outputPatchComputeMaps(const char *filename, int tag)
Definition: Node.C:1521
static int forceNeeded(int)
Definition: Output.C:460
BigReal alchLambda
Bool pairInteractionOn
void swap(Array< T > &s, Array< T > &t)
Definition: MsmMap.h:319
void split(int iStream, int numStreams)
Definition: Random.h:77
RequireReduction * amd_reduction
Definition: Controller.h:238
char outputFilename[NAMD_FILENAME_BUFFER_SIZE]
void multigratorTemperature(int step, int callNumber)
Definition: Controller.C:853
int64_t num_deg_freedom(int isInitialReport=0) const
Definition: Molecule.h:524
void gather_energies(IMDEnergies *energies)
Definition: IMDOutput.C:24
float Eimpr
Definition: imd.h:36
BigReal min_energy
Definition: Controller.h:94
BigReal heat
Definition: Controller.h:162
#define MOVETO(X)
Definition: Controller.C:560
float Etot
Definition: imd.h:29
Tensor pressure_slow
Definition: Controller.h:70
BigReal accelMDTalpha
SimpleBroadcastObject< BigReal > stochRescaleCoefficient
Definition: Broadcasts.h:77
void langevinPiston1(int)
Definition: Controller.C:987
BigReal alchWork
Definition: Controller.h:151
void outputTiEnergy(int step)
Definition: Controller.C:3750
static char * FEPTITLE2(int X)
Definition: Controller.h:357
BigReal electEnergy_f
Definition: Controller.h:114
int numCalcCrossterms
Definition: Molecule.h:626
ofstream_namd xstFile
Definition: Controller.h:252
Tensor langevinPiston_strainRate
Definition: Controller.h:33
BigReal langevinPistonDecay
double memusage_MB()
Definition: memusage.h:13
BigReal temperature
Definition: Controller.h:161
int pressureProfileCount
Definition: Controller.h:246
BigReal alchLambdaIDWS
std::vector< BigReal > multigratorNu
Definition: Controller.h:213
NamdState *const state
Definition: Controller.h:236
BigReal adaptTempDTavenum
Definition: Controller.h:316
BigReal adaptTempDTave
Definition: Controller.h:315
Lattice checkpoint_lattice
Definition: Controller.h:269
PressureProfileReduction * ppint
Definition: Controller.h:244
void rescale(Tensor factor)
Definition: Lattice.h:60
ofstream_namd adaptTempRestartFile
Definition: Controller.h:327
int adaptTempBin
Definition: Controller.h:319
BigReal berendsenPressureTarget
BigReal ljEnergy_ti_1
Definition: Controller.h:130
BigReal adaptTempDtMin
Definition: Controller.h:325
ControllerState checkpoint_state
Definition: Controller.h:270
Tensor pressure
Definition: Controller.h:167
Controller(NamdState *s)
Definition: Controller.C:128
BigReal stochRescalePeriod
Definition: Random.h:37
void require(void)
Definition: ReductionMgr.h:341
void writeExtendedSystemData(int step, ofstream_namd &file)
Definition: Controller.C:3643
ControllerBroadcasts * broadcast
Definition: Controller.h:251
BigReal computeAlchWork(const int step)
Definition: Controller.C:3879
float Edihe
Definition: imd.h:35
Tensor groupPressure_slow
Definition: Controller.h:75
BigReal * adaptTempPotEnergyVar
Definition: Controller.h:311
BigReal getEnergyTailCorr(const BigReal, const int)
BigReal groupPressure_avg
Definition: Controller.h:83
int numLonepairs
Number of lone pairs.
Definition: Molecule.h:581
BigReal drudeBondTempAvg
Definition: Controller.h:156
int64_t numDegFreedom
Definition: Controller.h:101
std::vector< BigReal > multigratorZeta
Definition: Controller.h:216
void stochRescaleVelocities(int)
Definition: Controller.C:1348
int64 numCalcFullExclusions
Definition: Molecule.h:628
int64_t num_group_deg_freedom() const
Definition: Molecule.h:511
int accelMDGcMDPrepSteps
BigReal rescaleTemp
void NAMD_bug(const char *err_msg)
Definition: common.C:129
BigReal recent_dEdl_elec_2
Definition: Controller.h:147
BigReal drudeBondTemp
Definition: Controller.h:155
BigReal bondedEnergyDiff_f
Definition: Controller.h:113
void correctMomentum(int step)
Definition: Controller.C:1255
BigReal recent_dEdl_elec_1
Definition: Controller.h:146
BigReal accelMDFirstStep
float Eangle
Definition: imd.h:34
BigReal getBondLambda(const BigReal)
void berendsenPressure(int)
Definition: Controller.C:950
BigReal adaptTempDtMax
Definition: Controller.h:326
zVector strainRate2
int numFixedRigidBonds
Definition: Molecule.h:606
int recent_TiNo
Definition: Controller.h:152
static char * ETITLE(int X)
Definition: Controller.C:1413
bool is_open() const
Definition: fstream_namd.h:30
void awaken(void)
Definition: Controller.h:45
void rebalance(Sequencer *seq, PatchID id)
BigReal getCurrentLambda2(const int)
SimpleBroadcastObject< Tensor > velocityRescaleTensor2
Definition: Broadcasts.h:72
int numFepInitial
Definition: Molecule.h:608
void langevinPiston2(int)
Definition: Controller.C:1142
BigReal kineticEnergyCentered
Definition: Controller.h:160
void printTiming(int)
Definition: Controller.C:2967
int computeChecksum
Definition: Controller.h:89
BigReal accelMDdV
Definition: Controller.h:50
BigReal electEnergySlow_ti_1
Definition: Controller.h:129
char adaptTempRestartFile[NAMD_FILENAME_BUFFER_SIZE]
int marginViolations
Definition: Controller.h:90
BigReal langevinTemp
SubmitReduction * submit_reduction
Definition: Controller.h:239
int numCalcDihedrals
Definition: Molecule.h:624
BigReal groLJEnergy
Definition: Controller.h:107
void recvCheckpointAck(checkpoint &cp)
Definition: Controller.C:4117
BigReal rescaleVelocities_sumTemps
Definition: Controller.h:174
SimpleBroadcastObject< int > scriptBarrier
Definition: Broadcasts.h:83
BigReal scriptArg1
BigReal x
Definition: Vector.h:66
int numFixedAtoms
Definition: Molecule.h:597
BigReal goTotalEnergy
Definition: Controller.h:111
float T
Definition: imd.h:28
BigReal x
Definition: Controller.C:590
int berendsenPressure_count
Definition: Controller.h:35
int numCalcImpropers
Definition: Molecule.h:625
SimpleBroadcastObject< BigReal > velocityRescaleFactor2
Definition: Broadcasts.h:73
RequireReduction * multigratorReduction
Definition: Controller.h:217
BigReal electEnergyPME_ti_1
Definition: Controller.h:141
Bool adaptTempAutoDt
Definition: Controller.h:324
int numAtoms
Definition: Molecule.h:557
void receivePressure(int step, int minimize=0)
Definition: Controller.C:1420
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t float const PatchPairRecord *__restrict__ const int *__restrict__ const int2 *__restrict__ const unsigned int *__restrict__ unsigned int *__restrict__ int *__restrict__ int *__restrict__ TileListStat *__restrict__ const BoundingBox *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ const int numPatches
BigReal dE
Definition: Controller.h:119
char alchOutFile[NAMD_FILENAME_BUFFER_SIZE]
#define END_OF_RUN
Definition: Output.h:26
int multigratorNoseHooverChainLength
BigReal multigratorTemperatureTarget
#define GET_VECTOR(O, R, A)
Definition: ReductionMgr.h:54
BigReal multigratorPressureRelaxationTime
BigReal groGaussEnergy
Definition: Controller.h:108
void NAMD_die(const char *err_msg)
Definition: common.C:85
Tensor groupPressure_normal
Definition: Controller.h:73
BigReal exp_dE_ByRT
Definition: Controller.h:118
void run(void)
Definition: Controller.C:268
static char * FORMAT(BigReal X)
Definition: ComputeQM.C:367
static LdbCoordinator * Object()
BigReal reassignIncr
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
BigReal getVirialTailCorr(const BigReal)
Tensor virial_amd
Definition: Controller.h:72
int32 tstep
Definition: imd.h:27
BigReal min_f_dot_f
Definition: Controller.h:95
void terminate(void)
Definition: Controller.C:4178
int ldbSteps
Definition: Controller.h:220
SimpleBroadcastObject< BigReal > velocityRescaleFactor
Definition: Broadcasts.h:68
zVector strainRate
void publish(int tag, const T &t)
SimpleBroadcastObject< BigReal > minimizeCoefficient
Definition: Broadcasts.h:78
Bool pressureProfileEwaldOn
float Epot
Definition: imd.h:30
BigReal volume(void) const
Definition: Lattice.h:277
void reassignVelocities(int)
Definition: Controller.C:1308
SimpleBroadcastObject< Vector > accelMDRescaleFactor
Definition: Broadcasts.h:85
int checkpoint_task
Definition: Controller.h:277
Random * random
Definition: Controller.h:234
void getData(int firsttimestep, int step, const Lattice &lattice, BigReal *total)
Definition: Controller.C:90
Tensor outer(const Vector &v1, const Vector &v2)
Definition: Tensor.h:241
void printEnergies(int step, int minimize)
Definition: Controller.C:3035
BigReal multigratorTemperatureRelaxationTime
BigReal accelMDdVAverage
Definition: Controller.h:301
BigReal xx
Definition: Tensor.h:17
BigReal adaptTempDt
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
float Evdw
Definition: imd.h:31
float Ebond
Definition: imd.h:33
void NAMD_backup_file(const char *filename, const char *extension)
Definition: common.C:167
BigReal accelMDGSigma0P
BigReal item(int i) const
Definition: ReductionMgr.h:340
char restartFilename[NAMD_FILENAME_BUFFER_SIZE]
int numConstraints
Definition: Molecule.h:589
virtual ~Controller(void)
Definition: Controller.C:248
BigReal dudx
Definition: Controller.C:590
BigReal multigatorCalcEnthalpy(BigReal potentialEnergy, int step, int minimize)
Definition: Controller.C:919
BigReal langevinPistonPeriod
static char * FEPTITLE(int X)
Definition: Controller.h:343
BigReal electEnergy_ti_2
Definition: Controller.h:131
BigReal zz
Definition: Tensor.h:19
int64 numCalcExclusions
Definition: Molecule.h:627
BigReal net_dEdl_elec_1
Definition: Controller.h:136
void outputFepEnergy(int step)
Definition: Controller.C:3676
BigReal length2(void) const
Definition: Vector.h:173
int pairlistWarnings
Definition: Controller.h:91
unsigned int randomSeed
BigReal alchTemp
BigReal initialTemp
CollectionMaster *const collection
Definition: Controller.h:249
long long int64
Definition: common.h:34
BigReal * adaptTempPotEnergyAveDen
Definition: Controller.h:308
int num_fixed_groups() const
Definition: Molecule.h:504
int pressureProfileAtomTypes
#define simParams
Definition: Output.C:127
BigReal * pressureProfileAverage
Definition: Controller.h:247
RequireReduction * min_reduction
Definition: Controller.h:60
BigReal recent_dEdl_lj_1
Definition: Controller.h:148
#define CALCULATE
Definition: Controller.C:553
char accelMDGRestartFile[NAMD_FILENAME_BUFFER_SIZE]
void printFepMessage(int)
Definition: Controller.C:1274
BigReal stochRescaleTimefactor
Definition: Controller.h:195
BigReal * adaptTempPotEnergyVarNum
Definition: Controller.h:309
void calc_accelMDG_force_factor(BigReal k, BigReal E, BigReal testV, Tensor vir_orig, BigReal *dV, BigReal *factor, Tensor *vir)
Definition: Controller.C:1765
BigReal electEnergyPME_ti_2
Definition: Controller.h:142
#define PRINT_BRACKET
BigReal cumAlchWork
Definition: Controller.h:140
Definition: Tensor.h:15
BigReal adaptTempCg
Definition: Controller.h:322
Tensor pressure_normal
Definition: Controller.h:68
Tensor strainRate_old
Definition: Controller.h:205
Tensor pressure_tavg
Definition: Controller.h:85
int rescaleVelocities_numTemps
Definition: Controller.h:175
BigReal y
Definition: Vector.h:66
Vector b() const
Definition: Lattice.h:253
BigReal electEnergy
Definition: Controller.h:104
void calc_accelMDG_E_k(int iE, int V_n, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV, BigReal *k0, BigReal *k, BigReal *E, int *iEused, char *warn)
Definition: Controller.C:1736
int getNumStepsToRun(void)
Tensor groupPressure
Definition: Controller.h:168
BigReal getCurrentLambda(const int)
BigReal ljEnergy_f
Definition: Controller.h:116
int slowFreq
Definition: Controller.h:80
void traceBarrier(int, int)
Definition: Controller.C:4150
Bool pressureProfileOn
int * adaptTempPotEnergySamples
Definition: Controller.h:312
BigReal yy
Definition: Tensor.h:18
PressureProfileReduction * ppnonbonded
Definition: Controller.h:243
void printDynamicsEnergies(int)
Definition: Controller.C:3023
void multigratorPressure(int step, int callNumber)
Definition: Controller.C:778
BigReal tCoupleTemp
#define LIMIT_SCALING(VAR, MIN, MAX, FLAG)
int multigratorPressureFreq
BigReal totalEnergy
Definition: Controller.h:103
BigReal kineticEnergyHalfstep
Definition: Controller.h:159
RequireReduction * willRequire(int setID, int size=-1)
Definition: ReductionMgr.C:526
BigReal electEnergySlow_ti_2
Definition: Controller.h:132
BigReal getVdwLambda(const BigReal)
void tcoupleVelocities(int)
Definition: Controller.C:1326
BigReal adaptTempBetaMax
Definition: Controller.h:318
void enqueueForces(int seq)
static Tensor symmetric(const Vector &v1, const Vector &v2)
Definition: Tensor.h:45
SimParameters *const simParams
Definition: Controller.h:235
Tensor controlPressure_normal
Definition: Controller.h:76
BigReal net_dEdl_elec_2
Definition: Controller.h:137
static char * TITITLE(int X)
Definition: Controller.h:364
int numDrudeAtoms
Number of Drude particles.
Definition: Molecule.h:582
void submit(void)
Definition: ReductionMgr.h:323
BigReal adaptTempAutoDt
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:83
void resumeAfterTraceBarrier(int)
Definition: Controller.C:4157
BigReal recent_dEdl_lj_2
Definition: Controller.h:149
char xstFilename[NAMD_FILENAME_BUFFER_SIZE]
#define CALLBACKLIST(LABEL, VALUE)
void sendCheckpointReq(int remote, const char *key, int task, Lattice &lat, ControllerState &cs)
Definition: Node.C:1280
BigReal accelMDGSigma0D
BigReal bondedEnergy_ti_2
Definition: Controller.h:127
valT values[SORTTILELISTSKERNEL_ITEMS_PER_THREAD]
int tavg_count
Definition: Controller.h:87
BigReal reassignHold
int b_p() const
Definition: Lattice.h:274
RequireReduction * reduction
Definition: Controller.h:237
PressureProfileReduction * ppbonded
Definition: Controller.h:242
BigReal * adaptTempPotEnergyAve
Definition: Controller.h:310
void writeTiEnergyData(int step, ofstream_namd &file)
Definition: Controller.C:3944
Lattice origLattice
Definition: Controller.h:281
int numCalcTholes
Definition: Molecule.h:631
gridSize x
ofstream_namd fepFile
Definition: Controller.h:258
void rebalanceLoad(int)
Definition: Controller.C:4126
SimpleBroadcastObject< Tensor > positionRescaleFactor2
Definition: Broadcasts.h:74
void writeFepEnergyData(int step, ofstream_namd &file)
Definition: Controller.C:3908
void open(const char *_fname, std::ios_base::openmode _mode=std::ios_base::out)
Definition: fstream_namd.C:11
int numRigidBonds
Definition: Molecule.h:605
BigReal net_dEdl_bond_1
Definition: Controller.h:134
void outputExtendedSystem(int step)
Definition: Controller.C:3968
SimpleBroadcastObject< Tensor > velocityRescaleTensor
Definition: Broadcasts.h:71
int a_p() const
Definition: Lattice.h:273
BigReal bondedEnergy_ti_1
Definition: Controller.h:126
Tensor positionRescaleFactor
Definition: Controller.h:206
Tensor langevinPiston_origStrainRate
Definition: Controller.h:204
BigReal goNonnativeEnergy
Definition: Controller.h:110
void(* namd_sighandler_t)(int)
Definition: Controller.C:426
BigReal langevinPistonTarget
Molecule * molecule
Definition: Node.h:176
int controlNumDegFreedom
Definition: Controller.h:169
Vector a() const
Definition: Lattice.h:252
int numCalcAngles
Definition: Molecule.h:623
Tensor berendsenPressure_avg
Definition: Controller.h:34
int stepInFullRun
Definition: Controller.h:102
int outputEnergiesPrecision
int avg_count
Definition: Controller.h:84
static char * FEPTITLE_BACK(int X)
Definition: Controller.h:350
BigReal multigratorXiT
Definition: Controller.h:210
Tensor pressure_nbond
Definition: Controller.h:69
BigReal multigratorXi
Definition: Controller.h:209
Vector c() const
Definition: Lattice.h:254
int multigratorTemperatureFreq
Tensor controlPressure_nbond
Definition: Controller.h:77
BigReal stochRescaleTemp
#define FORCE_OUTPUT
Definition: Output.h:28
BigReal kineticEnergy
Definition: Controller.h:158
void adaptTempUpdate(int step, int minimize=0)
Definition: Controller.C:2477
Tensor momentumSqrSum
Definition: Controller.h:211
BigReal getElecLambda(const BigReal)
double BigReal
Definition: common.h:114
BigReal min_f_dot_v
Definition: Controller.h:96
int c_p() const
Definition: Lattice.h:275
static void my_sigint_handler(int sig)
Definition: Controller.C:422
BigReal ljEnergy_ti_2
Definition: Controller.h:133
BigReal recent_dEdl_bond_1
Definition: Controller.h:144