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