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)
1389 {
1390  static char tmp_string[25];
1391  const double maxnum = 99999999999.9999;
1392  if ( X > maxnum ) X = maxnum;
1393  if ( X < -maxnum ) X = -maxnum;
1394  sprintf(tmp_string," %14.4f",X);
1395  return tmp_string;
1396 }
1397 
1398 static char *FORMAT(const char *X)
1399 {
1400  static char tmp_string[25];
1401  sprintf(tmp_string," %14s",X);
1402  return tmp_string;
1403 }
1404 
1405 static char *ETITLE(int X)
1406 {
1407  static char tmp_string[21];
1408  sprintf(tmp_string,"ENERGY: %7d",X);
1409  return tmp_string;
1410 }
1411 
1412 void Controller::receivePressure(int step, int minimize)
1413 {
1414  Node *node = Node::Object();
1415  Molecule *molecule = node->molecule;
1416  SimParameters *simParameters = node->simParameters;
1417  Lattice &lattice = state->lattice;
1418 
1419  reduction->require();
1420 
1421  Tensor virial_normal;
1422  Tensor virial_nbond;
1423  Tensor virial_slow;
1424 #ifdef ALTVIRIAL
1425  Tensor altVirial_normal;
1426  Tensor altVirial_nbond;
1427  Tensor altVirial_slow;
1428 #endif
1429  Tensor intVirial_normal;
1430  Tensor intVirial_nbond;
1431  Tensor intVirial_slow;
1432  Vector extForce_normal;
1433  Vector extForce_nbond;
1434  Vector extForce_slow;
1435 
1436 #if 1
1437  numDegFreedom = molecule->num_deg_freedom();
1438  int64_t numGroupDegFreedom = molecule->num_group_deg_freedom();
1439  int numFixedGroups = molecule->num_fixed_groups();
1440  int numFixedAtoms = molecule->num_fixed_atoms();
1441 #endif
1442 #if 0
1443  int numAtoms = molecule->numAtoms;
1444  numDegFreedom = 3 * numAtoms;
1445  int numGroupDegFreedom = 3 * molecule->numHydrogenGroups;
1446  int numFixedAtoms =
1447  ( simParameters->fixedAtomsOn ? molecule->numFixedAtoms : 0 );
1448  int numLonepairs = molecule->numLonepairs;
1449  int numFixedGroups = ( numFixedAtoms ? molecule->numFixedGroups : 0 );
1450  if ( numFixedAtoms ) numDegFreedom -= 3 * numFixedAtoms;
1451  if (numLonepairs) numDegFreedom -= 3 * numLonepairs;
1452  if ( numFixedGroups ) numGroupDegFreedom -= 3 * numFixedGroups;
1453  if ( ! ( numFixedAtoms || molecule->numConstraints
1454  || simParameters->comMove || simParameters->langevinOn ) ) {
1455  numDegFreedom -= 3;
1456  numGroupDegFreedom -= 3;
1457  }
1458  if (simParameters->pairInteractionOn) {
1459  // this doesn't attempt to deal with fixed atoms or constraints
1460  numDegFreedom = 3 * molecule->numFepInitial;
1461  }
1462  int numRigidBonds = molecule->numRigidBonds;
1463  int numFixedRigidBonds =
1464  ( simParameters->fixedAtomsOn ? molecule->numFixedRigidBonds : 0 );
1465 
1466  // numLonepairs is subtracted here because all lonepairs have a rigid bond
1467  // to oxygen, but all of the LP degrees of freedom are dealt with above
1468  numDegFreedom -= ( numRigidBonds - numFixedRigidBonds - numLonepairs);
1469 #endif
1470 
1473 
1474  BigReal groupKineticEnergyHalfstep = kineticEnergyHalfstep -
1476  BigReal groupKineticEnergyCentered = kineticEnergyCentered -
1478 
1479  BigReal atomTempHalfstep = 2.0 * kineticEnergyHalfstep
1480  / ( numDegFreedom * BOLTZMANN );
1481  BigReal atomTempCentered = 2.0 * kineticEnergyCentered
1482  / ( numDegFreedom * BOLTZMANN );
1483  BigReal groupTempHalfstep = 2.0 * groupKineticEnergyHalfstep
1484  / ( numGroupDegFreedom * BOLTZMANN );
1485  BigReal groupTempCentered = 2.0 * groupKineticEnergyCentered
1486  / ( numGroupDegFreedom * BOLTZMANN );
1487 
1488  /* test code for comparing different temperatures
1489  iout << "TEMPTEST: " << step << " " <<
1490  atomTempHalfstep << " " <<
1491  atomTempCentered << " " <<
1492  groupTempHalfstep << " " <<
1493  groupTempCentered << "\n" << endi;
1494  iout << "Number of degrees of freedom: " << numDegFreedom << " " <<
1495  numGroupDegFreedom << "\n" << endi;
1496  */
1497 
1498  GET_TENSOR(momentumSqrSum,reduction,REDUCTION_MOMENTUM_SQUARED);
1499 
1500  GET_TENSOR(virial_normal,reduction,REDUCTION_VIRIAL_NORMAL);
1501  GET_TENSOR(virial_nbond,reduction,REDUCTION_VIRIAL_NBOND);
1502  GET_TENSOR(virial_slow,reduction,REDUCTION_VIRIAL_SLOW);
1503 
1504 #ifdef ALTVIRIAL
1505  GET_TENSOR(altVirial_normal,reduction,REDUCTION_ALT_VIRIAL_NORMAL);
1506  GET_TENSOR(altVirial_nbond,reduction,REDUCTION_ALT_VIRIAL_NBOND);
1507  GET_TENSOR(altVirial_slow,reduction,REDUCTION_ALT_VIRIAL_SLOW);
1508 #endif
1509 
1510  GET_TENSOR(intVirial_normal,reduction,REDUCTION_INT_VIRIAL_NORMAL);
1511  GET_TENSOR(intVirial_nbond,reduction,REDUCTION_INT_VIRIAL_NBOND);
1512  GET_TENSOR(intVirial_slow,reduction,REDUCTION_INT_VIRIAL_SLOW);
1513 
1514  GET_VECTOR(extForce_normal,reduction,REDUCTION_EXT_FORCE_NORMAL);
1515  GET_VECTOR(extForce_nbond,reduction,REDUCTION_EXT_FORCE_NBOND);
1516  GET_VECTOR(extForce_slow,reduction,REDUCTION_EXT_FORCE_SLOW);
1517  // APH NOTE: These four lines are now done in calcPressure()
1518  // Vector extPosition = lattice.origin();
1519  // virial_normal -= outer(extForce_normal,extPosition);
1520  // virial_nbond -= outer(extForce_nbond,extPosition);
1521  // virial_slow -= outer(extForce_slow,extPosition);
1522 
1525 
1526  if (simParameters->drudeOn) {
1527  BigReal drudeComKE
1529  BigReal drudeBondKE
1531  int g_bond = 3 * molecule->numDrudeAtoms;
1532  int g_com = numDegFreedom - g_bond;
1533  temperature = 2.0 * drudeComKE / (g_com * BOLTZMANN);
1534  drudeBondTemp = (g_bond!=0 ? (2.*drudeBondKE/(g_bond*BOLTZMANN)) : 0.);
1535  }
1536 
1537  // Calculate number of degrees of freedom (controlNumDegFreedom)
1538  if ( simParameters->useGroupPressure )
1539  {
1540  controlNumDegFreedom = molecule->numHydrogenGroups - numFixedGroups;
1541  if ( ! ( numFixedAtoms || molecule->numConstraints
1542  || simParameters->comMove || simParameters->langevinOn ) ) {
1543  controlNumDegFreedom -= 1;
1544  }
1545  }
1546  else
1547  {
1549  }
1550  if (simParameters->fixCellDims) {
1551  if (simParameters->fixCellDimX) controlNumDegFreedom -= 1;
1552  if (simParameters->fixCellDimY) controlNumDegFreedom -= 1;
1553  if (simParameters->fixCellDimZ) controlNumDegFreedom -= 1;
1554  }
1555 
1556  // Calculate pressure tensors using virials
1557  calcPressure(step, minimize,
1558  virial_normal, virial_nbond, virial_slow,
1559  intVirial_normal, intVirial_nbond, intVirial_slow,
1560  extForce_normal, extForce_nbond, extForce_slow);
1561 
1562 #ifdef DEBUG_PRESSURE
1563  iout << iINFO << "Control pressure = " << controlPressure <<
1564  " = " << controlPressure_normal << " + " <<
1565  controlPressure_nbond << " + " << controlPressure_slow << "\n" << endi;
1566 #endif
1567  if ( simParams->accelMDOn && simParams->accelMDDebugOn && ! (step % simParameters->accelMDOutFreq) ) {
1568  iout << " PRESS " << trace(pressure)*PRESSUREFACTOR/3. << "\n"
1569  << " PNORMAL " << trace(pressure_normal)*PRESSUREFACTOR/3. << "\n"
1570  << " PNBOND " << trace(pressure_nbond)*PRESSUREFACTOR/3. << "\n"
1571  << " PSLOW " << trace(pressure_slow)*PRESSUREFACTOR/3. << "\n"
1572  << " PAMD " << trace(pressure_amd)*PRESSUREFACTOR/3. << "\n"
1573  << " GPRESS " << trace(groupPressure)*PRESSUREFACTOR/3. << "\n"
1574  << " GPNORMAL " << trace(groupPressure_normal)*PRESSUREFACTOR/3. << "\n"
1575  << " GPNBOND " << trace(groupPressure_nbond)*PRESSUREFACTOR/3. << "\n"
1576  << " GPSLOW " << trace(groupPressure_slow)*PRESSUREFACTOR/3. << "\n"
1577  << endi;
1578  }
1579 }
1580 
1581 //
1582 // Calculates all pressure tensors using virials
1583 //
1584 // Sets variables:
1585 // pressure, pressure_normal, pressure_nbond, pressure_slow
1586 // groupPressure, groupPressure_normal, groupPressure_nbond, groupPressure_slow
1587 // controlPressure, controlPressure_normal, controlPressure_nbond, controlPressure_slow
1588 // pressure_amd
1589 //
1590 void Controller::calcPressure(int step, int minimize,
1591  const Tensor& virial_normal_in, const Tensor& virial_nbond_in, const Tensor& virial_slow_in,
1592  const Tensor& intVirial_normal, const Tensor& intVirial_nbond, const Tensor& intVirial_slow,
1593  const Vector& extForce_normal, const Vector& extForce_nbond, const Vector& extForce_slow) {
1594 
1595  Tensor virial_normal = virial_normal_in;
1596  Tensor virial_nbond = virial_nbond_in;
1597  Tensor virial_slow = virial_slow_in;
1598 
1599  // Tensor tmp = virial_normal;
1600  // fprintf(stderr, "%1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf %1.2lf\n",
1601  // tmp.xx, tmp.xy, tmp.xz, tmp.yx, tmp.yy, tmp.yz, tmp.zx, tmp.zy, tmp.zz);
1602 
1603  Node *node = Node::Object();
1604  Molecule *molecule = node->molecule;
1605  SimParameters *simParameters = node->simParameters;
1606  Lattice &lattice = state->lattice;
1607 
1608  BigReal volume;
1609 
1610  Vector extPosition = lattice.origin();
1611  virial_normal -= outer(extForce_normal,extPosition);
1612  virial_nbond -= outer(extForce_nbond,extPosition);
1613  virial_slow -= outer(extForce_slow,extPosition);
1614 
1615  if ( (volume=lattice.volume()) != 0. )
1616  {
1617 
1618  if (simParameters->LJcorrection && volume) {
1619 #ifdef MEM_OPT_VERSION
1620  NAMD_bug("LJcorrection not supported in memory optimized build.");
1621 #else
1622  // Apply tail correction to pressure
1623  BigReal alchLambda = simParameters->getCurrentLambda(step);
1624  virial_normal += Tensor::identity(molecule->getVirialTailCorr(alchLambda) / volume);
1625 #endif
1626  }
1627 
1628  // kinetic energy component included in virials
1629  pressure_normal = virial_normal / volume;
1630  groupPressure_normal = ( virial_normal - intVirial_normal ) / volume;
1631 
1632  if (simParameters->accelMDOn) {
1633  pressure_amd = virial_amd / volume;
1636  }
1637 
1638  if ( minimize || ! ( step % nbondFreq ) )
1639  {
1640  pressure_nbond = virial_nbond / volume;
1641  groupPressure_nbond = ( virial_nbond - intVirial_nbond ) / volume;
1642  }
1643 
1644  if ( minimize || ! ( step % slowFreq ) )
1645  {
1646  pressure_slow = virial_slow / volume;
1647  groupPressure_slow = ( virial_slow - intVirial_slow ) / volume;
1648  }
1649 
1653  }
1654  else
1655  {
1656  pressure = Tensor();
1657  groupPressure = Tensor();
1658  }
1659 
1660  if ( simParameters->useGroupPressure )
1661  {
1666  }
1667  else
1668  {
1673  }
1674 
1675  if ( simParameters->useFlexibleCell ) {
1676  // use symmetric pressure to control rotation
1677  // controlPressure_normal = symmetric(controlPressure_normal);
1678  // controlPressure_nbond = symmetric(controlPressure_nbond);
1679  // controlPressure_slow = symmetric(controlPressure_slow);
1680  // controlPressure = symmetric(controlPressure);
1681  // only use on-diagonal components for now
1686  if ( simParameters->useConstantRatio ) {
1687 #define AVGXY(T) T.xy = T.yx = 0; T.xx = T.yy = 0.5 * ( T.xx + T.yy );\
1688  T.xz = T.zx = T.yz = T.zy = 0.5 * ( T.xz + T.yz );
1693 #undef AVGXY
1694  }
1695  } else {
1702  controlPressure =
1703  Tensor::identity(trace(controlPressure)/3.);
1704  }
1705 }
1706 
1708 (BigReal testV, int n,
1709  BigReal *Vmax, BigReal *Vmin, BigReal *Vavg, BigReal *M2, BigReal *sigmaV){
1710  BigReal delta;
1711 
1712  if(testV > *Vmax){
1713  *Vmax = testV;
1714  }
1715  else if(testV < *Vmin){
1716  *Vmin = testV;
1717  }
1718 
1719  //mean and std calculated by Online Algorithm
1720  delta = testV - *Vavg;
1721  *Vavg += delta / (BigReal)n;
1722  *M2 += delta * (testV - (*Vavg));
1723 
1724  *sigmaV = sqrt(*M2/n);
1725 }
1726 
1728 (int iE, int step, BigReal sigma0, BigReal Vmax, BigReal Vmin, BigReal Vavg, BigReal sigmaV,
1729  BigReal* k0, BigReal* k, BigReal* E, int *iEused, char *warn){
1730  switch(iE){
1731  case 2:
1732  *k0 = (1-sigma0/sigmaV) * (Vmax-Vmin) / (Vavg-Vmin);
1733  // if k0 <=0 OR k0 > 1, switch to iE=1 mode
1734  if(*k0 > 1.)
1735  *warn |= 1;
1736  else if(*k0 <= 0.)
1737  *warn |= 2;
1738  //else stay in iE=2 mode
1739  else{
1740  *E = Vmin + (Vmax-Vmin)/(*k0);
1741  *iEused = 2;
1742  break;
1743  }
1744  case 1:
1745  *E = Vmax;
1746  *k0 = (sigma0 / sigmaV) * (Vmax-Vmin) / (Vmax-Vavg);
1747  if(*k0 > 1.)
1748  *k0 = 1.;
1749  *iEused = 1;
1750  break;
1751  }
1752 
1753  *k = *k0 / (Vmax - Vmin);
1754 }
1755 
1757 (BigReal k, BigReal E, BigReal testV, Tensor vir_orig,
1758  BigReal *dV, BigReal *factor, Tensor *vir){
1759  BigReal VE = testV - E;
1760  //if V < E, apply boost
1761  if(VE < 0){
1762  *factor = k * VE;
1763  *vir += vir_orig * (*factor);
1764  *dV += (*factor) * VE/2.;
1765  *factor += 1.;
1766  }
1767  else{
1768  *factor = 1.;
1769  }
1770 }
1771 
1773 (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){
1774  FILE *rest;
1775  char timestepstr[20];
1776  char *filename;
1777  int baselen;
1778  char *restart_name;
1779  const char *bsuffix;
1780 
1781  if(lasttime || simParams->restartFrequency == 0){
1782  filename = simParams->outputFilename;
1783  bsuffix = ".BAK";
1784  }
1785  else{
1786  filename = simParams->restartFilename;
1787  bsuffix = ".old";
1788  }
1789 
1790  baselen = strlen(filename);
1791  restart_name = new char[baselen+26];
1792 
1793  strcpy(restart_name, filename);
1794  if ( simParams->restartSave && !lasttime) {
1795  sprintf(timestepstr,".%d",step_n);
1796  strcat(restart_name, timestepstr);
1797  bsuffix = ".BAK";
1798  }
1799  strcat(restart_name, ".gamd");
1800 
1801  if(write_topic){
1802  NAMD_backup_file(restart_name,bsuffix);
1803 
1804  rest = fopen(restart_name, "w");
1805  if(rest){
1806  fprintf(rest, "# NAMD accelMDG restart file\n"
1807  "# D/T Vn Vmax Vmin Vavg sigmaV M2 E k\n"
1808  "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
1809  type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
1810  fclose(rest);
1811  iout << "GAUSSIAN ACCELERATED MD RESTART DATA WRITTEN TO " << restart_name << "\n" << endi;
1812  }
1813  else
1814  iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
1815  }
1816  else{
1817  rest = fopen(restart_name, "a");
1818  if(rest){
1819  fprintf(rest, "%c %d %.17lg %.17lg %.17lg %.17lg %.17le %.17lg %.17lg\n",
1820  type, V_n-1, Vmax, Vmin, Vavg, sigmaV, M2, E, k);
1821  fclose(rest);
1822  }
1823  else
1824  iout << iWARN << "Cannot write accelMDG restart file\n" << endi;
1825  }
1826 
1827  delete[] restart_name;
1828 }
1829 
1830 
1831 void Controller::rescaleaccelMD(int step, int minimize)
1832 {
1833  if ( !simParams->accelMDOn ) return;
1834 
1836 
1837  // copy all to submit_reduction
1838  for ( int i=0; i<REDUCTION_MAX_RESERVED; ++i ) {
1840  }
1842 
1843  if (step == simParams->firstTimestep) {
1844  accelMDdVAverage = 0;
1845  accelMDdV = 0;
1846  }
1847 // if ( minimize || ((step < simParams->accelMDFirstStep ) || (step > simParams->accelMDLastStep ))) return;
1848  if ( minimize || (step < simParams->accelMDFirstStep ) || ( simParams->accelMDLastStep > 0 && step > simParams->accelMDLastStep )) return;
1849 
1850  Node *node = Node::Object();
1851  Molecule *molecule = node->molecule;
1852  Lattice &lattice = state->lattice;
1853 
1854  const BigReal accelMDE = simParams->accelMDE;
1855  const BigReal accelMDalpha = simParams->accelMDalpha;
1856  const BigReal accelMDTE = simParams->accelMDTE;
1857  const BigReal accelMDTalpha = simParams->accelMDTalpha;
1858  const int accelMDOutFreq = simParams->accelMDOutFreq;
1859 
1860  //GaMD
1861  static BigReal VmaxP, VminP, VavgP, M2P=0, sigmaVP;
1862  static BigReal VmaxD, VminD, VavgD, M2D=0, sigmaVD;
1863  static BigReal k0P, kP, EP;
1864  static BigReal k0D, kD, ED;
1865  static int V_n = 1, iEusedD, iEusedP;
1866  static char warnD, warnP;
1867  static int restfreq;
1868 
1871  const int ntcmd = simParams->accelMDFirstStep + simParams->accelMDGcMDSteps;
1872  const int ntebprep = ntcmd + simParams->accelMDGEquiPrepSteps;
1873  const int nteb = ntcmd + simParams->accelMDGEquiSteps;
1874  const int ntave = simParams->accelMDGStatWindow;
1875 
1876  BigReal volume;
1877  BigReal bondEnergy;
1878  BigReal angleEnergy;
1879  BigReal dihedralEnergy;
1880  BigReal improperEnergy;
1881  BigReal crosstermEnergy;
1882  BigReal boundaryEnergy;
1883  BigReal miscEnergy;
1884  BigReal amd_electEnergy;
1885  BigReal amd_ljEnergy;
1886  BigReal amd_ljEnergy_Corr = 0.;
1887  BigReal amd_electEnergySlow;
1888  BigReal amd_groLJEnergy;
1889  BigReal amd_groGaussEnergy;
1890  BigReal amd_goTotalEnergy;
1891  BigReal potentialEnergy;
1892  BigReal factor_dihe = 1;
1893  BigReal factor_tot = 1;
1894  BigReal testV;
1895  BigReal dV = 0.;
1896  Vector accelMDfactor;
1897  Tensor vir; //auto initialized to 0
1898  Tensor vir_dihe;
1899  Tensor vir_normal;
1900  Tensor vir_nbond;
1901  Tensor vir_slow;
1902 
1903  volume = lattice.volume();
1904 
1905  bondEnergy = amd_reduction->item(REDUCTION_BOND_ENERGY);
1906  angleEnergy = amd_reduction->item(REDUCTION_ANGLE_ENERGY);
1907  dihedralEnergy = amd_reduction->item(REDUCTION_DIHEDRAL_ENERGY);
1908  improperEnergy = amd_reduction->item(REDUCTION_IMPROPER_ENERGY);
1909  crosstermEnergy = amd_reduction->item(REDUCTION_CROSSTERM_ENERGY);
1910  boundaryEnergy = amd_reduction->item(REDUCTION_BC_ENERGY);
1911  miscEnergy = amd_reduction->item(REDUCTION_MISC_ENERGY);
1912 
1913  GET_TENSOR(vir_dihe,amd_reduction,REDUCTION_VIRIAL_AMD_DIHE);
1914  GET_TENSOR(vir_normal,amd_reduction,REDUCTION_VIRIAL_NORMAL);
1915  GET_TENSOR(vir_nbond,amd_reduction,REDUCTION_VIRIAL_NBOND);
1916  GET_TENSOR(vir_slow,amd_reduction,REDUCTION_VIRIAL_SLOW);
1917 
1918  if ( !( step % nbondFreq ) ) {
1919  amd_electEnergy = amd_reduction->item(REDUCTION_ELECT_ENERGY);
1920  amd_ljEnergy = amd_reduction->item(REDUCTION_LJ_ENERGY);
1921  amd_groLJEnergy = amd_reduction->item(REDUCTION_GRO_LJ_ENERGY);
1922  amd_groGaussEnergy = amd_reduction->item(REDUCTION_GRO_GAUSS_ENERGY);
1925  amd_goTotalEnergy = goNativeEnergy + goNonnativeEnergy;
1926  } else {
1927  amd_electEnergy = electEnergy;
1928  amd_ljEnergy = ljEnergy;
1929  amd_groLJEnergy = groLJEnergy;
1930  amd_groGaussEnergy = groGaussEnergy;
1931  amd_goTotalEnergy = goTotalEnergy;
1932  }
1933 
1934  if( simParams->LJcorrection && volume ) {
1935  // Obtain tail correction (copied from printEnergies())
1936  // This value is only printed for sanity check
1937  // accelMD currently does not 'boost' tail correction
1938  amd_ljEnergy_Corr = molecule->tail_corr_ener / volume;
1939  }
1940 
1941  if ( !( step % slowFreq ) ) {
1942  amd_electEnergySlow = amd_reduction->item(REDUCTION_ELECT_ENERGY_SLOW);
1943  } else {
1944  amd_electEnergySlow = electEnergySlow;
1945  }
1946 
1947  potentialEnergy = bondEnergy + angleEnergy + dihedralEnergy +
1948  improperEnergy + amd_electEnergy + amd_electEnergySlow + amd_ljEnergy +
1949  crosstermEnergy + boundaryEnergy + miscEnergy +
1950  amd_goTotalEnergy + amd_groLJEnergy + amd_groGaussEnergy;
1951 
1952  //GaMD
1953  if(simParams->accelMDG){
1954  // if first time running accelMDG module
1955  if(step == firststep) {
1956  iEusedD = iEusedP = simParams->accelMDGiE;
1957  warnD = warnP = '\0';
1958 
1959  //restart file reading
1961  FILE *rest = fopen(simParams->accelMDGRestartFile, "r");
1962  char line[256];
1963  int dihe_n=0, tot_n=0;
1964  if(!rest){
1965  sprintf(line, "Cannot open accelMDG restart file: %s", simParams->accelMDGRestartFile);
1966  NAMD_die(line);
1967  }
1968 
1969  while(fgets(line, 256, rest)){
1970  if(line[0] == 'D'){
1971  dihe_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
1972  &V_n, &VmaxD, &VminD, &VavgD, &sigmaVD, &M2D, &ED, &kD);
1973  }
1974  else if(line[0] == 'T'){
1975  tot_n = sscanf(line+1, " %d %la %la %la %la %la %la %la",
1976  &V_n, &VmaxP, &VminP, &VavgP, &sigmaVP, &M2P, &EP, &kP);
1977  }
1978  }
1979 
1980  fclose(rest);
1981 
1982  V_n++;
1983 
1984  //check dihe settings
1986  if(dihe_n !=8)
1987  NAMD_die("Invalid dihedral potential energy format in accelMDG restart file");
1988  k0D = kD * (VmaxD - VminD);
1989  iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: DIHED"
1990  << " Vmax " << VmaxD << " Vmin " << VminD
1991  << " Vavg " << VavgD << " sigmaV " << sigmaVD
1992  << " E " << ED << " k " << kD
1993  << "\n" << endi;
1994  }
1995 
1996  //check tot settings
1998  if(tot_n !=8)
1999  NAMD_die("Invalid total potential energy format in accelMDG restart file");
2000  k0P = kP * (VmaxP - VminP);
2001  iout << "GAUSSIAN ACCELERATED MD READ FROM RESTART FILE: TOTAL"
2002  << " Vmax " << VmaxP << " Vmin " << VminP
2003  << " Vavg " << VavgP << " sigmaV " << sigmaVP
2004  << " E " << EP << " k " << kP
2005  << "\n" << endi;
2006  }
2007 
2008  iEusedD = (ED == VmaxD) ? 1 : 2;
2009  iEusedP = (EP == VmaxP) ? 1 : 2;
2010  }
2011  //local restfreq follows NAMD restartfreq (default: 500)
2013  }
2014 
2015  //for dihedral potential
2017  testV = dihedralEnergy + crosstermEnergy;
2018 
2019  //write restart file every restartfreq steps
2020  if(step > firststep && step % restfreq == 0)
2021  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2022  true, false);
2023  //write restart file at the end of the simulation
2024  if( simParams->accelMDLastStep > 0 ){
2025  if( step == simParams->accelMDLastStep )
2026  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2027  true, true);
2028  }
2029  else if(step == simParams->N)
2030  write_accelMDG_rest_file(step, 'D', V_n, VmaxD, VminD, VavgD, sigmaVD, M2D, ED, kD,
2031  true, true);
2032 
2033  //conventional MD
2034  if(step < ntcmd){
2035  //very first step
2036  if(step == firststep && !simParams->accelMDGRestart){
2037  //initialize stat
2038  VmaxD = VminD = VavgD = testV;
2039  M2D = sigmaVD = 0.;
2040  }
2041  //first step after cmdprep
2042  else if(step == ntcmdprep && ntcmdprep != 0){
2043  //reset stat
2044  VmaxD = VminD = VavgD = testV;
2045  M2D = sigmaVD = 0.;
2046  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
2047  }
2048  //every ntave steps
2049  else if(ntave > 0 && step % ntave == 0){
2050  //update Vmax Vmin
2051  if(testV > VmaxD) VmaxD = testV;
2052  if(testV < VminD) VminD = testV;
2053  //reset avg and std
2054  VavgD = testV;
2055  M2D = sigmaVD = 0.;
2056  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
2057  }
2058  //normal steps
2059  else
2060  calc_accelMDG_mean_std(testV, V_n,
2061  &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2062 
2063  //last cmd step
2064  if(step == ntcmd - 1){
2066  VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2067  }
2068  }
2069  //equilibration
2070  else if(step < nteb){
2071  calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
2072  &dV, &factor_dihe, &vir);
2073 
2074  //first step after cmd
2075  if(step == ntcmd && simParams->accelMDGresetVaftercmd){
2076  //reset stat
2077  VmaxD = VminD = VavgD = testV;
2078  M2D = sigmaVD = 0.;
2079  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL STATISTICS\n" << endi;
2080  }
2081  //every ntave steps
2082  else if(ntave > 0 && step % ntave == 0){
2083  //update Vmax Vmin
2084  if(testV > VmaxD) VmaxD = testV;
2085  if(testV < VminD) VminD = testV;
2086  //reset avg and std
2087  VavgD = testV;
2088  M2D = sigmaVD = 0.;
2089  iout << "GAUSSIAN ACCELERATED MD: RESET DIHEDRAL POTENTIAL AVG AND STD\n" << endi;
2090  }
2091  else
2092  calc_accelMDG_mean_std(testV, V_n,
2093  &VmaxD, &VminD, &VavgD, &M2D, &sigmaVD) ;
2094 
2095  //steps after ebprep
2096  if(step >= ntebprep)
2097  if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2099  VmaxD, VminD, VavgD, sigmaVD, &k0D, &kD, &ED, &iEusedD, &warnD);
2100  }
2101  //production
2102  else{
2103  calc_accelMDG_force_factor(kD, ED, testV, vir_dihe,
2104  &dV, &factor_dihe, &vir);
2105  }
2106 
2107  }
2108  //for total potential
2110  Tensor vir_tot = vir_normal + vir_nbond + vir_slow;
2111  testV = potentialEnergy;
2112  if(simParams->accelMDdual){
2113  testV -= dihedralEnergy + crosstermEnergy;
2114  vir_tot -= vir_dihe;
2115  }
2116 
2117  //write restart file every restartfreq steps
2118  if(step > firststep && step % restfreq == 0)
2119  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2120  !simParams->accelMDdual, false);
2121  //write restart file at the end of simulation
2122  if( simParams->accelMDLastStep > 0 ){
2123  if( step == simParams->accelMDLastStep )
2124  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2125  !simParams->accelMDdual, true);
2126  }
2127  else if(step == simParams->N)
2128  write_accelMDG_rest_file(step, 'T', V_n, VmaxP, VminP, VavgP, sigmaVP, M2P, EP, kP,
2129  !simParams->accelMDdual, true);
2130 
2131  //conventional MD
2132  if(step < ntcmd){
2133  //very first step
2134  if(step == firststep && !simParams->accelMDGRestart){
2135  //initialize stat
2136  VmaxP = VminP = VavgP = testV;
2137  M2P = sigmaVP = 0.;
2138  }
2139  //first step after cmdprep
2140  else if(step == ntcmdprep && ntcmdprep != 0){
2141  //reset stat
2142  VmaxP = VminP = VavgP = testV;
2143  M2P = sigmaVP = 0.;
2144  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
2145  }
2146  //every ntave steps
2147  else if(ntave > 0 && step % ntave == 0){
2148  //update Vmax Vmin
2149  if(testV > VmaxP) VmaxP = testV;
2150  if(testV < VminP) VminP = testV;
2151  //reset avg and std
2152  VavgP = testV;
2153  M2P = sigmaVP = 0.;
2154  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
2155  }
2156  //normal steps
2157  else
2158  calc_accelMDG_mean_std(testV, V_n,
2159  &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2160  //last cmd step
2161  if(step == ntcmd - 1){
2163  VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2164  }
2165  }
2166  //equilibration
2167  else if(step < nteb){
2168  calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
2169  &dV, &factor_tot, &vir);
2170 
2171  //first step after cmd
2172  if(step == ntcmd && simParams->accelMDGresetVaftercmd){
2173  //reset stat
2174  VmaxP = VminP = VavgP = testV;
2175  M2P = sigmaVP = 0.;
2176  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL STATISTICS\n" << endi;
2177  }
2178  //every ntave steps
2179  else if(ntave > 0 && step % ntave == 0){
2180  //update Vmax Vmin
2181  if(testV > VmaxP) VmaxP = testV;
2182  if(testV < VminP) VminP = testV;
2183  //reset avg and std
2184  VavgP = testV;
2185  M2P = sigmaVP = 0.;
2186  iout << "GAUSSIAN ACCELERATED MD: RESET TOTAL POTENTIAL AVG AND STD\n" << endi;
2187  }
2188  else
2189  calc_accelMDG_mean_std(testV, V_n,
2190  &VmaxP, &VminP, &VavgP, &M2P, &sigmaVP);
2191 
2192  //steps after ebprep
2193  if(step >= ntebprep)
2194  if((ntave > 0 && step % ntave == 0) || ntave <= 0)
2196  VmaxP, VminP, VavgP, sigmaVP, &k0P, &kP, &EP, &iEusedP, &warnP);
2197  }
2198  //production
2199  else{
2200  calc_accelMDG_force_factor(kP, EP, testV, vir_tot,
2201  &dV, &factor_tot, &vir);
2202  }
2203 
2204  }
2205  accelMDdVAverage += dV;
2206 
2207  //first step after ntcmdprep
2208  if((ntcmdprep > 0 && step == ntcmdprep) ||
2209  (step < nteb && ntave > 0 && step % ntave == 0) ||
2210  (simParams->accelMDGresetVaftercmd && step == ntcmd)){
2211  V_n = 1;
2212  }
2213 
2214  if(step < nteb)
2215  V_n++;
2216 
2217  }
2218  //aMD
2219  else{
2220  if (simParams->accelMDdihe) {
2221 
2222  testV = dihedralEnergy + crosstermEnergy;
2223  if ( testV < accelMDE ) {
2224  factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2225  factor_dihe *= factor_dihe;
2226  vir = vir_dihe * (factor_dihe - 1.0);
2227  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2228  accelMDdVAverage += dV;
2229  }
2230 
2231  } else if (simParams->accelMDdual) {
2232 
2233  testV = dihedralEnergy + crosstermEnergy;
2234  if ( testV < accelMDE ) {
2235  factor_dihe = accelMDalpha/(accelMDalpha + accelMDE - testV);
2236  factor_dihe *= factor_dihe;
2237  vir = vir_dihe * (factor_dihe - 1.0) ;
2238  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2239  }
2240 
2241  testV = potentialEnergy - dihedralEnergy - crosstermEnergy;
2242  if ( testV < accelMDTE ) {
2243  factor_tot = accelMDTalpha/(accelMDTalpha + accelMDTE - testV);
2244  factor_tot *= factor_tot;
2245  vir += (vir_normal - vir_dihe + vir_nbond + vir_slow) * (factor_tot - 1.0);
2246  dV += (accelMDTE - testV)*(accelMDTE - testV)/(accelMDTalpha + accelMDTE - testV);
2247  }
2248  accelMDdVAverage += dV;
2249 
2250  } else {
2251 
2252  testV = potentialEnergy;
2253  if ( testV < accelMDE ) {
2254  factor_tot = accelMDalpha/(accelMDalpha + accelMDE - testV);
2255  factor_tot *= factor_tot;
2256  vir = (vir_normal + vir_nbond + vir_slow) * (factor_tot - 1.0);
2257  dV = (accelMDE - testV)*(accelMDE - testV)/(accelMDalpha + accelMDE - testV);
2258  accelMDdVAverage += dV;
2259  }
2260  }
2261  }
2262 
2263  accelMDfactor[0]=factor_dihe;
2264  accelMDfactor[1]=factor_tot;
2265  accelMDfactor[2]=1;
2266  broadcast->accelMDRescaleFactor.publish(step,accelMDfactor);
2267  virial_amd = vir;
2268  accelMDdV = dV;
2269 
2270  if ( factor_tot < 0.001 ) {
2271  iout << iWARN << "accelMD is using a very high boost potential, simulation may become unstable!"
2272  << "\n" << endi;
2273  }
2274 
2275  if ( ! (step % accelMDOutFreq) ) {
2276  if ( !(step == simParams->firstTimestep) ) {
2277  accelMDdVAverage = accelMDdVAverage/accelMDOutFreq;
2278  }
2279  iout << "ACCELERATED MD: STEP " << step
2280  << " dV " << dV
2281  << " dVAVG " << accelMDdVAverage
2282  << " BOND " << bondEnergy
2283  << " ANGLE " << angleEnergy
2284  << " DIHED " << dihedralEnergy+crosstermEnergy
2285  << " IMPRP " << improperEnergy
2286  << " ELECT " << amd_electEnergy+amd_electEnergySlow
2287  << " VDW " << amd_ljEnergy
2288  << " POTENTIAL " << potentialEnergy;
2289  if(amd_ljEnergy_Corr)
2290  iout << " LJcorr " << amd_ljEnergy_Corr;
2291  iout << "\n" << endi;
2292  if(simParams->accelMDG){
2294  iout << "GAUSSIAN ACCELERATED MD: DIHED iE " << iEusedD
2295  << " Vmax " << VmaxD << " Vmin " << VminD
2296  << " Vavg " << VavgD << " sigmaV " << sigmaVD
2297  << " E " << ED << " k0 " << k0D << " k " << kD
2298  << "\n" << endi;
2299  if(warnD & 1)
2300  iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 > 1, "
2301  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2302  if(warnD & 2)
2303  iout << "GAUSSIAN ACCELERATED MD: DIHED !!! WARNING: k0 <= 0, "
2304  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2306  iout << "GAUSSIAN ACCELERATED MD: TOTAL iE " << iEusedP
2307  << " Vmax " << VmaxP << " Vmin " << VminP
2308  << " Vavg " << VavgP << " sigmaV " << sigmaVP
2309  << " E " << EP << " k0 " << k0P << " k " << kP
2310  << "\n" << endi;
2311  if(warnP & 1)
2312  iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 > 1, "
2313  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2314  if(warnP & 2)
2315  iout << "GAUSSIAN ACCELERATED MD: TOTAL !!! WARNING: k0 <= 0, "
2316  << "SWITCHED TO iE = 1 MODE !!!\n" << endi;
2317  warnD = warnP = '\0';
2318  }
2319 
2320  accelMDdVAverage = 0;
2321 
2322  if (simParams->accelMDDebugOn) {
2323  Tensor p_normal;
2324  Tensor p_nbond;
2325  Tensor p_slow;
2326  Tensor p;
2327  if ( volume != 0. ) {
2328  p_normal = vir_normal/volume;
2329  p_nbond = vir_nbond/volume;
2330  p_slow = vir_slow/volume;
2331  p = vir/volume;
2332  }
2333  iout << " accelMD Scaling Factor: " << accelMDfactor << "\n"
2334  << " accelMD PNORMAL " << trace(p_normal)*PRESSUREFACTOR/3. << "\n"
2335  << " accelMD PNBOND " << trace(p_nbond)*PRESSUREFACTOR/3. << "\n"
2336  << " accelMD PSLOW " << trace(p_slow)*PRESSUREFACTOR/3. << "\n"
2337  << " accelMD PAMD " << trace(p)*PRESSUREFACTOR/3. << "\n"
2338  << endi;
2339  }
2340  }
2341 }
2342 
2344  if (!simParams->adaptTempOn) return;
2345  iout << iINFO << "INITIALISING ADAPTIVE TEMPERING\n" << endi;
2346  adaptTempDtMin = 0;
2347  adaptTempDtMax = 0;
2348  adaptTempAutoDt = false;
2349  if (simParams->adaptTempBins == 0) {
2350  iout << iINFO << "READING ADAPTIVE TEMPERING RESTART FILE\n" << endi;
2351  std::ifstream adaptTempRead(simParams->adaptTempInFile);
2352  if (adaptTempRead) {
2353  int readInt;
2354  BigReal readReal;
2355  bool readBool;
2356  // step
2357  adaptTempRead >> readInt;
2358  // Start with min and max temperatures
2359  adaptTempRead >> adaptTempT; // KELVIN
2360  adaptTempRead >> adaptTempBetaMin; // KELVIN
2361  adaptTempRead >> adaptTempBetaMax; // KELVIN
2362  adaptTempBetaMin = 1./adaptTempBetaMin; // KELVIN^-1
2363  adaptTempBetaMax = 1./adaptTempBetaMax; // KELVIN^-1
2364  // In case file is manually edited
2365  if (adaptTempBetaMin > adaptTempBetaMax){
2366  readReal = adaptTempBetaMax;
2367  adaptTempBetaMax = adaptTempBetaMin;
2368  adaptTempBetaMin = adaptTempBetaMax;
2369  }
2370  adaptTempRead >> adaptTempBins;
2371  adaptTempRead >> adaptTempCg;
2372  adaptTempRead >> adaptTempDt;
2380  adaptTempDBeta = (adaptTempBetaMax - adaptTempBetaMin)/(adaptTempBins);
2381  for(int j = 0; j < adaptTempBins; ++j) {
2382  adaptTempRead >> adaptTempPotEnergyAve[j];
2383  adaptTempRead >> adaptTempPotEnergyVar[j];
2384  adaptTempRead >> adaptTempPotEnergySamples[j];
2385  adaptTempRead >> adaptTempPotEnergyAveNum[j];
2386  adaptTempRead >> adaptTempPotEnergyVarNum[j];
2387  adaptTempRead >> adaptTempPotEnergyAveDen[j];
2388  adaptTempBetaN[j] = adaptTempBetaMin + j * adaptTempDBeta;
2389  }
2390  adaptTempRead.close();
2391  }
2392  else NAMD_die("Could not open ADAPTIVE TEMPERING restart file.\n");
2393  }
2394  else {
2409  if (simParams->langevinOn)
2411  else if (simParams->rescaleFreq > 0)
2413  for(int j = 0; j < adaptTempBins; ++j){
2414  adaptTempPotEnergyAveNum[j] = 0.;
2415  adaptTempPotEnergyAveDen[j] = 0.;
2417  adaptTempPotEnergyVarNum[j] = 0.;
2418  adaptTempPotEnergyVar[j] = 0.;
2419  adaptTempPotEnergyAve[j] = 0.;
2421  }
2422  }
2423  if (simParams->adaptTempAutoDt > 0.0) {
2424  adaptTempAutoDt = true;
2427  }
2428  adaptTempDTave = 0;
2429  adaptTempDTavenum = 0;
2430  iout << iINFO << "ADAPTIVE TEMPERING: TEMPERATURE RANGE: [" << 1./adaptTempBetaMax << "," << 1./adaptTempBetaMin << "] KELVIN\n";
2431  iout << iINFO << "ADAPTIVE TEMPERING: NUMBER OF BINS TO STORE POT. ENERGY: " << adaptTempBins << "\n";
2432  iout << iINFO << "ADAPTIVE TEMPERING: ADAPTIVE BIN AVERAGING PARAMETER: " << adaptTempCg << "\n";
2433  iout << iINFO << "ADAPTIVE TEMPERING: SCALING CONSTANT FOR LANGEVIN TEMPERATURE UPDATES: " << adaptTempDt << "\n";
2434  iout << iINFO << "ADAPTIVE TEMPERING: INITIAL TEMPERATURE: " << adaptTempT<< "\n";
2435  if (simParams->adaptTempRestartFreq > 0) {
2439  NAMD_die("Error opening restart file");
2440  }
2441 }
2442 
2445  adaptTempRestartFile.seekp(std::ios::beg);
2446  iout << "ADAPTEMP: WRITING RESTART FILE AT STEP " << step << "\n" << endi;
2447  adaptTempRestartFile << step << " ";
2448  // Start with min and max temperatures
2449  adaptTempRestartFile << adaptTempT<< " "; // KELVIN
2450  adaptTempRestartFile << 1./adaptTempBetaMin << " "; // KELVIN
2451  adaptTempRestartFile << 1./adaptTempBetaMax << " "; // KELVIN
2455  adaptTempRestartFile << "\n" ;
2456  for(int j = 0; j < adaptTempBins; ++j) {
2463  adaptTempRestartFile << "\n";
2464  }
2466  }
2467 }
2468 
2469 void Controller::adaptTempUpdate(int step, int minimize)
2470 {
2471  //Beta = 1./T
2472  if ( !simParams->adaptTempOn ) return;
2473  int j = 0;
2474  if (step == simParams->firstTimestep) {
2475  adaptTempInit(step);
2476  return;
2477  }
2478  if ( minimize || (step < simParams->adaptTempFirstStep ) ||
2479  ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) return;
2480  const int adaptTempOutFreq = simParams->adaptTempOutFreq;
2481  const bool adaptTempDebug = simParams->adaptTempDebug;
2482  //Calculate Current inverse temperature and bin
2483  BigReal adaptTempBeta = 1./adaptTempT;
2484  adaptTempBin = (int)floor((adaptTempBeta - adaptTempBetaMin)/adaptTempDBeta);
2485 
2486  if (adaptTempBin < 0 || adaptTempBin > adaptTempBins)
2487  iout << iWARN << " adaptTempBin out of range: adaptTempBin: " << adaptTempBin
2488  << " adaptTempBeta: " << adaptTempBeta
2489  << " adaptTempDBeta: " << adaptTempDBeta
2490  << " betaMin:" << adaptTempBetaMin
2491  << " betaMax: " << adaptTempBetaMax << "\n";
2494 
2495  BigReal potentialEnergy;
2496  BigReal potEnergyAverage;
2497  BigReal potEnergyVariance;
2498  potentialEnergy = totalEnergy-kineticEnergy;
2499 
2500  //calculate new bin average and variance using adaptive averaging
2503  adaptTempPotEnergyVarNum[adaptTempBin] = adaptTempPotEnergyVarNum[adaptTempBin]*gammaAve + potentialEnergy*potentialEnergy;
2504 
2507  potEnergyVariance -= potEnergyAverage*potEnergyAverage;
2508 
2509  adaptTempPotEnergyAve[adaptTempBin] = potEnergyAverage;
2510  adaptTempPotEnergyVar[adaptTempBin] = potEnergyVariance;
2511 
2512  // Weighted integral of <Delta E^2>_beta dbeta <= Eq 4 of JCP 132 244101
2513  // Integrals of Eqs 5 and 6 is done as piecewise assuming <Delta E^2>_beta
2514  // is constant for each bin. This is to estimate <E(beta)> where beta \in
2515  // (beta_i,beta_{i+1}) using Eq 2 of JCP 132 244101
2516  if ( ! ( step % simParams->adaptTempFreq ) ) {
2517  // If adaptTempBin not at the edge, calculate improved average:
2518  if (adaptTempBin > 0 && adaptTempBin < adaptTempBins-1) {
2519  // Get Averaging Limits:
2520  BigReal deltaBeta = 0.04*adaptTempBeta; //0.08 used in paper - make variable
2521  BigReal betaPlus;
2522  BigReal betaMinus;
2523  int nMinus =0;
2524  int nPlus = 0;
2525  if ( adaptTempBeta-adaptTempBetaMin < deltaBeta ) deltaBeta = adaptTempBeta-adaptTempBetaMin;
2526  if ( adaptTempBetaMax-adaptTempBeta < deltaBeta ) deltaBeta = adaptTempBetaMax-adaptTempBeta;
2527  betaMinus = adaptTempBeta - deltaBeta;
2528  betaPlus = adaptTempBeta + deltaBeta;
2529  BigReal betaMinus2 = betaMinus*betaMinus;
2530  BigReal betaPlus2 = betaPlus*betaPlus;
2531  nMinus = (int)floor((betaMinus-adaptTempBetaMin)/adaptTempDBeta);
2532  nPlus = (int)floor((betaPlus-adaptTempBetaMin)/adaptTempDBeta);
2533  // Variables for <E(beta)> estimate:
2534  BigReal potEnergyAve0 = 0.0;
2535  BigReal potEnergyAve1 = 0.0;
2536  // Integral terms
2537  BigReal A0 = 0;
2538  BigReal A1 = 0;
2539  BigReal A2 = 0;
2540  //A0 phi_s integral for beta_minus < beta < beta_{i+1}
2541  BigReal betaNp1 = adaptTempBetaN[adaptTempBin+1];
2542  BigReal DeltaE2Ave;
2543  BigReal DeltaE2AveSum = 0;
2544  BigReal betaj;
2545  for (j = nMinus+1; j <= adaptTempBin; ++j) {
2546  potEnergyAve0 += adaptTempPotEnergyAve[j];
2547  DeltaE2Ave = adaptTempPotEnergyVar[j];
2548  DeltaE2AveSum += DeltaE2Ave;
2549  A0 += j*DeltaE2Ave;
2550  }
2551  A0 *= adaptTempDBeta;
2552  A0 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaMinus)*DeltaE2AveSum;
2553  A0 *= adaptTempDBeta;
2554  betaj = adaptTempBetaN[nMinus+1]-betaMinus;
2555  betaj *= betaj;
2556  A0 += 0.5*betaj*adaptTempPotEnergyVar[nMinus];
2557  A0 /= (betaNp1-betaMinus);
2558 
2559  //A1 phi_s integral for beta_{i+1} < beta < beta_plus
2560  DeltaE2AveSum = 0;
2561  for (j = adaptTempBin+1; j < nPlus; j++) {
2562  potEnergyAve1 += adaptTempPotEnergyAve[j];
2563  DeltaE2Ave = adaptTempPotEnergyVar[j];
2564  DeltaE2AveSum += DeltaE2Ave;
2565  A1 += j*DeltaE2Ave;
2566  }
2567  A1 *= adaptTempDBeta;
2568  A1 += (adaptTempBetaMin+0.5*adaptTempDBeta-betaPlus)*DeltaE2AveSum;
2569  A1 *= adaptTempDBeta;
2570  if ( nPlus < adaptTempBins ) {
2571  betaj = betaPlus-adaptTempBetaN[nPlus];
2572  betaj *= betaj;
2573  A1 -= 0.5*betaj*adaptTempPotEnergyVar[nPlus];
2574  }
2575  A1 /= (betaPlus-betaNp1);
2576 
2577  //A2 phi_t integral for beta_i
2578  A2 = 0.5*adaptTempDBeta*potEnergyVariance;
2579 
2580  // Now calculate a+ and a-
2581  DeltaE2Ave = A0-A1;
2582  BigReal aplus = 0;
2583  BigReal aminus = 0;
2584  if (DeltaE2Ave != 0) {
2585  aplus = (A0-A2)/(A0-A1);
2586  if (aplus < 0) {
2587  aplus = 0;
2588  }
2589  if (aplus > 1) {
2590  aplus = 1;
2591  }
2592  aminus = 1-aplus;
2593  potEnergyAve0 *= adaptTempDBeta;
2594  potEnergyAve0 += adaptTempPotEnergyAve[nMinus]*(adaptTempBetaN[nMinus+1]-betaMinus);
2595  potEnergyAve0 /= (betaNp1-betaMinus);
2596  potEnergyAve1 *= adaptTempDBeta;
2597  if ( nPlus < adaptTempBins ) {
2598  potEnergyAve1 += adaptTempPotEnergyAve[nPlus]*(betaPlus-adaptTempBetaN[nPlus]);
2599  }
2600  potEnergyAve1 /= (betaPlus-betaNp1);
2601  potEnergyAverage = aminus*potEnergyAve0;
2602  potEnergyAverage += aplus*potEnergyAve1;
2603  }
2604  if (simParams->adaptTempDebug) {
2605  iout << "ADAPTEMP DEBUG:" << "\n"
2606  << " adaptTempBin: " << adaptTempBin << "\n"
2607  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
2608  << " adaptTempBeta: " << adaptTempBeta << "\n"
2609  << " adaptTemp: " << adaptTempT<< "\n"
2610  << " betaMin: " << adaptTempBetaMin << "\n"
2611  << " betaMax: " << adaptTempBetaMax << "\n"
2612  << " gammaAve: " << gammaAve << "\n"
2613  << " deltaBeta: " << deltaBeta << "\n"
2614  << " betaMinus: " << betaMinus << "\n"
2615  << " betaPlus: " << betaPlus << "\n"
2616  << " nMinus: " << nMinus << "\n"
2617  << " nPlus: " << nPlus << "\n"
2618  << " A0: " << A0 << "\n"
2619  << " A1: " << A1 << "\n"
2620  << " A2: " << A2 << "\n"
2621  << " a+: " << aplus << "\n"
2622  << " a-: " << aminus << "\n"
2623  << endi;
2624  }
2625  }
2626  else {
2627  if (simParams->adaptTempDebug) {
2628  iout << "ADAPTEMP DEBUG:" << "\n"
2629  << " adaptTempBin: " << adaptTempBin << "\n"
2630  << " Samples: " << adaptTempPotEnergySamples[adaptTempBin] << "\n"
2631  << " adaptTempBeta: " << adaptTempBeta << "\n"
2632  << " adaptTemp: " << adaptTempT<< "\n"
2633  << " betaMin: " << adaptTempBetaMin << "\n"
2634  << " betaMax: " << adaptTempBetaMax << "\n"
2635  << " gammaAve: " << gammaAve << "\n"
2636  << " deltaBeta: N/A\n"
2637  << " betaMinus: N/A\n"
2638  << " betaPlus: N/A\n"
2639  << " nMinus: N/A\n"
2640  << " nPlus: N/A\n"
2641  << " A0: N/A\n"
2642  << " A1: N/A\n"
2643  << " A2: N/A\n"
2644  << " a+: N/A\n"
2645  << " a-: N/A\n"
2646  << endi;
2647  }
2648  }
2649 
2650  //dT is new temperature
2651  BigReal dT = ((potentialEnergy-potEnergyAverage)/BOLTZMANN+adaptTempT)*adaptTempDt;
2652  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
2653  dT += adaptTempT;
2654 
2655  // Check if dT in [adaptTempTmin,adaptTempTmax]. If not try simpler estimate of mean
2656  // This helps sampling with poor statistics in the bins surrounding adaptTempBin.
2657  if ( dT > 1./adaptTempBetaMin || dT < 1./adaptTempBetaMax ) {
2659  dT += random->gaussian()*sqrt(2.*adaptTempDt)*adaptTempT;
2660  dT += adaptTempT;
2661  // Check again, if not then keep original adaptTempTor assign random.
2662  if ( dT > 1./adaptTempBetaMin ) {
2663  if (!simParams->adaptTempRandom) {
2664  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
2665  // << " K higher than adaptTempTmax."
2666  // << " Keeping temperature at "
2667  // << adaptTempT<< "\n"<< endi;
2668  dT = adaptTempT;
2669  }
2670  else {
2671  //iout << iWARN << "ADAPTEMP: " << step << " T= " << dT
2672  // << " K higher than adaptTempTmax."
2673  // << " Assigning random temperature in range\n" << endi;
2674  dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
2675  dT = 1./dT;
2676  }
2677  }
2678  else if ( dT < 1./adaptTempBetaMax ) {
2679  if (!simParams->adaptTempRandom) {
2680  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
2681  // << " K lower than adaptTempTmin."
2682  // << " Keeping temperature at " << adaptTempT<< "\n" << endi;
2683  dT = adaptTempT;
2684  }
2685  else {
2686  //iout << iWARN << "ADAPTEMP: " << step << " T= "<< dT
2687  // << " K lower than adaptTempTmin."
2688  // << " Assigning random temperature in range\n" << endi;
2689  dT = adaptTempBetaMin + random->uniform()*(adaptTempBetaMax-adaptTempBetaMin);
2690  dT = 1./dT;
2691  }
2692  }
2693  else if (adaptTempAutoDt) {
2694  //update temperature step size counter
2695  //FOR "TRUE" ADAPTIVE TEMPERING
2696  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
2697  if (adaptTempTdiff > 0) {
2698  adaptTempDTave += adaptTempTdiff;
2700 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
2701  }
2702  if(adaptTempDTavenum == 100){
2703  BigReal Frac;
2705  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
2706  Frac = adaptTempDTave/Frac;
2707  //if average temperature jump is > 3% of temperature range,
2708  //modify jump size to match 3%
2709  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
2710  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
2711  Frac = adaptTempDtMax/Frac;
2712  iout << "ADAPTEMP: Updating adaptTempDt to ";
2713  adaptTempDt *= Frac;
2714  iout << adaptTempDt << "\n" << endi;
2715  }
2716  adaptTempDTave = 0;
2717  adaptTempDTavenum = 0;
2718  }
2719  }
2720  }
2721  else if (adaptTempAutoDt) {
2722  //update temperature step size counter
2723  // FOR "TRUE" ADAPTIVE TEMPERING
2724  BigReal adaptTempTdiff = fabs(dT-adaptTempT);
2725  if (adaptTempTdiff > 0) {
2726  adaptTempDTave += adaptTempTdiff;
2728 // iout << "ADAPTEMP: adapTempTdiff = " << adaptTempTdiff << "\n";
2729  }
2730  if(adaptTempDTavenum == 100){
2731  BigReal Frac;
2733  Frac = 1./adaptTempBetaMin-1./adaptTempBetaMax;
2734  Frac = adaptTempDTave/Frac;
2735  //if average temperature jump is > 3% of temperature range,
2736  //modify jump size to match 3%
2737  iout << "ADAPTEMP: " << step << " FRAC " << Frac << "\n";
2738  if (Frac > adaptTempDtMax || Frac < adaptTempDtMin) {
2739  Frac = adaptTempDtMax/Frac;
2740  iout << "ADAPTEMP: Updating adaptTempDt to ";
2741  adaptTempDt *= Frac;
2742  iout << adaptTempDt << "\n" << endi;
2743  }
2744  adaptTempDTave = 0;
2745  adaptTempDTavenum = 0;
2746 
2747  }
2748 
2749  }
2750 
2751  adaptTempT = dT;
2753  }
2754  adaptTempWriteRestart(step);
2755  if ( ! (step % adaptTempOutFreq) ) {
2756  iout << "ADAPTEMP: STEP " << step
2757  << " TEMP " << adaptTempT
2758  << " ENERGY " << std::setprecision(10) << potentialEnergy
2759  << " ENERGYAVG " << std::setprecision(10) << potEnergyAverage
2760  << " ENERGYVAR " << std::setprecision(10) << potEnergyVariance;
2761  iout << "\n" << endi;
2762  }
2763 
2764 }
2765 
2766 
2767 void Controller::compareChecksums(int step, int forgiving) {
2768  Node *node = Node::Object();
2769  Molecule *molecule = node->molecule;
2770 
2771  // Some basic correctness checking
2772  BigReal checksum, checksum_b;
2773  char errmsg[256];
2774 
2775  checksum = reduction->item(REDUCTION_ATOM_CHECKSUM);
2776  if ( ((int)checksum) != molecule->numAtoms ) {
2777  sprintf(errmsg, "Bad global atom count! (%d vs %d)\n",
2778  (int)checksum, molecule->numAtoms);
2779  NAMD_bug(errmsg);
2780  }
2781 
2783  if ( ((int)checksum) && ((int)checksum) != computeChecksum ) {
2784  if ( ! computeChecksum ) {
2785  computesPartitioned = 0;
2786  } else if ( (int)checksum > computeChecksum && ! computesPartitioned ) {
2787  computesPartitioned = 1;
2788  } else {
2789  NAMD_bug("Bad global compute count!\n");
2790  }
2791  computeChecksum = ((int)checksum);
2792  }
2793 
2794  checksum_b = checksum = reduction->item(REDUCTION_BOND_CHECKSUM);
2795  if ( checksum_b && (((int)checksum) != molecule->numCalcBonds) ) {
2796  sprintf(errmsg, "Bad global bond count! (%d vs %d)\n",
2797  (int)checksum, molecule->numCalcBonds);
2798  if ( forgiving && (((int)checksum) < molecule->numCalcBonds) )
2799  iout << iWARN << errmsg << endi;
2800  else NAMD_bug(errmsg);
2801  }
2802 
2804  if ( checksum_b && (((int)checksum) != molecule->numCalcAngles) ) {
2805  sprintf(errmsg, "Bad global angle count! (%d vs %d)\n",
2806  (int)checksum, molecule->numCalcAngles);
2807  if ( forgiving && (((int)checksum) < molecule->numCalcAngles) )
2808  iout << iWARN << errmsg << endi;
2809  else NAMD_bug(errmsg);
2810  }
2811 
2813  if ( checksum_b && (((int)checksum) != molecule->numCalcDihedrals) ) {
2814  sprintf(errmsg, "Bad global dihedral count! (%d vs %d)\n",
2815  (int)checksum, molecule->numCalcDihedrals);
2816  if ( forgiving && (((int)checksum) < molecule->numCalcDihedrals) )
2817  iout << iWARN << errmsg << endi;
2818  else NAMD_bug(errmsg);
2819  }
2820 
2822  if ( checksum_b && (((int)checksum) != molecule->numCalcImpropers) ) {
2823  sprintf(errmsg, "Bad global improper count! (%d vs %d)\n",
2824  (int)checksum, molecule->numCalcImpropers);
2825  if ( forgiving && (((int)checksum) < molecule->numCalcImpropers) )
2826  iout << iWARN << errmsg << endi;
2827  else NAMD_bug(errmsg);
2828  }
2829 
2831  if ( checksum_b && (((int)checksum) != molecule->numCalcTholes) ) {
2832  sprintf(errmsg, "Bad global Thole count! (%d vs %d)\n",
2833  (int)checksum, molecule->numCalcTholes);
2834  if ( forgiving && (((int)checksum) < molecule->numCalcTholes) )
2835  iout << iWARN << errmsg << endi;
2836  else NAMD_bug(errmsg);
2837  }
2838 
2840  if ( checksum_b && (((int)checksum) != molecule->numCalcAnisos) ) {
2841  sprintf(errmsg, "Bad global Aniso count! (%d vs %d)\n",
2842  (int)checksum, molecule->numCalcAnisos);
2843  if ( forgiving && (((int)checksum) < molecule->numCalcAnisos) )
2844  iout << iWARN << errmsg << endi;
2845  else NAMD_bug(errmsg);
2846  }
2847 
2849  if ( checksum_b && (((int)checksum) != molecule->numCalcCrossterms) ) {
2850  sprintf(errmsg, "Bad global crossterm count! (%d vs %d)\n",
2851  (int)checksum, molecule->numCalcCrossterms);
2852  if ( forgiving && (((int)checksum) < molecule->numCalcCrossterms) )
2853  iout << iWARN << errmsg << endi;
2854  else NAMD_bug(errmsg);
2855  }
2856 
2858  if ( ((int64)checksum) > molecule->numCalcExclusions &&
2859  ( ! simParams->mollyOn || step % slowFreq ) ) {
2860  if ( forgiving )
2861  iout << iWARN << "High global exclusion count ("
2862  << ((int64)checksum) << " vs "
2863  << molecule->numCalcExclusions << "), possible error!\n"
2864  << iWARN << "This warning is not unusual during minimization.\n"
2865  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
2866  else {
2867  char errmsg[256];
2868  sprintf(errmsg, "High global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2869  (long long)checksum, (long long)molecule->numCalcExclusions);
2870  NAMD_bug(errmsg);
2871  }
2872  }
2873 
2874  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcExclusions &&
2876  if ( forgiving || ! simParams->fullElectFrequency ) {
2877  iout << iWARN << "Low global exclusion count! ("
2878  << ((int64)checksum) << " vs " << molecule->numCalcExclusions << ")";
2879  if ( forgiving ) {
2880  iout << "\n"
2881  << iWARN << "This warning is not unusual during minimization.\n"
2882  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
2883  } else {
2884  iout << " System unstable or pairlistdist or cutoff too small.\n";
2885  }
2886  iout << endi;
2887  } else {
2888  char errmsg[256];
2889  sprintf(errmsg, "Low global exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2890  (long long)checksum, (long long)molecule->numCalcExclusions);
2891  NAMD_bug(errmsg);
2892  }
2893  }
2894 
2895 #ifdef NAMD_CUDA
2897  if ( ((int64)checksum) > molecule->numCalcFullExclusions &&
2898  ( ! simParams->mollyOn || step % slowFreq ) ) {
2899  if ( forgiving )
2900  iout << iWARN << "High global CUDA exclusion count ("
2901  << ((int64)checksum) << " vs "
2902  << molecule->numCalcFullExclusions << "), possible error!\n"
2903  << iWARN << "This warning is not unusual during minimization.\n"
2904  << iWARN << "Decreasing pairlistdist or cutoff that is too close to periodic cell size may avoid this.\n" << endi;
2905  else {
2906  char errmsg[256];
2907  sprintf(errmsg, "High global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too close to periodic cell size.\n",
2908  (long long)checksum, (long long)molecule->numCalcFullExclusions);
2909  NAMD_bug(errmsg);
2910  }
2911  }
2912 
2913  if ( ((int64)checksum) && ((int64)checksum) < molecule->numCalcFullExclusions &&
2915  if ( forgiving || ! simParams->fullElectFrequency ) {
2916  iout << iWARN << "Low global CUDA exclusion count! ("
2917  << ((int64)checksum) << " vs " << molecule->numCalcFullExclusions << ")";
2918  if ( forgiving ) {
2919  iout << "\n"
2920  << iWARN << "This warning is not unusual during minimization.\n"
2921  << iWARN << "Increasing pairlistdist or cutoff may avoid this.\n";
2922  } else {
2923  iout << " System unstable or pairlistdist or cutoff too small.\n";
2924  }
2925  iout << endi;
2926  } else {
2927  char errmsg[256];
2928  sprintf(errmsg, "Low global CUDA exclusion count! (%lld vs %lld) System unstable or pairlistdist or cutoff too small.\n",
2929  (long long)checksum, (long long)molecule->numCalcFullExclusions);
2930  NAMD_bug(errmsg);
2931  }
2932  }
2933 #endif
2934 
2936  if ( ((int)checksum) && ! marginViolations ) {
2937  iout << iERROR << "Margin is too small for " << ((int)checksum) <<
2938  " atoms during timestep " << step << ".\n" << iERROR <<
2939  "Incorrect nonbonded forces and energies may be calculated!\n" << endi;
2940  }
2941  marginViolations += (int)checksum;
2942 
2944  if ( simParams->outputPairlists && ((int)checksum) && ! pairlistWarnings ) {
2945  iout << iINFO <<
2946  "Pairlistdist is too small for " << ((int)checksum) <<
2947  " computes during timestep " << step << ".\n" << endi;
2948  }
2949  if ( simParams->outputPairlists ) pairlistWarnings += (int)checksum;
2950 
2952  if ( checksum ) {
2953  if ( forgiving )
2954  iout << iWARN << "Stray PME grid charges ignored!\n" << endi;
2955  else NAMD_bug("Stray PME grid charges detected!\n");
2956  }
2957 }
2958 
2959 void Controller::printTiming(int step) {
2960 
2961  if ( simParams->outputTiming && ! ( step % simParams->outputTiming ) )
2962  {
2963  const double endWTime = CmiWallTimer() - firstWTime;
2964  const double endCTime = CmiTimer() - firstCTime;
2965 
2966  // fflush about once per minute
2967  if ( (((int)endWTime)>>6) != (((int)startWTime)>>6 ) ) fflush_count = 2;
2968 
2969  const double elapsedW =
2970  (endWTime - startWTime) / simParams->outputTiming;
2971  const double elapsedC =
2972  (endCTime - startCTime) / simParams->outputTiming;
2973 
2974  const double remainingW = elapsedW * (simParams->N - step);
2975  const double remainingW_hours = remainingW / 3600;
2976 
2977  startWTime = endWTime;
2978  startCTime = endCTime;
2979 
2980 #ifdef NAMD_CUDA
2981  if ( simParams->outputEnergies < 60 &&
2982  step < (simParams->firstTimestep + 10 * simParams->outputTiming) ) {
2983  iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n" << endi;
2984  }
2985 #endif
2986  if ( step >= (simParams->firstTimestep + simParams->outputTiming) ) {
2987  CmiPrintf("TIMING: %d CPU: %g, %g/step Wall: %g, %g/step"
2988  ", %g hours remaining, %f MB of memory in use.\n",
2989  step, endCTime, elapsedC, endWTime, elapsedW,
2990  remainingW_hours, memusage_MB());
2991  if ( fflush_count ) { --fflush_count; fflush(stdout); }
2992  }
2993  }
2994 }
2995 
2997 
2998  rescaleaccelMD(step,1);
2999  receivePressure(step,1);
3000 
3001  Node *node = Node::Object();
3002  Molecule *molecule = node->molecule;
3003  compareChecksums(step,1);
3004 
3005  printEnergies(step,1);
3006 
3012 
3013 }
3014 
3016 
3017  Node *node = Node::Object();
3018  Molecule *molecule = node->molecule;
3019  SimParameters *simParameters = node->simParameters;
3020  Lattice &lattice = state->lattice;
3021 
3022  compareChecksums(step);
3023 
3024  printEnergies(step,0);
3025 }
3026 
3027 void Controller::printEnergies(int step, int minimize)
3028 {
3029  Node *node = Node::Object();
3030  Molecule *molecule = node->molecule;
3031  SimParameters *simParameters = node->simParameters;
3032  Lattice &lattice = state->lattice;
3033 
3034  // Drude model ANISO energy is added into BOND energy
3035  // and THOLE energy is added into ELECT energy
3036 
3037  BigReal bondEnergy;
3038  BigReal angleEnergy;
3039  BigReal dihedralEnergy;
3040  BigReal improperEnergy;
3041  BigReal crosstermEnergy;
3042  BigReal boundaryEnergy;
3043  BigReal miscEnergy;
3044  BigReal potentialEnergy;
3045  BigReal flatEnergy;
3046  BigReal smoothEnergy;
3047  BigReal work;
3048 
3049  Vector momentum;
3050  Vector angularMomentum;
3051  BigReal volume = lattice.volume();
3052 
3053  bondEnergy = reduction->item(REDUCTION_BOND_ENERGY);
3054  angleEnergy = reduction->item(REDUCTION_ANGLE_ENERGY);
3055  dihedralEnergy = reduction->item(REDUCTION_DIHEDRAL_ENERGY);
3056  improperEnergy = reduction->item(REDUCTION_IMPROPER_ENERGY);
3057  crosstermEnergy = reduction->item(REDUCTION_CROSSTERM_ENERGY);
3058  boundaryEnergy = reduction->item(REDUCTION_BC_ENERGY);
3059  miscEnergy = reduction->item(REDUCTION_MISC_ENERGY);
3060 
3061  if ( minimize || ! ( step % nbondFreq ) )
3062  {
3065 
3066  // JLai
3069 
3073 
3074 //fepb
3078 
3085 //fepe
3086  }
3087 
3088  if ( minimize || ! ( step % slowFreq ) )
3089  {
3091 //fepb
3093 
3098 //fepe
3099  }
3100 
3101  if (simParameters->LJcorrection && volume) {
3102 #ifdef MEM_OPT_VERSION
3103  NAMD_bug("LJcorrection not supported in memory optimized build.");
3104 #else
3105  // Apply tail correction to energy.
3106  BigReal alchLambda = simParameters->getCurrentLambda(step);
3107  ljEnergy += molecule->getEnergyTailCorr(alchLambda, 0) / volume;
3108 
3109  if (simParameters->alchOn) {
3110  if (simParameters->alchFepOn) {
3111  BigReal alchLambda2 = simParameters->getCurrentLambda2(step);
3112  ljEnergy_f += molecule->getEnergyTailCorr(alchLambda2, 0) / volume;
3113  } else if (simParameters->alchThermIntOn) {
3114  ljEnergy_ti_1 += molecule->getEnergyTailCorr(1.0, 1) / volume;
3115  ljEnergy_ti_2 += molecule->getEnergyTailCorr(0.0, 1) / volume;
3116  }
3117  }
3118 #endif
3119  }
3120 
3121 //fepb BKR - Compute alchemical work if using dynamic lambda. This is here so
3122 // that the cumulative work can be given during a callback.
3123  if (simParameters->alchLambdaFreq > 0) {
3124  if (step <=
3125  simParameters->firstTimestep + simParameters->alchEquilSteps) {
3126  cumAlchWork = 0.0;
3127  }
3128  alchWork = computeAlchWork(step);
3129  cumAlchWork += alchWork;
3130  }
3131 //fepe
3132 
3133  momentum.x = reduction->item(REDUCTION_MOMENTUM_X);
3134  momentum.y = reduction->item(REDUCTION_MOMENTUM_Y);
3135  momentum.z = reduction->item(REDUCTION_MOMENTUM_Z);
3136  angularMomentum.x = reduction->item(REDUCTION_ANGULAR_MOMENTUM_X);
3137  angularMomentum.y = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Y);
3138  angularMomentum.z = reduction->item(REDUCTION_ANGULAR_MOMENTUM_Z);
3139 
3140  // Ported by JLai
3141  potentialEnergy = (bondEnergy + angleEnergy + dihedralEnergy
3142  + improperEnergy + electEnergy + electEnergySlow + ljEnergy
3143  + crosstermEnergy + boundaryEnergy + miscEnergy + goTotalEnergy
3145  // End of port
3146  totalEnergy = potentialEnergy + kineticEnergy;
3147  flatEnergy = (totalEnergy
3149  if ( !(step%slowFreq) ) {
3150  // only adjust based on most accurate energies
3152  if ( smooth2_avg == XXXBIGREAL ) smooth2_avg = s;
3153  if ( step != simParams->firstTimestep ) {
3154  smooth2_avg *= 0.9375;
3155  smooth2_avg += 0.0625 * s;
3156  }
3157  }
3158  smoothEnergy = (flatEnergy + smooth2_avg
3160 
3161  // Reset values for accumulated heat and work.
3162  if (step <= simParameters->firstTimestep &&
3163  (simParameters->stochRescaleOn && simParameters->stochRescaleHeat)) {
3164  heat = 0.0;
3166  }
3167  if ( simParameters->outputMomenta && ! minimize &&
3168  ! ( step % simParameters->outputMomenta ) )
3169  {
3170  iout << "MOMENTUM: " << step
3171  << " P: " << momentum
3172  << " L: " << angularMomentum
3173  << "\n" << endi;
3174  }
3175 
3176  if ( simParameters->outputPressure ) {
3179  tavg_count += 1;
3180  if ( minimize || ! ( step % simParameters->outputPressure ) ) {
3181  iout << "PRESSURE: " << step << " "
3182  << PRESSUREFACTOR * pressure << "\n"
3183  << "GPRESSURE: " << step << " "
3184  << PRESSUREFACTOR * groupPressure << "\n";
3185  if ( tavg_count > 1 ) iout << "PRESSAVG: " << step << " "
3186  << (PRESSUREFACTOR/tavg_count) * pressure_tavg << "\n"
3187  << "GPRESSAVG: " << step << " "
3189  iout << endi;
3190  pressure_tavg = 0;
3191  groupPressure_tavg = 0;
3192  tavg_count = 0;
3193  }
3194  }
3195 
3196  // pressure profile reductions
3197  if (pressureProfileSlabs) {
3198  const int freq = simParams->pressureProfileFreq;
3199  const int arraysize = 3*pressureProfileSlabs;
3200 
3201  BigReal *total = new BigReal[arraysize];
3202  memset(total, 0, arraysize*sizeof(BigReal));
3203  const int first = simParams->firstTimestep;
3204 
3205  if (ppbonded) ppbonded->getData(first, step, lattice, total);
3206  if (ppnonbonded) ppnonbonded->getData(first, step, lattice, total);
3207  if (ppint) ppint->getData(first, step, lattice, total);
3208  for (int i=0; i<arraysize; i++) pressureProfileAverage[i] += total[i];
3210 
3211  if (!(step % freq)) {
3212  // convert NAMD internal virial to pressure in units of bar
3214 
3215  iout << "PRESSUREPROFILE: " << step << " ";
3216  if (step == first) {
3217  // output pressure profile for this step
3218  for (int i=0; i<arraysize; i++) {
3219  iout << total[i] * scalefac << " ";
3220  }
3221  } else {
3222  // output pressure profile averaged over the last count steps.
3223  scalefac /= pressureProfileCount;
3224  for (int i=0; i<arraysize; i++)
3225  iout << pressureProfileAverage[i]*scalefac << " ";
3226  }
3227  iout << "\n" << endi;
3228 
3229  // Clear the average for the next block
3230  memset(pressureProfileAverage, 0, arraysize*sizeof(BigReal));
3232  }
3233  delete [] total;
3234  }
3235 
3236  if ( step != simParams->firstTimestep || stepInFullRun == 0 ) { // skip repeated first step
3237  if ( stepInFullRun % simParams->firstLdbStep == 0 ) {
3238  int benchPhase = stepInFullRun / simParams->firstLdbStep;
3239  if ( benchPhase > 0 && benchPhase < 10 ) {
3240 #ifdef NAMD_CUDA
3241  if ( simParams->outputEnergies < 60 ) {
3242  iout << iWARN << "Energy evaluation is expensive, increase outputEnergies to improve performance.\n";
3243  }
3244 #endif
3245  iout << iINFO;
3246  if ( benchPhase < 4 ) iout << "Initial time: ";
3247  else iout << "Benchmark time: ";
3248  iout << CkNumPes() << " CPUs ";
3249  {
3250  BigReal wallPerStep =
3251  (CmiWallTimer() - startBenchTime) / simParams->firstLdbStep;
3252  BigReal ns = simParams->dt / 1000000.0;
3253  BigReal days = 1.0 / (24.0 * 60.0 * 60.0);
3254  BigReal daysPerNano = wallPerStep * days / ns;
3255  iout << wallPerStep << " s/step " << daysPerNano << " days/ns ";
3256  iout << memusage_MB() << " MB memory\n" << endi;
3257  }
3258  }
3259  startBenchTime = CmiWallTimer();
3260  }
3261  ++stepInFullRun;
3262  }
3263 
3264  printTiming(step);
3265 
3266  Vector pairVDWForce, pairElectForce;
3267  if ( simParameters->pairInteractionOn ) {
3268  GET_VECTOR(pairVDWForce,reduction,REDUCTION_PAIR_VDW_FORCE);
3269  GET_VECTOR(pairElectForce,reduction,REDUCTION_PAIR_ELECT_FORCE);
3270  }
3271 
3272  // Compute cumulative nonequilibrium work (including shadow work).
3273  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3274  work = totalEnergy - totalEnergy0 - heat;
3275  }
3276 
3277  // callback to Tcl with whatever we can
3278 #ifdef NAMD_TCL
3279 #define CALLBACKDATA(LABEL,VALUE) \
3280  labels << (LABEL) << " "; values << (VALUE) << " ";
3281 #define CALLBACKLIST(LABEL,VALUE) \
3282  labels << (LABEL) << " "; values << "{" << (VALUE) << "} ";
3283  if (step == simParams->N && node->getScript() && node->getScript()->doCallback()) {
3284  std::ostringstream labels, values;
3285 #if CMK_BLUEGENEL
3286  // the normal version below gives a compiler error
3287  values.precision(16);
3288 #else
3289  values << std::setprecision(16);
3290 #endif
3291  CALLBACKDATA("TS",step);
3292  CALLBACKDATA("BOND",bondEnergy);
3293  CALLBACKDATA("ANGLE",angleEnergy);
3294  CALLBACKDATA("DIHED",dihedralEnergy);
3295  CALLBACKDATA("CROSS",crosstermEnergy);
3296  CALLBACKDATA("IMPRP",improperEnergy);
3298  CALLBACKDATA("VDW",ljEnergy);
3299  CALLBACKDATA("BOUNDARY",boundaryEnergy);
3300  CALLBACKDATA("MISC",miscEnergy);
3301  CALLBACKDATA("POTENTIAL",potentialEnergy);
3302  CALLBACKDATA("KINETIC",kineticEnergy);
3303  CALLBACKDATA("TOTAL",totalEnergy);
3304  CALLBACKDATA("TEMP",temperature);
3305  CALLBACKLIST("PRESSURE",pressure*PRESSUREFACTOR);
3306  CALLBACKLIST("GPRESSURE",groupPressure*PRESSUREFACTOR);
3307  CALLBACKDATA("VOLUME",lattice.volume());
3308  CALLBACKLIST("CELL_A",lattice.a());
3309  CALLBACKLIST("CELL_B",lattice.b());
3310  CALLBACKLIST("CELL_C",lattice.c());
3311  CALLBACKLIST("CELL_O",lattice.origin());
3312  labels << "PERIODIC "; values << "{" << lattice.a_p() << " "
3313  << lattice.b_p() << " " << lattice.c_p() << "} ";
3314  if (simParameters->drudeOn) {
3315  CALLBACKDATA("DRUDEBOND",drudeBondTemp);
3316  }
3317  if ( simParameters->pairInteractionOn ) {
3318  CALLBACKLIST("VDW_FORCE",pairVDWForce);
3319  CALLBACKLIST("ELECT_FORCE",pairElectForce);
3320  }
3321  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3322  CALLBACKLIST("HEAT",heat);
3323  CALLBACKLIST("WORK",work);
3324  }
3325  if (simParameters->alchOn) {
3326  if (simParameters->alchThermIntOn) {
3327  CALLBACKLIST("BOND1", bondedEnergy_ti_1);
3330  CALLBACKLIST("VDW1", ljEnergy_ti_1);
3331  CALLBACKLIST("BOND2", bondedEnergy_ti_2);
3334  CALLBACKLIST("VDW2", ljEnergy_ti_2);
3335  if (simParameters->alchLambdaFreq > 0) {
3336  CALLBACKLIST("CUMALCHWORK", cumAlchWork);
3337  }
3338  } else if (simParameters->alchFepOn) {
3339  CALLBACKLIST("BOND2", bondEnergy + angleEnergy + dihedralEnergy +
3340  improperEnergy + bondedEnergyDiff_f);
3342  CALLBACKLIST("VDW2", ljEnergy_f);
3343  }
3344  }
3345 
3346  labels << '\0'; values << '\0'; // insane but makes Linux work
3347  state->callback_labelstring = labels.str();
3348  state->callback_valuestring = values.str();
3349  // node->getScript()->doCallback(labelstring.c_str(),valuestring.c_str());
3350  }
3351 #undef CALLBACKDATA
3352 #endif
3353 
3355 
3356  temp_avg += temperature;
3357  pressure_avg += trace(pressure)/3.;
3358  groupPressure_avg += trace(groupPressure)/3.;
3359  avg_count += 1;
3360 
3362  ! (step % simParams->outputPairlists) ) {
3363  iout << iINFO << pairlistWarnings <<
3364  " pairlist warnings in past " << simParams->outputPairlists <<
3365  " steps.\n" << endi;
3366  pairlistWarnings = 0;
3367  }
3368 
3369  BigReal enthalpy;
3370  if (simParameters->multigratorOn && ((step % simParameters->outputEnergies) == 0)) {
3371  enthalpy = multigatorCalcEnthalpy(potentialEnergy, step, minimize);
3372  }
3373 
3374  // NO CALCULATIONS OR REDUCTIONS BEYOND THIS POINT!!!
3375  if ( ! minimize && step % simParameters->outputEnergies ) return;
3376  // ONLY OUTPUT SHOULD OCCUR BELOW THIS LINE!!!
3377 
3378  if (simParameters->IMDon && !(step % simParameters->IMDfreq)) {
3379  IMDEnergies energies;
3380  energies.tstep = step;
3381  energies.T = temp_avg/avg_count;
3382  energies.Etot = totalEnergy;
3383  energies.Epot = potentialEnergy;
3384  energies.Evdw = ljEnergy;
3385  energies.Eelec = electEnergy + electEnergySlow;
3386  energies.Ebond = bondEnergy;
3387  energies.Eangle = angleEnergy;
3388  energies.Edihe = dihedralEnergy + crosstermEnergy;
3389  energies.Eimpr = improperEnergy;
3390  Node::Object()->imd->gather_energies(&energies);
3391  }
3392 
3393  if ( marginViolations ) {
3394  iout << iERROR << marginViolations <<
3395  " margin violations detected since previous energy output.\n" << endi;
3396  }
3397  marginViolations = 0;
3398 
3399  if ( (step % (10 * (minimize?1:simParameters->outputEnergies) ) ) == 0 )
3400  {
3401  iout << "ETITLE: TS";
3402  iout << FORMAT("BOND");
3403  iout << FORMAT("ANGLE");
3404  iout << FORMAT("DIHED");
3405  if ( ! simParameters->mergeCrossterms ) iout << FORMAT("CROSS");
3406  iout << FORMAT("IMPRP");
3407  iout << " ";
3408  iout << FORMAT("ELECT");
3409  iout << FORMAT("VDW");
3410  iout << FORMAT("BOUNDARY");
3411  iout << FORMAT("MISC");
3412  iout << FORMAT("KINETIC");
3413  iout << " ";
3414  iout << FORMAT("TOTAL");
3415  iout << FORMAT("TEMP");
3416  iout << FORMAT("POTENTIAL");
3417  // iout << FORMAT("TOTAL2");
3418  iout << FORMAT("TOTAL3");
3419  iout << FORMAT("TEMPAVG");
3420  if ( volume != 0. ) {
3421  iout << " ";
3422  iout << FORMAT("PRESSURE");
3423  iout << FORMAT("GPRESSURE");
3424  iout << FORMAT("VOLUME");
3425  iout << FORMAT("PRESSAVG");
3426  iout << FORMAT("GPRESSAVG");
3427  }
3428  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3429  iout << " ";
3430  iout << FORMAT("HEAT");
3431  iout << FORMAT("WORK");
3432  }
3433  if (simParameters->drudeOn) {
3434  iout << " ";
3435  iout << FORMAT("DRUDEBOND");
3436  iout << FORMAT("DRBONDAVG");
3437  }
3438  // Ported by JLai
3439  if (simParameters->goGroPair) {
3440  iout << " ";
3441  iout << FORMAT("GRO_PAIR_LJ");
3442  iout << FORMAT("GRO_PAIR_GAUSS");
3443  }
3444 
3445  if (simParameters->goForcesOn) {
3446  iout << " ";
3447  iout << FORMAT("NATIVE");
3448  iout << FORMAT("NONNATIVE");
3449  //iout << FORMAT("REL_NATIVE");
3450  //iout << FORMAT("REL_NONNATIVE");
3451  iout << FORMAT("GOTOTAL");
3452  //iout << FORMAT("GOAVG");
3453  }
3454  // End of port -- JLai
3455 
3456  if (simParameters->alchOn) {
3457  if (simParameters->alchThermIntOn) {
3458  iout << "\nTITITLE: TS";
3459  iout << FORMAT("BOND1");
3460  iout << FORMAT("ELECT1");
3461  iout << FORMAT("VDW1");
3462  iout << FORMAT("BOND2");
3463  iout << " ";
3464  iout << FORMAT("ELECT2");
3465  iout << FORMAT("VDW2");
3466  if (simParameters->alchLambdaFreq > 0) {
3467  iout << FORMAT("LAMBDA");
3468  iout << FORMAT("ALCHWORK");
3469  iout << FORMAT("CUMALCHWORK");
3470  }
3471  } else if (simParameters->alchFepOn) {
3472  iout << "\nFEPTITLE: TS";
3473  iout << FORMAT("BOND2");
3474  iout << FORMAT("ELECT2");
3475  iout << FORMAT("VDW2");
3476  if (simParameters->alchLambdaFreq > 0) {
3477  iout << FORMAT("LAMBDA");
3478  }
3479  }
3480  }
3481 
3482  iout << "\n\n" << endi;
3483 
3484  if (simParameters->qmForcesOn) {
3485  iout << "QMETITLE: TS";
3486  iout << FORMAT("QMID");
3487  iout << FORMAT("ENERGY");
3488  if (simParameters->PMEOn) iout << FORMAT("PMECORRENERGY");
3489  iout << "\n\n" << endi;
3490  }
3491 
3492  }
3493 
3494  // N.B. HP's aCC compiler merges FORMAT calls in the same expression.
3495  // Need separate statements because data returned in static array.
3496  iout << ETITLE(step);
3497  iout << FORMAT(bondEnergy);
3498  iout << FORMAT(angleEnergy);
3499  if ( simParameters->mergeCrossterms ) {
3500  iout << FORMAT(dihedralEnergy+crosstermEnergy);
3501  } else {
3502  iout << FORMAT(dihedralEnergy);
3503  iout << FORMAT(crosstermEnergy);
3504  }
3505  iout << FORMAT(improperEnergy);
3506  iout << " ";
3508  iout << FORMAT(ljEnergy);
3509  iout << FORMAT(boundaryEnergy);
3510  iout << FORMAT(miscEnergy);
3511  iout << FORMAT(kineticEnergy);
3512  iout << " ";
3513  iout << FORMAT(totalEnergy);
3514  iout << FORMAT(temperature);
3515  iout << FORMAT(potentialEnergy);
3516  // iout << FORMAT(flatEnergy);
3517  iout << FORMAT(smoothEnergy);
3519  if ( volume != 0. )
3520  {
3521  iout << " ";
3522  iout << FORMAT(trace(pressure)*PRESSUREFACTOR/3.);
3523  iout << FORMAT(trace(groupPressure)*PRESSUREFACTOR/3.);
3524  iout << FORMAT(volume);
3527  }
3528  if (simParameters->stochRescaleOn && simParameters->stochRescaleHeat) {
3529  iout << " ";
3530  iout << FORMAT(heat);
3531  iout << FORMAT(work);
3532  }
3533  if (simParameters->drudeOn) {
3534  iout << " ";
3537  }
3538  // Ported by JLai
3539  if (simParameters->goGroPair) {
3540  iout << " ";
3541  iout << FORMAT(groLJEnergy);
3543  }
3544 
3545  if (simParameters->goForcesOn) {
3546  iout << " ";
3549  //iout << FORMAT(relgoNativeEnergy);
3550  //iout << FORMAT(relgoNonnativeEnergy);
3552  //iout << FORMAT("not implemented");
3553  } // End of port -- JLai
3554 
3555  if (simParameters->alchOn) {
3556  if (simParameters->alchThermIntOn) {
3557  iout << "\n";
3558  iout << TITITLE(step);
3564  iout << " ";
3568  if (simParameters->alchLambdaFreq > 0) {
3569  iout << FORMAT(simParameters->getCurrentLambda(step));
3570  iout << FORMAT(alchWork);
3571  iout << FORMAT(cumAlchWork);
3572  }
3573  } else if (simParameters->alchFepOn) {
3574  iout << "\n";
3575  iout << FEPTITLE2(step);
3576  iout << FORMAT(bondEnergy + angleEnergy + dihedralEnergy
3577  + improperEnergy + bondedEnergyDiff_f);
3579  iout << FORMAT(ljEnergy_f);
3580  if (simParameters->alchLambdaFreq > 0) {
3581  iout << FORMAT(simParameters->getCurrentLambda(step));
3582  }
3583  }
3584  }
3585 
3586  iout << "\n\n" << endi;
3587 
3588 #if(CMK_CCS_AVAILABLE && CMK_WEB_MODE)
3589  char webout[80];
3590  sprintf(webout,"%d %d %d %d",(int)totalEnergy,
3591  (int)(potentialEnergy),
3592  (int)kineticEnergy,(int)temperature);
3593  CApplicationDepositNode0Data(webout);
3594 #endif
3595 
3596  if (simParameters->pairInteractionOn) {
3597  iout << "PAIR INTERACTION:";
3598  iout << " STEP: " << step;
3599  iout << " VDW_FORCE: ";
3600  iout << FORMAT(pairVDWForce.x);
3601  iout << FORMAT(pairVDWForce.y);
3602  iout << FORMAT(pairVDWForce.z);
3603  iout << " ELECT_FORCE: ";
3604  iout << FORMAT(pairElectForce.x);
3605  iout << FORMAT(pairElectForce.y);
3606  iout << FORMAT(pairElectForce.z);
3607  iout << "\n" << endi;
3608  }
3609  drudeBondTempAvg = 0;
3610  temp_avg = 0;
3611  pressure_avg = 0;
3612  groupPressure_avg = 0;
3613  avg_count = 0;
3614 
3615  if ( fflush_count ) {
3616  --fflush_count;
3617  fflush(stdout);
3618  }
3619 }
3620 
3622  Lattice &lattice = state->lattice;
3623  file << "#$LABELS step";
3624  if ( lattice.a_p() ) file << " a_x a_y a_z";
3625  if ( lattice.b_p() ) file << " b_x b_y b_z";
3626  if ( lattice.c_p() ) file << " c_x c_y c_z";
3627  file << " o_x o_y o_z";
3628  if ( simParams->langevinPistonOn ) {
3629  file << " s_x s_y s_z s_u s_v s_w";
3630  }
3631  file << std::endl;
3632 }
3633 
3635  Lattice &lattice = state->lattice;
3636  file.precision(12);
3637  file << step;
3638  if ( lattice.a_p() ) file << " " << lattice.a().x << " " << lattice.a().y << " " << lattice.a().z;
3639  if ( lattice.b_p() ) file << " " << lattice.b().x << " " << lattice.b().y << " " << lattice.b().z;
3640  if ( lattice.c_p() ) file << " " << lattice.c().x << " " << lattice.c().y << " " << lattice.c().z;
3641  file << " " << lattice.origin().x << " " << lattice.origin().y << " " << lattice.origin().z;
3642  if ( simParams->langevinPistonOn ) {
3643  Vector strainRate = diagonal(langevinPiston_strainRate);
3644  Vector strainRate2 = off_diagonal(langevinPiston_strainRate);
3645  file << " " << strainRate.x;
3646  file << " " << strainRate.y;
3647  file << " " << strainRate.z;
3648  file << " " << strainRate2.x;
3649  file << " " << strainRate2.y;
3650  file << " " << strainRate2.z;
3651  }
3652  file << std::endl;
3653 }
3654 
3656 {
3657  if ( Output::coordinateNeeded(timestep) ){
3658  collection->enqueuePositions(timestep,state->lattice);
3659  }
3660  if ( Output::velocityNeeded(timestep) )
3661  collection->enqueueVelocities(timestep);
3662  if ( Output::forceNeeded(timestep) )
3663  collection->enqueueForces(timestep);
3664 }
3665 
3666 //Modifications for alchemical fep
3668  if (simParams->alchOn && simParams->alchFepOn) {
3669  const int stepInRun = step - simParams->firstTimestep;
3670  const int alchEquilSteps = simParams->alchEquilSteps;
3671  const BigReal alchLambda = simParams->alchLambda;
3672  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
3673 
3674  if (alchEnsembleAvg && (stepInRun == 0 || stepInRun == alchEquilSteps)) {
3675  FepNo = 0;
3676  exp_dE_ByRT = 0.0;
3677  net_dE = 0.0;
3678  }
3682 
3683  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. ||
3684  (step / simParams->alchIDWSFreq) % 2 == 1 )) {
3685  // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
3686  FepNo++;
3687  exp_dE_ByRT += exp(-dE/RT);
3688  net_dE += dE;
3689  }
3690 
3691  if (simParams->alchOutFreq) {
3692  if (stepInRun == 0) {
3693  if (!fepFile.is_open()) {
3696  iout << "OPENING FEP ENERGY OUTPUT FILE\n" << endi;
3697  if(alchEnsembleAvg){
3698  fepSum = 0.0;
3699  fepFile << "# STEP Elec "
3700  << "vdW dE dE_avg Temp dG\n"
3701  << "# l l+dl "
3702  << " l l+dl E(l+dl)-E(l)" << std::endl;
3703  }
3704  else{
3705  fepFile << "# STEP Elec "
3706  << "vdW dE Temp\n"
3707  << "# l l+dl "
3708  << " l l+dl E(l+dl)-E(l)" << std::endl;
3709  }
3710  }
3711  if(!step){
3712  fepFile << "#NEW FEP WINDOW: "
3713  << "LAMBDA SET TO " << alchLambda << " LAMBDA2 "
3714  << simParams->alchLambda2;
3715  if ( simParams->alchLambdaIDWS >= 0. ) {
3716  fepFile << " LAMBDA_IDWS " << simParams->alchLambdaIDWS;
3717  }
3718  fepFile << std::endl;
3719  }
3720  }
3721  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3722  fepFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3723  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3724  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3725  }
3726  if ((simParams->N) && ((step%simParams->alchOutFreq) == 0)) {
3727  writeFepEnergyData(step, fepFile);
3728  fepFile.flush();
3729  }
3730  if (alchEnsembleAvg && (step == simParams->N)) {
3731  fepSum = fepSum + dG;
3732  fepFile << "#Free energy change for lambda window [ " << alchLambda
3733  << " " << simParams->alchLambda2 << " ] is " << dG
3734  << " ; net change until now is " << fepSum << std::endl;
3735  fepFile.flush();
3736  }
3737  }
3738  }
3739 }
3740 
3743  const int stepInRun = step - simParams->firstTimestep;
3744  const int alchEquilSteps = simParams->alchEquilSteps;
3745  const int alchLambdaFreq = simParams->alchLambdaFreq;
3746 
3747  if (stepInRun == 0 || stepInRun == alchEquilSteps) {
3748  TiNo = 0;
3749  net_dEdl_bond_1 = 0;
3750  net_dEdl_bond_2 = 0;
3751  net_dEdl_elec_1 = 0;
3752  net_dEdl_elec_2 = 0;
3753  net_dEdl_lj_1 = 0;
3754  net_dEdl_lj_2 = 0;
3755  }
3756  TiNo++;
3765 
3766  // Don't accumulate block averages (you'll get a divide by zero!) or even
3767  // open the TI output if alchOutFreq is zero.
3768  if (simParams->alchOutFreq) {
3769  if (stepInRun == 0 || stepInRun == alchEquilSteps
3770  || (! ((step - 1) % simParams->alchOutFreq))) {
3771  // output of instantaneous dU/dl now replaced with running average
3772  // over last alchOutFreq steps (except for step 0)
3773  recent_TiNo = 0;
3774  recent_dEdl_bond_1 = 0;
3775  recent_dEdl_bond_2 = 0;
3776  recent_dEdl_elec_1 = 0;
3777  recent_dEdl_elec_2 = 0;
3778  recent_dEdl_lj_1 = 0;
3779  recent_dEdl_lj_2 = 0;
3780  recent_alchWork = 0;
3781  }
3782  recent_TiNo++;
3786  + electEnergyPME_ti_1);
3788  + electEnergyPME_ti_2);
3792 
3793  if (stepInRun == 0) {
3794  if (!tiFile.is_open()) {
3797  /* BKR - This has been rather drastically updated to better match
3798  stdout. This was necessary for several reasons:
3799  1) PME global scaling is obsolete (now removed)
3800  2) scaling of bonded terms was added
3801  3) alchemical work is now accumulated when switching is active
3802  */
3803  iout << "OPENING TI ENERGY OUTPUT FILE\n" << endi;
3804  tiFile << "#TITITLE: TS";
3805  tiFile << FORMAT("BOND1");
3806  tiFile << FORMAT("AVGBOND1");
3807  tiFile << FORMAT("ELECT1");
3808  tiFile << FORMAT("AVGELECT1");
3809  tiFile << " ";
3810  tiFile << FORMAT("VDW1");
3811  tiFile << FORMAT("AVGVDW1");
3812  tiFile << FORMAT("BOND2");
3813  tiFile << FORMAT("AVGBOND2");
3814  tiFile << FORMAT("ELECT2");
3815  tiFile << " ";
3816  tiFile << FORMAT("AVGELECT2");
3817  tiFile << FORMAT("VDW2");
3818  tiFile << FORMAT("AVGVDW2");
3819  if (alchLambdaFreq > 0) {
3820  tiFile << FORMAT("ALCHWORK");
3821  tiFile << FORMAT("CUMALCHWORK");
3822  }
3823  tiFile << std::endl;
3824  }
3825 
3826  if (alchLambdaFreq > 0) {
3827  tiFile << "#ALCHEMICAL SWITCHING ACTIVE "
3828  << simParams->alchLambda << " --> " << simParams->getCurrentLambda2(step)
3829  << "\n#LAMBDA SCHEDULE: "
3830  << "dL: " << simParams->getLambdaDelta()
3831  << " Freq: " << alchLambdaFreq;
3832  }
3833  else {
3834  const BigReal alchLambda = simParams->alchLambda;
3835  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
3836  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
3837  const BigReal elec_lambda_1 = simParams->getElecLambda(alchLambda);
3838  const BigReal elec_lambda_2 = simParams->getElecLambda(1-alchLambda);
3839  const BigReal vdw_lambda_1 = simParams->getVdwLambda(alchLambda);
3840  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-alchLambda);
3841  tiFile << "#NEW TI WINDOW: "
3842  << "LAMBDA " << alchLambda
3843  << "\n#PARTITION 1 SCALING: BOND " << bond_lambda_1
3844  << " VDW " << vdw_lambda_1 << " ELEC " << elec_lambda_1
3845  << "\n#PARTITION 2 SCALING: BOND " << bond_lambda_2
3846  << " VDW " << vdw_lambda_2 << " ELEC " << elec_lambda_2;
3847  }
3848  tiFile << "\n#CONSTANT TEMPERATURE: " << simParams->alchTemp << " K"
3849  << std::endl;
3850  }
3851 
3852  if ((alchEquilSteps) && (stepInRun == alchEquilSteps)) {
3853  tiFile << "#" << alchEquilSteps << " STEPS OF EQUILIBRATION AT "
3854  << "LAMBDA " << simParams->alchLambda << " COMPLETED\n"
3855  << "#STARTING COLLECTION OF ENSEMBLE AVERAGE" << std::endl;
3856  }
3857  if ((step%simParams->alchOutFreq) == 0) {
3858  writeTiEnergyData(step, tiFile);
3859  tiFile.flush();
3860  }
3861  }
3862  }
3863 }
3864 
3865 /*
3866  Work is accumulated whenever alchLambda changes. In the current scheme we
3867  always increment lambda _first_, then integrate in time. Therefore the work
3868  is wrt the "old" lambda before the increment.
3869 */
3871  // alchemical scaling factors for groups 1/2 at the previous lambda
3872  const BigReal oldLambda = simParams->getCurrentLambda(step-1);
3873  const BigReal bond_lambda_1 = simParams->getBondLambda(oldLambda);
3874  const BigReal bond_lambda_2 = simParams->getBondLambda(1-oldLambda);
3875  const BigReal elec_lambda_1 = simParams->getElecLambda(oldLambda);
3876  const BigReal elec_lambda_2 = simParams->getElecLambda(1-oldLambda);
3877  const BigReal vdw_lambda_1 = simParams->getVdwLambda(oldLambda);
3878  const BigReal vdw_lambda_2 = simParams->getVdwLambda(1-oldLambda);
3879  // alchemical scaling factors for groups 1/2 at the new lambda
3880  const BigReal alchLambda = simParams->getCurrentLambda(step);
3881  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda);
3882  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda);
3883  const BigReal elec_lambda_12 = simParams->getElecLambda(alchLambda);
3884  const BigReal elec_lambda_22 = simParams->getElecLambda(1-alchLambda);
3885  const BigReal vdw_lambda_12 = simParams->getVdwLambda(alchLambda);
3886  const BigReal vdw_lambda_22 = simParams->getVdwLambda(1-alchLambda);
3887 
3888  return ((bond_lambda_12 - bond_lambda_1)*bondedEnergy_ti_1 +
3889  (elec_lambda_12 - elec_lambda_1)*
3891  (vdw_lambda_12 - vdw_lambda_1)*ljEnergy_ti_1 +
3892  (bond_lambda_22 - bond_lambda_2)*bondedEnergy_ti_2 +
3893  (elec_lambda_22 - elec_lambda_2)*
3895  (vdw_lambda_22 - vdw_lambda_2)*ljEnergy_ti_2
3896  );
3897 }
3898 
3903 
3904  const bool alchEnsembleAvg = simParams->alchEnsembleAvg;
3905  const int stepInRun = step - simParams->firstTimestep;
3906  const BigReal alchLambda = simParams->alchLambda;
3907 
3908  if(stepInRun){
3909  if ( simParams->alchLambdaIDWS >= 0. &&
3910  (step / simParams->alchIDWSFreq) % 2 == 0 ) {
3911  // IDWS is active and we are on a "backward" timestep
3912  fepFile << FEPTITLE_BACK(step);
3913  } else {
3914  // "forward" timestep
3915  fepFile << FEPTITLE(step);
3916  }
3917  fepFile << FORMAT(eeng);
3918  fepFile << FORMAT(eeng_f);
3919  fepFile << FORMAT(ljEnergy);
3921  fepFile << FORMAT(dE);
3922  if(alchEnsembleAvg){
3923  BigReal dE_avg = net_dE / FepNo;
3924  fepFile << FORMAT(dE_avg);
3925  }
3927  if(alchEnsembleAvg){
3928  dG = -(RT * log(exp_dE_ByRT / FepNo));
3929  fepFile << FORMAT(dG);
3930  }
3931  fepFile << std::endl;
3932  }
3933 }
3934 
3936  tiFile << TITITLE(step);
3941  tiFile << " ";
3947  tiFile << " ";
3951  if (simParams->alchLambdaFreq > 0) {
3954  }
3955  tiFile << std::endl;
3956 }
3957 //fepe
3958 
3960 {
3961 
3962  if ( step >= 0 ) {
3963 
3964  // Write out eXtended System Trajectory (XST) file
3965  if ( simParams->xstFrequency &&
3966  ((step % simParams->xstFrequency) == 0) )
3967  {
3968  if ( ! xstFile.is_open() )
3969  {
3970  iout << "OPENING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
3973  while (!xstFile) {
3974  if ( errno == EINTR ) {
3975  CkPrintf("Warning: Interrupted system call opening XST trajectory file, retrying.\n");
3976  xstFile.clear();
3978  continue;
3979  }
3980  char err_msg[257];
3981  sprintf(err_msg, "Error opening XST trajectory file %s",simParams->xstFilename);
3982  NAMD_err(err_msg);
3983  }
3984  xstFile << "# NAMD extended system trajectory file" << std::endl;
3986  }
3988  xstFile.flush();
3989  }
3990 
3991  // Write out eXtended System Configuration (XSC) files
3992  // Output a restart file
3993  if ( simParams->restartFrequency &&
3994  ((step % simParams->restartFrequency) == 0) &&
3995  (step != simParams->firstTimestep) )
3996  {
3997  iout << "WRITING EXTENDED SYSTEM TO RESTART FILE AT STEP "
3998  << step << "\n" << endi;
3999  char fname[140];
4000  const char *bsuffix = ".old";
4001  strcpy(fname, simParams->restartFilename);
4002  if ( simParams->restartSave ) {
4003  char timestepstr[20];
4004  sprintf(timestepstr,".%d",step);
4005  strcat(fname, timestepstr);
4006  bsuffix = ".BAK";
4007  }
4008  strcat(fname, ".xsc");
4009  NAMD_backup_file(fname,bsuffix);
4010  ofstream_namd xscFile(fname);
4011  while (!xscFile) {
4012  if ( errno == EINTR ) {
4013  CkPrintf("Warning: Interrupted system call opening XSC restart file, retrying.\n");
4014  xscFile.clear();
4015  xscFile.open(fname);
4016  continue;
4017  }
4018  char err_msg[257];
4019  sprintf(err_msg, "Error opening XSC restart file %s",fname);
4020  NAMD_err(err_msg);
4021  }
4022  xscFile << "# NAMD extended system configuration restart file" << std::endl;
4023  writeExtendedSystemLabels(xscFile);
4024  writeExtendedSystemData(step,xscFile);
4025  if (!xscFile) {
4026  char err_msg[257];
4027  sprintf(err_msg, "Error writing XSC restart file %s",fname);
4028  NAMD_err(err_msg);
4029  }
4030  }
4031 
4032  }
4033 
4034  // Output final coordinates
4035  if (step == FILE_OUTPUT || step == END_OF_RUN)
4036  {
4037  int realstep = ( step == FILE_OUTPUT ?
4039  iout << "WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP "
4040  << realstep << "\n" << endi;
4041  static char fname[140];
4042  strcpy(fname, simParams->outputFilename);
4043  strcat(fname, ".xsc");
4044  NAMD_backup_file(fname);
4045  ofstream_namd xscFile(fname);
4046  while (!xscFile) {
4047  if ( errno == EINTR ) {
4048  CkPrintf("Warning: Interrupted system call opening XSC output file, retrying.\n");
4049  xscFile.clear();
4050  xscFile.open(fname);
4051  continue;
4052  }
4053  char err_msg[257];
4054  sprintf(err_msg, "Error opening XSC output file %s",fname);
4055  NAMD_err(err_msg);
4056  }
4057  xscFile << "# NAMD extended system configuration output file" << std::endl;
4058  writeExtendedSystemLabels(xscFile);
4059  writeExtendedSystemData(realstep,xscFile);
4060  if (!xscFile) {
4061  char err_msg[257];
4062  sprintf(err_msg, "Error writing XSC output file %s",fname);
4063  NAMD_err(err_msg);
4064  }
4065  }
4066 
4067  // Close trajectory file
4068  if (step == END_OF_RUN) {
4069  if ( xstFile.is_open() ) {
4070  xstFile.close();
4071  iout << "CLOSING EXTENDED SYSTEM TRAJECTORY FILE\n" << endi;
4072  }
4073  }
4074 
4075 }
4076 
4077 
4078 void Controller::recvCheckpointReq(const char *key, int task, checkpoint &cp) { // responding replica
4079  switch ( task ) {
4081  if ( ! checkpoints.count(key) ) {
4082  checkpoints[key] = new checkpoint;
4083  }
4084  *checkpoints[key] = cp;
4085  break;
4087  if ( ! checkpoints.count(key) ) {
4088  NAMD_die("Unable to load checkpoint, requested key was never stored.");
4089  }
4090  cp = *checkpoints[key];
4091  break;
4093  if ( ! checkpoints.count(key) ) {
4094  NAMD_die("Unable to swap checkpoint, requested key was never stored.");
4095  }
4096  std::swap(cp,*checkpoints[key]);
4097  break;
4099  if ( ! checkpoints.count(key) ) {
4100  NAMD_die("Unable to free checkpoint, requested key was never stored.");
4101  }
4102  delete checkpoints[key];
4103  checkpoints.erase(key);
4104  break;
4105  }
4106 }
4107 
4108 void Controller::recvCheckpointAck(checkpoint &cp) { // initiating replica
4110  state->lattice = cp.lattice;
4111  *(ControllerState*)this = cp.state;
4112  }
4113  CkpvAccess(_qd)->process();
4114 }
4115 
4116 
4118 {
4119  if ( ! ldbSteps ) {
4121  }
4122  if ( ! --ldbSteps ) {
4123  startBenchTime -= CmiWallTimer();
4124  Node::Object()->outputPatchComputeMaps("before_ldb", step);
4126  startBenchTime += CmiWallTimer();
4127  fflush_count = 3;
4128  }
4129 }
4130 
4131 void Controller::cycleBarrier(int doBarrier, int step) {
4132 #if USE_BARRIER
4133  if (doBarrier) {
4134  broadcast->cycleBarrier.publish(step,1);
4135  CkPrintf("Cycle time at sync Wall: %f CPU %f\n",
4136  CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4137  }
4138 #endif
4139 }
4140 
4141 void Controller::traceBarrier(int turnOnTrace, int step) {
4142  CkPrintf("Cycle time at trace sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4143  CProxy_Node nd(CkpvAccess(BOCclass_group).node);
4144  nd.traceBarrier(turnOnTrace, step);
4145  CthSuspend();
4146 }
4147 
4149  broadcast->traceBarrier.publish(step,1);
4150  CkPrintf("Cycle time at trace sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4151  awaken();
4152 }
4153 
4154 #ifdef MEASURE_NAMD_WITH_PAPI
4155 void Controller::papiMeasureBarrier(int turnOnMeasure, int step){
4156  CkPrintf("Cycle time at PAPI measurement sync (begin) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4157  CProxy_Node nd(CkpvAccess(BOCclass_group).node);
4158  nd.papiMeasureBarrier(turnOnMeasure, step);
4159  CthSuspend();
4160 }
4161 
4162 void Controller::resumeAfterPapiMeasureBarrier(int step){
4163  broadcast->papiMeasureBarrier.publish(step,1);
4164  CkPrintf("Cycle time at PAPI measurement sync (end) Wall at step %d: %f CPU %f\n", step, CmiWallTimer()-firstWTime,CmiTimer()-firstCTime);
4165  awaken();
4166 }
4167 #endif
4168 
4170  BackEnd::awaken();
4171  CthFree(thread);
4172  CthSuspend();
4173 }
4174 
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:603
Bool accelMDGresetVaftercmd
BigReal berendsenPressureCompressibility
BigReal net_dE
Definition: Controller.h:120
void enqueueCollections(int)
Definition: Controller.C:3655
#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:107
#define XXXBIGREAL
Definition: Controller.C:60
char restartFilename[128]
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:4078
void cycleBarrier(int, int)
Definition: Controller.C:4131
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:1831
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:1773
BigReal getLambdaDelta(void)
int numCalcBonds
Definition: Molecule.h:621
ControllerState state
Definition: Controller.h:274
BigReal adaptTempCgamma
int numCalcAnisos
Definition: Molecule.h:631
void compareChecksums(int, int=0)
Definition: Controller.C:2767
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:102
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:1590
int eventEndOfTimeStep
Definition: Node.C:285
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:45
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:2443
int numHydrogenGroups
Definition: Molecule.h:599
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:1708
Tensor controlPressure
Definition: Controller.h:170
void adaptTempInit(int step)
Definition: Controller.C:2343
void minimize()
Definition: Controller.C:594
BigReal totalEnergy0
Definition: Controller.h:164
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cudaTextureObject_t cudaTextureObject_t cudaTextureObject_t const int const 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__ float4 *__restrict__ float4 *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ float *__restrict__ const int numPatches
#define EVAL_MEASURE
Definition: Output.h:27
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:2996
void integrate(int)
Definition: Controller.C:429
BigReal multigratorPressureTarget
#define PRESSUREFACTOR
Definition: common.h:47
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:3621
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
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:108
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:87
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:497
BigReal tail_corr_ener
Definition: Molecule.h:470
#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:1512
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
void multigratorTemperature(int step, int callNumber)
Definition: Controller.C:853
int64_t num_deg_freedom(int isInitialReport=0) const
Definition: Molecule.h:523
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:3741
static char * FEPTITLE2(int X)
Definition: Controller.h:357
BigReal electEnergy_f
Definition: Controller.h:114
int numCalcCrossterms
Definition: Molecule.h:625
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
char outputFilename[128]
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
<