Sequencer.C

Go to the documentation of this file.
00001 
00007 /*****************************************************************************
00008  * $Source: /home/cvs/namd/cvsroot/namd2/src/Sequencer.C,v $
00009  * $Author: jim $
00010  * $Date: 2016/08/26 19:40:32 $
00011  * $Revision: 1.1230 $
00012  *****************************************************************************/
00013 
00014 //for gbis debugging; print net force on each atom
00015 #define PRINT_FORCES 0
00016 
00017 #include "InfoStream.h"
00018 #include "Node.h"
00019 #include "SimParameters.h"
00020 #include "Sequencer.h"
00021 #include "HomePatch.h"
00022 #include "ReductionMgr.h"
00023 #include "CollectionMgr.h"
00024 #include "BroadcastObject.h"
00025 #include "Output.h"
00026 #include "Controller.h"
00027 #include "Broadcasts.h"
00028 #include "Molecule.h"
00029 #include "NamdOneTools.h"
00030 #include "LdbCoordinator.h"
00031 #include "Thread.h"
00032 #include "Random.h"
00033 #include "PatchMap.inl"
00034 #include "ComputeMgr.h"
00035 #include "ComputeGlobal.h"
00036 
00037 #define MIN_DEBUG_LEVEL 4
00038 //#define DEBUGM
00039 #include "Debug.h"
00040 
00041 #if USE_HPM
00042 #define START_HPM_STEP  1000
00043 #define STOP_HPM_STEP   1500
00044 #endif
00045 
00046 Sequencer::Sequencer(HomePatch *p) :
00047         simParams(Node::Object()->simParameters),
00048         patch(p),
00049         collection(CollectionMgr::Object()),
00050         ldbSteps(0)
00051 {
00052     broadcast = new ControllerBroadcasts(& patch->ldObjHandle);
00053     reduction = ReductionMgr::Object()->willSubmit(
00054                   simParams->accelMDOn ? REDUCTIONS_AMD : REDUCTIONS_BASIC );
00055     min_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MINIMIZER,1);
00056     if (simParams->pressureProfileOn) {
00057       int ntypes = simParams->pressureProfileAtomTypes;
00058       int nslabs = simParams->pressureProfileSlabs;
00059       pressureProfileReduction = 
00060         ReductionMgr::Object()->willSubmit(
00061                 REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
00062     } else {
00063       pressureProfileReduction = NULL;
00064     }
00065     if (simParams->multigratorOn) {
00066       multigratorReduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
00067     } else {
00068       multigratorReduction = NULL;
00069     }
00070     ldbCoordinator = (LdbCoordinator::Object());
00071     random = new Random(simParams->randomSeed);
00072     random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
00073 
00074     rescaleVelocities_numTemps = 0;
00075     berendsenPressure_count = 0;
00076 //    patch->write_tip4_props();
00077 }
00078 
00079 Sequencer::~Sequencer(void)
00080 {
00081     delete broadcast;
00082     delete reduction;
00083     delete min_reduction;
00084     if (pressureProfileReduction) delete pressureProfileReduction;
00085     delete random;
00086     if (multigratorReduction) delete multigratorReduction;
00087 }
00088 
00089 // Invoked by thread
00090 void Sequencer::threadRun(Sequencer* arg)
00091 {
00092     LdbCoordinator::Object()->startWork(arg->patch->ldObjHandle);
00093     arg->algorithm();
00094 }
00095 
00096 // Invoked by Node::run() via HomePatch::runSequencer()
00097 void Sequencer::run(void)
00098 {
00099     // create a Thread and invoke it
00100     DebugM(4, "::run() - this = " << this << "\n" );
00101     thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
00102     CthSetStrategyDefault(thread);
00103     priority = PATCH_PRIORITY(patch->getPatchID());
00104     awaken();
00105 }
00106 
00107 void Sequencer::suspend(void)
00108 {
00109     LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
00110     CthSuspend();
00111     LdbCoordinator::Object()->startWork(patch->ldObjHandle);
00112 }
00113 
00114 // Defines sequence of operations on a patch.  e.g. when
00115 // to push out information for Compute objects to consume
00116 // when to migrate atoms, when to add forces to velocity update.
00117 void Sequencer::algorithm(void)
00118 {
00119   int scriptTask;
00120   int scriptSeq = 0;
00121   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
00122     switch ( scriptTask ) {
00123       case SCRIPT_OUTPUT:
00124         submitCollections(FILE_OUTPUT);
00125         break;
00126       case SCRIPT_FORCEOUTPUT:
00127         submitCollections(FORCE_OUTPUT);
00128         break;
00129       case SCRIPT_MEASURE:
00130         submitCollections(EVAL_MEASURE);
00131         break;
00132       case SCRIPT_REINITVELS:
00133         reinitVelocities();
00134         break;
00135       case SCRIPT_RESCALEVELS:
00136         rescaleVelocitiesByFactor(simParams->scriptArg1);
00137         break;
00138       case SCRIPT_RELOADCHARGES:
00139         reloadCharges();
00140         break;
00141       case SCRIPT_CHECKPOINT:
00142         patch->checkpoint();
00143         checkpoint_berendsenPressure_count = berendsenPressure_count;
00144         break;
00145       case SCRIPT_REVERT:
00146         patch->revert();
00147         berendsenPressure_count = checkpoint_berendsenPressure_count;
00148         pairlistsAreValid = 0;
00149         break;
00150       case SCRIPT_CHECKPOINT_STORE:
00151       case SCRIPT_CHECKPOINT_LOAD:
00152       case SCRIPT_CHECKPOINT_SWAP:
00153       case SCRIPT_CHECKPOINT_FREE:
00154         patch->exchangeCheckpoint(scriptTask,berendsenPressure_count);
00155         break;
00156       case SCRIPT_ATOMSENDRECV:
00157       case SCRIPT_ATOMSEND:
00158       case SCRIPT_ATOMRECV:
00159         patch->exchangeAtoms(scriptTask);
00160         break;
00161       case SCRIPT_MINIMIZE:
00162         minimize();
00163         break;
00164       case SCRIPT_RUN:
00165       case SCRIPT_CONTINUE:
00166         integrate(scriptTask);
00167         break;
00168       default:
00169         NAMD_bug("Unknown task in Sequencer::algorithm");
00170     }
00171   }
00172   submitCollections(END_OF_RUN);
00173   terminate();
00174 }
00175 
00176 
00177 extern int eventEndOfTimeStep;
00178 
00179 void Sequencer::integrate(int scriptTask) {
00180     char traceNote[24];
00181     char tracePrefix[20];
00182     sprintf(tracePrefix, "p:%d,s:",patch->patchID);
00183 //    patch->write_tip4_props();
00184 
00185     int &step = patch->flags.step;
00186     step = simParams->firstTimestep;
00187 
00188     // drag switches
00189     const Bool rotDragOn = simParams->rotDragOn;
00190     const Bool movDragOn = simParams->movDragOn;
00191 
00192     const int commOnly = simParams->commOnly;
00193 
00194     int &maxForceUsed = patch->flags.maxForceUsed;
00195     int &maxForceMerged = patch->flags.maxForceMerged;
00196     maxForceUsed = Results::normal;
00197     maxForceMerged = Results::normal;
00198 
00199     const int numberOfSteps = simParams->N;
00200     const int stepsPerCycle = simParams->stepsPerCycle;
00201     const BigReal timestep = simParams->dt;
00202 
00203     // what MTS method?
00204     const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00205 
00206     const int nonbondedFrequency = simParams->nonbondedFrequency;
00207     slowFreq = nonbondedFrequency;
00208     const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00209     int &doNonbonded = patch->flags.doNonbonded;
00210     doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
00211     if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00212     if ( doNonbonded ) maxForceUsed = Results::nbond;
00213 
00214     // Do we do full electrostatics?
00215     const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00216     const int fullElectFrequency = simParams->fullElectFrequency;
00217     if ( dofull ) slowFreq = fullElectFrequency;
00218     const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00219     int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00220     doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
00221     if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00222                                         maxForceMerged = Results::slow;
00223     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00224 
00225     const Bool accelMDOn = simParams->accelMDOn;
00226     const Bool accelMDdihe = simParams->accelMDdihe;
00227     const Bool accelMDdual = simParams->accelMDdual;
00228     if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
00229 
00230     // Is adaptive tempering on?
00231     const Bool adaptTempOn = simParams->adaptTempOn;
00232     adaptTempT = simParams->initialTemp;
00233     if (simParams->langevinOn)
00234         adaptTempT = simParams->langevinTemp;
00235     else if (simParams->rescaleFreq > 0)
00236         adaptTempT = simParams->rescaleTemp;
00237         
00238 
00239     int &doMolly = patch->flags.doMolly;
00240     doMolly = simParams->mollyOn && doFullElectrostatics;
00241     // BEGIN LA
00242     int &doLoweAndersen = patch->flags.doLoweAndersen;
00243     doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00244     // END LA
00245 
00246     int &doGBIS = patch->flags.doGBIS;
00247     doGBIS = simParams->GBISOn;
00248 
00249     int &doLCPO = patch->flags.doLCPO;
00250     doLCPO = simParams->LCPOOn;
00251 
00252     int zeroMomentum = simParams->zeroMomentum;
00253     
00254     // Do we need to return forces to TCL script or Colvar module?
00255     int doTcl = simParams->tclForcesOn;
00256         int doColvars = simParams->colvarsOn;
00257     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00258 
00259     // Bother to calculate energies?
00260     int &doEnergy = patch->flags.doEnergy;
00261     int energyFrequency = simParams->outputEnergies;
00262 
00263     const int reassignFreq = simParams->reassignFreq;
00264 
00265     int &doVirial = patch->flags.doVirial;
00266     doVirial = 1;
00267 
00268   if ( scriptTask == SCRIPT_RUN ) {
00269 
00270 //    printf("Doing initial rattle\n");
00271     rattle1(0.,0);  // enforce rigid bond constraints on initial positions
00272 
00273     if (simParams->lonepairs) {
00274       patch->atomMapper->registerIDsFullAtom(
00275                 patch->atom.begin(),patch->atom.end());
00276     }
00277 
00278     if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00279        reassignVelocities(timestep,step);
00280     }
00281 
00282     doEnergy = ! ( step % energyFrequency );
00283     if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00284     //Update energy every timestep for adaptive tempering
00285     if ( adaptTempOn ) doEnergy=1;
00286     runComputeObjects(1,step<numberOfSteps); // must migrate here!
00287     rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD 
00288     adaptTempUpdate(step); // update adaptive tempering temperature
00289 
00290     if ( staleForces || doTcl || doColvars ) {
00291       if ( doNonbonded ) saveForce(Results::nbond);
00292       if ( doFullElectrostatics ) saveForce(Results::slow);
00293     }
00294     if ( ! commOnly ) {
00295       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
00296     }
00297     minimizationQuenchVelocity();
00298     rattle1(-timestep,0);
00299     submitHalfstep(step);
00300     if ( ! commOnly ) {
00301       newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
00302     }
00303     rattle1(timestep,1);
00304     if (doTcl || doColvars)  // include constraint forces
00305       computeGlobal->saveTotalForces(patch);
00306     submitHalfstep(step);
00307     if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00308     if ( ! commOnly ) {
00309       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
00310     }
00311     submitReductions(step);
00312     if(traceIsOn()){
00313         traceUserEvent(eventEndOfTimeStep);
00314         sprintf(traceNote, "%s%d",tracePrefix,step); 
00315         traceUserSuppliedNote(traceNote);
00316     }
00317     rebalanceLoad(step);
00318 
00319   } // scriptTask == SCRIPT_RUN
00320 
00321   bool doMultigratorRattle = false;
00322 
00323     for ( ++step; step <= numberOfSteps; ++step )
00324     {
00325 
00326       rescaleVelocities(step);
00327       tcoupleVelocities(timestep,step);
00328       berendsenPressure(step);
00329 
00330       if ( ! commOnly ) {
00331         newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); 
00332       }
00333 
00334       // We do RATTLE here if multigrator thermostat was applied in the previous step
00335       if (doMultigratorRattle) rattle1(timestep, 1);
00336 
00337       /* reassignment based on half-step velocities
00338          if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00339          addVelocityToPosition(0.5*timestep);
00340          reassignVelocities(timestep,step);
00341          addVelocityToPosition(0.5*timestep);
00342          rattle1(0.,0);
00343          rattle1(-timestep,0);
00344          addVelocityToPosition(-1.0*timestep);
00345          rattle1(timestep,0);
00346          } */
00347 
00348       maximumMove(timestep);
00349       if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) {
00350         if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00351         // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
00352         langevinVelocities(timestep);
00353         langevinPiston(step);
00354         if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00355       } else {
00356         // If Langevin is not used, take full time step directly instread of two half steps      
00357         if ( ! commOnly ) addVelocityToPosition(timestep); 
00358       }
00359 
00360       // impose hard wall potential for Drude bond length
00361       hardWallDrude(timestep, 1);
00362 
00363       minimizationQuenchVelocity();
00364 
00365       doNonbonded = !(step%nonbondedFrequency);
00366       doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00367 
00368       if ( zeroMomentum && doFullElectrostatics )
00369         correctMomentum(step,slowstep);
00370 
00371       submitHalfstep(step);
00372 
00373       doMolly = simParams->mollyOn && doFullElectrostatics;
00374       // BEGIN LA
00375       doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00376       // END LA
00377 
00378       maxForceUsed = Results::normal;
00379       if ( doNonbonded ) maxForceUsed = Results::nbond;
00380       if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00381       if ( accelMDOn && (accelMDdihe || accelMDdual))  maxForceUsed = Results::amdf;
00382 
00383       // Migrate Atoms on stepsPerCycle
00384       doEnergy = ! ( step % energyFrequency );
00385       if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00386       if ( adaptTempOn ) doEnergy=1; 
00387 
00388       // Multigrator 
00389       if (simParams->multigratorOn) {
00390         doVirial = (!(step % energyFrequency) || ((simParams->outputPressure > 0) && !(step % simParams->outputPressure)) 
00391           || !(step % simParams->multigratorPressureFreq));
00392         doKineticEnergy = (!(step % energyFrequency) || !(step % simParams->multigratorTemperatureFreq));
00393         doMomenta = (simParams->outputMomenta > 0) && !(step % simParams->outputMomenta);
00394       } else {
00395         doVirial = 1;
00396         doKineticEnergy = 1;
00397         doMomenta = 1;
00398       }
00399 
00400       runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00401 
00402       rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
00403      
00404       if ( staleForces || doTcl || doColvars ) {
00405         if ( doNonbonded ) saveForce(Results::nbond);
00406         if ( doFullElectrostatics ) saveForce(Results::slow);
00407       }
00408 
00409       // reassignment based on full-step velocities
00410       if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00411         reassignVelocities(timestep,step);
00412         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00413         rattle1(-timestep,0);
00414       }
00415 
00416       if ( ! commOnly ) {
00417         langevinVelocitiesBBK1(timestep);
00418         newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00419         langevinVelocitiesBBK2(timestep);
00420       }
00421 
00422       // add drag to each atom's positions
00423       if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00424       if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00425 
00426       rattle1(timestep,1);
00427       if (doTcl || doColvars)  // include constraint forces
00428         computeGlobal->saveTotalForces(patch);
00429 
00430       submitHalfstep(step);
00431       if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00432 
00433       if ( ! commOnly ) {
00434         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00435       }
00436 
00437         // rattle2(timestep,step);
00438 
00439         submitReductions(step);
00440         submitCollections(step);
00441        //Update adaptive tempering temperature
00442         adaptTempUpdate(step);
00443 
00444       // Multigrator temperature and pressure steps
00445       multigratorTemperature(step, 1);
00446       multigratorPressure(step, 1);
00447       multigratorPressure(step, 2);
00448       multigratorTemperature(step, 2);
00449       doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
00450 
00451 #if CYCLE_BARRIER
00452         cycleBarrier(!((step+1) % stepsPerCycle), step);
00453 #elif PME_BARRIER
00454         cycleBarrier(doFullElectrostatics, step);
00455 #elif  STEP_BARRIER
00456         cycleBarrier(1, step);
00457 #endif
00458 
00459          if(Node::Object()->specialTracing || simParams->statsOn){
00460                  int bstep = simParams->traceStartStep;
00461                  int estep = bstep + simParams->numTraceSteps;
00462                  if(step == bstep || step == estep){
00463                          traceBarrier(step);
00464                  }                       
00465          }
00466 
00467 #ifdef MEASURE_NAMD_WITH_PAPI    
00468          if(simParams->papiMeasure) {
00469                  int bstep = simParams->papiMeasureStartStep;
00470                  int estep = bstep + simParams->numPapiMeasureSteps;
00471                  if(step == bstep || step==estep) {
00472                          papiMeasureBarrier(step);
00473                  }
00474          }
00475 #endif
00476           
00477         if(traceIsOn()){
00478             traceUserEvent(eventEndOfTimeStep);
00479             sprintf(traceNote, "%s%d",tracePrefix,step); 
00480             traceUserSuppliedNote(traceNote);
00481         }
00482         rebalanceLoad(step);
00483 
00484 #if PME_BARRIER
00485         // a step before PME
00486         cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00487 #endif
00488 
00489 #if USE_HPM
00490         if(step == START_HPM_STEP)
00491           (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
00492 
00493         if(step == STOP_HPM_STEP)
00494           (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
00495 #endif
00496 
00497     }
00498 }
00499 
00500 // add moving drag to each atom's position
00501 void Sequencer::addMovDragToPosition(BigReal timestep) {
00502   FullAtom *atom = patch->atom.begin();
00503   int numAtoms = patch->numAtoms;
00504   Molecule *molecule = Node::Object()->molecule;   // need its methods
00505   const BigReal movDragGlobVel = simParams->movDragGlobVel;
00506   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00507   Vector movDragVel, dragIncrement;
00508   for ( int i = 0; i < numAtoms; ++i )
00509   {
00510     // skip if fixed atom or zero drag attribute
00511     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00512          || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00513     molecule->get_movdrag_params(movDragVel, atom[i].id);
00514     dragIncrement = movDragGlobVel * movDragVel * dt;
00515     atom[i].position += dragIncrement;
00516   }
00517 }
00518 
00519 // add rotating drag to each atom's position
00520 void Sequencer::addRotDragToPosition(BigReal timestep) {
00521   FullAtom *atom = patch->atom.begin();
00522   int numAtoms = patch->numAtoms;
00523   Molecule *molecule = Node::Object()->molecule;   // need its methods
00524   const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00525   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00526   BigReal rotDragVel, dAngle;
00527   Vector atomRadius;
00528   Vector rotDragAxis, rotDragPivot, dragIncrement;
00529   for ( int i = 0; i < numAtoms; ++i )
00530   {
00531     // skip if fixed atom or zero drag attribute
00532     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00533          || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00534     molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00535     dAngle = rotDragGlobVel * rotDragVel * dt;
00536     rotDragAxis /= rotDragAxis.length();
00537     atomRadius = atom[i].position - rotDragPivot;
00538     dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00539     atom[i].position += dragIncrement;
00540   }
00541 }
00542 
00543 void Sequencer::minimize() {
00544   const int numberOfSteps = simParams->N;
00545   const int stepsPerCycle = simParams->stepsPerCycle;
00546   int &step = patch->flags.step;
00547   step = simParams->firstTimestep;
00548 
00549   int &maxForceUsed = patch->flags.maxForceUsed;
00550   int &maxForceMerged = patch->flags.maxForceMerged;
00551   maxForceUsed = Results::normal;
00552   maxForceMerged = Results::normal;
00553   int &doNonbonded = patch->flags.doNonbonded;
00554   doNonbonded = 1;
00555   maxForceUsed = Results::nbond;
00556   maxForceMerged = Results::nbond;
00557   const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00558   int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00559   doFullElectrostatics = dofull;
00560   if ( dofull ) {
00561     maxForceMerged = Results::slow;
00562     maxForceUsed = Results::slow;
00563   }
00564   int &doMolly = patch->flags.doMolly;
00565   doMolly = simParams->mollyOn && doFullElectrostatics;
00566   // BEGIN LA
00567   int &doLoweAndersen = patch->flags.doLoweAndersen;
00568   doLoweAndersen = 0;
00569   // END LA
00570 
00571   int &doGBIS = patch->flags.doGBIS;
00572   doGBIS = simParams->GBISOn;
00573 
00574     int &doLCPO = patch->flags.doLCPO;
00575     doLCPO = simParams->LCPOOn;
00576 
00577     int doTcl = simParams->tclForcesOn;
00578         int doColvars = simParams->colvarsOn;
00579     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00580 
00581   int &doEnergy = patch->flags.doEnergy;
00582   doEnergy = 1;
00583 
00584   patch->rattle1(0.,0,0);  // enforce rigid bond constraints on initial positions
00585 
00586   if (simParams->lonepairs) {
00587     patch->atomMapper->registerIDsFullAtom(
00588                 patch->atom.begin(),patch->atom.end());
00589   }
00590 
00591   runComputeObjects(1,step<numberOfSteps); // must migrate here!
00592 
00593   if ( doTcl || doColvars ) {
00594     if ( doNonbonded ) saveForce(Results::nbond);
00595     if ( doFullElectrostatics ) saveForce(Results::slow);
00596     computeGlobal->saveTotalForces(patch);
00597   }
00598 
00599   BigReal fmax2 = TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00600 
00601   submitMinimizeReductions(step,fmax2);
00602   rebalanceLoad(step);
00603 
00604   int downhill = 1;  // start out just fixing bad contacts
00605   int minSeq = 0;
00606   for ( ++step; step <= numberOfSteps; ++step ) {
00607    BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00608    if ( downhill ) {
00609     if ( c ) minimizeMoveDownhill(fmax2);
00610     else {
00611       downhill = 0;
00612       fmax2 *= 10000.;
00613     }
00614    }
00615    if ( ! downhill ) {
00616     if ( ! c ) {  // new direction
00617       c = broadcast->minimizeCoefficient.get(minSeq++);
00618       newMinimizeDirection(c);  // v = c * v + f
00619       c = broadcast->minimizeCoefficient.get(minSeq++);
00620     }  // same direction
00621     newMinimizePosition(c);  // x = x + c * v
00622    }
00623 
00624     runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00625     if ( doTcl || doColvars ) {
00626       if ( doNonbonded ) saveForce(Results::nbond);
00627       if ( doFullElectrostatics ) saveForce(Results::slow);
00628       computeGlobal->saveTotalForces(patch);
00629     }
00630     submitMinimizeReductions(step,fmax2);
00631     submitCollections(step, 1);  // write out zeros for velocities
00632     rebalanceLoad(step);
00633   }
00634   quenchVelocities();  // zero out bogus velocity
00635 }
00636 
00637 // x = x + 0.1 * unit(f) for large f
00638 void Sequencer::minimizeMoveDownhill(BigReal fmax2) {
00639 
00640   FullAtom *a = patch->atom.begin();
00641   Force *f1 = patch->f[Results::normal].begin();  // includes nbond and slow
00642   int numAtoms = patch->numAtoms;
00643 
00644   for ( int i = 0; i < numAtoms; ++i ) {
00645     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
00646     Force f = f1[i];
00647     if ( f.length2() > fmax2 ) {
00648       a[i].position += ( 0.1 * f.unit() );
00649       int hgs = a[i].hydrogenGroupSize;  // 0 if not parent
00650       for ( int j=1; j<hgs; ++j ) {
00651         a[++i].position += ( 0.1 * f.unit() );
00652       }
00653     }
00654   }
00655 
00656   patch->rattle1(0.,0,0);
00657 }
00658 
00659 // v = c * v + f
00660 void Sequencer::newMinimizeDirection(BigReal c) {
00661   FullAtom *a = patch->atom.begin();
00662   Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
00663   const bool fixedAtomsOn = simParams->fixedAtomsOn;
00664   const bool drudeHardWallOn = simParams->drudeHardWallOn;
00665   int numAtoms = patch->numAtoms;
00666   BigReal maxv2 = 0.;
00667 
00668   for ( int i = 0; i < numAtoms; ++i ) {
00669     a[i].velocity *= c;
00670     a[i].velocity += f1[i];
00671     if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00672       a[i].velocity = a[i-1].velocity;
00673     }
00674     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00675     BigReal v2 = a[i].velocity.length2();
00676     if ( v2 > maxv2 ) maxv2 = v2;
00677   }
00678 
00679   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial); }
00680 
00681   maxv2 = 0.;
00682   for ( int i = 0; i < numAtoms; ++i ) {
00683     if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00684       a[i].velocity = a[i-1].velocity;
00685       if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00686     }
00687     BigReal v2 = a[i].velocity.length2();
00688     if ( v2 > maxv2 ) maxv2 = v2;
00689   }
00690 
00691   min_reduction->max(0,maxv2);
00692   min_reduction->submit();
00693 
00694   // prevent hydrogens from being left behind
00695   BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00696   // int adjustCount = 0;
00697   int hgs;
00698   for ( int i = 0; i < numAtoms; i += hgs ) {
00699     hgs = a[i].hydrogenGroupSize;
00700     BigReal minChildVel = a[i].velocity.length2();
00701     if ( minChildVel < fmax2 ) continue;
00702     int adjustChildren = 1;
00703     for ( int j = i+1; j < (i+hgs); ++j ) {
00704       if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
00705     }
00706     if ( adjustChildren ) {
00707       // if ( hgs > 1 ) ++adjustCount;
00708       for ( int j = i+1; j < (i+hgs); ++j ) {
00709         if (a[i].mass < 0.01) continue;  // lone pair
00710         a[j].velocity = a[i].velocity;
00711       }
00712     }
00713   }
00714   // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
00715 
00716 }
00717 
00718 // x = x + c * v
00719 void Sequencer::newMinimizePosition(BigReal c) {
00720   FullAtom *a = patch->atom.begin();
00721   int numAtoms = patch->numAtoms;
00722 
00723   for ( int i = 0; i < numAtoms; ++i ) {
00724     a[i].position += c * a[i].velocity;
00725   }
00726 
00727   if ( simParams->drudeHardWallOn ) {
00728     for ( int i = 1; i < numAtoms; ++i ) {
00729       if ( (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00730         a[i].position -= a[i-1].position;
00731       }
00732     }
00733   }
00734 
00735   patch->rattle1(0.,0,0);
00736 
00737   if ( simParams->drudeHardWallOn ) {
00738     for ( int i = 1; i < numAtoms; ++i ) {
00739       if ( (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00740         a[i].position += a[i-1].position;
00741       }
00742     }
00743   }
00744 }
00745 
00746 // v = 0
00747 void Sequencer::quenchVelocities() {
00748   FullAtom *a = patch->atom.begin();
00749   int numAtoms = patch->numAtoms;
00750 
00751   for ( int i = 0; i < numAtoms; ++i ) {
00752     a[i].velocity = 0;
00753   }
00754 }
00755 
00756 void Sequencer::submitMomentum(int step) {
00757 
00758   FullAtom *a = patch->atom.begin();
00759   const int numAtoms = patch->numAtoms;
00760 
00761   Vector momentum = 0;
00762   BigReal mass = 0;
00763 if ( simParams->zeroMomentumAlt ) {
00764   for ( int i = 0; i < numAtoms; ++i ) {
00765     momentum += a[i].mass * a[i].velocity;
00766     mass += 1.;
00767   }
00768 } else {
00769   for ( int i = 0; i < numAtoms; ++i ) {
00770     momentum += a[i].mass * a[i].velocity;
00771     mass += a[i].mass;
00772   }
00773 }
00774 
00775   ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
00776   reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
00777 }
00778 
00779 void Sequencer::correctMomentum(int step, BigReal drifttime) {
00780 
00781   if ( simParams->fixedAtomsOn )
00782     NAMD_die("Cannot zero momentum when fixed atoms are present.");
00783 
00784   const Vector dv = broadcast->momentumCorrection.get(step);
00785   const Vector dx = dv * ( drifttime / TIMEFACTOR );
00786 
00787   FullAtom *a = patch->atom.begin();
00788   const int numAtoms = patch->numAtoms;
00789 
00790 if ( simParams->zeroMomentumAlt ) {
00791   for ( int i = 0; i < numAtoms; ++i ) {
00792     BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00793     a[i].velocity += dv * rmass;
00794     a[i].position += dx * rmass;
00795   }
00796 } else {
00797   for ( int i = 0; i < numAtoms; ++i ) {
00798     a[i].velocity += dv;
00799     a[i].position += dx;
00800   }
00801 }
00802 
00803 }
00804 
00805 // --------- For Multigrator ---------
00806 void Sequencer::scalePositionsVelocities(const Tensor& posScale, const Tensor& velScale) {
00807   FullAtom *a = patch->atom.begin();
00808   int numAtoms = patch->numAtoms;
00809   Position origin = patch->lattice.origin();
00810   if ( simParams->fixedAtomsOn ) {
00811     NAMD_bug("Sequencer::scalePositionsVelocities, fixed atoms not implemented");
00812   }
00813   if ( simParams->useGroupPressure ) {
00814     int hgs;
00815     for ( int i = 0; i < numAtoms; i += hgs ) {
00816       hgs = a[i].hydrogenGroupSize;
00817       Position pos_cm(0.0, 0.0, 0.0);
00818       Velocity vel_cm(0.0, 0.0, 0.0);
00819       BigReal m_cm = 0.0;
00820       for (int j=0;j < hgs;++j) {
00821         m_cm += a[i+j].mass;
00822         pos_cm += a[i+j].mass*a[i+j].position;
00823         vel_cm += a[i+j].mass*a[i+j].velocity;
00824       }
00825       pos_cm /= m_cm;
00826       vel_cm /= m_cm;
00827       pos_cm -= origin;
00828       Position dpos = posScale*pos_cm;
00829       Velocity dvel = velScale*vel_cm;
00830       for (int j=0;j < hgs;++j) {
00831         a[i+j].position += dpos;
00832         a[i+j].velocity += dvel;
00833       }
00834     }
00835   } else {
00836     for ( int i = 0; i < numAtoms; i++) {
00837       a[i].position += posScale*(a[i].position-origin);
00838       a[i].velocity = velScale*a[i].velocity;
00839     }
00840   }
00841 }
00842 
00843 void Sequencer::multigratorPressure(int step, int callNumber) {
00844 // Calculate new positions, momenta, and volume using positionRescaleFactor and 
00845 // velocityRescaleTensor values returned from Controller::multigratorPressureCalcScale()
00846   if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
00847     FullAtom *a = patch->atom.begin();
00848     int numAtoms = patch->numAtoms;
00849 
00850     // Receive scaling factors from Controller
00851     Tensor scaleTensor    = (callNumber == 1) ? broadcast->positionRescaleFactor.get(step) : broadcast->positionRescaleFactor2.get(step);
00852     Tensor velScaleTensor = (callNumber == 1) ? broadcast->velocityRescaleTensor.get(step) : broadcast->velocityRescaleTensor2.get(step);
00853     Tensor posScaleTensor = scaleTensor;
00854     posScaleTensor -= Tensor::identity();
00855     if (simParams->useGroupPressure) {
00856       velScaleTensor -= Tensor::identity();
00857     }
00858 
00859     // Scale volume
00860     patch->lattice.rescale(scaleTensor);
00861     // Scale positions and velocities
00862     scalePositionsVelocities(posScaleTensor, velScaleTensor);
00863 
00864     if (!patch->flags.doFullElectrostatics) NAMD_bug("Sequencer::multigratorPressure, doFullElectrostatics must be true");
00865 
00866     // Calculate new forces
00867     // NOTE: We should not need to migrate here since any migration should have happened in the
00868     // previous call to runComputeObjects inside the MD loop in Sequencer::integrate()
00869     const int numberOfSteps = simParams->N;
00870     const int stepsPerCycle = simParams->stepsPerCycle;
00871     runComputeObjects(0 , step<numberOfSteps, 1);
00872 
00873     reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
00874     reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
00875 
00876     // Virials etc.
00877     Tensor virialNormal;
00878     Tensor momentumSqrSum;
00879     BigReal kineticEnergy = 0;
00880     if ( simParams->pairInteractionOn ) {
00881       if ( simParams->pairInteractionSelf ) {
00882         for ( int i = 0; i < numAtoms; ++i ) {
00883           if ( a[i].partition != 1 ) continue;
00884           kineticEnergy += a[i].mass * a[i].velocity.length2();
00885           virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
00886         }
00887       }
00888     } else {
00889       for ( int i = 0; i < numAtoms; ++i ) {
00890         if (a[i].mass < 0.01) continue;
00891         kineticEnergy += a[i].mass * a[i].velocity.length2();
00892         virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
00893       }
00894     }
00895     if (!simParams->useGroupPressure) momentumSqrSum = virialNormal;
00896     kineticEnergy *= 0.5;
00897     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
00898     ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virialNormal);
00899 
00900     if ( simParams->fixedAtomsOn ) {
00901       Tensor fixVirialNormal;
00902       Tensor fixVirialNbond;
00903       Tensor fixVirialSlow;
00904       Vector fixForceNormal = 0;
00905       Vector fixForceNbond = 0;
00906       Vector fixForceSlow = 0;
00907 
00908       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
00909 
00910       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, fixVirialNormal);
00911       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, fixVirialNbond);
00912       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, fixVirialSlow);
00913       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NORMAL, fixForceNormal);
00914       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NBOND, fixForceNbond);
00915       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_SLOW, fixForceSlow);
00916     }
00917 
00918     // Internal virial and group momentum
00919     Tensor intVirialNormal;
00920     Tensor intVirialNormal2;
00921     Tensor intVirialNbond;
00922     Tensor intVirialSlow;
00923     int hgs;
00924     for ( int i = 0; i < numAtoms; i += hgs ) {
00925       hgs = a[i].hydrogenGroupSize;
00926       int j;
00927       BigReal m_cm = 0;
00928       Position x_cm(0,0,0);
00929       Velocity v_cm(0,0,0);
00930       for ( j = i; j < (i+hgs); ++j ) {
00931         m_cm += a[j].mass;
00932         x_cm += a[j].mass * a[j].position;
00933         v_cm += a[j].mass * a[j].velocity;
00934       }
00935       if (simParams->useGroupPressure) momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
00936       x_cm /= m_cm;
00937       v_cm /= m_cm;
00938       if (simParams->fixedAtomsOn) NAMD_bug("Sequencer::multigratorPressure, simParams->fixedAtomsOn not implemented yet");
00939       if ( simParams->pairInteractionOn ) {
00940         if ( simParams->pairInteractionSelf ) {
00941           NAMD_bug("Sequencer::multigratorPressure, this part needs to be implemented correctly");
00942           for ( j = i; j < (i+hgs); ++j ) {
00943             if ( a[j].partition != 1 ) continue;
00944             BigReal mass = a[j].mass;
00945             Vector v = a[j].velocity;
00946             Vector dv = v - v_cm;
00947             intVirialNormal2.outerAdd (mass, v, dv);
00948             Vector dx = a[j].position - x_cm;
00949             intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
00950             intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
00951             intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
00952           }
00953         }
00954       } else {
00955         for ( j = i; j < (i+hgs); ++j ) {
00956           BigReal mass = a[j].mass;
00957           Vector v = a[j].velocity;
00958           Vector dv = v - v_cm;
00959           intVirialNormal2.outerAdd(mass, v, dv);
00960           Vector dx = a[j].position - x_cm;
00961           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
00962           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
00963           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
00964         }
00965       }
00966     }
00967 
00968     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
00969     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal2);
00970     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NBOND, intVirialNbond);
00971     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_SLOW, intVirialSlow);
00972     ADD_TENSOR_OBJECT(reduction, REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
00973 
00974     reduction->submit();
00975   }
00976 }
00977 
00978 void Sequencer::scaleVelocities(const BigReal velScale) {
00979   FullAtom *a = patch->atom.begin();
00980   int numAtoms = patch->numAtoms;
00981   for ( int i = 0; i < numAtoms; i++) {
00982     a[i].velocity *= velScale;
00983   }
00984 }
00985 
00986 BigReal Sequencer::calcKineticEnergy() {
00987   FullAtom *a = patch->atom.begin();
00988   int numAtoms = patch->numAtoms;
00989   BigReal kineticEnergy = 0.0;
00990   if ( simParams->pairInteractionOn ) {
00991     if ( simParams->pairInteractionSelf ) {
00992       for (int i = 0; i < numAtoms; ++i ) {
00993         if ( a[i].partition != 1 ) continue;
00994         kineticEnergy += a[i].mass * a[i].velocity.length2();
00995       }
00996     }
00997   } else {
00998     for (int i = 0; i < numAtoms; ++i ) {
00999       kineticEnergy += a[i].mass * a[i].velocity.length2();
01000     }
01001   }
01002   kineticEnergy *= 0.5;
01003   return kineticEnergy;
01004 }
01005 
01006 void Sequencer::multigratorTemperature(int step, int callNumber) {
01007   if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
01008     // Scale velocities
01009     BigReal velScale = (callNumber == 1) ? broadcast->velocityRescaleFactor.get(step) : broadcast->velocityRescaleFactor2.get(step);
01010     scaleVelocities(velScale);
01011     // Calculate new kineticEnergy
01012     BigReal kineticEnergy = calcKineticEnergy();
01013     multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY) += kineticEnergy;
01014     if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
01015       // If this is a pressure cycle, calculate new momentum squared sum
01016       FullAtom *a = patch->atom.begin();
01017       int numAtoms = patch->numAtoms;
01018       Tensor momentumSqrSum;
01019       if (simParams->useGroupPressure) {
01020         int hgs;
01021         for ( int i = 0; i < numAtoms; i += hgs ) {
01022           hgs = a[i].hydrogenGroupSize;
01023           int j;
01024           BigReal m_cm = 0;
01025           Position x_cm(0,0,0);
01026           Velocity v_cm(0,0,0);
01027           for ( j = i; j < (i+hgs); ++j ) {
01028             m_cm += a[j].mass;
01029             x_cm += a[j].mass * a[j].position;
01030             v_cm += a[j].mass * a[j].velocity;
01031           }
01032           momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
01033         }
01034       } else {
01035         for ( int i = 0; i < numAtoms; i++) {
01036           momentumSqrSum.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01037         }
01038       }
01039       ADD_TENSOR_OBJECT(multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
01040     }
01041     // Submit reductions (kineticEnergy and, if applicable, momentumSqrSum)
01042     multigratorReduction->submit();
01043 
01044   }
01045 }
01046 // --------- End Multigrator ---------
01047 
01048 void Sequencer::newtonianVelocities(BigReal stepscale, const BigReal timestep, 
01049                                     const BigReal nbondstep, 
01050                                     const BigReal slowstep, 
01051                                     const int staleForces, 
01052                                     const int doNonbonded,
01053                                     const int doFullElectrostatics)
01054 {
01055   // Deterministic velocity update, account for multigrator
01056   if (staleForces || (doNonbonded && doFullElectrostatics)) {
01057     addForceToMomentum3(stepscale*timestep, Results::normal, 0,
01058                         stepscale*nbondstep, Results::nbond, staleForces,
01059                         stepscale*slowstep, Results::slow, staleForces);
01060   } else {
01061     addForceToMomentum(stepscale*timestep);
01062     if (staleForces || doNonbonded)
01063       addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces);
01064     if (staleForces || doFullElectrostatics)
01065       addForceToMomentum(stepscale*slowstep, Results::slow, staleForces);
01066   }
01067 }
01068 
01069 void Sequencer::langevinVelocities(BigReal dt_fs)
01070 {
01071 // This routine is used for the BAOAB integrator,
01072 // Ornstein-Uhlenbeck exact solve for the O-part.
01073 // See B. Leimkuhler and C. Matthews, AMRX (2012)
01074 // Routine originally written by JPhillips, with fresh errors by CMatthews June2012
01075 
01076   if ( simParams->langevinOn && simParams->langevin_useBAOAB )
01077   {
01078     FullAtom *a = patch->atom.begin();
01079     int numAtoms = patch->numAtoms;
01080     Molecule *molecule = Node::Object()->molecule;
01081     BigReal dt = dt_fs * 0.001;  // convert to ps
01082     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
01083     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
01084     {
01085         kbT = BOLTZMANN*adaptTempT;
01086     }
01087 
01088     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01089     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01090 
01091     for ( int i = 0; i < numAtoms; ++i )
01092     {
01093       BigReal dt_gamma = dt * a[i].langevinParam;
01094       if ( ! dt_gamma ) continue;
01095 
01096       BigReal f1 = exp( -dt_gamma );
01097       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT * 
01098                          ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
01099       a[i].velocity *= f1;
01100       a[i].velocity += f2 * random->gaussian_vector();
01101     }
01102   }
01103 }
01104 
01105 void Sequencer::langevinVelocitiesBBK1(BigReal dt_fs)
01106 {
01107   if ( simParams->langevinOn && !simParams->langevin_useBAOAB )
01108   {
01109     FullAtom *a = patch->atom.begin();
01110     int numAtoms = patch->numAtoms;
01111     Molecule *molecule = Node::Object()->molecule;
01112     BigReal dt = dt_fs * 0.001;  // convert to ps
01113     int i;
01114 
01115     if (simParams->drudeOn) {
01116       for (i = 0;  i < numAtoms;  i++) {
01117 
01118         if (i < numAtoms-1 &&
01119             a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01120           //printf("*** Found Drude particle %d\n", a[i+1].id);
01121           // i+1 is a Drude particle with parent i
01122 
01123           // convert from Cartesian coordinates to (COM,bond) coordinates
01124           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
01125           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
01126           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
01127           BigReal dt_gamma;
01128 
01129           // use Langevin damping factor i for v_com
01130           dt_gamma = dt * a[i].langevinParam;
01131           if (dt_gamma != 0.0) {
01132             v_com *= ( 1. - 0.5 * dt_gamma );
01133           }
01134 
01135           // use Langevin damping factor i+1 for v_bnd
01136           dt_gamma = dt * a[i+1].langevinParam;
01137           if (dt_gamma != 0.0) {
01138             v_bnd *= ( 1. - 0.5 * dt_gamma );
01139           }
01140 
01141           // convert back
01142           a[i].velocity = v_com - m * v_bnd;
01143           a[i+1].velocity = v_bnd + a[i].velocity;
01144 
01145           i++;  // +1 from loop, we've updated both particles
01146         }
01147         else {
01148           BigReal dt_gamma = dt * a[i].langevinParam;
01149           if ( ! dt_gamma ) continue;
01150 
01151           a[i].velocity *= ( 1. - 0.5 * dt_gamma );
01152         }
01153 
01154       } // end for
01155     } // end if drudeOn
01156     else {
01157 
01158       for ( i = 0; i < numAtoms; ++i )
01159       {
01160         BigReal dt_gamma = dt * a[i].langevinParam;
01161         if ( ! dt_gamma ) continue;
01162 
01163         a[i].velocity *= ( 1. - 0.5 * dt_gamma );
01164       }
01165 
01166     } // end else
01167 
01168   } // end if langevinOn
01169 }
01170 
01171 
01172 void Sequencer::langevinVelocitiesBBK2(BigReal dt_fs)
01173 {
01174   if ( simParams->langevinOn && !simParams->langevin_useBAOAB ) 
01175   {
01176     rattle1(dt_fs,1);  // conserve momentum if gammas differ
01177 
01178     FullAtom *a = patch->atom.begin();
01179     int numAtoms = patch->numAtoms;
01180     Molecule *molecule = Node::Object()->molecule;
01181     BigReal dt = dt_fs * 0.001;  // convert to ps
01182     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
01183     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
01184     {
01185         kbT = BOLTZMANN*adaptTempT;
01186     }
01187     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01188     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01189     int i;
01190 
01191     if (simParams->drudeOn) {
01192       BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp);  // drude bond Temp
01193 
01194       for (i = 0;  i < numAtoms;  i++) {
01195 
01196         if (i < numAtoms-1 &&
01197             a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01198           //printf("*** Found Drude particle %d\n", a[i+1].id);
01199           // i+1 is a Drude particle with parent i
01200 
01201           // convert from Cartesian coordinates to (COM,bond) coordinates
01202           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
01203           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
01204           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
01205           BigReal dt_gamma;
01206 
01207           // use Langevin damping factor i for v_com
01208           dt_gamma = dt * a[i].langevinParam;
01209           if (dt_gamma != 0.0) {
01210             BigReal mass = a[i].mass + a[i+1].mass;
01211             v_com += random->gaussian_vector() *
01212               sqrt( 2 * dt_gamma * kbT *
01213                   ( a[i].partition ? tempFactor : 1.0 ) / mass );
01214             v_com /= ( 1. + 0.5 * dt_gamma );
01215           }
01216 
01217           // use Langevin damping factor i+1 for v_bnd
01218           dt_gamma = dt * a[i+1].langevinParam;
01219           if (dt_gamma != 0.0) {
01220             BigReal mass = a[i+1].mass * (1. - m);
01221             v_bnd += random->gaussian_vector() *
01222               sqrt( 2 * dt_gamma * kbT_bnd *
01223                   ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
01224             v_bnd /= ( 1. + 0.5 * dt_gamma );
01225           }
01226 
01227           // convert back
01228           a[i].velocity = v_com - m * v_bnd;
01229           a[i+1].velocity = v_bnd + a[i].velocity;
01230 
01231           i++;  // +1 from loop, we've updated both particles
01232         }
01233         else { 
01234           BigReal dt_gamma = dt * a[i].langevinParam;
01235           if ( ! dt_gamma ) continue;
01236 
01237           a[i].velocity += random->gaussian_vector() *
01238             sqrt( 2 * dt_gamma * kbT *
01239                 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
01240           a[i].velocity /= ( 1. + 0.5 * dt_gamma );
01241         }
01242 
01243       } // end for
01244     } // end if drudeOn
01245     else {
01246 
01247       for ( i = 0; i < numAtoms; ++i )
01248       {
01249         BigReal dt_gamma = dt * a[i].langevinParam;
01250         if ( ! dt_gamma ) continue;
01251 
01252         a[i].velocity += random->gaussian_vector() *
01253           sqrt( 2 * dt_gamma * kbT *
01254               ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
01255         a[i].velocity /= ( 1. + 0.5 * dt_gamma );
01256       }
01257 
01258     } // end else
01259 
01260   } // end if langevinOn
01261 }
01262 
01263 
01264 void Sequencer::berendsenPressure(int step)
01265 {
01266   if ( simParams->berendsenPressureOn ) {
01267   berendsenPressure_count += 1;
01268   const int freq = simParams->berendsenPressureFreq;
01269   if ( ! (berendsenPressure_count % freq ) ) {
01270    berendsenPressure_count = 0;
01271    FullAtom *a = patch->atom.begin();
01272    int numAtoms = patch->numAtoms;
01273    Tensor factor = broadcast->positionRescaleFactor.get(step);
01274    patch->lattice.rescale(factor);
01275    if ( simParams->useGroupPressure )
01276    {
01277     int hgs;
01278     for ( int i = 0; i < numAtoms; i += hgs ) {
01279       int j;
01280       hgs = a[i].hydrogenGroupSize;
01281       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
01282         for ( j = i; j < (i+hgs); ++j ) {
01283           a[j].position = patch->lattice.apply_transform(
01284                                 a[j].fixedPosition,a[j].transform);
01285         }
01286         continue;
01287       }
01288       BigReal m_cm = 0;
01289       Position x_cm(0,0,0);
01290       for ( j = i; j < (i+hgs); ++j ) {
01291         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01292         m_cm += a[j].mass;
01293         x_cm += a[j].mass * a[j].position;
01294       }
01295       x_cm /= m_cm;
01296       Position new_x_cm = x_cm;
01297       patch->lattice.rescale(new_x_cm,factor);
01298       Position delta_x_cm = new_x_cm - x_cm;
01299       for ( j = i; j < (i+hgs); ++j ) {
01300         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01301           a[j].position = patch->lattice.apply_transform(
01302                                 a[j].fixedPosition,a[j].transform);
01303           continue;
01304         }
01305         a[j].position += delta_x_cm;
01306       }
01307     }
01308    }
01309    else
01310    {
01311     for ( int i = 0; i < numAtoms; ++i )
01312     {
01313       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01314         a[i].position = patch->lattice.apply_transform(
01315                                 a[i].fixedPosition,a[i].transform);
01316         continue;
01317       }
01318       patch->lattice.rescale(a[i].position,factor);
01319     }
01320    }
01321   }
01322   } else {
01323     berendsenPressure_count = 0;
01324   }
01325 }
01326 
01327 void Sequencer::langevinPiston(int step)
01328 {
01329   if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
01330   {
01331    FullAtom *a = patch->atom.begin();
01332    int numAtoms = patch->numAtoms;
01333    Tensor factor = broadcast->positionRescaleFactor.get(step);
01334    // JCP FIX THIS!!!
01335    Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
01336    patch->lattice.rescale(factor);
01337    Molecule *mol = Node::Object()->molecule;
01338    if ( simParams->useGroupPressure )
01339    {
01340     int hgs;
01341     for ( int i = 0; i < numAtoms; i += hgs ) {
01342       int j;
01343       hgs = a[i].hydrogenGroupSize;
01344       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
01345         for ( j = i; j < (i+hgs); ++j ) {
01346           a[j].position = patch->lattice.apply_transform(
01347                                 a[j].fixedPosition,a[j].transform);
01348         }
01349         continue;
01350       }
01351       BigReal m_cm = 0;
01352       Position x_cm(0,0,0);
01353       Velocity v_cm(0,0,0);
01354       for ( j = i; j < (i+hgs); ++j ) {
01355         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01356         m_cm += a[j].mass;
01357         x_cm += a[j].mass * a[j].position;
01358         v_cm += a[j].mass * a[j].velocity;
01359       }
01360       x_cm /= m_cm;
01361       Position new_x_cm = x_cm;
01362       patch->lattice.rescale(new_x_cm,factor);
01363       Position delta_x_cm = new_x_cm - x_cm;
01364       v_cm /= m_cm;
01365       Velocity delta_v_cm;
01366       delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
01367       delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
01368       delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
01369       for ( j = i; j < (i+hgs); ++j ) {
01370         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01371           a[j].position = patch->lattice.apply_transform(
01372                                 a[j].fixedPosition,a[j].transform);
01373           continue;
01374         }
01375         if ( mol->is_atom_exPressure(a[j].id) ) continue;
01376         a[j].position += delta_x_cm;
01377         a[j].velocity += delta_v_cm;
01378       }
01379     }
01380    }
01381    else
01382    {
01383     for ( int i = 0; i < numAtoms; ++i )
01384     {
01385       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01386         a[i].position = patch->lattice.apply_transform(
01387                                 a[i].fixedPosition,a[i].transform);
01388         continue;
01389       }
01390       if ( mol->is_atom_exPressure(a[i].id) ) continue;
01391       patch->lattice.rescale(a[i].position,factor);
01392       a[i].velocity.x *= velFactor.x;
01393       a[i].velocity.y *= velFactor.y;
01394       a[i].velocity.z *= velFactor.z;
01395     }
01396    }
01397   }
01398 }
01399 
01400 void Sequencer::rescaleVelocities(int step)
01401 {
01402   const int rescaleFreq = simParams->rescaleFreq;
01403   if ( rescaleFreq > 0 ) {
01404     FullAtom *a = patch->atom.begin();
01405     int numAtoms = patch->numAtoms;
01406     ++rescaleVelocities_numTemps;
01407     if ( rescaleVelocities_numTemps == rescaleFreq ) {
01408       BigReal factor = broadcast->velocityRescaleFactor.get(step);
01409       for ( int i = 0; i < numAtoms; ++i )
01410       {
01411         a[i].velocity *= factor;
01412       }
01413       rescaleVelocities_numTemps = 0;
01414     }
01415   }
01416 }
01417 
01418 void Sequencer::rescaleaccelMD (int step, int doNonbonded, int doFullElectrostatics)
01419 {
01420    if (!simParams->accelMDOn) return;
01421    if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
01422 
01423    Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
01424    const BigReal factor_dihe = accelMDfactor[0];
01425    const BigReal factor_tot  = accelMDfactor[1];
01426    const int numAtoms = patch->numAtoms;
01427 
01428    if (simParams->accelMDdihe && factor_tot <1 )
01429        NAMD_die("accelMD broadcasting error!\n");
01430    if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
01431        NAMD_die("accelMD broadcasting error!\n");
01432 
01433    if (simParams->accelMDdihe && factor_dihe < 1) {
01434         for (int i = 0; i < numAtoms; ++i)
01435           if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01436               patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);         
01437    }
01438 
01439    if ( !simParams->accelMDdihe && factor_tot < 1) {
01440         for (int i = 0; i < numAtoms; ++i)
01441             patch->f[Results::normal][i] *= factor_tot;
01442         if (doNonbonded) {
01443             for (int i = 0; i < numAtoms; ++i)
01444                  patch->f[Results::nbond][i] *= factor_tot;
01445         }
01446         if (doFullElectrostatics) {
01447             for (int i = 0; i < numAtoms; ++i)
01448                  patch->f[Results::slow][i] *= factor_tot;
01449         }
01450    }
01451 
01452    if (simParams->accelMDdual && factor_dihe < 1) {
01453         for (int i = 0; i < numAtoms; ++i)
01454            if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01455                patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
01456    }
01457 
01458 }
01459 
01460 void Sequencer::adaptTempUpdate(int step)
01461 {
01462    //check if adaptive tempering is enabled and in the right timestep range
01463    if (!simParams->adaptTempOn) return;
01464    if ( (step < simParams->adaptTempFirstStep ) || 
01465      ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) {
01466         if (simParams->langevinOn) // restore langevin temperature
01467             adaptTempT = simParams->langevinTemp;
01468         return;
01469    }
01470    // Get Updated Temperature
01471    if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
01472     adaptTempT = broadcast->adaptTemperature.get(step);
01473 }
01474 
01475 void Sequencer::reassignVelocities(BigReal timestep, int step)
01476 {
01477   const int reassignFreq = simParams->reassignFreq;
01478   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01479     FullAtom *a = patch->atom.begin();
01480     int numAtoms = patch->numAtoms;
01481     BigReal newTemp = simParams->reassignTemp;
01482     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01483     if ( simParams->reassignIncr > 0.0 ) {
01484       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01485         newTemp = simParams->reassignHold;
01486     } else {
01487       if ( newTemp < simParams->reassignHold )
01488         newTemp = simParams->reassignHold;
01489     }
01490     BigReal kbT = BOLTZMANN * newTemp;
01491 
01492     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01493     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01494 
01495     for ( int i = 0; i < numAtoms; ++i )
01496     {
01497       a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
01498         sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
01499           * random->gaussian_vector() );
01500     }
01501   } else {
01502     NAMD_bug("Sequencer::reassignVelocities called improperly!");
01503   }
01504 }
01505 
01506 void Sequencer::reinitVelocities(void)
01507 {
01508   FullAtom *a = patch->atom.begin();
01509   int numAtoms = patch->numAtoms;
01510   BigReal newTemp = simParams->initialTemp;
01511   BigReal kbT = BOLTZMANN * newTemp;
01512 
01513   int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01514   BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01515 
01516   for ( int i = 0; i < numAtoms; ++i )
01517   {
01518     a[i].velocity = ( ( (simParams->fixedAtomsOn && a[i].atomFixed) || a[i].mass <= 0.) ? Vector(0,0,0) :
01519       sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
01520         * random->gaussian_vector() );
01521     if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass >= 0.001 ) {
01522       a[i+1].velocity = a[i].velocity;  // zero is good enough
01523       ++i;
01524     }
01525   }
01526 }
01527 
01528 void Sequencer::rescaleVelocitiesByFactor(BigReal factor)
01529 {
01530   FullAtom *a = patch->atom.begin();
01531   int numAtoms = patch->numAtoms;
01532   for ( int i = 0; i < numAtoms; ++i )
01533   {
01534     a[i].velocity *= factor;
01535   }
01536 }
01537 
01538 void Sequencer::reloadCharges()
01539 {
01540   FullAtom *a = patch->atom.begin();
01541   int numAtoms = patch->numAtoms;
01542   Molecule *molecule = Node::Object()->molecule;
01543   for ( int i = 0; i < numAtoms; ++i )
01544   {
01545     a[i].charge = molecule->atomcharge(a[i].id);
01546   }
01547 }
01548 
01549 void Sequencer::tcoupleVelocities(BigReal dt_fs, int step)
01550 {
01551   if ( simParams->tCoupleOn )
01552   {
01553     FullAtom *a = patch->atom.begin();
01554     int numAtoms = patch->numAtoms;
01555     BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
01556     Molecule *molecule = Node::Object()->molecule;
01557     BigReal dt = dt_fs * 0.001;  // convert to ps
01558     coefficient *= dt;
01559     for ( int i = 0; i < numAtoms; ++i )
01560     {
01561       BigReal f1 = exp( coefficient * a[i].langevinParam );
01562       a[i].velocity *= f1;
01563     }
01564   }
01565 }
01566 
01567 void Sequencer::saveForce(const int ftag)
01568 {
01569   patch->saveForce(ftag);
01570 }
01571 
01572 void Sequencer::addForceToMomentum(BigReal dt, const int ftag, const int useSaved)
01573 {
01574 #if CMK_BLUEGENEL
01575   CmiNetworkProgressAfter (0);
01576 #endif
01577   patch->addForceToMomentum(dt,ftag,useSaved);
01578 }
01579 
01580 void Sequencer::addForceToMomentum3(const BigReal timestep1, const int ftag1, const int useSaved1,
01581         const BigReal timestep2, const int ftag2, const int useSaved2,
01582         const BigReal timestep3, const int ftag3, const int useSaved3) {
01583 #if CMK_BLUEGENEL
01584   CmiNetworkProgressAfter (0);
01585 #endif
01586   patch->addForceToMomentum3(timestep1,ftag1,useSaved1, timestep2,ftag2,useSaved2, timestep3,ftag3,useSaved3);  
01587 }
01588 
01589 void Sequencer::addVelocityToPosition(BigReal dt)
01590 {
01591 #if CMK_BLUEGENEL
01592   CmiNetworkProgressAfter (0);
01593 #endif
01594   patch->addVelocityToPosition(dt);
01595 }
01596 
01597 void Sequencer::hardWallDrude(BigReal dt, int pressure)
01598 {
01599   if ( simParams->drudeHardWallOn ) {
01600     Tensor virial;
01601     Tensor *vp = ( pressure ? &virial : 0 );
01602     if ( patch->hardWallDrude(dt, vp, pressureProfileReduction) ) {
01603       iout << iERROR << "Constraint failure in HardWallDrude(); "
01604         << "simulation may become unstable.\n" << endi;
01605       Node::Object()->enableEarlyExit();
01606       terminate();
01607     }
01608     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01609   }
01610 }
01611 
01612 void Sequencer::rattle1(BigReal dt, int pressure)
01613 {
01614   if ( simParams->rigidBonds != RIGID_NONE ) {
01615     Tensor virial;
01616     Tensor *vp = ( pressure ? &virial : 0 );
01617     if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
01618       iout << iERROR << 
01619         "Constraint failure; simulation has become unstable.\n" << endi;
01620       Node::Object()->enableEarlyExit();
01621       terminate();
01622     }
01623     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01624   }
01625 }
01626 
01627 // void Sequencer::rattle2(BigReal dt, int step)
01628 // {
01629 //   if ( simParams->rigidBonds != RIGID_NONE ) {
01630 //     Tensor virial;
01631 //     patch->rattle2(dt, &virial);
01632 //     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01633 //     // we need to add to alt and int virial because not included in forces
01634 // #ifdef ALTVIRIAL
01635 //     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01636 // #endif
01637 //     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,virial);
01638 //   }
01639 // }
01640 
01641 void Sequencer::maximumMove(BigReal timestep)
01642 {
01643   FullAtom *a = patch->atom.begin();
01644   int numAtoms = patch->numAtoms;
01645   if ( simParams->maximumMove ) {
01646     const BigReal dt = timestep / TIMEFACTOR;
01647     const BigReal maxvel = simParams->maximumMove / dt;
01648     const BigReal maxvel2 = maxvel * maxvel;
01649     for ( int i=0; i<numAtoms; ++i ) {
01650       if ( a[i].velocity.length2() > maxvel2 ) {
01651         a[i].velocity *= ( maxvel / a[i].velocity.length() );
01652       }
01653     }
01654   } else {
01655     const BigReal dt = timestep / TIMEFACTOR;
01656     const BigReal maxvel = simParams->cutoff / dt;
01657     const BigReal maxvel2 = maxvel * maxvel;
01658     int killme = 0;
01659     for ( int i=0; i<numAtoms; ++i ) {
01660       killme = killme || ( a[i].velocity.length2() > maxvel2 );
01661     }
01662     if ( killme ) {
01663       killme = 0;
01664       for ( int i=0; i<numAtoms; ++i ) {
01665         if ( a[i].velocity.length2() > maxvel2 ) {
01666           ++killme;
01667           iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
01668             << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
01669             << ( PDBVELFACTOR * maxvel ) << ", atom "
01670             << i << " of " << numAtoms << " on patch "
01671             << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
01672         }
01673       }
01674       iout << iERROR << 
01675         "Atoms moving too fast; simulation has become unstable ("
01676         << killme << " atoms on patch " << patch->patchID
01677         << " pe " << CkMyPe() << ").\n" << endi;
01678       Node::Object()->enableEarlyExit();
01679       terminate();
01680     }
01681   }
01682 }
01683 
01684 void Sequencer::minimizationQuenchVelocity(void)
01685 {
01686   if ( simParams->minimizeOn ) {
01687     FullAtom *a = patch->atom.begin();
01688     int numAtoms = patch->numAtoms;
01689     for ( int i=0; i<numAtoms; ++i ) {
01690       a[i].velocity = 0.;
01691     }
01692   }
01693 }
01694 
01695 void Sequencer::submitHalfstep(int step)
01696 {
01697   // velocity-dependent quantities *** ONLY ***
01698   // positions are not at half-step when called
01699   FullAtom *a = patch->atom.begin();
01700   int numAtoms = patch->numAtoms;
01701 
01702 #if CMK_BLUEGENEL
01703   CmiNetworkProgressAfter (0);
01704 #endif
01705 
01706   // For non-Multigrator doKineticEnergy = 1 always
01707   Tensor momentumSqrSum;
01708   if (doKineticEnergy || patch->flags.doVirial)
01709   {
01710     BigReal kineticEnergy = 0;
01711     Tensor virial;
01712     if ( simParams->pairInteractionOn ) {
01713       if ( simParams->pairInteractionSelf ) {
01714         for ( int i = 0; i < numAtoms; ++i ) {
01715           if ( a[i].partition != 1 ) continue;
01716           kineticEnergy += a[i].mass * a[i].velocity.length2();
01717           virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01718         }
01719       }
01720     } else {
01721       for ( int i = 0; i < numAtoms; ++i ) {
01722         if (a[i].mass < 0.01) continue;
01723         kineticEnergy += a[i].mass * a[i].velocity.length2();
01724         virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01725       }
01726     }
01727     
01728     if (simParams->multigratorOn && !simParams->useGroupPressure) {
01729       momentumSqrSum = virial;
01730     }
01731     kineticEnergy *= 0.5 * 0.5;
01732     reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
01733     virial *= 0.5;
01734     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01735 #ifdef ALTVIRIAL
01736     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01737 #endif
01738   }
01739  
01740   if (pressureProfileReduction) {
01741     int nslabs = simParams->pressureProfileSlabs;
01742     const Lattice &lattice = patch->lattice;
01743     BigReal idz = nslabs/lattice.c().z;
01744     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01745     int useGroupPressure = simParams->useGroupPressure;
01746 
01747     // Compute kinetic energy partition, possibly subtracting off
01748     // internal kinetic energy if group pressure is enabled.
01749     // Since the regular pressure is 1/2 mvv and the internal kinetic
01750     // term that is subtracted off for the group pressure is
01751     // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
01752     // 1/2 m * v * v_cm.  The factor of 1/2 is because submitHalfstep
01753     // gets called twice per timestep.
01754     int hgs;
01755     for (int i=0; i<numAtoms; i += hgs) {
01756       int j, ppoffset;
01757       hgs = a[i].hydrogenGroupSize;
01758       int partition = a[i].partition;
01759 
01760       BigReal m_cm = 0;
01761       Velocity v_cm(0,0,0);
01762       for (j=i; j< i+hgs; ++j) {
01763         m_cm += a[j].mass;
01764         v_cm += a[j].mass * a[j].velocity;
01765       }
01766       v_cm /= m_cm;
01767       for (j=i; j < i+hgs; ++j) {
01768         BigReal mass = a[j].mass;
01769         if (! (useGroupPressure && j != i)) {
01770           BigReal z = a[j].position.z;
01771           int slab = (int)floor((z-zmin)*idz);
01772           if (slab < 0) slab += nslabs;
01773           else if (slab >= nslabs) slab -= nslabs;
01774           ppoffset = 3*(slab + partition*nslabs);
01775         }
01776         BigReal wxx, wyy, wzz;
01777         if (useGroupPressure) {
01778           wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
01779           wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
01780           wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
01781         } else {
01782           wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
01783           wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
01784           wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
01785         }
01786         pressureProfileReduction->item(ppoffset  ) += wxx;
01787         pressureProfileReduction->item(ppoffset+1) += wyy;
01788         pressureProfileReduction->item(ppoffset+2) += wzz;
01789       }
01790     }
01791   } 
01792 
01793   // For non-Multigrator doKineticEnergy = 1 always
01794   if (doKineticEnergy || patch->flags.doVirial)
01795   {
01796     BigReal intKineticEnergy = 0;
01797     Tensor intVirialNormal;
01798 
01799     int hgs;
01800     for ( int i = 0; i < numAtoms; i += hgs ) {
01801 
01802 #if CMK_BLUEGENEL
01803       CmiNetworkProgress ();
01804 #endif
01805 
01806       hgs = a[i].hydrogenGroupSize;
01807       int j;
01808       BigReal m_cm = 0;
01809       Velocity v_cm(0,0,0);
01810       for ( j = i; j < (i+hgs); ++j ) {
01811         m_cm += a[j].mass;
01812         v_cm += a[j].mass * a[j].velocity;
01813       }
01814       if (simParams->multigratorOn && simParams->useGroupPressure) {
01815         momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
01816       }
01817       v_cm /= m_cm;
01818       if ( simParams->pairInteractionOn ) {
01819         if ( simParams->pairInteractionSelf ) {
01820           for ( j = i; j < (i+hgs); ++j ) {
01821             if ( a[j].partition != 1 ) continue;
01822             BigReal mass = a[j].mass;
01823             Vector v = a[j].velocity;
01824             Vector dv = v - v_cm;
01825             intKineticEnergy += mass * (v * dv);
01826             intVirialNormal.outerAdd (mass, v, dv);
01827           }
01828         }
01829       } else {
01830         for ( j = i; j < (i+hgs); ++j ) {
01831           BigReal mass = a[j].mass;
01832           Vector v = a[j].velocity;
01833           Vector dv = v - v_cm;
01834           intKineticEnergy += mass * (v * dv);
01835           intVirialNormal.outerAdd(mass, v, dv);
01836         }
01837       }
01838     }
01839 
01840     intKineticEnergy *= 0.5 * 0.5;
01841     reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
01842     intVirialNormal *= 0.5;
01843     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01844     if ( simParams->multigratorOn) {
01845       momentumSqrSum *= 0.5;
01846       ADD_TENSOR_OBJECT(reduction,REDUCTION_MOMENTUM_SQUARED,momentumSqrSum);
01847     }
01848   }
01849 
01850 }
01851 
01852 void Sequencer::calcFixVirial(Tensor& fixVirialNormal, Tensor& fixVirialNbond, Tensor& fixVirialSlow,
01853   Vector& fixForceNormal, Vector& fixForceNbond, Vector& fixForceSlow) {
01854 
01855   FullAtom *a = patch->atom.begin();
01856   int numAtoms = patch->numAtoms;
01857 
01858   for ( int j = 0; j < numAtoms; j++ ) {
01859     if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01860       Vector dx = a[j].fixedPosition;
01861       // all negative because fixed atoms cancels these forces
01862       fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
01863       fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
01864       fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
01865       fixForceNormal -= patch->f[Results::normal][j];
01866       fixForceNbond -= patch->f[Results::nbond][j];
01867       fixForceSlow -= patch->f[Results::slow][j];
01868     }
01869   }
01870 }
01871 
01872 void Sequencer::submitReductions(int step)
01873 {
01874   FullAtom *a = patch->atom.begin();
01875   int numAtoms = patch->numAtoms;
01876 
01877 #if CMK_BLUEGENEL
01878   CmiNetworkProgressAfter(0);
01879 #endif
01880 
01881   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01882   reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
01883 
01884   // For non-Multigrator doKineticEnergy = 1 always
01885   if (doKineticEnergy || doMomenta || patch->flags.doVirial)
01886   {
01887     BigReal kineticEnergy = 0;
01888     Vector momentum = 0;
01889     Vector angularMomentum = 0;
01890     Vector o = patch->lattice.origin();
01891     int i;
01892     if ( simParams->pairInteractionOn ) {
01893       if ( simParams->pairInteractionSelf ) {
01894         for (i = 0; i < numAtoms; ++i ) {
01895           if ( a[i].partition != 1 ) continue;
01896           kineticEnergy += a[i].mass * a[i].velocity.length2();
01897           momentum += a[i].mass * a[i].velocity;
01898           angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01899         }
01900       }
01901     } else {
01902       for (i = 0; i < numAtoms; ++i ) {
01903         kineticEnergy += a[i].mass * a[i].velocity.length2();
01904         momentum += a[i].mass * a[i].velocity;
01905         angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01906       }
01907       if (simParams->drudeOn) {
01908         BigReal drudeComKE = 0.;
01909         BigReal drudeBondKE = 0.;
01910 
01911         for (i = 0;  i < numAtoms;  i++) {
01912           if (i < numAtoms-1 &&
01913               a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01914             // i+1 is a Drude particle with parent i
01915 
01916             // convert from Cartesian coordinates to (COM,bond) coordinates
01917             BigReal m_com = (a[i].mass + a[i+1].mass);  // mass of COM
01918             BigReal m = a[i+1].mass / m_com;  // mass ratio
01919             BigReal m_bond = a[i+1].mass * (1. - m);  // mass of bond
01920             Vector v_bond = a[i+1].velocity - a[i].velocity;  // vel of bond
01921             Vector v_com = a[i].velocity + m * v_bond;  // vel of COM
01922 
01923             drudeComKE += m_com * v_com.length2();
01924             drudeBondKE += m_bond * v_bond.length2();
01925 
01926             i++;  // +1 from loop, we've updated both particles
01927           }
01928           else {
01929             drudeComKE += a[i].mass * a[i].velocity.length2();
01930           }
01931         } // end for
01932 
01933         drudeComKE *= 0.5;
01934         drudeBondKE *= 0.5;
01935         reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY)
01936           += drudeComKE;
01937         reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY)
01938           += drudeBondKE;
01939       } // end drudeOn
01940 
01941     } // end else
01942 
01943     kineticEnergy *= 0.5;
01944     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
01945     ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
01946     ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);  
01947   }
01948 
01949 #ifdef ALTVIRIAL
01950   // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
01951   {
01952     Tensor altVirial;
01953     for ( int i = 0; i < numAtoms; ++i ) {
01954       altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
01955     }
01956     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
01957   }
01958   {
01959     Tensor altVirial;
01960     for ( int i = 0; i < numAtoms; ++i ) {
01961       altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
01962     }
01963     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
01964   }
01965   {
01966     Tensor altVirial;
01967     for ( int i = 0; i < numAtoms; ++i ) {
01968       altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
01969     }
01970     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
01971   }
01972 #endif
01973 
01974   // For non-Multigrator doKineticEnergy = 1 always
01975   if (doKineticEnergy || patch->flags.doVirial)
01976   {
01977     BigReal intKineticEnergy = 0;
01978     Tensor intVirialNormal;
01979     Tensor intVirialNbond;
01980     Tensor intVirialSlow;
01981 
01982     int hgs;
01983     for ( int i = 0; i < numAtoms; i += hgs ) {
01984 #if CMK_BLUEGENEL
01985       CmiNetworkProgress();
01986 #endif
01987       hgs = a[i].hydrogenGroupSize;
01988       int j;
01989       BigReal m_cm = 0;
01990       Position x_cm(0,0,0);
01991       Velocity v_cm(0,0,0);
01992       for ( j = i; j < (i+hgs); ++j ) {
01993         m_cm += a[j].mass;
01994         x_cm += a[j].mass * a[j].position;
01995         v_cm += a[j].mass * a[j].velocity;
01996       }
01997       x_cm /= m_cm;
01998       v_cm /= m_cm;
01999       int fixedAtomsOn = simParams->fixedAtomsOn;
02000       if ( simParams->pairInteractionOn ) {
02001         int pairInteractionSelf = simParams->pairInteractionSelf;
02002         for ( j = i; j < (i+hgs); ++j ) {
02003           if ( a[j].partition != 1 &&
02004                ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
02005           // net force treated as zero for fixed atoms
02006           if ( fixedAtomsOn && a[j].atomFixed ) continue;
02007           BigReal mass = a[j].mass;
02008           Vector v = a[j].velocity;
02009           Vector dv = v - v_cm;
02010           intKineticEnergy += mass * (v * dv);
02011           Vector dx = a[j].position - x_cm;
02012           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02013           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02014           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02015         }
02016       } else {
02017         for ( j = i; j < (i+hgs); ++j ) {
02018           // net force treated as zero for fixed atoms
02019           if ( fixedAtomsOn && a[j].atomFixed ) continue;
02020           BigReal mass = a[j].mass;
02021           Vector v = a[j].velocity;
02022           Vector dv = v - v_cm;
02023           intKineticEnergy += mass * (v * dv);
02024           Vector dx = a[j].position - x_cm;
02025           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02026           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02027           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02028         }
02029       }
02030     }
02031 
02032     intKineticEnergy *= 0.5;
02033     reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
02034     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02035     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
02036     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
02037   }
02038 
02039   if (pressureProfileReduction && simParams->useGroupPressure) {
02040     // subtract off internal virial term, calculated as for intVirial.
02041     int nslabs = simParams->pressureProfileSlabs;
02042     const Lattice &lattice = patch->lattice;
02043     BigReal idz = nslabs/lattice.c().z;
02044     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
02045     int useGroupPressure = simParams->useGroupPressure;
02046 
02047     int hgs;
02048     for (int i=0; i<numAtoms; i += hgs) {
02049       int j;
02050       hgs = a[i].hydrogenGroupSize;
02051       BigReal m_cm = 0;
02052       Position x_cm(0,0,0);
02053       for (j=i; j< i+hgs; ++j) {
02054         m_cm += a[j].mass;
02055         x_cm += a[j].mass * a[j].position;
02056       }
02057       x_cm /= m_cm;
02058       
02059       BigReal z = a[i].position.z;
02060       int slab = (int)floor((z-zmin)*idz);
02061       if (slab < 0) slab += nslabs;
02062       else if (slab >= nslabs) slab -= nslabs;
02063       int partition = a[i].partition;
02064       int ppoffset = 3*(slab + nslabs*partition);
02065       for (j=i; j < i+hgs; ++j) {
02066         BigReal mass = a[j].mass;
02067         Vector dx = a[j].position - x_cm;
02068         const Vector &fnormal = patch->f[Results::normal][j];
02069         const Vector &fnbond  = patch->f[Results::nbond][j];
02070         const Vector &fslow   = patch->f[Results::slow][j];
02071         BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
02072         BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
02073         BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
02074         pressureProfileReduction->item(ppoffset  ) -= wxx;
02075         pressureProfileReduction->item(ppoffset+1) -= wyy;
02076         pressureProfileReduction->item(ppoffset+2) -= wzz;
02077       }
02078     }
02079   }
02080 
02081   // For non-Multigrator doVirial = 1 always
02082   if (patch->flags.doVirial)
02083   {
02084     if ( simParams->fixedAtomsOn ) {
02085       Tensor fixVirialNormal;
02086       Tensor fixVirialNbond;
02087       Tensor fixVirialSlow;
02088       Vector fixForceNormal = 0;
02089       Vector fixForceNbond = 0;
02090       Vector fixForceSlow = 0;
02091 
02092       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
02093 
02094       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
02095       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
02096       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
02097       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
02098       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
02099       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
02100     }
02101   }
02102 
02103   reduction->submit();
02104   if (pressureProfileReduction) pressureProfileReduction->submit();
02105 }
02106 
02107 void Sequencer::submitMinimizeReductions(int step, BigReal fmax2)
02108 {
02109   FullAtom *a = patch->atom.begin();
02110   Force *f1 = patch->f[Results::normal].begin();
02111   Force *f2 = patch->f[Results::nbond].begin();
02112   Force *f3 = patch->f[Results::slow].begin();
02113   const bool fixedAtomsOn = simParams->fixedAtomsOn;
02114   const bool drudeHardWallOn = simParams->drudeHardWallOn;
02115   const double drudeBondLen = simParams->drudeBondLen;
02116   const double drudeBondLen2 = drudeBondLen * drudeBondLen;
02117   const double drudeStep = 0.1/(TIMEFACTOR*TIMEFACTOR);
02118   const double drudeMove = 0.01;
02119   const double drudeStep2 = drudeStep * drudeStep;
02120   const double drudeMove2 = drudeMove * drudeMove;
02121   int numAtoms = patch->numAtoms;
02122 
02123   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
02124 
02125   for ( int i = 0; i < numAtoms; ++i ) {
02126     f1[i] += f2[i] + f3[i];  // add all forces
02127     if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
02128       if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
02129         if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
02130           a[i].position += drudeMove * f1[i].unit();
02131         } else {
02132           a[i].position += drudeStep * f1[i];
02133         }
02134         if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
02135           a[i].position = a[i-1].position + drudeBondLen * (a[i].position - a[i-1].position).unit();
02136         }
02137       }
02138       Vector netf = f1[i-1] + f1[i];
02139       if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
02140       f1[i-1] = netf;
02141       f1[i] = 0.;
02142     }
02143     if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
02144   }
02145 
02146   f2 = f3 = 0;  // included in f1
02147 
02148   BigReal maxv2 = 0.;
02149 
02150   for ( int i = 0; i < numAtoms; ++i ) {
02151     BigReal v2 = a[i].velocity.length2();
02152     if ( v2 > 0. ) {
02153       if ( v2 > maxv2 ) maxv2 = v2;
02154     } else {
02155       v2 = f1[i].length2();
02156       if ( v2 > maxv2 ) maxv2 = v2;
02157     }
02158   }
02159 
02160   if ( fmax2 > 10. * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR )
02161   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial, true /* forces */); }
02162 
02163   BigReal fdotf = 0;
02164   BigReal fdotv = 0;
02165   BigReal vdotv = 0;
02166   int numHuge = 0;
02167   for ( int i = 0; i < numAtoms; ++i ) {
02168     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
02169     if ( drudeHardWallOn && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) continue; // drude particle
02170     Force f = f1[i];
02171     BigReal ff = f * f;
02172     if ( ff > fmax2 ) {
02173       if (simParams->printBadContacts) {
02174         CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
02175       }
02176       ++numHuge;
02177       // pad scaling so minimizeMoveDownhill() doesn't miss them
02178       BigReal fmult = 1.01 * sqrt(fmax2/ff);
02179       f *= fmult;  ff = f * f;
02180       f1[i] *= fmult;
02181     }
02182     fdotf += ff;
02183     fdotv += f * a[i].velocity;
02184     vdotv += a[i].velocity * a[i].velocity;
02185   }
02186 
02187   reduction->item(REDUCTION_MIN_F_DOT_F) += fdotf;
02188   reduction->item(REDUCTION_MIN_F_DOT_V) += fdotv;
02189   reduction->item(REDUCTION_MIN_V_DOT_V) += vdotv;
02190   reduction->item(REDUCTION_MIN_HUGE_COUNT) += numHuge;
02191 
02192   {
02193     Tensor intVirialNormal;
02194     Tensor intVirialNbond;
02195     Tensor intVirialSlow;
02196 
02197     int hgs;
02198     for ( int i = 0; i < numAtoms; i += hgs ) {
02199       hgs = a[i].hydrogenGroupSize;
02200       int j;
02201       BigReal m_cm = 0;
02202       Position x_cm(0,0,0);
02203       for ( j = i; j < (i+hgs); ++j ) {
02204         m_cm += a[j].mass;
02205         x_cm += a[j].mass * a[j].position;
02206       }
02207       x_cm /= m_cm;
02208       for ( j = i; j < (i+hgs); ++j ) {
02209         BigReal mass = a[j].mass;
02210         // net force treated as zero for fixed atoms
02211         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
02212         Vector dx = a[j].position - x_cm;
02213         intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02214         intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02215         intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02216       }
02217     }
02218 
02219     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02220     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
02221     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
02222   }
02223 
02224   if ( simParams->fixedAtomsOn ) {
02225     Tensor fixVirialNormal;
02226     Tensor fixVirialNbond;
02227     Tensor fixVirialSlow;
02228     Vector fixForceNormal = 0;
02229     Vector fixForceNbond = 0;
02230     Vector fixForceSlow = 0;
02231 
02232     calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
02233 
02234     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
02235     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
02236     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
02237     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
02238     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
02239     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
02240   }
02241 
02242   reduction->submit();
02243 }
02244 
02245 void Sequencer::submitCollections(int step, int zeroVel)
02246 {
02247   int prec = Output::coordinateNeeded(step);
02248   if ( prec )
02249     collection->submitPositions(step,patch->atom,patch->lattice,prec);
02250   if ( Output::velocityNeeded(step) )
02251     collection->submitVelocities(step,zeroVel,patch->atom);
02252   if ( Output::forceNeeded(step) ) {
02253     int maxForceUsed = patch->flags.maxForceUsed;
02254     if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
02255     collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
02256   }
02257 }
02258 
02259 void Sequencer::runComputeObjects(int migration, int pairlists, int pressureStep)
02260 {
02261   if ( migration ) pairlistsAreValid = 0;
02262 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
02263   if ( pairlistsAreValid &&
02264        ( patch->flags.doFullElectrostatics || ! simParams->fullElectFrequency )
02265                          && ( pairlistsAge > (
02266 #else
02267   if ( pairlistsAreValid && ( pairlistsAge > (
02268 #endif
02269          (simParams->stepsPerCycle - 1) / simParams->pairlistsPerCycle ) ) ) {
02270     pairlistsAreValid = 0;
02271   }
02272   if ( ! simParams->usePairlists ) pairlists = 0;
02273   patch->flags.usePairlists = pairlists || pairlistsAreValid;
02274   patch->flags.savePairlists =
02275         pairlists && ! pairlistsAreValid;
02276 
02277   if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
02278 
02279   patch->positionsReady(migration);  // updates flags.sequence
02280 
02281   int seq = patch->flags.sequence;
02282   int basePriority = ( (seq & 0xffff) << 15 )
02283                      + PATCH_PRIORITY(patch->getPatchID());
02284   if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
02285     priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
02286     suspend(); // until all deposit boxes close
02287     patch->gbisComputeAfterP1();
02288     priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
02289     suspend();
02290     patch->gbisComputeAfterP2();
02291     priority = basePriority + COMPUTE_HOME_PRIORITY;
02292     suspend();
02293   } else {
02294     priority = basePriority + COMPUTE_HOME_PRIORITY;
02295     suspend(); // until all deposit boxes close
02296   }
02297 
02298   if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
02299     pairlistsAreValid = 1;
02300     pairlistsAge = 0;
02301   }
02302   // For multigrator, do not age pairlist during pressure step
02303   // NOTE: for non-multigrator pressureStep = 0 always
02304   if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
02305 
02306   if (simParams->lonepairs) {
02307     {
02308       Tensor virial;
02309       patch->redistrib_lonepair_forces(Results::normal, &virial);
02310       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02311     }
02312     if (patch->flags.doNonbonded) {
02313       Tensor virial;
02314       patch->redistrib_lonepair_forces(Results::nbond, &virial);
02315       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02316     }
02317     if (patch->flags.doFullElectrostatics) {
02318       Tensor virial;
02319       patch->redistrib_lonepair_forces(Results::slow, &virial);
02320       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02321     }
02322   } else if (simParams->watmodel == WAT_TIP4) {
02323     {
02324       Tensor virial;
02325       patch->redistrib_tip4p_forces(Results::normal, &virial);
02326       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02327     }
02328     if (patch->flags.doNonbonded) {
02329       Tensor virial;
02330       patch->redistrib_tip4p_forces(Results::nbond, &virial);
02331       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02332     }
02333     if (patch->flags.doFullElectrostatics) {
02334       Tensor virial;
02335       patch->redistrib_tip4p_forces(Results::slow, &virial);
02336       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02337     }
02338   } else if (simParams->watmodel == WAT_SWM4) {
02339     {
02340       Tensor virial;
02341       patch->redistrib_swm4_forces(Results::normal, &virial);
02342       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02343     }
02344     if (patch->flags.doNonbonded) {
02345       Tensor virial;
02346       patch->redistrib_swm4_forces(Results::nbond, &virial);
02347       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02348     }
02349     if (patch->flags.doFullElectrostatics) {
02350       Tensor virial;
02351       patch->redistrib_swm4_forces(Results::slow, &virial);
02352       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02353     }
02354   }
02355 
02356   if ( patch->flags.doMolly ) {
02357     Tensor virial;
02358     patch->mollyMollify(&virial);
02359     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
02360   }
02361 
02362 
02363   // BEGIN LA
02364   if (patch->flags.doLoweAndersen) {
02365       patch->loweAndersenFinish();
02366   }
02367   // END LA
02368 #ifdef NAMD_CUDA_XXX
02369   int numAtoms = patch->numAtoms;
02370   FullAtom *a = patch->atom.begin();
02371   for ( int i=0; i<numAtoms; ++i ) {
02372     CkPrintf("%d %g %g %g\n", a[i].id,
02373         patch->f[Results::normal][i].x +
02374         patch->f[Results::nbond][i].x +
02375         patch->f[Results::slow][i].x,
02376         patch->f[Results::normal][i].y + 
02377         patch->f[Results::nbond][i].y +
02378         patch->f[Results::slow][i].y,
02379         patch->f[Results::normal][i].z +
02380         patch->f[Results::nbond][i].z +
02381         patch->f[Results::slow][i].z);
02382     CkPrintf("%d %g %g %g\n", a[i].id,
02383         patch->f[Results::normal][i].x,
02384         patch->f[Results::nbond][i].x,
02385         patch->f[Results::slow][i].x);
02386     CkPrintf("%d %g %g %g\n", a[i].id,
02387         patch->f[Results::normal][i].y,
02388         patch->f[Results::nbond][i].y,
02389         patch->f[Results::slow][i].y);
02390     CkPrintf("%d %g %g %g\n", a[i].id,
02391         patch->f[Results::normal][i].z,
02392         patch->f[Results::nbond][i].z,
02393         patch->f[Results::slow][i].z);
02394   }
02395 #endif
02396 
02397 #if PRINT_FORCES
02398   int numAtoms = patch->numAtoms;
02399   FullAtom *a = patch->atom.begin();
02400   for ( int i=0; i<numAtoms; ++i ) {
02401     float fxNo = patch->f[Results::normal][i].x;
02402     float fxNb = patch->f[Results::nbond][i].x;
02403     float fxSl = patch->f[Results::slow][i].x;
02404     float fyNo = patch->f[Results::normal][i].y;
02405     float fyNb = patch->f[Results::nbond][i].y;
02406     float fySl = patch->f[Results::slow][i].y;
02407     float fzNo = patch->f[Results::normal][i].z;
02408     float fzNb = patch->f[Results::nbond][i].z;
02409     float fzSl = patch->f[Results::slow][i].z;
02410     float fx = fxNo+fxNb+fxSl;
02411     float fy = fyNo+fyNb+fySl;
02412     float fz = fzNo+fzNb+fzSl;
02413 
02414                 float f = sqrt(fx*fx+fy*fy+fz*fz);
02415     int id = patch->pExt[i].id;
02416     int seq = patch->flags.sequence;
02417     float x = patch->p[i].position.x;
02418     float y = patch->p[i].position.y;
02419     float z = patch->p[i].position.z;
02420     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
02421     CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
02422     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
02423 //fxNo,fyNo,fzNo,
02424 //fxNb,fyNb,fzNb,
02425 //fxSl,fySl,fzSl,
02426 fx,fy,fz
02427 );
02428         }
02429 #endif
02430 }
02431 
02432 void Sequencer::rebalanceLoad(int timestep) {
02433   if ( ! ldbSteps ) {
02434     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02435   }
02436   if ( ! --ldbSteps ) {
02437     patch->submitLoadStats(timestep);
02438     ldbCoordinator->rebalance(this,patch->getPatchID());
02439     pairlistsAreValid = 0;
02440   }
02441 }
02442 
02443 void Sequencer::cycleBarrier(int doBarrier, int step) {
02444 #if USE_BARRIER
02445         if (doBarrier)
02446           broadcast->cycleBarrier.get(step);
02447 #endif
02448 }
02449 
02450 void Sequencer::traceBarrier(int step){
02451         broadcast->traceBarrier.get(step);
02452 }
02453 
02454 #ifdef MEASURE_NAMD_WITH_PAPI
02455 void Sequencer::papiMeasureBarrier(int step){
02456         broadcast->papiMeasureBarrier.get(step);
02457 }
02458 #endif
02459 
02460 void Sequencer::terminate() {
02461   LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
02462   CthFree(thread);
02463   CthSuspend();
02464 }
02465 

Generated on Thu Sep 21 01:17:14 2017 for NAMD by  doxygen 1.4.7