Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

Sequencer Class Reference

#include <Sequencer.h>

List of all members.

Public Member Functions

 Sequencer (HomePatch *p)
virtual ~Sequencer (void)
void run (void)
void awaken (void)
void suspend (void)

Protected Member Functions

virtual void algorithm (void)
void integrate ()
void minimize ()
void runComputeObjects (int migration=1, int pairlists=0)
void submitReductions (int)
void submitHalfstep (int)
void submitMinimizeReductions (int)
void submitCollections (int step, int zeroVel=0)
void submitMomentum (int step)
void correctMomentum (int step, BigReal drifttime)
void saveForce (const int ftag=Results::normal)
void addForceToMomentum (BigReal, const int ftag=Results::normal, const int useSaved=0, const int pressure=0)
void addVelocityToPosition (BigReal)
void addRotDragToPosition (BigReal)
void addMovDragToPosition (BigReal)
void newMinimizeDirection (BigReal)
void newMinimizePosition (BigReal)
void quenchVelocities ()
void rattle1 (BigReal, int)
void rattle2 (BigReal, int)
void maximumMove (BigReal)
void minimizationQuenchVelocity (void)
void reloadCharges ()
void rescaleVelocities (int)
void reassignVelocities (BigReal, int)
void reinitVelocities (void)
void rescaleVelocitiesByFactor (BigReal)
void tcoupleVelocities (BigReal, int)
void berendsenPressure (int)
void langevinPiston (int)
void langevinVelocities (BigReal)
void langevinVelocitiesBBK1 (BigReal)
void langevinVelocitiesBBK2 (BigReal)
void cycleBarrier (int, int)
void terminate (void)
void rebalanceLoad (int timestep)

Protected Attributes

int pairlistsAreValid
int pairlistsAge
int rescaleVelocities_numTemps
int berendsenPressure_count
int checkpoint_berendsenPressure_count
int slowFreq
Random * random
SimParameters *const simParams
HomePatch *const patch
SubmitReduction * reduction
SubmitReduction * pressureProfileReduction
CollectionMgr *const collection
ControllerBroadcasts * broadcast
int ldbSteps

Friends

class HomePatch

Constructor & Destructor Documentation

Sequencer::Sequencer (  HomePatch *  p  ) 
 

Definition at line 38 of file Sequencer.C.

References berendsenPressure_count, broadcast, Patch::getPatchID(), PatchMap::numPatches(), PatchMap::Object(), LdbCoordinator::Object(), ReductionMgr::Object(), patch, SimParameters::pressureProfileOn, pressureProfileReduction, random, SimParameters::randomSeed, reduction, REDUCTIONS_BASIC, REDUCTIONS_PPROF_INTERNAL, rescaleVelocities_numTemps, simParams, simParams, Random::split(), and ReductionMgr::willSubmit().

00038                                  :
00039         simParams(Node::Object()->simParameters),
00040         patch(p),
00041         collection(CollectionMgr::Object()),
00042         ldbSteps(0)
00043 {
00044     broadcast = new ControllerBroadcasts;
00045     reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00046     if (simParams->pressureProfileOn) {
00047       pressureProfileReduction = 
00048         ReductionMgr::Object()->willSubmit(REDUCTIONS_PPROF_INTERNAL);
00049     } else {
00050       pressureProfileReduction = NULL;
00051     }
00052     ldbCoordinator = (LdbCoordinator::Object());
00053     random = new Random(simParams->randomSeed);
00054     random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
00055 
00056     rescaleVelocities_numTemps = 0;
00057     berendsenPressure_count = 0;
00058 }

Sequencer::~Sequencer (  void   )  [virtual]
 

Definition at line 60 of file Sequencer.C.

00061 {
00062     delete broadcast;
00063     delete reduction;
00064     delete pressureProfileReduction;
00065     delete random;
00066 }

Member Function Documentation

void Sequencer::addForceToMomentum (  BigReal  ,
const int  ftag = Results::normal,
const int  useSaved = 0,
const int  pressure = 0
)  [protected]
 

Definition at line 844 of file Sequencer.C.

References ADD_TENSOR_OBJECT, HomePatch::addForceToMomentum(), patch, reduction, REDUCTION_VIRIAL_NORMAL, simParams, and SimParameters::watmodel.

Referenced by integrate().

00846 {
00847 #if CMK_VERSION_BLUEGENE
00848   CmiNetworkProgressAfter (0);
00849 #endif
00850   Tensor virial;
00851   Tensor *vp = ( pressure ? &virial : 0 );
00852   patch->addForceToMomentum(dt,ftag,useSaved,vp);
00853   if (simParams->watmodel == WAT_TIP4) ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
00854 }

