NAMD
Sequencer.C
Go to the documentation of this file.
1 
7 /*****************************************************************************
8  * $Source: /home/cvs/namd/cvsroot/namd2/src/Sequencer.C,v $
9  * $Author: jim $
10  * $Date: 2016/08/26 19:40:32 $
11  * $Revision: 1.1230 $
12  *****************************************************************************/
13 
14 // The UPPER_BOUND macro is used to eliminate all of the per atom
15 // computation done for the numerical integration in Sequencer::integrate()
16 // other than the actual force computation and atom migration.
17 // The idea is to "turn off" the integration for doing performance
18 // profiling in order to get an upper bound on the speedup available
19 // by moving the integration parts to the GPU.
20 //
21 // Define it in the Make.config file, i.e. CXXOPTS += -DUPPER_BOUND
22 // or simply uncomment the line below.
23 //
24 //#define UPPER_BOUND
25 
26 //for gbis debugging; print net force on each atom
27 #define PRINT_FORCES 0
28 
29 #include "InfoStream.h"
30 #include "Node.h"
31 #include "SimParameters.h"
32 #include "Sequencer.h"
33 #include "HomePatch.h"
34 #include "ReductionMgr.h"
35 #include "CollectionMgr.h"
36 #include "BroadcastObject.h"
37 #include "Output.h"
38 #include "Controller.h"
39 #include "Broadcasts.h"
40 #include "Molecule.h"
41 #include "NamdOneTools.h"
42 #include "LdbCoordinator.h"
43 #include "Thread.h"
44 #include "Random.h"
45 #include "PatchMap.inl"
46 #include "ComputeMgr.h"
47 #include "ComputeGlobal.h"
48 #include "NamdEventsProfiling.h"
49 
50 #define MIN_DEBUG_LEVEL 4
51 //#define DEBUGM
52 //
53 // Define NL_DEBUG below to activate D_*() macros in integrate_SOA()
54 // for debugging.
55 //
56 //#define NL_DEBUG
57 #include "Debug.h"
58 
59 #if USE_HPM
60 #define START_HPM_STEP 1000
61 #define STOP_HPM_STEP 1500
62 #endif
63 
64 #define SPECIAL_PATCH_ID 91
65 
67  simParams(Node::Object()->simParameters),
68  patch(p),
69  collection(CollectionMgr::Object()),
70  ldbSteps(0)
71 {
77  int ntypes = simParams->pressureProfileAtomTypes;
78  int nslabs = simParams->pressureProfileSlabs;
81  REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
82  } else {
84  }
85  if (simParams->multigratorOn) {
87  } else {
88  multigratorReduction = NULL;
89  }
90  ldbCoordinator = (LdbCoordinator::Object());
93 
94  // Is soluteScaling enabled?
96  // If so, we must "manually" perform charge scaling on startup because
97  // Sequencer will not get a scripting task for initial charge scaling.
98  // Subsequent rescalings will take place through a scripting task.
100  }
101 
103  stochRescale_count = 0;
105 // patch->write_tip4_props();
106 }
107 
109 {
110  delete broadcast;
111  delete reduction;
112  delete min_reduction;
114  delete random;
116 }
117 
118 // Invoked by thread
119 void Sequencer::threadRun(Sequencer* arg)
120 {
122  arg->algorithm();
123 }
124 
125 // Invoked by Node::run() via HomePatch::runSequencer()
126 void Sequencer::run(void)
127 {
128  // create a Thread and invoke it
129  DebugM(4, "::run() - this = " << this << "\n" );
130  thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
131  CthSetStrategyDefault(thread);
132  priority = PATCH_PRIORITY(patch->getPatchID());
133  awaken();
134 }
135 
137 {
139  CthSuspend();
141 }
142 
143 // Defines sequence of operations on a patch. e.g. when
144 // to push out information for Compute objects to consume
145 // when to migrate atoms, when to add forces to velocity update.
147 {
148  int scriptTask;
149  int scriptSeq = 0;
150  // Blocking receive for the script barrier.
151  while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
152  switch ( scriptTask ) {
153  case SCRIPT_OUTPUT:
155  break;
156  case SCRIPT_FORCEOUTPUT:
158  break;
159  case SCRIPT_MEASURE:
161  break;
162  case SCRIPT_REINITVELS:
164  break;
165  case SCRIPT_RESCALEVELS:
167  break;
170  break;
172  reloadCharges();
173  break;
174  case SCRIPT_CHECKPOINT:
175  patch->checkpoint();
177  break;
178  case SCRIPT_REVERT:
179  patch->revert();
181  pairlistsAreValid = 0;
182  break;
188  break;
189  case SCRIPT_ATOMSENDRECV:
190  case SCRIPT_ATOMSEND:
191  case SCRIPT_ATOMRECV:
192  patch->exchangeAtoms(scriptTask);
193  break;
194  case SCRIPT_MINIMIZE:
195  minimize();
196  break;
197  case SCRIPT_RUN:
198  case SCRIPT_CONTINUE:
199  //
200  // DJH: Call a cleaned up version of integrate().
201  //
202  // We could test for simulation options and call a more basic version
203  // of integrate() where we can avoid performing most tests.
204  //
205  integrate(scriptTask);
206  break;
207  default:
208  NAMD_bug("Unknown task in Sequencer::algorithm");
209  }
210  }
212  terminate();
213 }
214 
215 
216 extern int eventEndOfTimeStep;
217 
218 void Sequencer::integrate(int scriptTask) {
219  char traceNote[24];
220  char tracePrefix[20];
221  sprintf(tracePrefix, "p:%d,s:",patch->patchID);
222 // patch->write_tip4_props();
223 
224  //
225  // DJH: Copy all data into SOA (structure of arrays)
226  // from AOS (array of structures) data structure.
227  //
228  //patch->copy_all_to_SOA();
229 
230 #ifdef TIMER_COLLECTION
231  TimerSet& t = patch->timerSet;
232 #endif
233  TIMER_INIT_WIDTH(t, KICK, simParams->timerBinWidth);
234  TIMER_INIT_WIDTH(t, MAXMOVE, simParams->timerBinWidth);
235  TIMER_INIT_WIDTH(t, DRIFT, simParams->timerBinWidth);
236  TIMER_INIT_WIDTH(t, PISTON, simParams->timerBinWidth);
237  TIMER_INIT_WIDTH(t, SUBMITHALF, simParams->timerBinWidth);
238  TIMER_INIT_WIDTH(t, VELBBK1, simParams->timerBinWidth);
239  TIMER_INIT_WIDTH(t, VELBBK2, simParams->timerBinWidth);
240  TIMER_INIT_WIDTH(t, RATTLE1, simParams->timerBinWidth);
241  TIMER_INIT_WIDTH(t, SUBMITFULL, simParams->timerBinWidth);
242  TIMER_INIT_WIDTH(t, SUBMITCOLLECT, simParams->timerBinWidth);
243 
244  int &step = patch->flags.step;
245  step = simParams->firstTimestep;
246 
247  // drag switches
248  const Bool rotDragOn = simParams->rotDragOn;
249  const Bool movDragOn = simParams->movDragOn;
250 
251  const int commOnly = simParams->commOnly;
252 
253  int &maxForceUsed = patch->flags.maxForceUsed;
254  int &maxForceMerged = patch->flags.maxForceMerged;
255  maxForceUsed = Results::normal;
256  maxForceMerged = Results::normal;
257 
258  const int numberOfSteps = simParams->N;
259  const int stepsPerCycle = simParams->stepsPerCycle;
260  const BigReal timestep = simParams->dt;
261 
262  // what MTS method?
263  const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
264 
265  const int nonbondedFrequency = simParams->nonbondedFrequency;
266  slowFreq = nonbondedFrequency;
267  const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
268  int &doNonbonded = patch->flags.doNonbonded;
269  doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
270  if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
271  if ( doNonbonded ) maxForceUsed = Results::nbond;
272 
273  // Do we do full electrostatics?
274  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
275  const int fullElectFrequency = simParams->fullElectFrequency;
276  if ( dofull ) slowFreq = fullElectFrequency;
277  const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
278  int &doFullElectrostatics = patch->flags.doFullElectrostatics;
279  doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
280  if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
281  maxForceMerged = Results::slow;
282  if ( doFullElectrostatics ) maxForceUsed = Results::slow;
283 
284 //#ifndef UPPER_BOUND
285  const Bool accelMDOn = simParams->accelMDOn;
286  const Bool accelMDdihe = simParams->accelMDdihe;
287  const Bool accelMDdual = simParams->accelMDdual;
288  if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
289 
290  // Is adaptive tempering on?
291  const Bool adaptTempOn = simParams->adaptTempOn;
293  if (simParams->langevinOn)
295  else if (simParams->rescaleFreq > 0)
297 
298 
299  int &doMolly = patch->flags.doMolly;
300  doMolly = simParams->mollyOn && doFullElectrostatics;
301  // BEGIN LA
302  int &doLoweAndersen = patch->flags.doLoweAndersen;
303  doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
304  // END LA
305 
306  int &doGBIS = patch->flags.doGBIS;
307  doGBIS = simParams->GBISOn;
308 
309  int &doLCPO = patch->flags.doLCPO;
310  doLCPO = simParams->LCPOOn;
311 
312  int zeroMomentum = simParams->zeroMomentum;
313 
314  // Do we need to return forces to TCL script or Colvar module?
315  int doTcl = simParams->tclForcesOn;
316  int doColvars = simParams->colvarsOn;
317 //#endif
319 
320  // Bother to calculate energies?
321  int &doEnergy = patch->flags.doEnergy;
322  int energyFrequency = simParams->outputEnergies;
323 #ifndef UPPER_BOUND
324  const int reassignFreq = simParams->reassignFreq;
325 #endif
326 
327  int &doVirial = patch->flags.doVirial;
328  doVirial = 1;
329 
330  if ( scriptTask == SCRIPT_RUN ) {
331 
332 #ifndef UPPER_BOUND
333 // printf("Doing initial rattle\n");
334 #ifndef UPPER_BOUND
335 D_MSG("rattle1()");
336  TIMER_START(t, RATTLE1);
337  rattle1(0.,0); // enforce rigid bond constraints on initial positions
338  TIMER_STOP(t, RATTLE1);
339 #endif
340 
343  patch->atom.begin(),patch->atom.end());
344  }
345 
346  if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
347  reassignVelocities(timestep,step);
348  }
349 #endif
350 
351  doEnergy = ! ( step % energyFrequency );
352 #ifndef UPPER_BOUND
353  if ( accelMDOn && !accelMDdihe ) doEnergy=1;
354  //Update energy every timestep for adaptive tempering
355  if ( adaptTempOn ) doEnergy=1;
356 #endif
357 D_MSG("runComputeObjects()");
358  runComputeObjects(1,step<numberOfSteps); // must migrate here!
359 #ifndef UPPER_BOUND
360  rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
361  adaptTempUpdate(step); // update adaptive tempering temperature
362 #endif
363 
364 #ifndef UPPER_BOUND
365  if ( staleForces || doTcl || doColvars ) {
366  if ( doNonbonded ) saveForce(Results::nbond);
367  if ( doFullElectrostatics ) saveForce(Results::slow);
368  }
369  if ( ! commOnly ) {
370 D_MSG("newtonianVelocities()");
371  TIMER_START(t, KICK);
372  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
373  TIMER_STOP(t, KICK);
374  }
376 #ifndef UPPER_BOUND
377 D_MSG("rattle1()");
378  TIMER_START(t, RATTLE1);
379  rattle1(-timestep,0);
380  TIMER_STOP(t, RATTLE1);
381 #endif
382 D_MSG("submitHalfstep()");
383  TIMER_START(t, SUBMITHALF);
384  submitHalfstep(step);
385  TIMER_STOP(t, SUBMITHALF);
386  if ( ! commOnly ) {
387 D_MSG("newtonianVelocities()");
388  TIMER_START(t, KICK);
389  newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
390  TIMER_STOP(t, KICK);
391  }
392 D_MSG("rattle1()");
393  TIMER_START(t, RATTLE1);
394  rattle1(timestep,1);
395  TIMER_STOP(t, RATTLE1);
396  if (doTcl || doColvars) // include constraint forces
397  computeGlobal->saveTotalForces(patch);
398 D_MSG("submitHalfstep()");
399  TIMER_START(t, SUBMITHALF);
400  submitHalfstep(step);
401  TIMER_STOP(t, SUBMITHALF);
402  if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
403  if ( ! commOnly ) {
404 D_MSG("newtonianVelocities()");
405  TIMER_START(t, KICK);
406  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
407  TIMER_STOP(t, KICK);
408  }
409 #endif
410 D_MSG("submitReductions()");
411  TIMER_START(t, SUBMITFULL);
412  submitReductions(step);
413  TIMER_STOP(t, SUBMITFULL);
414 #ifndef UPPER_BOUND
415  if(0){ // if(traceIsOn()){
416  traceUserEvent(eventEndOfTimeStep);
417  sprintf(traceNote, "%s%d",tracePrefix,step);
418  traceUserSuppliedNote(traceNote);
419  }
420 #endif
421  rebalanceLoad(step);
422 
423  } // scriptTask == SCRIPT_RUN
424 
425 #ifndef UPPER_BOUND
426  bool doMultigratorRattle = false;
427 #endif
428 
429  //
430  // DJH: There are a lot of mod operations below and elsewhere to
431  // test step number against the frequency of something happening.
432  // Mod and integer division are expensive!
433  // Might be better to replace with counters and test equality.
434  //
435 #if 0
436  for(int i = 0; i < NamdProfileEvent::EventsCount; i++)
437  CkPrintf("-------------- [%d] %s -------------\n", i, NamdProfileEventStr[i]);
438 #endif
439 
440 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
441  int& eon = patch->flags.event_on;
442  int epid = (simParams->beginEventPatchID <= patch->getPatchID()
443  && patch->getPatchID() <= simParams->endEventPatchID);
444  int beginStep = simParams->beginEventStep;
445  int endStep = simParams->endEventStep;
446  bool controlProfiling = patch->getPatchID() == 0;
447 #endif
448 
449  for ( ++step; step <= numberOfSteps; ++step )
450  {
451 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
452  eon = epid && (beginStep < step && step <= endStep);
453 
454  if (controlProfiling && step == beginStep) {
456  }
457  if (controlProfiling && step == endStep) {
459  }
460  char buf[32];
461  sprintf(buf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::INTEGRATE_1], patch->getPatchID());
462  NAMD_EVENT_START_EX(eon, NamdProfileEvent::INTEGRATE_1, buf);
463 #endif
464 #ifndef UPPER_BOUND
465 
466  rescaleVelocities(step);
467  tcoupleVelocities(timestep,step);
468  stochRescaleVelocities(timestep,step);
469  berendsenPressure(step);
470 
471  if ( ! commOnly ) {
472  TIMER_START(t, KICK);
473  newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
474  TIMER_STOP(t, KICK);
475  }
476 
477  // We do RATTLE here if multigrator thermostat was applied in the previous step
478  if (doMultigratorRattle) rattle1(timestep, 1);
479 
480  /* reassignment based on half-step velocities
481  if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
482  addVelocityToPosition(0.5*timestep);
483  reassignVelocities(timestep,step);
484  addVelocityToPosition(0.5*timestep);
485  rattle1(0.,0);
486  rattle1(-timestep,0);
487  addVelocityToPosition(-1.0*timestep);
488  rattle1(timestep,0);
489  } */
490 
491  TIMER_START(t, MAXMOVE);
492  maximumMove(timestep);
493  TIMER_STOP(t, MAXMOVE);
494 
495  NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_1); // integrate 1
496 
498  if ( ! commOnly ) {
499  TIMER_START(t, DRIFT);
500  addVelocityToPosition(0.5*timestep);
501  TIMER_STOP(t, DRIFT);
502  }
503  // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
504  langevinVelocities(timestep);
505 
506  // There is a blocking receive inside of langevinPiston()
507  // that might suspend the current thread of execution,
508  // so split profiling around this conditional block.
509  langevinPiston(step);
510 
511  if ( ! commOnly ) {
512  TIMER_START(t, DRIFT);
513  addVelocityToPosition(0.5*timestep);
514  TIMER_STOP(t, DRIFT);
515  }
516  } else {
517  // If Langevin is not used, take full time step directly instread of two half steps
518  if ( ! commOnly ) {
519  TIMER_START(t, DRIFT);
520  addVelocityToPosition(timestep);
521  TIMER_STOP(t, DRIFT);
522  }
523  }
524 
525  NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_2);
526 
527  // impose hard wall potential for Drude bond length
528  hardWallDrude(timestep, 1);
529 
531 #endif // UPPER_BOUND
532 
533  doNonbonded = !(step%nonbondedFrequency);
534  doFullElectrostatics = (dofull && !(step%fullElectFrequency));
535 
536 #ifndef UPPER_BOUND
537  if ( zeroMomentum && doFullElectrostatics ) {
538  // There is a blocking receive inside of correctMomentum().
539  correctMomentum(step,slowstep);
540  }
541 
542  // There are NO sends in submitHalfstep() just local summation
543  // into the Reduction struct.
544  TIMER_START(t, SUBMITHALF);
545  submitHalfstep(step);
546  TIMER_STOP(t, SUBMITHALF);
547 
548  doMolly = simParams->mollyOn && doFullElectrostatics;
549  // BEGIN LA
550  doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
551  // END LA
552 
553  maxForceUsed = Results::normal;
554  if ( doNonbonded ) maxForceUsed = Results::nbond;
555  if ( doFullElectrostatics ) maxForceUsed = Results::slow;
556  if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
557 
558  // Migrate Atoms on stepsPerCycle
559  doEnergy = ! ( step % energyFrequency );
560  if ( accelMDOn && !accelMDdihe ) doEnergy=1;
561  if ( adaptTempOn ) doEnergy=1;
562 
563  // Multigrator
564  if (simParams->multigratorOn) {
565  doVirial = (!(step % energyFrequency) || ((simParams->outputPressure > 0) && !(step % simParams->outputPressure))
566  || !(step % simParams->multigratorPressureFreq));
567  doKineticEnergy = (!(step % energyFrequency) || !(step % simParams->multigratorTemperatureFreq));
569  } else {
570  doVirial = 1;
571  doKineticEnergy = 1;
572  doMomenta = 1;
573  }
574 #endif
575  NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_2); // integrate 2
576 
577  // The current thread of execution will suspend in runComputeObjects().
578  runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
579 
580  NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_3);
581 
582 #ifndef UPPER_BOUND
583  rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
584 
585  if ( staleForces || doTcl || doColvars ) {
586  if ( doNonbonded ) saveForce(Results::nbond);
587  if ( doFullElectrostatics ) saveForce(Results::slow);
588  }
589 
590  // reassignment based on full-step velocities
591  if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
592  reassignVelocities(timestep,step);
593  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
594  rattle1(-timestep,0);
595  }
596 
597  if ( ! commOnly ) {
598  TIMER_START(t, VELBBK1);
599  langevinVelocitiesBBK1(timestep);
600  TIMER_STOP(t, VELBBK1);
601  TIMER_START(t, KICK);
602  newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
603  TIMER_STOP(t, KICK);
604  TIMER_START(t, VELBBK2);
605  langevinVelocitiesBBK2(timestep);
606  TIMER_STOP(t, VELBBK2);
607  }
608 
609  // add drag to each atom's positions
610  if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
611  if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
612 
613  TIMER_START(t, RATTLE1);
614  rattle1(timestep,1);
615  TIMER_STOP(t, RATTLE1);
616  if (doTcl || doColvars) // include constraint forces
617  computeGlobal->saveTotalForces(patch);
618 
619  TIMER_START(t, SUBMITHALF);
620  submitHalfstep(step);
621  TIMER_STOP(t, SUBMITHALF);
622  if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
623 
624  if ( ! commOnly ) {
625  TIMER_START(t, KICK);
626  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
627  TIMER_STOP(t, KICK);
628  }
629 
630  // rattle2(timestep,step);
631 #endif
632 
633  TIMER_START(t, SUBMITFULL);
634  submitReductions(step);
635  TIMER_STOP(t, SUBMITFULL);
636  TIMER_START(t, SUBMITCOLLECT);
637  submitCollections(step);
638  TIMER_STOP(t, SUBMITCOLLECT);
639 #ifndef UPPER_BOUND
640  //Update adaptive tempering temperature
641  adaptTempUpdate(step);
642 
643  // Multigrator temperature and pressure steps
644  multigratorTemperature(step, 1);
645  multigratorPressure(step, 1);
646  multigratorPressure(step, 2);
647  multigratorTemperature(step, 2);
648  doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
649 
650  NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_3); // integrate 3
651 #endif
652 
653 #if CYCLE_BARRIER
654  cycleBarrier(!((step+1) % stepsPerCycle), step);
655 #elif PME_BARRIER
656  cycleBarrier(doFullElectrostatics, step);
657 #elif STEP_BARRIER
658  cycleBarrier(1, step);
659 #endif
660 
661 #ifndef UPPER_BOUND
662  if(Node::Object()->specialTracing || simParams->statsOn){
663  int bstep = simParams->traceStartStep;
664  int estep = bstep + simParams->numTraceSteps;
665  if(step == bstep || step == estep){
666  traceBarrier(step);
667  }
668  }
669 
670 #ifdef MEASURE_NAMD_WITH_PAPI
671  if(simParams->papiMeasure) {
672  int bstep = simParams->papiMeasureStartStep;
673  int estep = bstep + simParams->numPapiMeasureSteps;
674  if(step == bstep || step==estep) {
675  papiMeasureBarrier(step);
676  }
677  }
678 #endif
679 
680  if(0){ // if(traceIsOn()){
681  traceUserEvent(eventEndOfTimeStep);
682  sprintf(traceNote, "%s%d",tracePrefix,step);
683  traceUserSuppliedNote(traceNote);
684  }
685 #endif // UPPER_BOUND
686  rebalanceLoad(step);
687 
688 #if PME_BARRIER
689  // a step before PME
690  cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
691 #endif
692 
693 #if USE_HPM
694  if(step == START_HPM_STEP)
695  (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
696 
697  if(step == STOP_HPM_STEP)
698  (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
699 #endif
700 
701  }
702 
703  TIMER_DONE(t);
704 #ifdef TIMER_COLLECTION
705  if (patch->patchID == SPECIAL_PATCH_ID) {
706  printf("Timer collection reporting in microseconds for "
707  "Patch %d\n", patch->patchID);
708  TIMER_REPORT(t);
709  }
710 #endif // TIMER_COLLECTION
711  //
712  // DJH: Copy updates of SOA back into AOS.
713  //
714  //patch->copy_updates_to_AOS();
715 }
716 
717 // add moving drag to each atom's position
719  FullAtom *atom = patch->atom.begin();
720  int numAtoms = patch->numAtoms;
721  Molecule *molecule = Node::Object()->molecule; // need its methods
722  const BigReal movDragGlobVel = simParams->movDragGlobVel;
723  const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
724  Vector movDragVel, dragIncrement;
725  for ( int i = 0; i < numAtoms; ++i )
726  {
727  // skip if fixed atom or zero drag attribute
728  if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
729  || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
730  molecule->get_movdrag_params(movDragVel, atom[i].id);
731  dragIncrement = movDragGlobVel * movDragVel * dt;
732  atom[i].position += dragIncrement;
733  }
734 }
735 
736 // add rotating drag to each atom's position
738  FullAtom *atom = patch->atom.begin();
739  int numAtoms = patch->numAtoms;
740  Molecule *molecule = Node::Object()->molecule; // need its methods
741  const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
742  const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
743  BigReal rotDragVel, dAngle;
744  Vector atomRadius;
745  Vector rotDragAxis, rotDragPivot, dragIncrement;
746  for ( int i = 0; i < numAtoms; ++i )
747  {
748  // skip if fixed atom or zero drag attribute
749  if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
750  || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
751  molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
752  dAngle = rotDragGlobVel * rotDragVel * dt;
753  rotDragAxis /= rotDragAxis.length();
754  atomRadius = atom[i].position - rotDragPivot;
755  dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
756  atom[i].position += dragIncrement;
757  }
758 }
759 
761  //
762  // DJH: Copy all data into SOA (structure of arrays)
763  // from AOS (array of structures) data structure.
764  //
765  //patch->copy_all_to_SOA();
766 
767  const int numberOfSteps = simParams->N;
768  const int stepsPerCycle = simParams->stepsPerCycle;
769  int &step = patch->flags.step;
770  step = simParams->firstTimestep;
771 
772  int &maxForceUsed = patch->flags.maxForceUsed;
773  int &maxForceMerged = patch->flags.maxForceMerged;
774  maxForceUsed = Results::normal;
775  maxForceMerged = Results::normal;
776  int &doNonbonded = patch->flags.doNonbonded;
777  doNonbonded = 1;
778  maxForceUsed = Results::nbond;
779  maxForceMerged = Results::nbond;
780  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
781  int &doFullElectrostatics = patch->flags.doFullElectrostatics;
782  doFullElectrostatics = dofull;
783  if ( dofull ) {
784  maxForceMerged = Results::slow;
785  maxForceUsed = Results::slow;
786  }
787  int &doMolly = patch->flags.doMolly;
788  doMolly = simParams->mollyOn && doFullElectrostatics;
789  // BEGIN LA
790  int &doLoweAndersen = patch->flags.doLoweAndersen;
791  doLoweAndersen = 0;
792  // END LA
793 
794  int &doGBIS = patch->flags.doGBIS;
795  doGBIS = simParams->GBISOn;
796 
797  int &doLCPO = patch->flags.doLCPO;
798  doLCPO = simParams->LCPOOn;
799 
800  int doTcl = simParams->tclForcesOn;
801  int doColvars = simParams->colvarsOn;
803 
804  int &doEnergy = patch->flags.doEnergy;
805  doEnergy = 1;
806 
807  patch->rattle1(0.,0,0); // enforce rigid bond constraints on initial positions
808 
811  patch->atom.begin(),patch->atom.end());
812  }
813 
814  runComputeObjects(1,step<numberOfSteps); // must migrate here!
815 
816  if ( doTcl || doColvars ) {
817  if ( doNonbonded ) saveForce(Results::nbond);
818  if ( doFullElectrostatics ) saveForce(Results::slow);
819  computeGlobal->saveTotalForces(patch);
820  }
821 
823 
824  submitMinimizeReductions(step,fmax2);
825  rebalanceLoad(step);
826 
827  int downhill = 1; // start out just fixing bad contacts
828  int minSeq = 0;
829  for ( ++step; step <= numberOfSteps; ++step ) {
830 
831  // Blocking receive for the minimization coefficient.
832  BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
833 
834  if ( downhill ) {
835  if ( c ) minimizeMoveDownhill(fmax2);
836  else {
837  downhill = 0;
838  fmax2 *= 10000.;
839  }
840  }
841  if ( ! downhill ) {
842  if ( ! c ) { // new direction
843 
844  // Blocking receive for the minimization coefficient.
845  c = broadcast->minimizeCoefficient.get(minSeq++);
846 
847  newMinimizeDirection(c); // v = c * v + f
848 
849  // Blocking receive for the minimization coefficient.
850  c = broadcast->minimizeCoefficient.get(minSeq++);
851 
852  } // same direction
853  newMinimizePosition(c); // x = x + c * v
854  }
855 
856  runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
857  if ( doTcl || doColvars ) {
858  if ( doNonbonded ) saveForce(Results::nbond);
859  if ( doFullElectrostatics ) saveForce(Results::slow);
860  computeGlobal->saveTotalForces(patch);
861  }
862  submitMinimizeReductions(step,fmax2);
863  submitCollections(step, 1); // write out zeros for velocities
864  rebalanceLoad(step);
865  }
866  quenchVelocities(); // zero out bogus velocity
867 
868  //
869  // DJH: Copy updates of SOA back into AOS.
870  //
871  //patch->copy_updates_to_AOS();
872 }
873 
874 // x = x + 0.1 * unit(f) for large f
876 
877  FullAtom *a = patch->atom.begin();
878  Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
879  int numAtoms = patch->numAtoms;
880 
881  for ( int i = 0; i < numAtoms; ++i ) {
882  if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
883  Force f = f1[i];
884  if ( f.length2() > fmax2 ) {
885  a[i].position += ( 0.1 * f.unit() );
886  int hgs = a[i].hydrogenGroupSize; // 0 if not parent
887  for ( int j=1; j<hgs; ++j ) {
888  a[++i].position += ( 0.1 * f.unit() );
889  }
890  }
891  }
892 
893  patch->rattle1(0.,0,0);
894 }
895 
896 // v = c * v + f
898  FullAtom *a = patch->atom.begin();
899  Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
900  const bool fixedAtomsOn = simParams->fixedAtomsOn;
901  const bool drudeHardWallOn = simParams->drudeHardWallOn;
902  int numAtoms = patch->numAtoms;
903  BigReal maxv2 = 0.;
904 
905  for ( int i = 0; i < numAtoms; ++i ) {
906  a[i].velocity *= c;
907  a[i].velocity += f1[i];
908  if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
909  a[i].velocity = a[i-1].velocity;
910  }
911  if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
912  BigReal v2 = a[i].velocity.length2();
913  if ( v2 > maxv2 ) maxv2 = v2;
914  }
915 
916  { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial); }
917 
918  maxv2 = 0.;
919  for ( int i = 0; i < numAtoms; ++i ) {
920  if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
921  a[i].velocity = a[i-1].velocity;
922  }
923  if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
924  BigReal v2 = a[i].velocity.length2();
925  if ( v2 > maxv2 ) maxv2 = v2;
926  }
927 
928  min_reduction->max(0,maxv2);
930 
931  // prevent hydrogens from being left behind
932  BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
933  // int adjustCount = 0;
934  int hgs;
935  for ( int i = 0; i < numAtoms; i += hgs ) {
936  hgs = a[i].hydrogenGroupSize;
937  BigReal minChildVel = a[i].velocity.length2();
938  if ( minChildVel < fmax2 ) continue;
939  int adjustChildren = 1;
940  for ( int j = i+1; j < (i+hgs); ++j ) {
941  if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
942  }
943  if ( adjustChildren ) {
944  // if ( hgs > 1 ) ++adjustCount;
945  for ( int j = i+1; j < (i+hgs); ++j ) {
946  if (a[i].mass < 0.01) continue; // lone pair
947  a[j].velocity = a[i].velocity;
948  }
949  }
950  }
951  // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
952 
953 }
954 
955 // x = x + c * v
957  FullAtom *a = patch->atom.begin();
958  int numAtoms = patch->numAtoms;
959 
960  for ( int i = 0; i < numAtoms; ++i ) {
961  a[i].position += c * a[i].velocity;
962  }
963 
964  if ( simParams->drudeHardWallOn ) {
965  for ( int i = 1; i < numAtoms; ++i ) {
966  if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
967  a[i].position -= a[i-1].position;
968  }
969  }
970  }
971 
972  patch->rattle1(0.,0,0);
973 
974  if ( simParams->drudeHardWallOn ) {
975  for ( int i = 1; i < numAtoms; ++i ) {
976  if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
977  a[i].position += a[i-1].position;
978  }
979  }
980  }
981 }
982 
983 // v = 0
985  FullAtom *a = patch->atom.begin();
986  int numAtoms = patch->numAtoms;
987 
988  for ( int i = 0; i < numAtoms; ++i ) {
989  a[i].velocity = 0;
990  }
991 }
992 
994 
995  FullAtom *a = patch->atom.begin();
996  const int numAtoms = patch->numAtoms;
997 
998  Vector momentum = 0;
999  BigReal mass = 0;
1000 if ( simParams->zeroMomentumAlt ) {
1001  for ( int i = 0; i < numAtoms; ++i ) {
1002  momentum += a[i].mass * a[i].velocity;
1003  mass += 1.;
1004  }
1005 } else {
1006  for ( int i = 0; i < numAtoms; ++i ) {
1007  momentum += a[i].mass * a[i].velocity;
1008  mass += a[i].mass;
1009  }
1010 }
1011 
1012  ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
1014 }
1015 
1017 
1018  //
1019  // DJH: This test should be done in SimParameters.
1020  //
1021  if ( simParams->fixedAtomsOn )
1022  NAMD_die("Cannot zero momentum when fixed atoms are present.");
1023 
1024  // Blocking receive for the momentum correction vector.
1025  const Vector dv = broadcast->momentumCorrection.get(step);
1026 
1027  const Vector dx = dv * ( drifttime / TIMEFACTOR );
1028 
1029  FullAtom *a = patch->atom.begin();
1030  const int numAtoms = patch->numAtoms;
1031 
1032 if ( simParams->zeroMomentumAlt ) {
1033  for ( int i = 0; i < numAtoms; ++i ) {
1034  a[i].velocity += dv * a[i].recipMass;
1035  a[i].position += dx * a[i].recipMass;
1036  }
1037 } else {
1038  for ( int i = 0; i < numAtoms; ++i ) {
1039  a[i].velocity += dv;
1040  a[i].position += dx;
1041  }
1042 }
1043 
1044 }
1045 
1046 // --------- For Multigrator ---------
1047 void Sequencer::scalePositionsVelocities(const Tensor& posScale, const Tensor& velScale) {
1048  FullAtom *a = patch->atom.begin();
1049  int numAtoms = patch->numAtoms;
1050  Position origin = patch->lattice.origin();
1051  if ( simParams->fixedAtomsOn ) {
1052  NAMD_bug("Sequencer::scalePositionsVelocities, fixed atoms not implemented");
1053  }
1054  if ( simParams->useGroupPressure ) {
1055  int hgs;
1056  for ( int i = 0; i < numAtoms; i += hgs ) {
1057  hgs = a[i].hydrogenGroupSize;
1058  Position pos_cm(0.0, 0.0, 0.0);
1059  Velocity vel_cm(0.0, 0.0, 0.0);
1060  BigReal m_cm = 0.0;
1061  for (int j=0;j < hgs;++j) {
1062  m_cm += a[i+j].mass;
1063  pos_cm += a[i+j].mass*a[i+j].position;
1064  vel_cm += a[i+j].mass*a[i+j].velocity;
1065  }
1066  pos_cm /= m_cm;
1067  vel_cm /= m_cm;
1068  pos_cm -= origin;
1069  Position dpos = posScale*pos_cm;
1070  Velocity dvel = velScale*vel_cm;
1071  for (int j=0;j < hgs;++j) {
1072  a[i+j].position += dpos;
1073  a[i+j].velocity += dvel;
1074  }
1075  }
1076  } else {
1077  for ( int i = 0; i < numAtoms; i++) {
1078  a[i].position += posScale*(a[i].position-origin);
1079  a[i].velocity = velScale*a[i].velocity;
1080  }
1081  }
1082 }
1083 
1084 void Sequencer::multigratorPressure(int step, int callNumber) {
1085 // Calculate new positions, momenta, and volume using positionRescaleFactor and
1086 // velocityRescaleTensor values returned from Controller::multigratorPressureCalcScale()
1088  FullAtom *a = patch->atom.begin();
1089  int numAtoms = patch->numAtoms;
1090 
1091  // Blocking receive (get) scaling factors from Controller
1092  Tensor scaleTensor = (callNumber == 1) ? broadcast->positionRescaleFactor.get(step) : broadcast->positionRescaleFactor2.get(step);
1093  Tensor velScaleTensor = (callNumber == 1) ? broadcast->velocityRescaleTensor.get(step) : broadcast->velocityRescaleTensor2.get(step);
1094  Tensor posScaleTensor = scaleTensor;
1095  posScaleTensor -= Tensor::identity();
1096  if (simParams->useGroupPressure) {
1097  velScaleTensor -= Tensor::identity();
1098  }
1099 
1100  // Scale volume
1101  patch->lattice.rescale(scaleTensor);
1102  // Scale positions and velocities
1103  scalePositionsVelocities(posScaleTensor, velScaleTensor);
1104 
1105  if (!patch->flags.doFullElectrostatics) NAMD_bug("Sequencer::multigratorPressure, doFullElectrostatics must be true");
1106 
1107  // Calculate new forces
1108  // NOTE: We should not need to migrate here since any migration should have happened in the
1109  // previous call to runComputeObjects inside the MD loop in Sequencer::integrate()
1110  const int numberOfSteps = simParams->N;
1111  const int stepsPerCycle = simParams->stepsPerCycle;
1112  runComputeObjects(0 , step<numberOfSteps, 1);
1113 
1114  reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
1116 
1117  // Virials etc.
1118  Tensor virialNormal;
1119  Tensor momentumSqrSum;
1120  BigReal kineticEnergy = 0;
1121  if ( simParams->pairInteractionOn ) {
1122  if ( simParams->pairInteractionSelf ) {
1123  for ( int i = 0; i < numAtoms; ++i ) {
1124  if ( a[i].partition != 1 ) continue;
1125  kineticEnergy += a[i].mass * a[i].velocity.length2();
1126  virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1127  }
1128  }
1129  } else {
1130  for ( int i = 0; i < numAtoms; ++i ) {
1131  if (a[i].mass < 0.01) continue;
1132  kineticEnergy += a[i].mass * a[i].velocity.length2();
1133  virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1134  }
1135  }
1136  if (!simParams->useGroupPressure) momentumSqrSum = virialNormal;
1137  kineticEnergy *= 0.5;
1139  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virialNormal);
1140 
1141  if ( simParams->fixedAtomsOn ) {
1142  Tensor fixVirialNormal;
1143  Tensor fixVirialNbond;
1144  Tensor fixVirialSlow;
1145  Vector fixForceNormal = 0;
1146  Vector fixForceNbond = 0;
1147  Vector fixForceSlow = 0;
1148 
1149  calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
1150 
1151  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, fixVirialNormal);
1152  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, fixVirialNbond);
1153  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, fixVirialSlow);
1154  ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NORMAL, fixForceNormal);
1155  ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NBOND, fixForceNbond);
1156  ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_SLOW, fixForceSlow);
1157  }
1158 
1159  // Internal virial and group momentum
1160  Tensor intVirialNormal;
1161  Tensor intVirialNormal2;
1162  Tensor intVirialNbond;
1163  Tensor intVirialSlow;
1164  int hgs;
1165  for ( int i = 0; i < numAtoms; i += hgs ) {
1166  hgs = a[i].hydrogenGroupSize;
1167  int j;
1168  BigReal m_cm = 0;
1169  Position x_cm(0,0,0);
1170  Velocity v_cm(0,0,0);
1171  for ( j = i; j < (i+hgs); ++j ) {
1172  m_cm += a[j].mass;
1173  x_cm += a[j].mass * a[j].position;
1174  v_cm += a[j].mass * a[j].velocity;
1175  }
1176  if (simParams->useGroupPressure) momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
1177  x_cm /= m_cm;
1178  v_cm /= m_cm;
1179  if (simParams->fixedAtomsOn) NAMD_bug("Sequencer::multigratorPressure, simParams->fixedAtomsOn not implemented yet");
1180  if ( simParams->pairInteractionOn ) {
1181  if ( simParams->pairInteractionSelf ) {
1182  NAMD_bug("Sequencer::multigratorPressure, this part needs to be implemented correctly");
1183  for ( j = i; j < (i+hgs); ++j ) {
1184  if ( a[j].partition != 1 ) continue;
1185  BigReal mass = a[j].mass;
1186  Vector v = a[j].velocity;
1187  Vector dv = v - v_cm;
1188  intVirialNormal2.outerAdd (mass, v, dv);
1189  Vector dx = a[j].position - x_cm;
1190  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
1191  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
1192  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
1193  }
1194  }
1195  } else {
1196  for ( j = i; j < (i+hgs); ++j ) {
1197  BigReal mass = a[j].mass;
1198  Vector v = a[j].velocity;
1199  Vector dv = v - v_cm;
1200  intVirialNormal2.outerAdd(mass, v, dv);
1201  Vector dx = a[j].position - x_cm;
1202  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
1203  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
1204  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
1205  }
1206  }
1207  }
1208 
1209  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
1210  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal2);
1211  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NBOND, intVirialNbond);
1212  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_SLOW, intVirialSlow);
1213  ADD_TENSOR_OBJECT(reduction, REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
1214 
1215  reduction->submit();
1216  }
1217 }
1218 
1219 void Sequencer::scaleVelocities(const BigReal velScale) {
1220  FullAtom *a = patch->atom.begin();
1221  int numAtoms = patch->numAtoms;
1222  for ( int i = 0; i < numAtoms; i++) {
1223  a[i].velocity *= velScale;
1224  }
1225 }
1226 
1228  FullAtom *a = patch->atom.begin();
1229  int numAtoms = patch->numAtoms;
1230  BigReal kineticEnergy = 0.0;
1231  if ( simParams->pairInteractionOn ) {
1232  if ( simParams->pairInteractionSelf ) {
1233  for (int i = 0; i < numAtoms; ++i ) {
1234  if ( a[i].partition != 1 ) continue;
1235  kineticEnergy += a[i].mass * a[i].velocity.length2();
1236  }
1237  }
1238  } else {
1239  for (int i = 0; i < numAtoms; ++i ) {
1240  kineticEnergy += a[i].mass * a[i].velocity.length2();
1241  }
1242  }
1243  kineticEnergy *= 0.5;
1244  return kineticEnergy;
1245 }
1246 
1247 void Sequencer::multigratorTemperature(int step, int callNumber) {
1249  // Blocking receive (get) velocity scaling factor.
1250  BigReal velScale = (callNumber == 1) ? broadcast->velocityRescaleFactor.get(step) : broadcast->velocityRescaleFactor2.get(step);
1251  scaleVelocities(velScale);
1252  // Calculate new kineticEnergy
1253  BigReal kineticEnergy = calcKineticEnergy();
1255  if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
1256  // If this is a pressure cycle, calculate new momentum squared sum
1257  FullAtom *a = patch->atom.begin();
1258  int numAtoms = patch->numAtoms;
1259  Tensor momentumSqrSum;
1260  if (simParams->useGroupPressure) {
1261  int hgs;
1262  for ( int i = 0; i < numAtoms; i += hgs ) {
1263  hgs = a[i].hydrogenGroupSize;
1264  int j;
1265  BigReal m_cm = 0;
1266  Position x_cm(0,0,0);
1267  Velocity v_cm(0,0,0);
1268  for ( j = i; j < (i+hgs); ++j ) {
1269  m_cm += a[j].mass;
1270  x_cm += a[j].mass * a[j].position;
1271  v_cm += a[j].mass * a[j].velocity;
1272  }
1273  momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
1274  }
1275  } else {
1276  for ( int i = 0; i < numAtoms; i++) {
1277  momentumSqrSum.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1278  }
1279  }
1280  ADD_TENSOR_OBJECT(multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
1281  }
1282  // Submit reductions (kineticEnergy and, if applicable, momentumSqrSum)
1284 
1285  }
1286 }
1287 // --------- End Multigrator ---------
1288 
1289 //
1290 // DJH: Calls one or more addForceToMomentum which in turn calls HomePatch
1291 // versions. We should inline to reduce the number of function calls.
1292 //
1293 void Sequencer::newtonianVelocities(BigReal stepscale, const BigReal timestep,
1294  const BigReal nbondstep,
1295  const BigReal slowstep,
1296  const int staleForces,
1297  const int doNonbonded,
1298  const int doFullElectrostatics)
1299 {
1300  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1301  NamdProfileEvent::NEWTONIAN_VELOCITIES);
1302 
1303  // Deterministic velocity update, account for multigrator
1304  if (staleForces || (doNonbonded && doFullElectrostatics)) {
1305  addForceToMomentum3(stepscale*timestep, Results::normal, 0,
1306  stepscale*nbondstep, Results::nbond, staleForces,
1307  stepscale*slowstep, Results::slow, staleForces);
1308  } else {
1309  addForceToMomentum(stepscale*timestep);
1310  if (staleForces || doNonbonded)
1311  addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces);
1312  if (staleForces || doFullElectrostatics)
1313  addForceToMomentum(stepscale*slowstep, Results::slow, staleForces);
1314  }
1315 }
1316 
1318 {
1319 // This routine is used for the BAOAB integrator,
1320 // Ornstein-Uhlenbeck exact solve for the O-part.
1321 // See B. Leimkuhler and C. Matthews, AMRX (2012)
1322 // Routine originally written by JPhillips, with fresh errors by CMatthews June2012
1323 
1325  {
1326  FullAtom *a = patch->atom.begin();
1327  int numAtoms = patch->numAtoms;
1328  Molecule *molecule = Node::Object()->molecule;
1329  BigReal dt = dt_fs * 0.001; // convert to ps
1332  {
1333  kbT = BOLTZMANN*adaptTempT;
1334  }
1335 
1336  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1337  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1338 
1339  for ( int i = 0; i < numAtoms; ++i )
1340  {
1341  BigReal dt_gamma = dt * a[i].langevinParam;
1342  if ( ! dt_gamma ) continue;
1343 
1344  BigReal f1 = exp( -dt_gamma );
1345  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
1346  ( a[i].partition ? tempFactor : 1.0 ) *
1347  a[i].recipMass );
1348  a[i].velocity *= f1;
1349  a[i].velocity += f2 * random->gaussian_vector();
1350  }
1351  }
1352 }
1353 
1355 {
1356  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1357  NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1);
1359  {
1360  FullAtom *a = patch->atom.begin();
1361  int numAtoms = patch->numAtoms;
1362  Molecule *molecule = Node::Object()->molecule;
1363  BigReal dt = dt_fs * 0.001; // convert to ps
1364  int i;
1365 
1366  if (simParams->drudeOn) {
1367  for (i = 0; i < numAtoms; i++) {
1368 
1369  if (i < numAtoms-1 &&
1370  a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
1371  //printf("*** Found Drude particle %d\n", a[i+1].id);
1372  // i+1 is a Drude particle with parent i
1373 
1374  // convert from Cartesian coordinates to (COM,bond) coordinates
1375  BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
1376  Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
1377  Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
1378  BigReal dt_gamma;
1379 
1380  // use Langevin damping factor i for v_com
1381  dt_gamma = dt * a[i].langevinParam;
1382  if (dt_gamma != 0.0) {
1383  v_com *= ( 1. - 0.5 * dt_gamma );
1384  }
1385 
1386  // use Langevin damping factor i+1 for v_bnd
1387  dt_gamma = dt * a[i+1].langevinParam;
1388  if (dt_gamma != 0.0) {
1389  v_bnd *= ( 1. - 0.5 * dt_gamma );
1390  }
1391 
1392  // convert back
1393  a[i].velocity = v_com - m * v_bnd;
1394  a[i+1].velocity = v_bnd + a[i].velocity;
1395 
1396  i++; // +1 from loop, we've updated both particles
1397  }
1398  else {
1399  BigReal dt_gamma = dt * a[i].langevinParam;
1400  if ( ! dt_gamma ) continue;
1401 
1402  a[i].velocity *= ( 1. - 0.5 * dt_gamma );
1403  }
1404 
1405  } // end for
1406  } // end if drudeOn
1407  else {
1408 
1409  //
1410  // DJH: The conditional inside loop prevents vectorization and doesn't
1411  // avoid much work since addition and multiplication are cheap.
1412  //
1413  for ( i = 0; i < numAtoms; ++i )
1414  {
1415  BigReal dt_gamma = dt * a[i].langevinParam;
1416  if ( ! dt_gamma ) continue;
1417 
1418  a[i].velocity *= ( 1. - 0.5 * dt_gamma );
1419  }
1420 
1421  } // end else
1422 
1423  } // end if langevinOn
1424 }
1425 
1426 
1428 {
1429  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1430  NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2);
1432  {
1433  //
1434  // DJH: This call is expensive. Avoid calling when gammas don't differ.
1435  // Set flag in SimParameters and make this call conditional.
1436  //
1437  TIMER_START(patch->timerSet, RATTLE1);
1438  rattle1(dt_fs,1); // conserve momentum if gammas differ
1439  TIMER_STOP(patch->timerSet, RATTLE1);
1440 
1441  FullAtom *a = patch->atom.begin();
1442  int numAtoms = patch->numAtoms;
1443  Molecule *molecule = Node::Object()->molecule;
1444  BigReal dt = dt_fs * 0.001; // convert to ps
1447  {
1448  kbT = BOLTZMANN*adaptTempT;
1449  }
1450  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1451  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1452  int i;
1453 
1454  if (simParams->drudeOn) {
1455  BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp); // drude bond Temp
1456 
1457  for (i = 0; i < numAtoms; i++) {
1458 
1459  if (i < numAtoms-1 &&
1460  a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
1461  //printf("*** Found Drude particle %d\n", a[i+1].id);
1462  // i+1 is a Drude particle with parent i
1463 
1464  // convert from Cartesian coordinates to (COM,bond) coordinates
1465  BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
1466  Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
1467  Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
1468  BigReal dt_gamma;
1469 
1470  // use Langevin damping factor i for v_com
1471  dt_gamma = dt * a[i].langevinParam;
1472  if (dt_gamma != 0.0) {
1473  BigReal mass = a[i].mass + a[i+1].mass;
1474  v_com += random->gaussian_vector() *
1475  sqrt( 2 * dt_gamma * kbT *
1476  ( a[i].partition ? tempFactor : 1.0 ) / mass );
1477  v_com /= ( 1. + 0.5 * dt_gamma );
1478  }
1479 
1480  // use Langevin damping factor i+1 for v_bnd
1481  dt_gamma = dt * a[i+1].langevinParam;
1482  if (dt_gamma != 0.0) {
1483  BigReal mass = a[i+1].mass * (1. - m);
1484  v_bnd += random->gaussian_vector() *
1485  sqrt( 2 * dt_gamma * kbT_bnd *
1486  ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
1487  v_bnd /= ( 1. + 0.5 * dt_gamma );
1488  }
1489 
1490  // convert back
1491  a[i].velocity = v_com - m * v_bnd;
1492  a[i+1].velocity = v_bnd + a[i].velocity;
1493 
1494  i++; // +1 from loop, we've updated both particles
1495  }
1496  else {
1497  BigReal dt_gamma = dt * a[i].langevinParam;
1498  if ( ! dt_gamma ) continue;
1499 
1500  a[i].velocity += random->gaussian_vector() *
1501  sqrt( 2 * dt_gamma * kbT *
1502  ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
1503  a[i].velocity /= ( 1. + 0.5 * dt_gamma );
1504  }
1505 
1506  } // end for
1507  } // end if drudeOn
1508  else {
1509 
1510  //
1511  // DJH: For case using same gamma (the Langevin parameter),
1512  // no partitions (e.g. FEP), and no adaptive tempering (adaptTempMD),
1513  // we can precompute constants. Then by lifting the RNG from the
1514  // loop (filling up an array of random numbers), we can vectorize
1515  // loop and simplify arithmetic to just addition and multiplication.
1516  //
1517  for ( i = 0; i < numAtoms; ++i )
1518  {
1519  BigReal dt_gamma = dt * a[i].langevinParam;
1520  if ( ! dt_gamma ) continue;
1521 
1522  a[i].velocity += random->gaussian_vector() *
1523  sqrt( 2 * dt_gamma * kbT *
1524  ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
1525  a[i].velocity /= ( 1. + 0.5 * dt_gamma );
1526  }
1527 
1528  } // end else
1529 
1530  } // end if langevinOn
1531 }
1532 
1533 
1535 {
1536  if ( simParams->berendsenPressureOn ) {
1538  const int freq = simParams->berendsenPressureFreq;
1539  if ( ! (berendsenPressure_count % freq ) ) {
1541  FullAtom *a = patch->atom.begin();
1542  int numAtoms = patch->numAtoms;
1543  // Blocking receive for the updated lattice scaling factor.
1544  Tensor factor = broadcast->positionRescaleFactor.get(step);
1545  patch->lattice.rescale(factor);
1546  if ( simParams->useGroupPressure )
1547  {
1548  int hgs;
1549  for ( int i = 0; i < numAtoms; i += hgs ) {
1550  int j;
1551  hgs = a[i].hydrogenGroupSize;
1552  if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
1553  for ( j = i; j < (i+hgs); ++j ) {
1555  a[j].fixedPosition,a[j].transform);
1556  }
1557  continue;
1558  }
1559  BigReal m_cm = 0;
1560  Position x_cm(0,0,0);
1561  for ( j = i; j < (i+hgs); ++j ) {
1562  if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
1563  m_cm += a[j].mass;
1564  x_cm += a[j].mass * a[j].position;
1565  }
1566  x_cm /= m_cm;
1567  Position new_x_cm = x_cm;
1568  patch->lattice.rescale(new_x_cm,factor);
1569  Position delta_x_cm = new_x_cm - x_cm;
1570  for ( j = i; j < (i+hgs); ++j ) {
1571  if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
1573  a[j].fixedPosition,a[j].transform);
1574  continue;
1575  }
1576  a[j].position += delta_x_cm;
1577  }
1578  }
1579  }
1580  else
1581  {
1582  for ( int i = 0; i < numAtoms; ++i )
1583  {
1584  if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
1586  a[i].fixedPosition,a[i].transform);
1587  continue;
1588  }
1589  patch->lattice.rescale(a[i].position,factor);
1590  }
1591  }
1592  }
1593  } else {
1595  }
1596 }
1597 
1599 {
1600  if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
1601  {
1602  //
1603  // DJH: Loops below simplify if we lift out special cases of fixed atoms
1604  // and pressure excluded atoms and make them their own branch.
1605  //
1606  FullAtom *a = patch->atom.begin();
1607  int numAtoms = patch->numAtoms;
1608  // Blocking receive for the updated lattice scaling factor.
1609  Tensor factor = broadcast->positionRescaleFactor.get(step);
1610  TIMER_START(patch->timerSet, PISTON);
1611  // JCP FIX THIS!!!
1612  Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
1613  patch->lattice.rescale(factor);
1614  Molecule *mol = Node::Object()->molecule;
1615  if ( simParams->useGroupPressure )
1616  {
1617  int hgs;
1618  for ( int i = 0; i < numAtoms; i += hgs ) {
1619  int j;
1620  hgs = a[i].hydrogenGroupSize;
1621  if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
1622  for ( j = i; j < (i+hgs); ++j ) {
1624  a[j].fixedPosition,a[j].transform);
1625  }
1626  continue;
1627  }
1628  BigReal m_cm = 0;
1629  Position x_cm(0,0,0);
1630  Velocity v_cm(0,0,0);
1631  for ( j = i; j < (i+hgs); ++j ) {
1632  if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
1633  m_cm += a[j].mass;
1634  x_cm += a[j].mass * a[j].position;
1635  v_cm += a[j].mass * a[j].velocity;
1636  }
1637  x_cm /= m_cm;
1638  Position new_x_cm = x_cm;
1639  patch->lattice.rescale(new_x_cm,factor);
1640  Position delta_x_cm = new_x_cm - x_cm;
1641  v_cm /= m_cm;
1642  Velocity delta_v_cm;
1643  delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
1644  delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
1645  delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
1646  for ( j = i; j < (i+hgs); ++j ) {
1647  if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
1649  a[j].fixedPosition,a[j].transform);
1650  continue;
1651  }
1652  if ( mol->is_atom_exPressure(a[j].id) ) continue;
1653  a[j].position += delta_x_cm;
1654  a[j].velocity += delta_v_cm;
1655  }
1656  }
1657  }
1658  else
1659  {
1660  for ( int i = 0; i < numAtoms; ++i )
1661  {
1662  if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
1664  a[i].fixedPosition,a[i].transform);
1665  continue;
1666  }
1667  if ( mol->is_atom_exPressure(a[i].id) ) continue;
1668  patch->lattice.rescale(a[i].position,factor);
1669  a[i].velocity.x *= velFactor.x;
1670  a[i].velocity.y *= velFactor.y;
1671  a[i].velocity.z *= velFactor.z;
1672  }
1673  }
1674  TIMER_STOP(patch->timerSet, PISTON);
1675  }
1676 }
1677 
1679 {
1680  const int rescaleFreq = simParams->rescaleFreq;
1681  if ( rescaleFreq > 0 ) {
1682  FullAtom *a = patch->atom.begin();
1683  int numAtoms = patch->numAtoms;
1685  if ( rescaleVelocities_numTemps == rescaleFreq ) {
1686  // Blocking receive for the velcity scaling factor.
1687  BigReal factor = broadcast->velocityRescaleFactor.get(step);
1688  for ( int i = 0; i < numAtoms; ++i )
1689  {
1690  a[i].velocity *= factor;
1691  }
1693  }
1694  }
1695 }
1696 
1697 void Sequencer::rescaleaccelMD (int step, int doNonbonded, int doFullElectrostatics)
1698 {
1699  if (!simParams->accelMDOn) return;
1700  if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
1701 
1702  // Blocking receive for the Accelerated MD scaling factors.
1703  Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
1704  const BigReal factor_dihe = accelMDfactor[0];
1705  const BigReal factor_tot = accelMDfactor[1];
1706  const int numAtoms = patch->numAtoms;
1707 
1708  if (simParams->accelMDdihe && factor_tot <1 )
1709  NAMD_die("accelMD broadcasting error!\n");
1710  if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
1711  NAMD_die("accelMD broadcasting error!\n");
1712 
1713  if (simParams->accelMDdihe && factor_dihe < 1) {
1714  for (int i = 0; i < numAtoms; ++i)
1715  if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
1716  patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);
1717  }
1718 
1719  if ( !simParams->accelMDdihe && factor_tot < 1) {
1720  for (int i = 0; i < numAtoms; ++i)
1721  patch->f[Results::normal][i] *= factor_tot;
1722  if (doNonbonded) {
1723  for (int i = 0; i < numAtoms; ++i)
1724  patch->f[Results::nbond][i] *= factor_tot;
1725  }
1726  if (doFullElectrostatics) {
1727  for (int i = 0; i < numAtoms; ++i)
1728  patch->f[Results::slow][i] *= factor_tot;
1729  }
1730  }
1731 
1732  if (simParams->accelMDdual && factor_dihe < 1) {
1733  for (int i = 0; i < numAtoms; ++i)
1734  if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
1735  patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
1736  }
1737 
1738 }
1739 
1741 {
1742  //check if adaptive tempering is enabled and in the right timestep range
1743  if (!simParams->adaptTempOn) return;
1744  if ( (step < simParams->adaptTempFirstStep ) ||
1746  if (simParams->langevinOn) // restore langevin temperature
1748  return;
1749  }
1750  // Get Updated Temperature
1751  if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
1752  // Blocking receive for the updated adaptive tempering temperature.
1754 }
1755 
1757 {
1758  const int reassignFreq = simParams->reassignFreq;
1759  if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1760  FullAtom *a = patch->atom.begin();
1761  int numAtoms = patch->numAtoms;
1762  BigReal newTemp = simParams->reassignTemp;
1763  newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
1764  if ( simParams->reassignIncr > 0.0 ) {
1765  if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
1766  newTemp = simParams->reassignHold;
1767  } else {
1768  if ( newTemp < simParams->reassignHold )
1769  newTemp = simParams->reassignHold;
1770  }
1771  BigReal kbT = BOLTZMANN * newTemp;
1772 
1773  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1774  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1775 
1776  for ( int i = 0; i < numAtoms; ++i )
1777  {
1778  a[i].velocity = ( ( simParams->fixedAtomsOn &&
1779  a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
1780  sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) *
1781  random->gaussian_vector() );
1782  }
1783  } else {
1784  NAMD_bug("Sequencer::reassignVelocities called improperly!");
1785  }
1786 }
1787 
1789 {
1790  FullAtom *a = patch->atom.begin();
1791  int numAtoms = patch->numAtoms;
1792  BigReal newTemp = simParams->initialTemp;
1793  BigReal kbT = BOLTZMANN * newTemp;
1794 
1795  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1796  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1797 
1798  for ( int i = 0; i < numAtoms; ++i )
1799  {
1800  a[i].velocity = ( ( (simParams->fixedAtomsOn && a[i].atomFixed) ||
1801  a[i].mass <= 0.) ? Vector(0,0,0) :
1802  sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) *
1803  random->gaussian_vector() );
1804  if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass > 0.05 ) {
1805  a[i+1].velocity = a[i].velocity; // zero is good enough
1806  ++i;
1807  }
1808  }
1809 }
1810 
1812 {
1813  FullAtom *a = patch->atom.begin();
1814  int numAtoms = patch->numAtoms;
1815  for ( int i = 0; i < numAtoms; ++i )
1816  {
1817  a[i].velocity *= factor;
1818  }
1819 }
1820 
1822 {
1823  FullAtom *a = patch->atom.begin();
1824  int numAtoms = patch->numAtoms;
1825  Molecule *molecule = Node::Object()->molecule;
1826  for ( int i = 0; i < numAtoms; ++i )
1827  {
1828  a[i].charge = molecule->atomcharge(a[i].id);
1829  }
1830 }
1831 
1832 // REST2 solute charge scaling
1834 {
1835  FullAtom *a = patch->atom.begin();
1836  int numAtoms = patch->numAtoms;
1837  Molecule *molecule = Node::Object()->molecule;
1838  BigReal sqrt_factor = sqrt(factor);
1839  // apply scaling to the original charge (stored in molecule)
1840  // of just the marked solute atoms
1841  for ( int i = 0; i < numAtoms; ++i ) {
1842  if (molecule->get_ss_type(a[i].id)) {
1843  a[i].charge = sqrt_factor * molecule->atomcharge(a[i].id);
1844  }
1845  }
1846 }
1847 
1849 {
1850  if ( simParams->tCoupleOn )
1851  {
1852  FullAtom *a = patch->atom.begin();
1853  int numAtoms = patch->numAtoms;
1854  // Blocking receive for the temperature coupling coefficient.
1855  BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
1856  Molecule *molecule = Node::Object()->molecule;
1857  BigReal dt = dt_fs * 0.001; // convert to ps
1858  coefficient *= dt;
1859  for ( int i = 0; i < numAtoms; ++i )
1860  {
1861  BigReal f1 = exp( coefficient * a[i].langevinParam );
1862  a[i].velocity *= f1;
1863  }
1864  }
1865 }
1866 
1873 {
1874  if ( simParams->stochRescaleOn )
1875  {
1876  FullAtom *a = patch->atom.begin();
1877  int numAtoms = patch->numAtoms;
1880  {
1881  // Blocking receive for the temperature coupling coefficient.
1882  BigReal coefficient = broadcast->stochRescaleCoefficient.get(step);
1883  DebugM(4, "stochastically rescaling velocities at step " << step << " by " << coefficient << "\n");
1884  for ( int i = 0; i < numAtoms; ++i )
1885  {
1886  a[i].velocity *= coefficient;
1887  }
1888  stochRescale_count = 0;
1889  }
1890  }
1891 }
1892 
1893 void Sequencer::saveForce(const int ftag)
1894 {
1895  patch->saveForce(ftag);
1896 }
1897 
1898 //
1899 // DJH: Need to change division by TIMEFACTOR into multiplication by
1900 // reciprocal of TIMEFACTOR. Done several times for each iteration of
1901 // the integrate() loop.
1902 //
1903 
1905  BigReal timestep, const int ftag, const int useSaved
1906  ) {
1907  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1908  NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
1909 #if CMK_BLUEGENEL
1910  CmiNetworkProgressAfter (0);
1911 #endif
1912  const BigReal dt = timestep / TIMEFACTOR;
1913  FullAtom *atom_arr = patch->atom.begin();
1914  ForceList *f_use = (useSaved ? patch->f_saved : patch->f);
1915  const Force *force_arr = f_use[ftag].const_begin();
1916  patch->addForceToMomentum(atom_arr, force_arr, dt, patch->numAtoms);
1917 }
1918 
1920  const BigReal timestep1, const int ftag1, const int useSaved1,
1921  const BigReal timestep2, const int ftag2, const int useSaved2,
1922  const BigReal timestep3, const int ftag3, const int useSaved3
1923  ) {
1924  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1925  NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
1926 #if CMK_BLUEGENEL
1927  CmiNetworkProgressAfter (0);
1928 #endif
1929  const BigReal dt1 = timestep1 / TIMEFACTOR;
1930  const BigReal dt2 = timestep2 / TIMEFACTOR;
1931  const BigReal dt3 = timestep3 / TIMEFACTOR;
1932  ForceList *f_use1 = (useSaved1 ? patch->f_saved : patch->f);
1933  ForceList *f_use2 = (useSaved2 ? patch->f_saved : patch->f);
1934  ForceList *f_use3 = (useSaved3 ? patch->f_saved : patch->f);
1935  FullAtom *atom_arr = patch->atom.begin();
1936  const Force *force_arr1 = f_use1[ftag1].const_begin();
1937  const Force *force_arr2 = f_use2[ftag2].const_begin();
1938  const Force *force_arr3 = f_use3[ftag3].const_begin();
1939  patch->addForceToMomentum3 (atom_arr, force_arr1, force_arr2, force_arr3,
1940  dt1, dt2, dt3, patch->numAtoms);
1941 }
1942 
1944 {
1945  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1946  NamdProfileEvent::ADD_VELOCITY_TO_POSITION);
1947 #if CMK_BLUEGENEL
1948  CmiNetworkProgressAfter (0);
1949 #endif
1950  const BigReal dt = timestep / TIMEFACTOR;
1951  FullAtom *atom_arr = patch->atom.begin();
1952  patch->addVelocityToPosition(atom_arr, dt, patch->numAtoms);
1953 }
1954 
1956 {
1957  if ( simParams->drudeHardWallOn ) {
1958  Tensor virial;
1959  Tensor *vp = ( pressure ? &virial : 0 );
1960  if ( patch->hardWallDrude(dt, vp, pressureProfileReduction) ) {
1961  iout << iERROR << "Constraint failure in HardWallDrude(); "
1962  << "simulation may become unstable.\n" << endi;
1964  terminate();
1965  }
1966  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
1967  }
1968 }
1969 
1970 void Sequencer::rattle1(BigReal dt, int pressure)
1971 {
1972  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::RATTLE1);
1973  if ( simParams->rigidBonds != RIGID_NONE ) {
1974  Tensor virial;
1975  Tensor *vp = ( pressure ? &virial : 0 );
1976  if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
1977  iout << iERROR <<
1978  "Constraint failure; simulation has become unstable.\n" << endi;
1980  terminate();
1981  }
1982 #if 0
1983  printf("virial = %g %g %g %g %g %g %g %g %g\n",
1984  virial.xx, virial.xy, virial.xz,
1985  virial.yx, virial.yy, virial.yz,
1986  virial.zx, virial.zy, virial.zz);
1987 #endif
1988  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
1989 #if 0
1990  {
1991  const FullAtom *a = patch->atom.const_begin();
1992  for (int n=0; n < patch->numAtoms; n++) {
1993  printf("pos[%d] = %g %g %g\n", n,
1994  a[n].position.x, a[n].position.y, a[n].position.z);
1995  }
1996  for (int n=0; n < patch->numAtoms; n++) {
1997  printf("vel[%d] = %g %g %g\n", n,
1998  a[n].velocity.x, a[n].velocity.y, a[n].velocity.z);
1999  }
2000  if (pressure) {
2001  for (int n=0; n < patch->numAtoms; n++) {
2002  printf("force[%d] = %g %g %g\n", n,
2003  patch->f[Results::normal][n].x,
2004  patch->f[Results::normal][n].y,
2005  patch->f[Results::normal][n].z);
2006  }
2007  }
2008  }
2009 #endif
2010  }
2011 }
2012 
2013 // void Sequencer::rattle2(BigReal dt, int step)
2014 // {
2015 // if ( simParams->rigidBonds != RIGID_NONE ) {
2016 // Tensor virial;
2017 // patch->rattle2(dt, &virial);
2018 // ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
2019 // // we need to add to alt and int virial because not included in forces
2020 // #ifdef ALTVIRIAL
2021 // ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
2022 // #endif
2023 // ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,virial);
2024 // }
2025 // }
2026 
2028 {
2029  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::MAXIMUM_MOVE);
2030 
2031  FullAtom *a = patch->atom.begin();
2032  int numAtoms = patch->numAtoms;
2033  if ( simParams->maximumMove ) {
2034  const BigReal dt = timestep / TIMEFACTOR;
2035  const BigReal maxvel = simParams->maximumMove / dt;
2036  const BigReal maxvel2 = maxvel * maxvel;
2037  for ( int i=0; i<numAtoms; ++i ) {
2038  if ( a[i].velocity.length2() > maxvel2 ) {
2039  a[i].velocity *= ( maxvel / a[i].velocity.length() );
2040  }
2041  }
2042  } else {
2043  const BigReal dt = timestep / TIMEFACTOR;
2044  const BigReal maxvel = simParams->cutoff / dt;
2045  const BigReal maxvel2 = maxvel * maxvel;
2046  int killme = 0;
2047  for ( int i=0; i<numAtoms; ++i ) {
2048  killme = killme || ( a[i].velocity.length2() > maxvel2 );
2049  }
2050  if ( killme ) {
2051  killme = 0;
2052  for ( int i=0; i<numAtoms; ++i ) {
2053  if ( a[i].velocity.length2() > maxvel2 ) {
2054  ++killme;
2055  iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
2056  << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
2057  << ( PDBVELFACTOR * maxvel ) << ", atom "
2058  << i << " of " << numAtoms << " on patch "
2059  << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
2060  }
2061  }
2062  iout << iERROR <<
2063  "Atoms moving too fast; simulation has become unstable ("
2064  << killme << " atoms on patch " << patch->patchID
2065  << " pe " << CkMyPe() << ").\n" << endi;
2067  terminate();
2068  }
2069  }
2070 }
2071 
2073 {
2074  if ( simParams->minimizeOn ) {
2075  FullAtom *a = patch->atom.begin();
2076  int numAtoms = patch->numAtoms;
2077  for ( int i=0; i<numAtoms; ++i ) {
2078  a[i].velocity = 0.;
2079  }
2080  }
2081 }
2082 
2084 {
2085  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::SUBMIT_HALFSTEP);
2086 
2087  // velocity-dependent quantities *** ONLY ***
2088  // positions are not at half-step when called
2089  FullAtom *a = patch->atom.begin();
2090  int numAtoms = patch->numAtoms;
2091 
2092 #if CMK_BLUEGENEL
2093  CmiNetworkProgressAfter (0);
2094 #endif
2095 
2096  // For non-Multigrator doKineticEnergy = 1 always
2097  Tensor momentumSqrSum;
2099  {
2100  BigReal kineticEnergy = 0;
2101  Tensor virial;
2102  if ( simParams->pairInteractionOn ) {
2103  if ( simParams->pairInteractionSelf ) {
2104  for ( int i = 0; i < numAtoms; ++i ) {
2105  if ( a[i].partition != 1 ) continue;
2106  kineticEnergy += a[i].mass * a[i].velocity.length2();
2107  virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
2108  }
2109  }
2110  } else {
2111  for ( int i = 0; i < numAtoms; ++i ) {
2112  if (a[i].mass < 0.01) continue;
2113  kineticEnergy += a[i].mass * a[i].velocity.length2();
2114  virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
2115  }
2116  }
2117 
2119  momentumSqrSum = virial;
2120  }
2121  kineticEnergy *= 0.5 * 0.5;
2123  virial *= 0.5;
2124  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
2125 #ifdef ALTVIRIAL
2126  ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
2127 #endif
2128  }
2129 
2131  int nslabs = simParams->pressureProfileSlabs;
2132  const Lattice &lattice = patch->lattice;
2133  BigReal idz = nslabs/lattice.c().z;
2134  BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
2135  int useGroupPressure = simParams->useGroupPressure;
2136 
2137  // Compute kinetic energy partition, possibly subtracting off
2138  // internal kinetic energy if group pressure is enabled.
2139  // Since the regular pressure is 1/2 mvv and the internal kinetic
2140  // term that is subtracted off for the group pressure is
2141  // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
2142  // 1/2 m * v * v_cm. The factor of 1/2 is because submitHalfstep
2143  // gets called twice per timestep.
2144  int hgs;
2145  for (int i=0; i<numAtoms; i += hgs) {
2146  int j, ppoffset;
2147  hgs = a[i].hydrogenGroupSize;
2148  int partition = a[i].partition;
2149 
2150  BigReal m_cm = 0;
2151  Velocity v_cm(0,0,0);
2152  for (j=i; j< i+hgs; ++j) {
2153  m_cm += a[j].mass;
2154  v_cm += a[j].mass * a[j].velocity;
2155  }
2156  v_cm /= m_cm;
2157  for (j=i; j < i+hgs; ++j) {
2158  BigReal mass = a[j].mass;
2159  if (! (useGroupPressure && j != i)) {
2160  BigReal z = a[j].position.z;
2161  int slab = (int)floor((z-zmin)*idz);
2162  if (slab < 0) slab += nslabs;
2163  else if (slab >= nslabs) slab -= nslabs;
2164  ppoffset = 3*(slab + partition*nslabs);
2165  }
2166  BigReal wxx, wyy, wzz;
2167  if (useGroupPressure) {
2168  wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
2169  wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
2170  wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
2171  } else {
2172  wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
2173  wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
2174  wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
2175  }
2176  pressureProfileReduction->item(ppoffset ) += wxx;
2177  pressureProfileReduction->item(ppoffset+1) += wyy;
2178  pressureProfileReduction->item(ppoffset+2) += wzz;
2179  }
2180  }
2181  }
2182 
2183  // For non-Multigrator doKineticEnergy = 1 always
2185  {
2186  BigReal intKineticEnergy = 0;
2187  Tensor intVirialNormal;
2188 
2189  int hgs;
2190  for ( int i = 0; i < numAtoms; i += hgs ) {
2191 
2192 #if CMK_BLUEGENEL
2193  CmiNetworkProgress ();
2194 #endif
2195 
2196  hgs = a[i].hydrogenGroupSize;
2197  int j;
2198  BigReal m_cm = 0;
2199  Velocity v_cm(0,0,0);
2200  for ( j = i; j < (i+hgs); ++j ) {
2201  m_cm += a[j].mass;
2202  v_cm += a[j].mass * a[j].velocity;
2203  }
2205  momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
2206  }
2207  v_cm /= m_cm;
2208  if ( simParams->pairInteractionOn ) {
2209  if ( simParams->pairInteractionSelf ) {
2210  for ( j = i; j < (i+hgs); ++j ) {
2211  if ( a[j].partition != 1 ) continue;
2212  BigReal mass = a[j].mass;
2213  Vector v = a[j].velocity;
2214  Vector dv = v - v_cm;
2215  intKineticEnergy += mass * (v * dv);
2216  intVirialNormal.outerAdd (mass, v, dv);
2217  }
2218  }
2219  } else {
2220  for ( j = i; j < (i+hgs); ++j ) {
2221  BigReal mass = a[j].mass;
2222  Vector v = a[j].velocity;
2223  Vector dv = v - v_cm;
2224  intKineticEnergy += mass * (v * dv);
2225  intVirialNormal.outerAdd(mass, v, dv);
2226  }
2227  }
2228  }
2229 
2230  intKineticEnergy *= 0.5 * 0.5;
2232  intVirialNormal *= 0.5;
2233  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
2234  if ( simParams->multigratorOn) {
2235  momentumSqrSum *= 0.5;
2236  ADD_TENSOR_OBJECT(reduction,REDUCTION_MOMENTUM_SQUARED,momentumSqrSum);
2237  }
2238  }
2239 
2240 }
2241 
2242 void Sequencer::calcFixVirial(Tensor& fixVirialNormal, Tensor& fixVirialNbond, Tensor& fixVirialSlow,
2243  Vector& fixForceNormal, Vector& fixForceNbond, Vector& fixForceSlow) {
2244 
2245  FullAtom *a = patch->atom.begin();
2246  int numAtoms = patch->numAtoms;
2247 
2248  for ( int j = 0; j < numAtoms; j++ ) {
2249  if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
2250  Vector dx = a[j].fixedPosition;
2251  // all negative because fixed atoms cancels these forces
2252  fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
2253  fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
2254  fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
2255  fixForceNormal -= patch->f[Results::normal][j];
2256  fixForceNbond -= patch->f[Results::nbond][j];
2257  fixForceSlow -= patch->f[Results::slow][j];
2258  }
2259  }
2260 }
2261 
2263 {
2264 #ifndef UPPER_BOUND
2265  NAMD_EVENT_RANGE_2(patch->flags.event_on,
2266  NamdProfileEvent::SUBMIT_REDUCTIONS);
2267  FullAtom *a = patch->atom.begin();
2268 #endif
2269  int numAtoms = patch->numAtoms;
2270 
2271 #if CMK_BLUEGENEL
2272  CmiNetworkProgressAfter(0);
2273 #endif
2274 
2275  reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
2277 
2278 #ifndef UPPER_BOUND
2279  // For non-Multigrator doKineticEnergy = 1 always
2281  {
2282  BigReal kineticEnergy = 0;
2283  Vector momentum = 0;
2284  Vector angularMomentum = 0;
2285  Vector o = patch->lattice.origin();
2286  int i;
2287  if ( simParams->pairInteractionOn ) {
2288  if ( simParams->pairInteractionSelf ) {
2289  for (i = 0; i < numAtoms; ++i ) {
2290  if ( a[i].partition != 1 ) continue;
2291  kineticEnergy += a[i].mass * a[i].velocity.length2();
2292  momentum += a[i].mass * a[i].velocity;
2293  angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
2294  }
2295  }
2296  } else {
2297  for (i = 0; i < numAtoms; ++i ) {
2298  kineticEnergy += a[i].mass * a[i].velocity.length2();
2299  momentum += a[i].mass * a[i].velocity;
2300  angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
2301  }
2302  if (simParams->drudeOn) {
2303  BigReal drudeComKE = 0.;
2304  BigReal drudeBondKE = 0.;
2305 
2306  for (i = 0; i < numAtoms; i++) {
2307  if (i < numAtoms-1 &&
2308  a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
2309  // i+1 is a Drude particle with parent i
2310 
2311  // convert from Cartesian coordinates to (COM,bond) coordinates
2312  BigReal m_com = (a[i].mass + a[i+1].mass); // mass of COM
2313  BigReal m = a[i+1].mass / m_com; // mass ratio
2314  BigReal m_bond = a[i+1].mass * (1. - m); // mass of bond
2315  Vector v_bond = a[i+1].velocity - a[i].velocity; // vel of bond
2316  Vector v_com = a[i].velocity + m * v_bond; // vel of COM
2317 
2318  drudeComKE += m_com * v_com.length2();
2319  drudeBondKE += m_bond * v_bond.length2();
2320 
2321  i++; // +1 from loop, we've updated both particles
2322  }
2323  else {
2324  drudeComKE += a[i].mass * a[i].velocity.length2();
2325  }
2326  } // end for
2327 
2328  drudeComKE *= 0.5;
2329  drudeBondKE *= 0.5;
2331  += drudeComKE;
2333  += drudeBondKE;
2334  } // end drudeOn
2335 
2336  } // end else
2337 
2338  kineticEnergy *= 0.5;
2340  ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
2341  ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
2342  }
2343 
2344 #ifdef ALTVIRIAL
2345  // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
2346  {
2347  Tensor altVirial;
2348  for ( int i = 0; i < numAtoms; ++i ) {
2349  altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
2350  }
2351  ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
2352  }
2353  {
2354  Tensor altVirial;
2355  for ( int i = 0; i < numAtoms; ++i ) {
2356  altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
2357  }
2358  ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
2359  }
2360  {
2361  Tensor altVirial;
2362  for ( int i = 0; i < numAtoms; ++i ) {
2363  altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
2364  }
2365  ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
2366  }
2367 #endif
2368 
2369  // For non-Multigrator doKineticEnergy = 1 always
2371  {
2372  BigReal intKineticEnergy = 0;
2373  Tensor intVirialNormal;
2374  Tensor intVirialNbond;
2375  Tensor intVirialSlow;
2376 
2377  int hgs;
2378  for ( int i = 0; i < numAtoms; i += hgs ) {
2379 #if CMK_BLUEGENEL
2380  CmiNetworkProgress();
2381 #endif
2382  hgs = a[i].hydrogenGroupSize;
2383  int j;
2384  BigReal m_cm = 0;
2385  Position x_cm(0,0,0);
2386  Velocity v_cm(0,0,0);
2387  for ( j = i; j < (i+hgs); ++j ) {
2388  m_cm += a[j].mass;
2389  x_cm += a[j].mass * a[j].position;
2390  v_cm += a[j].mass * a[j].velocity;
2391  }
2392  x_cm /= m_cm;
2393  v_cm /= m_cm;
2394  int fixedAtomsOn = simParams->fixedAtomsOn;
2395  if ( simParams->pairInteractionOn ) {
2396  int pairInteractionSelf = simParams->pairInteractionSelf;
2397  for ( j = i; j < (i+hgs); ++j ) {
2398  if ( a[j].partition != 1 &&
2399  ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
2400  // net force treated as zero for fixed atoms
2401  if ( fixedAtomsOn && a[j].atomFixed ) continue;
2402  BigReal mass = a[j].mass;
2403  Vector v = a[j].velocity;
2404  Vector dv = v - v_cm;
2405  intKineticEnergy += mass * (v * dv);
2406  Vector dx = a[j].position - x_cm;
2407  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
2408  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
2409  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
2410  }
2411  } else {
2412  for ( j = i; j < (i+hgs); ++j ) {
2413  // net force treated as zero for fixed atoms
2414  if ( fixedAtomsOn && a[j].atomFixed ) continue;
2415  BigReal mass = a[j].mass;
2416  Vector v = a[j].velocity;
2417  Vector dv = v - v_cm;
2418  intKineticEnergy += mass * (v * dv);
2419  Vector dx = a[j].position - x_cm;
2420  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
2421  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
2422  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
2423  }
2424  }
2425  }
2426 
2427  intKineticEnergy *= 0.5;
2429  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
2430  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
2431  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
2432  }
2433 
2435  // subtract off internal virial term, calculated as for intVirial.
2436  int nslabs = simParams->pressureProfileSlabs;
2437  const Lattice &lattice = patch->lattice;
2438  BigReal idz = nslabs/lattice.c().z;
2439  BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
2440  int useGroupPressure = simParams->useGroupPressure;
2441 
2442  int hgs;
2443  for (int i=0; i<numAtoms; i += hgs) {
2444  int j;
2445  hgs = a[i].hydrogenGroupSize;
2446  BigReal m_cm = 0;
2447  Position x_cm(0,0,0);
2448  for (j=i; j< i+hgs; ++j) {
2449  m_cm += a[j].mass;
2450  x_cm += a[j].mass * a[j].position;
2451  }
2452  x_cm /= m_cm;
2453 
2454  BigReal z = a[i].position.z;
2455  int slab = (int)floor((z-zmin)*idz);
2456  if (slab < 0) slab += nslabs;
2457  else if (slab >= nslabs) slab -= nslabs;
2458  int partition = a[i].partition;
2459  int ppoffset = 3*(slab + nslabs*partition);
2460  for (j=i; j < i+hgs; ++j) {
2461  BigReal mass = a[j].mass;
2462  Vector dx = a[j].position - x_cm;
2463  const Vector &fnormal = patch->f[Results::normal][j];
2464  const Vector &fnbond = patch->f[Results::nbond][j];
2465  const Vector &fslow = patch->f[Results::slow][j];
2466  BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
2467  BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
2468  BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
2469  pressureProfileReduction->item(ppoffset ) -= wxx;
2470  pressureProfileReduction->item(ppoffset+1) -= wyy;
2471  pressureProfileReduction->item(ppoffset+2) -= wzz;
2472  }
2473  }
2474  }
2475 
2476  // For non-Multigrator doVirial = 1 always
2477  if (patch->flags.doVirial)
2478  {
2479  if ( simParams->fixedAtomsOn ) {
2480  Tensor fixVirialNormal;
2481  Tensor fixVirialNbond;
2482  Tensor fixVirialSlow;
2483  Vector fixForceNormal = 0;
2484  Vector fixForceNbond = 0;
2485  Vector fixForceSlow = 0;
2486 
2487  calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
2488 
2489  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
2490  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
2491  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
2492  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
2493  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
2494  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
2495  }
2496  }
2497 #endif // UPPER_BOUND
2498 
2499  reduction->submit();
2500 #ifndef UPPER_BOUND
2502 #endif
2503 }
2504 
2506 {
2507  FullAtom *a = patch->atom.begin();
2508  Force *f1 = patch->f[Results::normal].begin();
2509  Force *f2 = patch->f[Results::nbond].begin();
2510  Force *f3 = patch->f[Results::slow].begin();
2511  const bool fixedAtomsOn = simParams->fixedAtomsOn;
2512  const bool drudeHardWallOn = simParams->drudeHardWallOn;
2513  const double drudeBondLen = simParams->drudeBondLen;
2514  const double drudeBondLen2 = drudeBondLen * drudeBondLen;
2515  const double drudeStep = 0.1/(TIMEFACTOR*TIMEFACTOR);
2516  const double drudeMove = 0.01;
2517  const double drudeStep2 = drudeStep * drudeStep;
2518  const double drudeMove2 = drudeMove * drudeMove;
2519  int numAtoms = patch->numAtoms;
2520 
2521  reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
2522 
2523  for ( int i = 0; i < numAtoms; ++i ) {
2524  f1[i] += f2[i] + f3[i]; // add all forces
2525  if ( drudeHardWallOn && i && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) { // drude particle
2526  if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
2527  if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
2528  a[i].position += drudeMove * f1[i].unit();
2529  } else {
2530  a[i].position += drudeStep * f1[i];
2531  }
2532  if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
2533  a[i].position = a[i-1].position + drudeBondLen * (a[i].position - a[i-1].position).unit();
2534  }
2535  }
2536  Vector netf = f1[i-1] + f1[i];
2537  if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
2538  f1[i-1] = netf;
2539  f1[i] = 0.;
2540  }
2541  if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
2542  }
2543 
2544  f2 = f3 = 0; // included in f1
2545 
2546  BigReal maxv2 = 0.;
2547 
2548  for ( int i = 0; i < numAtoms; ++i ) {
2549  BigReal v2 = a[i].velocity.length2();
2550  if ( v2 > 0. ) {
2551  if ( v2 > maxv2 ) maxv2 = v2;
2552  } else {
2553  v2 = f1[i].length2();
2554  if ( v2 > maxv2 ) maxv2 = v2;
2555  }
2556  }
2557 
2558  if ( fmax2 > 10. * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR )
2559  { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial, true /* forces */); }
2560 
2561  BigReal fdotf = 0;
2562  BigReal fdotv = 0;
2563  BigReal vdotv = 0;
2564  int numHuge = 0;
2565  for ( int i = 0; i < numAtoms; ++i ) {
2566  if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
2567  if ( drudeHardWallOn && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) continue; // drude particle
2568  Force f = f1[i];
2569  BigReal ff = f * f;
2570  if ( ff > fmax2 ) {
2571  if (simParams->printBadContacts) {
2572  CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
2573  }
2574  ++numHuge;
2575  // pad scaling so minimizeMoveDownhill() doesn't miss them
2576  BigReal fmult = 1.01 * sqrt(fmax2/ff);
2577  f *= fmult; ff = f * f;
2578  f1[i] *= fmult;
2579  }
2580  fdotf += ff;
2581  fdotv += f * a[i].velocity;
2582  vdotv += a[i].velocity * a[i].velocity;
2583  }
2584 
2589 
2590  {
2591  Tensor intVirialNormal;
2592  Tensor intVirialNbond;
2593  Tensor intVirialSlow;
2594 
2595  int hgs;
2596  for ( int i = 0; i < numAtoms; i += hgs ) {
2597  hgs = a[i].hydrogenGroupSize;
2598  int j;
2599  BigReal m_cm = 0;
2600  Position x_cm(0,0,0);
2601  for ( j = i; j < (i+hgs); ++j ) {
2602  m_cm += a[j].mass;
2603  x_cm += a[j].mass * a[j].position;
2604  }
2605  x_cm /= m_cm;
2606  for ( j = i; j < (i+hgs); ++j ) {
2607  BigReal mass = a[j].mass;
2608  // net force treated as zero for fixed atoms
2609  if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
2610  Vector dx = a[j].position - x_cm;
2611  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
2612  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
2613  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
2614  }
2615  }
2616 
2617  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
2618  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
2619  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
2620  }
2621 
2622  if ( simParams->fixedAtomsOn ) {
2623  Tensor fixVirialNormal;
2624  Tensor fixVirialNbond;
2625  Tensor fixVirialSlow;
2626  Vector fixForceNormal = 0;
2627  Vector fixForceNbond = 0;
2628  Vector fixForceSlow = 0;
2629 
2630  calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
2631 
2632  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
2633  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
2634  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
2635  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
2636  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
2637  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
2638  }
2639 
2640  reduction->submit();
2641 }
2642 
2643 void Sequencer::submitCollections(int step, int zeroVel)
2644 {
2645  //
2646  // DJH: Copy updates of SOA back into AOS.
2647  // Do we need to update everything or is it safe to just update
2648  // positions and velocities separately, as needed?
2649  //
2650  //patch->copy_updates_to_AOS();
2651 
2652  NAMD_EVENT_RANGE_2(patch->flags.event_on,
2653  NamdProfileEvent::SUBMIT_COLLECTIONS);
2654  int prec = Output::coordinateNeeded(step);
2655  if ( prec ) {
2656  collection->submitPositions(step,patch->atom,patch->lattice,prec);
2657  }
2658  if ( Output::velocityNeeded(step) ) {
2659  collection->submitVelocities(step,zeroVel,patch->atom);
2660  }
2661  if ( Output::forceNeeded(step) ) {
2662  int maxForceUsed = patch->flags.maxForceUsed;
2663  if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
2664  collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
2665  }
2666 }
2667 
2668 void Sequencer::runComputeObjects(int migration, int pairlists, int pressureStep)
2669 {
2670  if ( migration ) pairlistsAreValid = 0;
2671 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
2672  if ( pairlistsAreValid &&
2674  && ( pairlistsAge > (
2675 #else
2676  if ( pairlistsAreValid && ( pairlistsAge > (
2677 #endif
2679  pairlistsAreValid = 0;
2680  }
2681  if ( ! simParams->usePairlists ) pairlists = 0;
2682  patch->flags.usePairlists = pairlists || pairlistsAreValid;
2683  patch->flags.savePairlists =
2684  pairlists && ! pairlistsAreValid;
2685 
2686  if ( simParams->singleTopology ) patch->reposition_all_alchpairs();
2687  if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
2688 
2689  //
2690  // DJH: Copy updates of SOA back into AOS.
2691  // The positionsReady() routine starts force computation and atom migration.
2692  //
2693  // We could reduce amount of copying here by checking migration status
2694  // and copying velocities only when migrating. Some types of simulation
2695  // always require velocities, such as Lowe-Anderson.
2696  //
2697  //patch->copy_updates_to_AOS();
2698 
2699  patch->positionsReady(migration); // updates flags.sequence
2700 
2701  int seq = patch->flags.sequence;
2702  int basePriority = ( (seq & 0xffff) << 15 )
2703  + PATCH_PRIORITY(patch->getPatchID());
2704  if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
2705  priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
2706  suspend(); // until all deposit boxes close
2707  patch->gbisComputeAfterP1();
2708  priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
2709  suspend();
2710  patch->gbisComputeAfterP2();
2711  priority = basePriority + COMPUTE_HOME_PRIORITY;
2712  suspend();
2713  } else {
2714  priority = basePriority + COMPUTE_HOME_PRIORITY;
2715  suspend(); // until all deposit boxes close
2716  }
2717 
2718  //
2719  // DJH: Copy all data into SOA from AOS.
2720  //
2721  // We need everything copied after atom migration.
2722  // When doing force computation without atom migration,
2723  // all data except forces will already be up-to-date in SOA
2724  // (except maybe for some special types of simulation).
2725  //
2726  //patch->copy_all_to_SOA();
2727 
2728  //
2729  // DJH: Copy forces to SOA.
2730  // Force available after suspend() has returned.
2731  //
2732  //patch->copy_forces_to_SOA();
2733 
2734  if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
2735  pairlistsAreValid = 1;
2736  pairlistsAge = 0;
2737  }
2738  // For multigrator, do not age pairlist during pressure step
2739  // NOTE: for non-multigrator pressureStep = 0 always
2740  if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
2741 
2742  if (simParams->lonepairs) {
2743  {
2744  Tensor virial;
2745  patch->redistrib_lonepair_forces(Results::normal, &virial);
2746  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
2747  }
2748  if (patch->flags.doNonbonded) {
2749  Tensor virial;
2750  patch->redistrib_lonepair_forces(Results::nbond, &virial);
2751  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
2752  }
2753  if (patch->flags.doFullElectrostatics) {
2754  Tensor virial;
2755  patch->redistrib_lonepair_forces(Results::slow, &virial);
2756  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
2757  }
2758  } else if (simParams->watmodel == WAT_TIP4) {
2759  {
2760  Tensor virial;
2761  patch->redistrib_tip4p_forces(Results::normal, &virial);
2762  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
2763  }
2764  if (patch->flags.doNonbonded) {
2765  Tensor virial;
2766  patch->redistrib_tip4p_forces(Results::nbond, &virial);
2767  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
2768  }
2769  if (patch->flags.doFullElectrostatics) {
2770  Tensor virial;
2771  patch->redistrib_tip4p_forces(Results::slow, &virial);
2772  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
2773  }
2774  } else if (simParams->watmodel == WAT_SWM4) {
2775  {
2776  Tensor virial;
2777  patch->redistrib_swm4_forces(Results::normal, &virial);
2778  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
2779  }
2780  if (patch->flags.doNonbonded) {
2781  Tensor virial;
2782  patch->redistrib_swm4_forces(Results::nbond, &virial);
2783  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
2784  }
2785  if (patch->flags.doFullElectrostatics) {
2786  Tensor virial;
2787  patch->redistrib_swm4_forces(Results::slow, &virial);
2788  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
2789  }
2790  }
2791 
2792  if (simParams->singleTopology) {
2793  patch->redistrib_alchpair_forces(Results::normal);
2794  if (patch->flags.doNonbonded) {
2795  patch->redistrib_alchpair_forces(Results::nbond);
2796  }
2797  if (patch->flags.doFullElectrostatics) {
2798  patch->redistrib_alchpair_forces(Results::slow);
2799  }
2800  }
2801 
2802  if ( patch->flags.doMolly ) {
2803  Tensor virial;
2804  patch->mollyMollify(&virial);
2805  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
2806  }
2807 
2808 
2809  // BEGIN LA
2810  if (patch->flags.doLoweAndersen) {
2811  patch->loweAndersenFinish();
2812  }
2813  // END LA
2814 #ifdef NAMD_CUDA_XXX
2815  int numAtoms = patch->numAtoms;
2816  FullAtom *a = patch->atom.begin();
2817  for ( int i=0; i<numAtoms; ++i ) {
2818  CkPrintf("%d %g %g %g\n", a[i].id,
2819  patch->f[Results::normal][i].x +
2820  patch->f[Results::nbond][i].x +
2821  patch->f[Results::slow][i].x,
2822  patch->f[Results::normal][i].y +
2823  patch->f[Results::nbond][i].y +
2824  patch->f[Results::slow][i].y,
2825  patch->f[Results::normal][i].z +
2826  patch->f[Results::nbond][i].z +
2827  patch->f[Results::slow][i].z);
2828  CkPrintf("%d %g %g %g\n", a[i].id,
2829  patch->f[Results::normal][i].x,
2830  patch->f[Results::nbond][i].x,
2831  patch->f[Results::slow][i].x);
2832  CkPrintf("%d %g %g %g\n", a[i].id,
2833  patch->f[Results::normal][i].y,
2834  patch->f[Results::nbond][i].y,
2835  patch->f[Results::slow][i].y);
2836  CkPrintf("%d %g %g %g\n", a[i].id,
2837  patch->f[Results::normal][i].z,
2838  patch->f[Results::nbond][i].z,
2839  patch->f[Results::slow][i].z);
2840  }
2841 #endif
2842 
2843 #if PRINT_FORCES
2844  int numAtoms = patch->numAtoms;
2845  FullAtom *a = patch->atom.begin();
2846  for ( int i=0; i<numAtoms; ++i ) {
2847  float fxNo = patch->f[Results::normal][i].x;
2848  float fxNb = patch->f[Results::nbond][i].x;
2849  float fxSl = patch->f[Results::slow][i].x;
2850  float fyNo = patch->f[Results::normal][i].y;
2851  float fyNb = patch->f[Results::nbond][i].y;
2852  float fySl = patch->f[Results::slow][i].y;
2853  float fzNo = patch->f[Results::normal][i].z;
2854  float fzNb = patch->f[Results::nbond][i].z;
2855  float fzSl = patch->f[Results::slow][i].z;
2856  float fx = fxNo+fxNb+fxSl;
2857  float fy = fyNo+fyNb+fySl;
2858  float fz = fzNo+fzNb+fzSl;
2859 
2860  float f = sqrt(fx*fx+fy*fy+fz*fz);
2861  int id = patch->pExt[i].id;
2862  int seq = patch->flags.sequence;
2863  float x = patch->p[i].position.x;
2864  float y = patch->p[i].position.y;
2865  float z = patch->p[i].position.z;
2866  //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
2867  CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
2868  //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
2869 //fxNo,fyNo,fzNo,
2870 //fxNb,fyNb,fzNb,
2871 //fxSl,fySl,fzSl,
2872 fx,fy,fz
2873 );
2874  }
2875 #endif
2876 }
2877 
2878 void Sequencer::rebalanceLoad(int timestep) {
2879  if ( ! ldbSteps ) {
2881  }
2882  if ( ! --ldbSteps ) {
2883  patch->submitLoadStats(timestep);
2884  ldbCoordinator->rebalance(this,patch->getPatchID());
2885  pairlistsAreValid = 0;
2886  }
2887 }
2888 
2889 void Sequencer::cycleBarrier(int doBarrier, int step) {
2890 #if USE_BARRIER
2891  if (doBarrier)
2892  // Blocking receive for the cycle barrier.
2893  broadcast->cycleBarrier.get(step);
2894 #endif
2895 }
2896 
2898  // Blocking receive for the trace barrier.
2899  broadcast->traceBarrier.get(step);
2900 }
2901 
2902 #ifdef MEASURE_NAMD_WITH_PAPI
2903 void Sequencer::papiMeasureBarrier(int step){
2904  // Blocking receive for the PAPI measure barrier.
2905  broadcast->papiMeasureBarrier.get(step);
2906 }
2907 #endif
2908 
2911  CthFree(thread);
2912  CthSuspend();
2913 }
2914 
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
SubmitReduction * multigratorReduction
Definition: Sequencer.h:119
Vector gaussian_vector(void)
Definition: Random.h:208
void rescaleVelocities(int)
Definition: Sequencer.C:1678
int doKineticEnergy
Definition: Sequencer.h:120
void minimizationQuenchVelocity(void)
Definition: Sequencer.C:2072
unsigned char partition
Definition: NamdTypes.h:56
void max(int i, BigReal v)
Definition: ReductionMgr.h:315
BigReal zy
Definition: Tensor.h:19
Bool berendsenPressureOn
Real langevinParam
Definition: NamdTypes.h:110
void tcoupleVelocities(BigReal, int)
Definition: Sequencer.C:1848
void addMovDragToPosition(BigReal)
Definition: Sequencer.C:718
void terminate(void)
Definition: Sequencer.C:2909
BigReal soluteScalingFactorCharge
virtual void algorithm(void)
Definition: Sequencer.C:146
#define NAMD_EVENT_STOP(eon, id)
SubmitReduction * pressureProfileReduction
Definition: Sequencer.h:135
void addVelocityToPosition(FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms)
Definition: HomePatch.C:2000
void minimize_rattle2(const BigReal, Tensor *virial, bool forces=false)
Definition: HomePatch.C:2962
void integrate(int)
Definition: Sequencer.C:218
static int velocityNeeded(int)
Definition: Output.C:365
void scaleVelocities(const BigReal velScale)
Definition: Sequencer.C:1219
#define GB1_COMPUTE_HOME_PRIORITY
Definition: Priorities.h:56
void addVelocityToPosition(BigReal)
Definition: Sequencer.C:1943
SubmitReduction * reduction
Definition: Sequencer.h:134
BigReal xz
Definition: Tensor.h:17
SubmitReduction * min_reduction
Definition: Sequencer.h:39
SimpleBroadcastObject< int > traceBarrier
Definition: Broadcasts.h:84
BigReal accelMDLastStep
int eventEndOfTimeStep
Definition: Node.C:285
void maximumMove(BigReal)
Definition: Sequencer.C:2027
int marginViolations
Definition: HomePatch.h:333
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
Bool is_atom_movdragged(int atomnum) const
Definition: Molecule.h:1214
#define BOLTZMANN
Definition: common.h:47
static int coordinateNeeded(int)
Definition: Output.C:198
Definition: Node.h:78
unsigned int atomFixed
Definition: NamdTypes.h:96
void cycleBarrier(int, int)
Definition: Sequencer.C:2889
#define FILE_OUTPUT
Definition: Output.h:25
Position fixedPosition
Definition: NamdTypes.h:102
Lattice & lattice
Definition: Patch.h:126
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:41
SimpleBroadcastObject< Vector > momentumCorrection
Definition: Broadcasts.h:79
void addRotDragToPosition(BigReal)
Definition: Sequencer.C:737
static PatchMap * Object()
Definition: PatchMap.h:27
void saveForce(const int ftag=Results::normal)
Definition: Sequencer.C:1893
void registerIDsFullAtom(const FullAtom *begin, const FullAtom *end)
Definition: AtomMap.C:50
#define EVAL_MEASURE
Definition: Output.h:27
Definition: Vector.h:64
void langevinVelocitiesBBK2(BigReal)
Definition: Sequencer.C:1427
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
int slowFreq
Definition: Sequencer.h:107
void newMinimizeDirection(BigReal)
Definition: Sequencer.C:897
void newMinimizePosition(BigReal)
Definition: Sequencer.C:956
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2381
void addForceToMomentum(FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms)
Definition: HomePatch.C:1932
BigReal reassignTemp
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
void startWork(const LDObjHandle &handle)
void langevinVelocitiesBBK1(BigReal)
Definition: Sequencer.C:1354
BigReal z
Definition: Vector.h:66
char const *const NamdProfileEventStr[]
Position position
Definition: NamdTypes.h:53
BigReal rotDragGlobVel
BigReal yz
Definition: Tensor.h:18
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
Definition: Tensor.h:255
int berendsenPressureFreq
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
void rattle1(BigReal, int)
Definition: Sequencer.C:1970
void saveTotalForces(HomePatch *)
if(ComputeNonbondedUtil::goMethod==2)
void submitVelocities(int seq, int zero, FullAtomList &a)
SimpleBroadcastObject< BigReal > adaptTemperature
Definition: Broadcasts.h:86
void rebalanceLoad(int timestep)
Definition: Sequencer.C:2878
void submitHalfstep(int)
Definition: Sequencer.C:2083
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:87
int doLoweAndersen
Definition: PatchTypes.h:26
Velocity velocity
Definition: NamdTypes.h:101
int pressureProfileSlabs
void minimizeMoveDownhill(BigReal fmax2)
Definition: Sequencer.C:875
void addForceToMomentum(BigReal, const int ftag=Results::normal, const int useSaved=0)
Definition: Sequencer.C:1904
BigReal length(void) const
Definition: Vector.h:169
void pauseWork(const LDObjHandle &handle)
void langevinPiston(int)
Definition: Sequencer.C:1598
Vector origin() const
Definition: Lattice.h:262
SimpleBroadcastObject< BigReal > tcoupleCoefficient
Definition: Broadcasts.h:76
static int forceNeeded(int)
Definition: Output.C:460
void exchangeCheckpoint(int scriptTask, int &bpc)
Definition: HomePatch.C:3477
AtomMapper * atomMapper
Definition: Patch.h:152
unsigned int groupFixed
Definition: NamdTypes.h:97
Bool pairInteractionOn
LDObjHandle ldObjHandle
Definition: HomePatch.h:479
void split(int iStream, int numStreams)
Definition: Random.h:77
void addForceToMomentum3(const BigReal timestep1, const int ftag1, const int useSaved1, const BigReal timestep2, const int ftag2, const int useSaved2, const BigReal timestep3, const int ftag3, const int useSaved3)
Definition: Sequencer.C:1919
Flags flags
Definition: Patch.h:127
SimpleBroadcastObject< BigReal > stochRescaleCoefficient
Definition: Broadcasts.h:77
void revert(void)
Definition: HomePatch.C:3459
void submitCollections(int step, int zeroVel=0)
Definition: Sequencer.C:2643
Charge charge
Definition: NamdTypes.h:54
BigReal calcKineticEnergy()
Definition: Sequencer.C:1227
#define SEQ_STK_SZ
Definition: Thread.h:11
void adaptTempUpdate(int)
Definition: Sequencer.C:1740
Bool langevin_useBAOAB
void rescale(Tensor factor)
Definition: Lattice.h:60
#define TIMER_START(T, TYPE)
Definition: HomePatch.h:265
#define NAIVE
Definition: SimParameters.h:49
void calcFixVirial(Tensor &fixVirialNormal, Tensor &fixVirialNbond, Tensor &fixVirialSlow, Vector &fixForceNormal, Vector &fixForceNbond, Vector &fixForceSlow)
Definition: Sequencer.C:2242
Definition: Random.h:37
Bool is_atom_rotdragged(int atomnum) const
Definition: Molecule.h:1230
#define NAMD_PROFILE_START()
void awaken(void)
Definition: Sequencer.h:29
#define NAMD_EVENT_START(eon, id)
int pairlistsAge
Definition: Sequencer.h:43
#define COMPUTE_HOME_PRIORITY
Definition: Priorities.h:76
BigReal rescaleTemp
void NAMD_bug(const char *err_msg)
Definition: common.C:123
#define TIMER_REPORT(T)
Definition: HomePatch.h:268
void multigratorPressure(int step, int callNumber)
Definition: Sequencer.C:1084
int doEnergy
Definition: PatchTypes.h:20
void berendsenPressure(int)
Definition: Sequencer.C:1534
void submitMomentum(int step)
Definition: Sequencer.C:993
int doFullElectrostatics
Definition: PatchTypes.h:23
BigReal yx
Definition: Tensor.h:18
Bool adaptTempLangevin
iterator end(void)
Definition: ResizeArray.h:37
int rescaleVelocities_numTemps
Definition: Sequencer.h:87
#define WAT_SWM4
Definition: common.h:192
void submitLoadStats(int timestep)
Definition: HomePatch.C:3619
void runComputeObjects(int migration=1, int pairlists=0, int pressureStep=0)
Definition: Sequencer.C:2668
void rebalance(Sequencer *seq, PatchID id)
void rescaleaccelMD(int, int, int)
Definition: Sequencer.C:1697
SimpleBroadcastObject< Tensor > velocityRescaleTensor2
Definition: Broadcasts.h:72
gridSize z
int Bool
Definition: common.h:133
unsigned char get_ss_type(int anum) const
Definition: Molecule.h:1351
BigReal drudeBondLen
Sequencer(HomePatch *p)
Definition: Sequencer.C:66
BigReal langevinTemp
int ldbSteps
Definition: Sequencer.h:140
int numAtoms
Definition: Patch.h:144
MTSChoices MTSAlgorithm
#define NAMD_EVENT_RANGE_2(eon, id)
void run(void)
Definition: Sequencer.C:126
SimpleBroadcastObject< int > scriptBarrier
Definition: Broadcasts.h:83
BigReal scriptArg1
BigReal x
Definition: Vector.h:66
void scalePositionsVelocities(const Tensor &posScale, const Tensor &velScale)
Definition: Sequencer.C:1047
BigReal adaptTempT
Definition: Sequencer.h:82
int maxForceUsed
Definition: PatchTypes.h:31
SimpleBroadcastObject< BigReal > velocityRescaleFactor2
Definition: Broadcasts.h:73
int sequence
Definition: PatchTypes.h:18
#define D_MSG(t)
Definition: Debug.h:149
#define END_OF_RUN
Definition: Output.h:26
void traceBarrier(int)
Definition: Sequencer.C:2897
int doNonbonded
Definition: PatchTypes.h:22
void NAMD_die(const char *err_msg)
Definition: common.C:83
static LdbCoordinator * Object()
BigReal reassignIncr
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
#define TIMER_INIT_WIDTH(T, TYPE, WIDTH)
Definition: HomePatch.h:264
void submitPositions(int seq, FullAtomList &a, Lattice l, int prec)
int berendsenPressure_count
Definition: Sequencer.h:104
SimpleBroadcastObject< BigReal > velocityRescaleFactor
Definition: Broadcasts.h:68
SimpleBroadcastObject< BigReal > minimizeCoefficient
Definition: Broadcasts.h:78
void reassignVelocities(BigReal, int)
Definition: Sequencer.C:1756
SimpleBroadcastObject< Vector > accelMDRescaleFactor
Definition: Broadcasts.h:85
int hardWallDrude(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2023
ComputeGlobal * computeGlobalObject
Definition: ComputeMgr.h:97
void saveForce(const int ftag=Results::normal)
Definition: HomePatch.C:1335
Random * random
Definition: Sequencer.h:131
BigReal xx
Definition: Tensor.h:17
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
void langevinVelocities(BigReal)
Definition: Sequencer.C:1317
void submitForces(int seq, FullAtomList &a, int maxForceUsed, ForceList *f)
void hardWallDrude(BigReal, int)
Definition: Sequencer.C:1955
void checkpoint(void)
Definition: HomePatch.C:3449
BigReal zz
Definition: Tensor.h:19
#define TIMER_STOP(T, TYPE)
Definition: HomePatch.h:266
const Lattice * lattice
Definition: GlobalMaster.h:141
PatchID getPatchID()
Definition: Patch.h:114
void suspend(void)
Definition: Sequencer.C:136
void multigratorTemperature(int step, int callNumber)
Definition: Sequencer.C:1247
BigReal length2(void) const
Definition: Vector.h:173
unsigned int randomSeed
BigReal initialTemp
void reinitVelocities(void)
Definition: Sequencer.C:1788
int pressureProfileAtomTypes
int checkpoint_berendsenPressure_count
Definition: Sequencer.h:105
#define simParams
Definition: Output.C:127
int numPatches(void) const
Definition: PatchMap.h:59
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
#define NAMD_EVENT_START_EX(eon, id, str)
Position apply_transform(Position data, const Transform &t) const
Definition: Lattice.h:132
CollectionMgr *const collection
Definition: Sequencer.h:137
void stochRescaleVelocities(BigReal, int)
Definition: Sequencer.C:1872
const PatchID patchID
Definition: Patch.h:143
#define GB2_COMPUTE_HOME_PRIORITY
Definition: Priorities.h:64
Definition: Tensor.h:15
void get_movdrag_params(Vector &v, int atomnum) const
Definition: Molecule.h:1321
BigReal xy
Definition: Tensor.h:17
#define NAMD_PROFILE_STOP()
int doVirial
Definition: PatchTypes.h:21
BigReal y
Definition: Vector.h:66
void get_rotdrag_params(BigReal &v, Vector &a, Vector &p, int atomnum) const
Definition: Molecule.h:1327
virtual ~Sequencer(void)
Definition: Sequencer.C:108
BigReal movDragGlobVel
int getNumStepsToRun(void)
int doLCPO
Definition: PatchTypes.h:29
void newtonianVelocities(BigReal, const BigReal, const BigReal, const BigReal, const int, const int, const int)
Definition: Sequencer.C:1293
void rescaleSoluteCharges(BigReal)
Definition: Sequencer.C:1833
Real atomcharge(int anum) const
Definition: Molecule.h:1048
Mass mass
Definition: NamdTypes.h:108
Bool pressureProfileOn
void submitMinimizeReductions(int, BigReal fmax2)
Definition: Sequencer.C:2505
#define ADD_VECTOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:27
BigReal yy
Definition: Tensor.h:18
int doMomenta
Definition: Sequencer.h:121
#define TIMER_DONE(T)
Definition: HomePatch.h:267
#define PDBVELFACTOR
Definition: common.h:50
const_iterator const_begin(void) const
Definition: ResizeArray.h:39
#define TIMEFACTOR
Definition: common.h:48
Bool pairInteractionSelf
int multigratorPressureFreq
gridSize y
#define WAT_TIP4
Definition: common.h:191
BigReal maximumMove
#define SPECIAL_PATCH_ID
Definition: Sequencer.C:64
void correctMomentum(int step, BigReal drifttime)
Definition: Sequencer.C:1016
int pairlistsAreValid
Definition: Sequencer.h:42
int doGBIS
Definition: PatchTypes.h:28
int stochRescale_count
Definition: Sequencer.h:100
void submit(void)
Definition: ReductionMgr.h:323
Bool is_atom_exPressure(int atomnum) const
Definition: Molecule.h:1447
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void addForceToMomentum3(FullAtom *__restrict atom_arr, const Force *__restrict force_arr1, const Force *__restrict force_arr2, const Force *__restrict force_arr3, const BigReal dt1, const BigReal dt2, const BigReal dt3, int num_atoms)
Definition: HomePatch.C:1961
ComputeMgr * computeMgr
Definition: Node.h:169
int maxForceMerged
Definition: PatchTypes.h:32
BigReal reassignHold
SubmitReduction * reduction
Used to submit restraint energy as MISC.
void quenchVelocities()
Definition: Sequencer.C:984
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
gridSize x
void enableEarlyExit(void)
Definition: Node.C:1356
void submitReductions(int)
Definition: Sequencer.C:2262
SimpleBroadcastObject< Tensor > positionRescaleFactor2
Definition: Broadcasts.h:74
#define RIGID_NONE
Definition: SimParameters.h:77
SimParameters *const simParams
Definition: Sequencer.h:132
SimpleBroadcastObject< Tensor > velocityRescaleTensor
Definition: Broadcasts.h:71
BigReal zx
Definition: Tensor.h:19
CompAtomExtList pExt
Definition: Patch.h:174
Molecule * molecule
Definition: Node.h:176
void rescaleVelocitiesByFactor(BigReal)
Definition: Sequencer.C:1811
int doMolly
Definition: PatchTypes.h:24
Vector unit(void) const
Definition: Vector.h:182
Vector c() const
Definition: Lattice.h:254
int multigratorTemperatureFreq
void reloadCharges()
Definition: Sequencer.C:1821
#define FORCE_OUTPUT
Definition: Output.h:28
double BigReal
Definition: common.h:114
void minimize()
Definition: Sequencer.C:760
int step
Definition: PatchTypes.h:16
#define PATCH_PRIORITY(PID)
Definition: Priorities.h:25
for(int i=0;i< n1;++i)
iterator begin(void)
Definition: ResizeArray.h:36
BigReal drudeTemp
void exchangeAtoms(int scriptTask)
Definition: HomePatch.C:3573