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

Generated on Mon Nov 20 01:17:14 2017 for NAMD by  doxygen 1.4.7