void Sequencer::addMovDragToPosition (  BigReal   )  [protected]
 

Definition at line 346 of file Sequencer.C.

References HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), BigReal, SimParameters::fixedAtomsOn, Molecule::get_movdrag_params(), Molecule::is_atom_movdragged(), Node::molecule, SimParameters::movDragGlobVel, Patch::numAtoms, Node::Object(), patch, CompAtom::position, and simParams.

Referenced by integrate().

00346                                                      {
00347   FullAtom *atom = patch->atom.begin();
00348   int numAtoms = patch->numAtoms;
00349   Molecule *molecule = Node::Object()->molecule;   // need its methods
00350   const BigReal movDragGlobVel = simParams->movDragGlobVel;
00351   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00352   Vector movDragVel, dragIncrement;
00353   for ( int i = 0; i < numAtoms; ++i )
00354   {
00355     // skip if fixed atom or zero drag attribute
00356     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00357          || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00358     molecule->get_movdrag_params(movDragVel, atom[i].id);
00359     dragIncrement = movDragGlobVel * movDragVel * dt;
00360     atom[i].position += dragIncrement;
00361   }
00362 }

void Sequencer::addRotDragToPosition (  BigReal   )  [protected]
 

Definition at line 365 of file Sequencer.C.

References HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), BigReal, SimParameters::fixedAtomsOn, Molecule::get_rotdrag_params(), Molecule::is_atom_rotdragged(), Vector::length(), Node::molecule, Patch::numAtoms, Node::Object(), patch, CompAtom::position, SimParameters::rotDragGlobVel, and simParams.

Referenced by integrate().

00365                                                      {
00366   FullAtom *atom = patch->atom.begin();
00367   int numAtoms = patch->numAtoms;
00368   Molecule *molecule = Node::Object()->molecule;   // need its methods
00369   const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00370   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00371   BigReal rotDragVel, dAngle;
00372   Vector atomRadius;
00373   Vector rotDragAxis, rotDragPivot, dragIncrement;
00374   for ( int i = 0; i < numAtoms; ++i )
00375   {
00376     // skip if fixed atom or zero drag attribute
00377     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00378          || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00379     molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00380     dAngle = rotDragGlobVel * rotDragVel * dt;
00381     rotDragAxis /= rotDragAxis.length();
00382     atomRadius = atom[i].position - rotDragPivot;
00383     dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00384     atom[i].position += dragIncrement;
00385   }
00386 }

void Sequencer::addVelocityToPosition (  BigReal   )  [protected]
 

Definition at line 856 of file Sequencer.C.

References HomePatch::addVelocityToPosition(), and patch.

Referenced by integrate().

00857 {
00858 #if CMK_VERSION_BLUEGENE
00859   CmiNetworkProgressAfter (0);
00860 #endif
00861   patch->addVelocityToPosition(dt);
00862 }

void Sequencer::algorithm (  void   )  [protected, virtual]
 

Definition at line 87 of file Sequencer.C.

References berendsenPressure_count, broadcast, HomePatch::checkpoint(), checkpoint_berendsenPressure_count, END_OF_RUN, EVAL_MEASURE, FILE_OUTPUT, SimpleBroadcastObject< T >::get(), integrate(), minimize(), pairlistsAreValid, patch, reinitVelocities(), reloadCharges(), rescaleVelocitiesByFactor(), HomePatch::revert(), SCRIPT_CHECKPOINT, SCRIPT_MEASURE, SCRIPT_MINIMIZE, SCRIPT_OUTPUT, SCRIPT_REINITVELS, SCRIPT_RELOADCHARGES, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, simParams, submitCollections(), and terminate().

00088 {
00089   int scriptTask;
00090   int scriptSeq = 0;
00091   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
00092     switch ( scriptTask ) {
00093       case SCRIPT_OUTPUT:
00094         submitCollections(FILE_OUTPUT);
00095         break;
00096       case SCRIPT_MEASURE:
00097         submitCollections(EVAL_MEASURE);
00098         break;
00099       case SCRIPT_REINITVELS:
00100         reinitVelocities();
00101         break;
00102       case SCRIPT_RESCALEVELS:
00103         rescaleVelocitiesByFactor(simParams->scriptArg1);
00104         break;
00105       case SCRIPT_RELOADCHARGES:
00106         reloadCharges();
00107         break;
00108       case SCRIPT_CHECKPOINT:
00109         patch->checkpoint();
00110         checkpoint_berendsenPressure_count = berendsenPressure_count;
00111         break;
00112       case SCRIPT_REVERT:
00113         patch->revert();
00114         berendsenPressure_count = checkpoint_berendsenPressure_count;
00115         pairlistsAreValid = 0;
00116         break;
00117       case SCRIPT_MINIMIZE:
00118         minimize();
00119         break;
00120       case SCRIPT_RUN:
00121         integrate();
00122         break;
00123     }
00124   }
00125   submitCollections(END_OF_RUN);
00126   terminate();
00127 }

