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

Generated on Mon Jul 16 01:17:15 2018 for NAMD by  doxygen 1.4.7