void Sequencer::awaken (  void   )  [inline]
 

Definition at line 28 of file Sequencer.h.

Referenced by LdbCoordinator::awakenSequencers(), HomePatch::boxClosed(), HomePatch::depositMigration(), and run().

00028 { CthAwaken(thread); }

void Sequencer::berendsenPressure (  int   )  [protected]
 

Definition at line 596 of file Sequencer.C.

References Lattice::apply_transform(), HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), berendsenPressure_count, SimParameters::berendsenPressureFreq, SimParameters::berendsenPressureOn, BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), CompAtom::groupFixed, CompAtom::hydrogenGroupSize, Patch::lattice, FullAtom::mass, Patch::numAtoms, patch, Position, CompAtom::position, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, and SimParameters::useGroupPressure.

Referenced by integrate().

00597 {
00598   if ( simParams->berendsenPressureOn ) {
00599   berendsenPressure_count += 1;
00600   const int freq = simParams->berendsenPressureFreq;
00601   if ( ! (berendsenPressure_count % freq ) ) {
00602    berendsenPressure_count = 0;
00603    FullAtom *a = patch->atom.begin();
00604    int numAtoms = patch->numAtoms;
00605    Tensor factor = broadcast->positionRescaleFactor.get(step);
00606    patch->lattice.rescale(factor);
00607    if ( simParams->useGroupPressure )
00608    {
00609     int hgs;
00610     for ( int i = 0; i < numAtoms; i += hgs ) {
00611       int j;
00612       hgs = a[i].hydrogenGroupSize;
00613       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00614         for ( j = i; j < (i+hgs); ++j ) {
00615           a[j].position = patch->lattice.apply_transform(
00616                                 a[j].fixedPosition,a[j].transform);
00617         }
00618         continue;
00619       }
00620       BigReal m_cm = 0;
00621       Position x_cm(0,0,0);
00622       for ( j = i; j < (i+hgs); ++j ) {
00623         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00624         m_cm += a[j].mass;
00625         x_cm += a[j].mass * a[j].position;
00626       }
00627       x_cm /= m_cm;
00628       Position new_x_cm = x_cm;
00629       patch->lattice.rescale(new_x_cm,factor);
00630       Position delta_x_cm = new_x_cm - x_cm;
00631       for ( j = i; j < (i+hgs); ++j ) {
00632         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00633           a[j].position = patch->lattice.apply_transform(
00634                                 a[j].fixedPosition,a[j].transform);
00635           continue;
00636         }
00637         a[j].position += delta_x_cm;
00638       }
00639     }
00640    }
00641    else
00642    {
00643     for ( int i = 0; i < numAtoms; ++i )
00644     {
00645       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00646         a[i].position = patch->lattice.apply_transform(
00647                                 a[i].fixedPosition,a[i].transform);
00648         continue;
00649       }
00650       patch->lattice.rescale(a[i].position,factor);
00651     }
00652    }
00653   }
00654   } else {
00655     berendsenPressure_count = 0;
00656   }
00657 }

void Sequencer::correctMomentum (  int  step,
BigReal  drifttime
)  [protected]
 

Definition at line 495 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), FullAtom::mass, ControllerBroadcasts::momentumCorrection, NAMD_die(), Patch::numAtoms, patch, CompAtom::position, simParams, FullAtom::velocity, and SimParameters::zeroMomentumAlt.

Referenced by integrate().

00495                                                            {
00496 
00497   if ( simParams->fixedAtomsOn )
00498     NAMD_die("Cannot zero momentum when fixed atoms are present.");
00499 
00500   const Vector dv = broadcast->momentumCorrection.get(step);
00501   const Vector dx = dv * ( drifttime / TIMEFACTOR );
00502 
00503   FullAtom *a = patch->atom.begin();
00504   const int numAtoms = patch->numAtoms;
00505 
00506 if ( simParams->zeroMomentumAlt ) {
00507   for ( int i = 0; i < numAtoms; ++i ) {
00508     BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00509     a[i].velocity += dv * rmass;
00510     a[i].position += dx * rmass;
00511   }
00512 } else {
00513   for ( int i = 0; i < numAtoms; ++i ) {
00514     a[i].velocity += dv;
00515     a[i].position += dx;
00516   }
00517 }
00518 
00519 }

void Sequencer::cycleBarrier (  int  ,
int 
)  [protected]
 

Definition at line 1433 of file Sequencer.C.

References broadcast.

Referenced by integrate().

01433                                                     {
01434 #if USE_BARRIER
01435         if (doBarrier)
01436           broadcast->cycleBarrier.get(step);
01437 #endif
01438 }

void Sequencer::integrate (   )  [protected]
 

Definition at line 130 of file Sequencer.C.

References addForceToMomentum(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), berendsenPressure(), BigReal, Bool, SimParameters::commOnly, ComputeMgr::computeGlobalObject, Node::computeMgr, correctMomentum(), cycleBarrier(), Flags::doEnergy, Flags::doFullElectrostatics, Flags::doMolly, Flags::doNonbonded, SimParameters::dt, SimParameters::firstTimestep, Patch::flags, SimParameters::FMAOn, SimParameters::fullDirectOn, SimParameters::fullElectFrequency, langevinPiston(), langevinVelocitiesBBK1(), langevinVelocitiesBBK2(), Flags::maxForceMerged, Flags::maxForceUsed, maximumMove(), minimizationQuenchVelocity(), SimParameters::mollyOn, SimParameters::movDragOn, SimParameters::MTSAlgorithm, SimParameters::N, SimParameters::nonbondedFrequency, Node::Object(), SimParameters::outputEnergies, patch, SimParameters::PMEOn, rattle1(), SimParameters::reassignFreq, reassignVelocities(), rebalanceLoad(), rescaleVelocities(), SimParameters::rotDragOn, runComputeObjects(), saveForce(), ComputeGlobal::saveTotalForces(), simParams, slowFreq, Flags::step, SimParameters::stepsPerCycle, submitCollections(), submitHalfstep(), submitMomentum(), submitReductions(), SimParameters::tclForcesOn, tcoupleVelocities(), and SimParameters::zeroMomentum.

Referenced by algorithm().

00130                           {
00131 
00132     int &step = patch->flags.step;
00133     step = simParams->firstTimestep;
00134 
00135     // drag switches
00136     const Bool rotDragOn = simParams->rotDragOn;
00137     const Bool movDragOn = simParams->movDragOn;
00138 
00139     const int commOnly = simParams->commOnly;
00140 
00141     int &maxForceUsed = patch->flags.maxForceUsed;
00142     int &maxForceMerged = patch->flags.maxForceMerged;
00143     maxForceUsed = Results::normal;
00144     maxForceMerged = Results::normal;
00145 
00146     const int numberOfSteps = simParams->N;
00147     const int stepsPerCycle = simParams->stepsPerCycle;
00148     const BigReal timestep = simParams->dt;
00149 
00150     // what MTS method?
00151     const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00152 
00153     const int nonbondedFrequency = simParams->nonbondedFrequency;
00154     slowFreq = nonbondedFrequency;
00155     const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00156     int &doNonbonded = patch->flags.doNonbonded;
00157     doNonbonded = !(step%nonbondedFrequency);
00158     if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00159     if ( doNonbonded ) maxForceUsed = Results::nbond;
00160 
00161     // Do we do full electrostatics?
00162     const int dofull = ( simParams->fullDirectOn ||
00163                         simParams->FMAOn || simParams->PMEOn );
00164     const int fullElectFrequency = simParams->fullElectFrequency;
00165     if ( dofull ) slowFreq = fullElectFrequency;
00166     const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00167     int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00168     doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00169     if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00170                                         maxForceMerged = Results::slow;
00171     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00172     int &doMolly = patch->flags.doMolly;
00173     doMolly = simParams->mollyOn && doFullElectrostatics;
00174 
00175     int zeroMomentum = simParams->zeroMomentum;
00176     
00177     // Do we need to return forces to TCL script?
00178     int doTcl = simParams->tclForcesOn;
00179     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00180 
00181     // Bother to calculate energies?
00182     int &doEnergy = patch->flags.doEnergy;
00183     int energyFrequency = simParams->outputEnergies;
00184 
00185     rattle1(0.,0);  // enforce rigid bond constraints on initial positions
00186 
00187     const int reassignFreq = simParams->reassignFreq;
00188     if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00189        reassignVelocities(timestep,step);
00190     }
00191 
00192     doEnergy = ! ( step % energyFrequency );
00193     runComputeObjects(1,step<numberOfSteps); // must migrate here!
00194     if ( staleForces || doTcl ) {
00195       if ( doNonbonded ) saveForce(Results::nbond);
00196       if ( doFullElectrostatics ) saveForce(Results::slow);
00197     }
00198     if ( ! commOnly ) {
00199       addForceToMomentum(-0.5*timestep);
00200       if (staleForces || doNonbonded)
00201                 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00202       if (staleForces || doFullElectrostatics)
00203                 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces);
00204     }
00205     minimizationQuenchVelocity();
00206     rattle1(-timestep,0);
00207     submitHalfstep(step);
00208     if ( ! commOnly ) {
00209       addForceToMomentum(timestep);
00210       if (staleForces || doNonbonded)
00211                 addForceToMomentum(nbondstep,Results::nbond,staleForces,0);
00212       if (staleForces || doFullElectrostatics)
00213                 addForceToMomentum(slowstep,Results::slow,staleForces,0);
00214     }
00215     rattle1(timestep,1);
00216     if (doTcl)  // include constraint forces
00217       computeGlobal->saveTotalForces(patch);
00218     submitHalfstep(step);
00219     if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00220     if ( ! commOnly ) {
00221       addForceToMomentum(-0.5*timestep);
00222       if (staleForces || doNonbonded)
00223                 addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00224       if (staleForces || doFullElectrostatics)
00225                 addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00226     }
00227     submitReductions(step);
00228     rebalanceLoad(step);
00229 
00230     for ( ++step; step <= numberOfSteps; ++step )
00231     {
00232 
00233       rescaleVelocities(step);
00234       tcoupleVelocities(timestep,step);
00235       berendsenPressure(step);
00236 
00237       if ( ! commOnly ) {
00238         addForceToMomentum(0.5*timestep);
00239         if (staleForces || doNonbonded)
00240           addForceToMomentum(0.5*nbondstep,Results::nbond,staleForces,1);
00241         if (staleForces || doFullElectrostatics)
00242           addForceToMomentum(0.5*slowstep,Results::slow,staleForces,1);
00243       }
00244 
00245       /* reassignment based on half-step velocities
00246          if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00247          addVelocityToPosition(0.5*timestep);
00248          reassignVelocities(timestep,step);
00249          addVelocityToPosition(0.5*timestep);
00250          rattle1(0.,0);
00251          rattle1(-timestep,0);
00252          addVelocityToPosition(-1.0*timestep);
00253          rattle1(timestep,0);
00254          } */
00255 
00256       maximumMove(timestep);
00257       if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00258       langevinPiston(step);
00259       if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00260 
00261       minimizationQuenchVelocity();
00262 
00263       doNonbonded = !(step%nonbondedFrequency);
00264       doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00265 
00266       if ( zeroMomentum && doFullElectrostatics )
00267         correctMomentum(step,slowstep);
00268 
00269       submitHalfstep(step);
00270 
00271       doMolly = simParams->mollyOn && doFullElectrostatics;
00272 
00273       maxForceUsed = Results::normal;
00274       if ( doNonbonded ) maxForceUsed = Results::nbond;
00275       if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00276 
00277       // Migrate Atoms on stepsPerCycle
00278       doEnergy = ! ( step % energyFrequency );
00279       runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00280       if ( staleForces || doTcl ) {
00281         if ( doNonbonded ) saveForce(Results::nbond);
00282         if ( doFullElectrostatics ) saveForce(Results::slow);
00283       }
00284 
00285       // reassignment based on full-step velocities
00286       if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00287         reassignVelocities(timestep,step);
00288         addForceToMomentum(-0.5*timestep);
00289         if (staleForces || doNonbonded)
00290           addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,0);
00291         if (staleForces || doFullElectrostatics)
00292           addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,0);
00293         rattle1(-timestep,0);
00294       }
00295 
00296       if ( ! commOnly ) {
00297         langevinVelocitiesBBK1(timestep);
00298         addForceToMomentum(timestep);
00299         if (staleForces || doNonbonded)
00300           addForceToMomentum(nbondstep,Results::nbond,staleForces,1);
00301         if (staleForces || doFullElectrostatics)
00302           addForceToMomentum(slowstep,Results::slow,staleForces,1);
00303         langevinVelocitiesBBK2(timestep);
00304       }
00305 
00306       // add drag to each atom's positions
00307       if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00308       if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00309 
00310       rattle1(timestep,1);
00311       if (doTcl)  // include constraint forces
00312         computeGlobal->saveTotalForces(patch);
00313 
00314       submitHalfstep(step);
00315       if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00316 
00317       if ( ! commOnly ) {
00318         addForceToMomentum(-0.5*timestep);
00319         if (staleForces || doNonbonded)
00320           addForceToMomentum(-0.5*nbondstep,Results::nbond,staleForces,1);
00321         if (staleForces || doFullElectrostatics)
00322           addForceToMomentum(-0.5*slowstep,Results::slow,staleForces,1);
00323       }
00324 
00325         // rattle2(timestep,step);
00326 
00327         submitReductions(step);
00328         submitCollections(step);
00329 
00330 #if CYCLE_BARRIER
00331         cycleBarrier(!((step+1) % stepsPerCycle), step);
00332 #elif PME_BARRIER
00333         cycleBarrier(doFullElectrostatics, step);
00334 #endif
00335 
00336         rebalanceLoad(step);
00337 
00338 #if PME_BARRIER
00339         // a step before PME
00340         cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00341 #endif
00342     }
00343 }

void Sequencer::langevinPiston (  int   )  [protected]
 

Definition at line 659 of file Sequencer.C.

References Lattice::apply_transform(), HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), BigReal, broadcast, SimParameters::fixedAtomsOn, SimpleBroadcastObject< T >::get(), CompAtom::groupFixed, CompAtom::hydrogenGroupSize, Molecule::is_atom_exPressure(), SimParameters::langevinPistonOn, Patch::lattice, FullAtom::mass, Node::molecule, Patch::numAtoms, Node::Object(), patch, Position, CompAtom::position, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, slowFreq, SimParameters::useGroupPressure, FullAtom::velocity, Velocity, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by integrate().

00660 {
00661   if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
00662   {
00663    FullAtom *a = patch->atom.begin();
00664    int numAtoms = patch->numAtoms;
00665    Tensor factor = broadcast->positionRescaleFactor.get(step);
00666    // JCP FIX THIS!!!
00667    Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
00668    patch->lattice.rescale(factor);
00669    Molecule *mol = Node::Object()->molecule;
00670    if ( simParams->useGroupPressure )
00671    {
00672     int hgs;
00673     for ( int i = 0; i < numAtoms; i += hgs ) {
00674       int j;
00675       hgs = a[i].hydrogenGroupSize;
00676       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
00677         for ( j = i; j < (i+hgs); ++j ) {
00678           a[j].position = patch->lattice.apply_transform(
00679                                 a[j].fixedPosition,a[j].transform);
00680         }
00681         continue;
00682       }
00683       BigReal m_cm = 0;
00684       Position x_cm(0,0,0);
00685       Velocity v_cm(0,0,0);
00686       for ( j = i; j < (i+hgs); ++j ) {
00687         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
00688         m_cm += a[j].mass;
00689         x_cm += a[j].mass * a[j].position;
00690         v_cm += a[j].mass * a[j].velocity;
00691       }
00692       x_cm /= m_cm;
00693       Position new_x_cm = x_cm;
00694       patch->lattice.rescale(new_x_cm,factor);
00695       Position delta_x_cm = new_x_cm - x_cm;
00696       v_cm /= m_cm;
00697       Velocity delta_v_cm;
00698       delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
00699       delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
00700       delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
00701       for ( j = i; j < (i+hgs); ++j ) {
00702         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
00703           a[j].position = patch->lattice.apply_transform(
00704                                 a[j].fixedPosition,a[j].transform);
00705           continue;
00706         }
00707         if ( mol->is_atom_exPressure(a[j].id) ) continue;
00708         a[j].position += delta_x_cm;
00709         a[j].velocity += delta_v_cm;
00710       }
00711     }
00712    }
00713    else
00714    {
00715     for ( int i = 0; i < numAtoms; ++i )
00716     {
00717       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
00718         a[i].position = patch->lattice.apply_transform(
00719                                 a[i].fixedPosition,a[i].transform);
00720         continue;
00721       }
00722       if ( mol->is_atom_exPressure(a[i].id) ) continue;
00723       patch->lattice.rescale(a[i].position,factor);
00724       a[i].velocity.x *= velFactor.x;
00725       a[i].velocity.y *= velFactor.y;
00726       a[i].velocity.z *= velFactor.z;
00727     }
00728    }
00729   }
00730 }

void Sequencer::langevinVelocities (  BigReal   )  [protected]
 

Definition at line 521 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, Random::gaussian_vector(), CompAtom::id, Molecule::langevin_param(), SimParameters::langevinOn, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, Node::molecule, Patch::numAtoms, Node::Object(), patch, random, simParams, and FullAtom::velocity.

00522 {
00523   if ( simParams->langevinOn )
00524   {
00525     FullAtom *a = patch->atom.begin();
00526     int numAtoms = patch->numAtoms;
00527     Molecule *molecule = Node::Object()->molecule;
00528     BigReal dt = dt_fs * 0.001;  // convert to ps
00529     BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00530 
00531     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00532     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00533 
00534     for ( int i = 0; i < numAtoms; ++i )
00535     {
00536       int aid = a[i].id;
00537       BigReal f1 = exp( -1. * dt * molecule->langevin_param(aid) );
00538       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT * 
00539                          ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00540 
00541       a[i].velocity *= f1;
00542       a[i].velocity += f2 * random->gaussian_vector();
00543     }
00544   }
00545 }

void Sequencer::langevinVelocitiesBBK1 (  BigReal   )  [protected]
 

Definition at line 547 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, CompAtom::id, Molecule::langevin_param(), SimParameters::langevinOn, Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, and FullAtom::velocity.

Referenced by integrate().

00548 {
00549   if ( simParams->langevinOn )
00550   {
00551     FullAtom *a = patch->atom.begin();
00552     int numAtoms = patch->numAtoms;
00553     Molecule *molecule = Node::Object()->molecule;
00554     BigReal dt = dt_fs * 0.001;  // convert to ps
00555     for ( int i = 0; i < numAtoms; ++i )
00556     {
00557       int aid = a[i].id;
00558       BigReal dt_gamma = dt * molecule->langevin_param(aid);
00559       if ( ! dt_gamma ) continue;
00560 
00561       a[i].velocity *= ( 1. - 0.5 * dt_gamma );
00562     }
00563   }
00564 }

void Sequencer::langevinVelocitiesBBK2 (  BigReal   )  [protected]
 

Definition at line 567 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, BOLTZMAN, Random::gaussian_vector(), CompAtom::id, Molecule::langevin_param(), SimParameters::langevinOn, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, Node::molecule, Patch::numAtoms, Node::Object(), patch, random, rattle1(), simParams, and FullAtom::velocity.

Referenced by integrate().

00568 {
00569   if ( simParams->langevinOn )
00570   {
00571     rattle1(dt_fs,1);  // conserve momentum if gammas differ
00572 
00573     FullAtom *a = patch->atom.begin();
00574     int numAtoms = patch->numAtoms;
00575     Molecule *molecule = Node::Object()->molecule;
00576     BigReal dt = dt_fs * 0.001;  // convert to ps
00577     BigReal kbT = BOLTZMAN*(simParams->langevinTemp);
00578 
00579     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
00580     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
00581 
00582     for ( int i = 0; i < numAtoms; ++i )
00583     {
00584       int aid = a[i].id;
00585       BigReal dt_gamma = dt * molecule->langevin_param(aid);
00586       if ( ! dt_gamma ) continue;
00587 
00588       a[i].velocity += random->gaussian_vector() *
00589              sqrt( 2 * dt_gamma * kbT *
00590                    ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
00591       a[i].velocity /= ( 1. + 0.5 * dt_gamma );
00592     }
00593   }
00594 }

void Sequencer::maximumMove (  BigReal   )  [protected]
 

Definition at line 893 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), BigReal, SimParameters::cutoff, Node::enableEarlyExit(), CompAtom::id, iERROR(), iout, Vector::length(), Vector::length2(), SimParameters::maximumMove, Patch::numAtoms, Node::Object(), patch, PDBVELFACTOR, simParams, terminate(), and FullAtom::velocity.

Referenced by integrate().

00894 {
00895   FullAtom *a = patch->atom.begin();
00896   int numAtoms = patch->numAtoms;
00897   if ( simParams->maximumMove ) {
00898     const BigReal dt = timestep / TIMEFACTOR;
00899     const BigReal maxvel = simParams->maximumMove / dt;
00900     const BigReal maxvel2 = maxvel * maxvel;
00901     for ( int i=0; i<numAtoms; ++i ) {
00902       if ( a[i].velocity.length2() > maxvel2 ) {
00903         a[i].velocity *= ( maxvel / a[i].velocity.length() );
00904       }
00905     }
00906   } else {
00907     const BigReal dt = timestep / TIMEFACTOR;
00908     const BigReal maxvel = simParams->cutoff / dt;
00909     const BigReal maxvel2 = maxvel * maxvel;
00910     int killme = 0;
00911     for ( int i=0; i<numAtoms; ++i ) {
00912       killme = killme || ( a[i].velocity.length2() > maxvel2 );
00913     }
00914     if ( killme ) {
00915       for ( int i=0; i<numAtoms; ++i ) {
00916         if ( a[i].velocity.length2() > maxvel2 ) {
00917           iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
00918             << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
00919             << ( PDBVELFACTOR * maxvel ) << ")\n" << endi;
00920         }
00921       }
00922       iout << iERROR << 
00923         "Atoms moving too fast; simulation has become unstable.\n" << endi;
00924       Node::Object()->enableEarlyExit();
00925       terminate();
00926     }
00927   }
00928 }

void Sequencer::minimizationQuenchVelocity (  void   )  [protected]
 

Definition at line 930 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::minimizeOn, Patch::numAtoms, patch, simParams, and FullAtom::velocity.

Referenced by integrate().

00931 {
00932   if ( simParams->minimizeOn ) {
00933     FullAtom *a = patch->atom.begin();
00934     int numAtoms = patch->numAtoms;
00935     for ( int i=0; i<numAtoms; ++i ) {
00936       a[i].velocity = 0.;
00937     }
00938   }
00939 }

void Sequencer::minimize (   )  [protected]
 

Definition at line 388 of file Sequencer.C.

References BigReal, broadcast, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doMolly, Flags::doNonbonded, SimParameters::firstTimestep, Patch::flags, SimParameters::FMAOn, SimParameters::fullDirectOn, SimpleBroadcastObject< T >::get(), Flags::maxForceMerged, Flags::maxForceUsed, ControllerBroadcasts::minimizeCoefficient, SimParameters::mollyOn, SimParameters::N, newMinimizeDirection(), newMinimizePosition(), patch, SimParameters::PMEOn, quenchVelocities(), rebalanceLoad(), runComputeObjects(), simParams, Flags::step, submitCollections(), and submitMinimizeReductions().

Referenced by algorithm().

00388                          {
00389   const int numberOfSteps = simParams->N;
00390   int &step = patch->flags.step;
00391   step = simParams->firstTimestep;
00392 
00393   int &maxForceUsed = patch->flags.maxForceUsed;
00394   int &maxForceMerged = patch->flags.maxForceMerged;
00395   maxForceUsed = Results::normal;
00396   maxForceMerged = Results::normal;
00397   int &doNonbonded = patch->flags.doNonbonded;
00398   doNonbonded = 1;
00399   maxForceUsed = Results::nbond;
00400   maxForceMerged = Results::nbond;
00401   const int dofull = ( simParams->fullDirectOn ||
00402                         simParams->FMAOn || simParams->PMEOn );
00403   int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00404   doFullElectrostatics = dofull;
00405   if ( dofull ) {
00406     maxForceMerged = Results::slow;
00407     maxForceUsed = Results::slow;
00408   }
00409   int &doMolly = patch->flags.doMolly;
00410   doMolly = simParams->mollyOn && doFullElectrostatics;
00411   int &doEnergy = patch->flags.doEnergy;
00412   doEnergy = 1;
00413 
00414   runComputeObjects(1,0); // must migrate here!
00415 
00416   submitMinimizeReductions(step);
00417   rebalanceLoad(step);
00418 
00419   int minSeq = 0;
00420   for ( ++step; step <= numberOfSteps; ++step ) {
00421     BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00422     if ( ! c ) {  // new direction
00423       c = broadcast->minimizeCoefficient.get(minSeq++);
00424       newMinimizeDirection(c);  // v = c * v + f
00425       c = broadcast->minimizeCoefficient.get(minSeq++);
00426     }  // same direction
00427     newMinimizePosition(c);  // x = x + c * v
00428 
00429     runComputeObjects(1,0);
00430     submitMinimizeReductions(step);
00431     submitCollections(step, 1);  // write out zeros for velocities
00432     rebalanceLoad(step);
00433   }
00434   quenchVelocities();  // zero out bogus velocity
00435 }

void Sequencer::newMinimizeDirection (  BigReal   )  [protected]
 

Definition at line 438 of file Sequencer.C.

References HomePatch::atom, CompAtom::atomFixed, ResizeArray< Elem >::begin(), Patch::f, SimParameters::fixedAtomsOn, Force, Patch::numAtoms,