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 (int)
void minimize ()
void runComputeObjects (int migration=1, int pairlists=0, int pressureStep=0)
void calcFixVirial (Tensor &fixVirialNormal, Tensor &fixVirialNbond, Tensor &fixVirialSlow, Vector &fixForceNormal, Vector &fixForceNbond, Vector &fixForceSlow)
void submitReductions (int)
void submitHalfstep (int)
void submitMinimizeReductions (int, BigReal fmax2)
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)
void addForceToMomentum3 (const BigReal timestep1, const int ftag1, const int useSaved1, const BigReal timestep2, const int ftag2, const int useSaved2, const BigReal timestep3, const int ftag3, const int useSaved3)
void addVelocityToPosition (BigReal)
void addRotDragToPosition (BigReal)
void addMovDragToPosition (BigReal)
void minimizeMoveDownhill (BigReal fmax2)
void newMinimizeDirection (BigReal)
void newMinimizePosition (BigReal)
void quenchVelocities ()
void hardWallDrude (BigReal, int)
void rattle1 (BigReal, int)
void maximumMove (BigReal)
void minimizationQuenchVelocity (void)
void reloadCharges ()
void adaptTempUpdate (int)
void rescaleVelocities (int)
void rescaleaccelMD (int, int, int)
void reassignVelocities (BigReal, int)
void reinitVelocities (void)
void rescaleVelocitiesByFactor (BigReal)
void tcoupleVelocities (BigReal, int)
void berendsenPressure (int)
void langevinPiston (int)
void newtonianVelocities (BigReal, const BigReal, const BigReal, const BigReal, const int, const int, const int)
void langevinVelocities (BigReal)
void langevinVelocitiesBBK1 (BigReal)
void langevinVelocitiesBBK2 (BigReal)
void scalePositionsVelocities (const Tensor &posScale, const Tensor &velScale)
void multigratorPressure (int step, int callNumber)
void scaleVelocities (const BigReal velScale)
BigReal calcKineticEnergy ()
void multigratorTemperature (int step, int callNumber)
void cycleBarrier (int, int)
void traceBarrier (int)
void terminate (void)
void rebalanceLoad (int timestep)

Protected Attributes

SubmitReductionmin_reduction
int pairlistsAreValid
int pairlistsAge
BigReal adaptTempT
int rescaleVelocities_numTemps
int berendsenPressure_count
int checkpoint_berendsenPressure_count
int slowFreq
SubmitReductionmultigratorReduction
int doKineticEnergy
int doMomenta
Randomrandom
SimParameters *const simParams
HomePatch *const patch
SubmitReductionreduction
SubmitReductionpressureProfileReduction
CollectionMgr *const collection
ControllerBroadcastsbroadcast
int ldbSteps

Friends

class HomePatch


Detailed Description

Definition at line 22 of file Sequencer.h.


Constructor & Destructor Documentation

Sequencer::Sequencer ( HomePatch p  ) 

Definition at line 46 of file Sequencer.C.

References SimParameters::accelMDOn, berendsenPressure_count, broadcast, Patch::getPatchID(), HomePatch::ldObjHandle, min_reduction, MULTIGRATOR_REDUCTION_MAX_RESERVED, SimParameters::multigratorOn, multigratorReduction, PatchMap::numPatches(), PatchMap::Object(), LdbCoordinator::Object(), ReductionMgr::Object(), patch, SimParameters::pressureProfileAtomTypes, SimParameters::pressureProfileOn, pressureProfileReduction, SimParameters::pressureProfileSlabs, random, SimParameters::randomSeed, reduction, REDUCTIONS_AMD, REDUCTIONS_BASIC, REDUCTIONS_MINIMIZER, REDUCTIONS_MULTIGRATOR, REDUCTIONS_PPROF_INTERNAL, rescaleVelocities_numTemps, simParams, Random::split(), and ReductionMgr::willSubmit().

00046                                  :
00047         simParams(Node::Object()->simParameters),
00048         patch(p),
00049         collection(CollectionMgr::Object()),
00050         ldbSteps(0)
00051 {
00052     broadcast = new ControllerBroadcasts(& patch->ldObjHandle);
00053     reduction = ReductionMgr::Object()->willSubmit(
00054                   simParams->accelMDOn ? REDUCTIONS_AMD : REDUCTIONS_BASIC );
00055     min_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MINIMIZER,1);
00056     if (simParams->pressureProfileOn) {
00057       int ntypes = simParams->pressureProfileAtomTypes;
00058       int nslabs = simParams->pressureProfileSlabs;
00059       pressureProfileReduction = 
00060         ReductionMgr::Object()->willSubmit(
00061                 REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
00062     } else {
00063       pressureProfileReduction = NULL;
00064     }
00065     if (simParams->multigratorOn) {
00066       multigratorReduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
00067     } else {
00068       multigratorReduction = NULL;
00069     }
00070     ldbCoordinator = (LdbCoordinator::Object());
00071     random = new Random(simParams->randomSeed);
00072     random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
00073 
00074     rescaleVelocities_numTemps = 0;
00075     berendsenPressure_count = 0;
00076 //    patch->write_tip4_props();
00077 }

Sequencer::~Sequencer ( void   )  [virtual]

Definition at line 79 of file Sequencer.C.

References broadcast, min_reduction, multigratorReduction, pressureProfileReduction, random, and reduction.

00080 {
00081     delete broadcast;
00082     delete reduction;
00083     delete min_reduction;
00084     if (pressureProfileReduction) delete pressureProfileReduction;
00085     delete random;
00086     if (multigratorReduction) delete multigratorReduction;
00087 }


Member Function Documentation

void Sequencer::adaptTempUpdate ( int   )  [protected]

Definition at line 1460 of file Sequencer.C.

References ControllerBroadcasts::adaptTemperature, SimParameters::adaptTempFreq, SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, adaptTempT, broadcast, SimParameters::firstTimestep, SimParameters::langevinOn, SimParameters::langevinTemp, and simParams.

Referenced by integrate().

01461 {
01462    //check if adaptive tempering is enabled and in the right timestep range
01463    if (!simParams->adaptTempOn) return;
01464    if ( (step < simParams->adaptTempFirstStep ) || 
01465      ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) {
01466         if (simParams->langevinOn) // restore langevin temperature
01467             adaptTempT = simParams->langevinTemp;
01468         return;
01469    }
01470    // Get Updated Temperature
01471    if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
01472     adaptTempT = broadcast->adaptTemperature.get(step);
01473 }

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

Definition at line 1572 of file Sequencer.C.

References HomePatch::addForceToMomentum(), and patch.

Referenced by newtonianVelocities().

01573 {
01574 #if CMK_BLUEGENEL
01575   CmiNetworkProgressAfter (0);
01576 #endif
01577   patch->addForceToMomentum(dt,ftag,useSaved);
01578 }

void Sequencer::addForceToMomentum3 ( const BigReal  timestep1,
const int  ftag1,
const int  useSaved1,
const BigReal  timestep2,
const int  ftag2,
const int  useSaved2,
const BigReal  timestep3,
const int  ftag3,
const int  useSaved3 
) [protected]

Definition at line 1580 of file Sequencer.C.

References HomePatch::addForceToMomentum3(), and patch.

Referenced by newtonianVelocities().

01582                                                                        {
01583 #if CMK_BLUEGENEL
01584   CmiNetworkProgressAfter (0);
01585 #endif
01586   patch->addForceToMomentum3(timestep1,ftag1,useSaved1, timestep2,ftag2,useSaved2, timestep3,ftag3,useSaved3);  
01587 }

void Sequencer::addMovDragToPosition ( BigReal   )  [protected]

Definition at line 501 of file Sequencer.C.

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

Referenced by integrate().

00501                                                      {
00502   FullAtom *atom = patch->atom.begin();
00503   int numAtoms = patch->numAtoms;
00504   Molecule *molecule = Node::Object()->molecule;   // need its methods
00505   const BigReal movDragGlobVel = simParams->movDragGlobVel;
00506   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00507   Vector movDragVel, dragIncrement;
00508   for ( int i = 0; i < numAtoms; ++i )
00509   {
00510     // skip if fixed atom or zero drag attribute
00511     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00512          || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00513     molecule->get_movdrag_params(movDragVel, atom[i].id);
00514     dragIncrement = movDragGlobVel * movDragVel * dt;
00515     atom[i].position += dragIncrement;
00516   }
00517 }

void Sequencer::addRotDragToPosition ( BigReal   )  [protected]

Definition at line 520 of file Sequencer.C.

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

Referenced by integrate().

00520                                                      {
00521   FullAtom *atom = patch->atom.begin();
00522   int numAtoms = patch->numAtoms;
00523   Molecule *molecule = Node::Object()->molecule;   // need its methods
00524   const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00525   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00526   BigReal rotDragVel, dAngle;
00527   Vector atomRadius;
00528   Vector rotDragAxis, rotDragPivot, dragIncrement;
00529   for ( int i = 0; i < numAtoms; ++i )
00530   {
00531     // skip if fixed atom or zero drag attribute
00532     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00533          || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00534     molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00535     dAngle = rotDragGlobVel * rotDragVel * dt;
00536     rotDragAxis /= rotDragAxis.length();
00537     atomRadius = atom[i].position - rotDragPivot;
00538     dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00539     atom[i].position += dragIncrement;
00540   }
00541 }

void Sequencer::addVelocityToPosition ( BigReal   )  [protected]

Definition at line 1589 of file Sequencer.C.

References HomePatch::addVelocityToPosition(), and patch.

Referenced by integrate().

01590 {
01591 #if CMK_BLUEGENEL
01592   CmiNetworkProgressAfter (0);
01593 #endif
01594   patch->addVelocityToPosition(dt);
01595 }

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

Definition at line 117 of file Sequencer.C.

References berendsenPressure_count, broadcast, HomePatch::checkpoint(), checkpoint_berendsenPressure_count, END_OF_RUN, EVAL_MEASURE, HomePatch::exchangeAtoms(), HomePatch::exchangeCheckpoint(), FILE_OUTPUT, FORCE_OUTPUT, SimpleBroadcastObject< T >::get(), integrate(), minimize(), NAMD_bug(), pairlistsAreValid, patch, reinitVelocities(), reloadCharges(), rescaleVelocitiesByFactor(), HomePatch::revert(), SCRIPT_ATOMRECV, SCRIPT_ATOMSEND, SCRIPT_ATOMSENDRECV, SCRIPT_CHECKPOINT, SCRIPT_CHECKPOINT_FREE, SCRIPT_CHECKPOINT_LOAD, SCRIPT_CHECKPOINT_STORE, SCRIPT_CHECKPOINT_SWAP, SCRIPT_CONTINUE, SCRIPT_END, SCRIPT_FORCEOUTPUT, SCRIPT_MEASURE, SCRIPT_MINIMIZE, SCRIPT_OUTPUT, SCRIPT_REINITVELS, SCRIPT_RELOADCHARGES, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, simParams, submitCollections(), and terminate().

00118 {
00119   int scriptTask;
00120   int scriptSeq = 0;
00121   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
00122     switch ( scriptTask ) {
00123       case SCRIPT_OUTPUT:
00124         submitCollections(FILE_OUTPUT);
00125         break;
00126       case SCRIPT_FORCEOUTPUT:
00127         submitCollections(FORCE_OUTPUT);
00128         break;
00129       case SCRIPT_MEASURE:
00130         submitCollections(EVAL_MEASURE);
00131         break;
00132       case SCRIPT_REINITVELS:
00133         reinitVelocities();
00134         break;
00135       case SCRIPT_RESCALEVELS:
00136         rescaleVelocitiesByFactor(simParams->scriptArg1);
00137         break;
00138       case SCRIPT_RELOADCHARGES:
00139         reloadCharges();
00140         break;
00141       case SCRIPT_CHECKPOINT:
00142         patch->checkpoint();
00143         checkpoint_berendsenPressure_count = berendsenPressure_count;
00144         break;
00145       case SCRIPT_REVERT:
00146         patch->revert();
00147         berendsenPressure_count = checkpoint_berendsenPressure_count;
00148         pairlistsAreValid = 0;
00149         break;
00150       case SCRIPT_CHECKPOINT_STORE:
00151       case SCRIPT_CHECKPOINT_LOAD:
00152       case SCRIPT_CHECKPOINT_SWAP:
00153       case SCRIPT_CHECKPOINT_FREE:
00154         patch->exchangeCheckpoint(scriptTask,berendsenPressure_count);
00155         break;
00156       case SCRIPT_ATOMSENDRECV:
00157       case SCRIPT_ATOMSEND:
00158       case SCRIPT_ATOMRECV:
00159         patch->exchangeAtoms(scriptTask);
00160         break;
00161       case SCRIPT_MINIMIZE:
00162         minimize();
00163         break;
00164       case SCRIPT_RUN:
00165       case SCRIPT_CONTINUE:
00166         integrate(scriptTask);
00167         break;
00168       default:
00169         NAMD_bug("Unknown task in Sequencer::algorithm");
00170     }
00171   }
00172   submitCollections(END_OF_RUN);
00173   terminate();
00174 }

void Sequencer::awaken ( void   )  [inline]

Definition at line 29 of file Sequencer.h.

References PRIORITY_SIZE.

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

00029                       {
00030       CthAwakenPrio(thread, CK_QUEUEING_IFIFO, PRIORITY_SIZE, &priority);
00031     }

void Sequencer::berendsenPressure ( int   )  [protected]

Definition at line 1264 of file Sequencer.C.

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

Referenced by integrate().

01265 {
01266   if ( simParams->berendsenPressureOn ) {
01267   berendsenPressure_count += 1;
01268   const int freq = simParams->berendsenPressureFreq;
01269   if ( ! (berendsenPressure_count % freq ) ) {
01270    berendsenPressure_count = 0;
01271    FullAtom *a = patch->atom.begin();
01272    int numAtoms = patch->numAtoms;
01273    Tensor factor = broadcast->positionRescaleFactor.get(step);
01274    patch->lattice.rescale(factor);
01275    if ( simParams->useGroupPressure )
01276    {
01277     int hgs;
01278     for ( int i = 0; i < numAtoms; i += hgs ) {
01279       int j;
01280       hgs = a[i].hydrogenGroupSize;
01281       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
01282         for ( j = i; j < (i+hgs); ++j ) {
01283           a[j].position = patch->lattice.apply_transform(
01284                                 a[j].fixedPosition,a[j].transform);
01285         }
01286         continue;
01287       }
01288       BigReal m_cm = 0;
01289       Position x_cm(0,0,0);
01290       for ( j = i; j < (i+hgs); ++j ) {
01291         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01292         m_cm += a[j].mass;
01293         x_cm += a[j].mass * a[j].position;
01294       }
01295       x_cm /= m_cm;
01296       Position new_x_cm = x_cm;
01297       patch->lattice.rescale(new_x_cm,factor);
01298       Position delta_x_cm = new_x_cm - x_cm;
01299       for ( j = i; j < (i+hgs); ++j ) {
01300         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01301           a[j].position = patch->lattice.apply_transform(
01302                                 a[j].fixedPosition,a[j].transform);
01303           continue;
01304         }
01305         a[j].position += delta_x_cm;
01306       }
01307     }
01308    }
01309    else
01310    {
01311     for ( int i = 0; i < numAtoms; ++i )
01312     {
01313       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01314         a[i].position = patch->lattice.apply_transform(
01315                                 a[i].fixedPosition,a[i].transform);
01316         continue;
01317       }
01318       patch->lattice.rescale(a[i].position,factor);
01319     }
01320    }
01321   }
01322   } else {
01323     berendsenPressure_count = 0;
01324   }
01325 }

void Sequencer::calcFixVirial ( Tensor fixVirialNormal,
Tensor fixVirialNbond,
Tensor fixVirialSlow,
Vector fixForceNormal,
Vector fixForceNbond,
Vector fixForceSlow 
) [protected]

Definition at line 1852 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::f, SimParameters::fixedAtomsOn, FullAtom::fixedPosition, j, Results::nbond, Results::normal, Patch::numAtoms, Tensor::outerAdd(), patch, simParams, and Results::slow.

Referenced by multigratorPressure().

01853                                                                        {
01854 
01855   FullAtom *a = patch->atom.begin();
01856   int numAtoms = patch->numAtoms;
01857 
01858   for ( int j = 0; j < numAtoms; j++ ) {
01859     if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01860       Vector dx = a[j].fixedPosition;
01861       // all negative because fixed atoms cancels these forces
01862       fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
01863       fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
01864       fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
01865       fixForceNormal -= patch->f[Results::normal][j];
01866       fixForceNbond -= patch->f[Results::nbond][j];
01867       fixForceSlow -= patch->f[Results::slow][j];
01868     }
01869   }
01870 }

BigReal Sequencer::calcKineticEnergy (  )  [protected]

Definition at line 986 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), Vector::length2(), FullAtom::mass, Patch::numAtoms, SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, partition(), patch, simParams, and FullAtom::velocity.

Referenced by multigratorTemperature().

00986                                      {
00987   FullAtom *a = patch->atom.begin();
00988   int numAtoms = patch->numAtoms;
00989   BigReal kineticEnergy = 0.0;
00990   if ( simParams->pairInteractionOn ) {
00991     if ( simParams->pairInteractionSelf ) {
00992       for (int i = 0; i < numAtoms; ++i ) {
00993         if ( a[i].partition != 1 ) continue;
00994         kineticEnergy += a[i].mass * a[i].velocity.length2();
00995       }
00996     }
00997   } else {
00998     for (int i = 0; i < numAtoms; ++i ) {
00999       kineticEnergy += a[i].mass * a[i].velocity.length2();
01000     }
01001   }
01002   kineticEnergy *= 0.5;
01003   return kineticEnergy;
01004 }

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

Definition at line 779 of file Sequencer.C.

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

Referenced by integrate().

00779                                                            {
00780 
00781   if ( simParams->fixedAtomsOn )
00782     NAMD_die("Cannot zero momentum when fixed atoms are present.");
00783 
00784   const Vector dv = broadcast->momentumCorrection.get(step);
00785   const Vector dx = dv * ( drifttime / TIMEFACTOR );
00786 
00787   FullAtom *a = patch->atom.begin();
00788   const int numAtoms = patch->numAtoms;
00789 
00790 if ( simParams->zeroMomentumAlt ) {
00791   for ( int i = 0; i < numAtoms; ++i ) {
00792     BigReal rmass = (a[i].mass > 0. ? 1. / a[i].mass : 0.);
00793     a[i].velocity += dv * rmass;
00794     a[i].position += dx * rmass;
00795   }
00796 } else {
00797   for ( int i = 0; i < numAtoms; ++i ) {
00798     a[i].velocity += dv;
00799     a[i].position += dx;
00800   }
00801 }
00802 
00803 }

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

Definition at line 2443 of file Sequencer.C.

References broadcast.

Referenced by integrate().

02443                                                     {
02444 #if USE_BARRIER
02445         if (doBarrier)
02446           broadcast->cycleBarrier.get(step);
02447 #endif
02448 }

void Sequencer::hardWallDrude ( BigReal  ,
int   
) [protected]

Definition at line 1597 of file Sequencer.C.

References ADD_TENSOR_OBJECT, SimParameters::drudeHardWallOn, Node::enableEarlyExit(), endi(), HomePatch::hardWallDrude(), iERROR(), iout, Node::Object(), patch, pressureProfileReduction, reduction, REDUCTION_VIRIAL_NORMAL, simParams, and terminate().

Referenced by integrate().

01598 {
01599   if ( simParams->drudeHardWallOn ) {
01600     Tensor virial;
01601     Tensor *vp = ( pressure ? &virial : 0 );
01602     if ( patch->hardWallDrude(dt, vp, pressureProfileReduction) ) {
01603       iout << iERROR << "Constraint failure in HardWallDrude(); "
01604         << "simulation may become unstable.\n" << endi;
01605       Node::Object()->enableEarlyExit();
01606       terminate();
01607     }
01608     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01609   }
01610 }

void Sequencer::integrate ( int   )  [protected]

Definition at line 179 of file Sequencer.C.

References SimParameters::accelMDdihe, SimParameters::accelMDdual, SimParameters::accelMDOn, SimParameters::adaptTempOn, adaptTempT, adaptTempUpdate(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), Results::amdf, HomePatch::atom, Patch::atomMapper, ResizeArray< Elem >::begin(), berendsenPressure(), SimParameters::colvarsOn, SimParameters::commOnly, ComputeMgr::computeGlobalObject, Node::computeMgr, correctMomentum(), cycleBarrier(), Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, doKineticEnergy, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, doMomenta, Flags::doNonbonded, Flags::doVirial, SimParameters::dt, ResizeArray< Elem >::end(), SimParameters::firstTimestep, Patch::flags, SimParameters::fullElectFrequency, SimParameters::GBISOn, hardWallDrude(), SimParameters::initialTemp, SimParameters::langevin_useBAOAB, SimParameters::langevinOn, langevinPiston(), SimParameters::langevinPistonOn, SimParameters::langevinTemp, langevinVelocities(), langevinVelocitiesBBK1(), langevinVelocitiesBBK2(), SimParameters::LCPOOn, SimParameters::lonepairs, SimParameters::loweAndersenOn, Flags::maxForceMerged, Flags::maxForceUsed, maximumMove(), minimizationQuenchVelocity(), SimParameters::mollyOn, SimParameters::movDragOn, SimParameters::MTSAlgorithm, SimParameters::multigratorOn, multigratorPressure(), SimParameters::multigratorPressureFreq, multigratorTemperature(), SimParameters::multigratorTemperatureFreq, SimParameters::N, NAIVE, Results::nbond, newtonianVelocities(), SimParameters::nonbondedFrequency, Results::normal, SimParameters::numTraceSteps, Node::Object(), SimParameters::outputEnergies, SimParameters::outputMomenta, SimParameters::outputPressure, patch, Patch::patchID, rattle1(), SimParameters::reassignFreq, reassignVelocities(), rebalanceLoad(), AtomMapper::registerIDsFullAtom(), rescaleaccelMD(), SimParameters::rescaleFreq, SimParameters::rescaleTemp, rescaleVelocities(), SimParameters::rotDragOn, runComputeObjects(), saveForce(), ComputeGlobal::saveTotalForces(), SCRIPT_RUN, simParams, Results::slow, slowFreq, SimParameters::statsOn, Flags::step, GlobalMaster::step, SimParameters::stepsPerCycle, submitCollections(), submitHalfstep(), submitMomentum(), submitReductions(), SimParameters::tclForcesOn, tcoupleVelocities(), traceBarrier(), SimParameters::traceStartStep, and SimParameters::zeroMomentum.

Referenced by algorithm().

00179                                         {
00180     char traceNote[24];
00181     char tracePrefix[20];
00182     sprintf(tracePrefix, "p:%d,s:",patch->patchID);
00183 //    patch->write_tip4_props();
00184 
00185     int &step = patch->flags.step;
00186     step = simParams->firstTimestep;
00187 
00188     // drag switches
00189     const Bool rotDragOn = simParams->rotDragOn;
00190     const Bool movDragOn = simParams->movDragOn;
00191 
00192     const int commOnly = simParams->commOnly;
00193 
00194     int &maxForceUsed = patch->flags.maxForceUsed;
00195     int &maxForceMerged = patch->flags.maxForceMerged;
00196     maxForceUsed = Results::normal;
00197     maxForceMerged = Results::normal;
00198 
00199     const int numberOfSteps = simParams->N;
00200     const int stepsPerCycle = simParams->stepsPerCycle;
00201     const BigReal timestep = simParams->dt;
00202 
00203     // what MTS method?
00204     const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00205 
00206     const int nonbondedFrequency = simParams->nonbondedFrequency;
00207     slowFreq = nonbondedFrequency;
00208     const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00209     int &doNonbonded = patch->flags.doNonbonded;
00210     doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
00211     if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00212     if ( doNonbonded ) maxForceUsed = Results::nbond;
00213 
00214     // Do we do full electrostatics?
00215     const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00216     const int fullElectFrequency = simParams->fullElectFrequency;
00217     if ( dofull ) slowFreq = fullElectFrequency;
00218     const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00219     int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00220     doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
00221     if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00222                                         maxForceMerged = Results::slow;
00223     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00224 
00225     const Bool accelMDOn = simParams->accelMDOn;
00226     const Bool accelMDdihe = simParams->accelMDdihe;
00227     const Bool accelMDdual = simParams->accelMDdual;
00228     if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
00229 
00230     // Is adaptive tempering on?
00231     const Bool adaptTempOn = simParams->adaptTempOn;
00232     adaptTempT = simParams->initialTemp;
00233     if (simParams->langevinOn)
00234         adaptTempT = simParams->langevinTemp;
00235     else if (simParams->rescaleFreq > 0)
00236         adaptTempT = simParams->rescaleTemp;
00237         
00238 
00239     int &doMolly = patch->flags.doMolly;
00240     doMolly = simParams->mollyOn && doFullElectrostatics;
00241     // BEGIN LA
00242     int &doLoweAndersen = patch->flags.doLoweAndersen;
00243     doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00244     // END LA
00245 
00246     int &doGBIS = patch->flags.doGBIS;
00247     doGBIS = simParams->GBISOn;
00248 
00249     int &doLCPO = patch->flags.doLCPO;
00250     doLCPO = simParams->LCPOOn;
00251 
00252     int zeroMomentum = simParams->zeroMomentum;
00253     
00254     // Do we need to return forces to TCL script or Colvar module?
00255     int doTcl = simParams->tclForcesOn;
00256         int doColvars = simParams->colvarsOn;
00257     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00258 
00259     // Bother to calculate energies?
00260     int &doEnergy = patch->flags.doEnergy;
00261     int energyFrequency = simParams->outputEnergies;
00262 
00263     const int reassignFreq = simParams->reassignFreq;
00264 
00265     int &doVirial = patch->flags.doVirial;
00266     doVirial = 1;
00267 
00268   if ( scriptTask == SCRIPT_RUN ) {
00269 
00270 //    printf("Doing initial rattle\n");
00271     rattle1(0.,0);  // enforce rigid bond constraints on initial positions
00272 
00273     if (simParams->lonepairs) {
00274       patch->atomMapper->registerIDsFullAtom(
00275                 patch->atom.begin(),patch->atom.end());
00276     }
00277 
00278     if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00279        reassignVelocities(timestep,step);
00280     }
00281 
00282     doEnergy = ! ( step % energyFrequency );
00283     if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00284     //Update energy every timestep for adaptive tempering
00285     if ( adaptTempOn ) doEnergy=1;
00286     runComputeObjects(1,step<numberOfSteps); // must migrate here!
00287     rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD 
00288     adaptTempUpdate(step); // update adaptive tempering temperature
00289 
00290     if ( staleForces || doTcl || doColvars ) {
00291       if ( doNonbonded ) saveForce(Results::nbond);
00292       if ( doFullElectrostatics ) saveForce(Results::slow);
00293     }
00294     if ( ! commOnly ) {
00295       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
00296     }
00297     minimizationQuenchVelocity();
00298     rattle1(-timestep,0);
00299     submitHalfstep(step);
00300     if ( ! commOnly ) {
00301       newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
00302     }
00303     rattle1(timestep,1);
00304     if (doTcl || doColvars)  // include constraint forces
00305       computeGlobal->saveTotalForces(patch);
00306     submitHalfstep(step);
00307     if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00308     if ( ! commOnly ) {
00309       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
00310     }
00311     submitReductions(step);
00312     if(traceIsOn()){
00313         traceUserEvent(eventEndOfTimeStep);
00314         sprintf(traceNote, "%s%d",tracePrefix,step); 
00315         traceUserSuppliedNote(traceNote);
00316     }
00317     rebalanceLoad(step);
00318 
00319   } // scriptTask == SCRIPT_RUN
00320 
00321   bool doMultigratorRattle = false;
00322 
00323     for ( ++step; step <= numberOfSteps; ++step )
00324     {
00325 
00326       rescaleVelocities(step);
00327       tcoupleVelocities(timestep,step);
00328       berendsenPressure(step);
00329 
00330       if ( ! commOnly ) {
00331         newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); 
00332       }
00333 
00334       // We do RATTLE here if multigrator thermostat was applied in the previous step
00335       if (doMultigratorRattle) rattle1(timestep, 1);
00336 
00337       /* reassignment based on half-step velocities
00338          if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00339          addVelocityToPosition(0.5*timestep);
00340          reassignVelocities(timestep,step);
00341          addVelocityToPosition(0.5*timestep);
00342          rattle1(0.,0);
00343          rattle1(-timestep,0);
00344          addVelocityToPosition(-1.0*timestep);
00345          rattle1(timestep,0);
00346          } */
00347 
00348       maximumMove(timestep);
00349       if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) {
00350         if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00351         // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
00352         langevinVelocities(timestep);
00353         langevinPiston(step);
00354         if ( ! commOnly ) addVelocityToPosition(0.5*timestep);
00355       } else {
00356         // If Langevin is not used, take full time step directly instread of two half steps      
00357         if ( ! commOnly ) addVelocityToPosition(timestep); 
00358       }
00359 
00360       // impose hard wall potential for Drude bond length
00361       hardWallDrude(timestep, 1);
00362 
00363       minimizationQuenchVelocity();
00364 
00365       doNonbonded = !(step%nonbondedFrequency);
00366       doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00367 
00368       if ( zeroMomentum && doFullElectrostatics )
00369         correctMomentum(step,slowstep);
00370 
00371       submitHalfstep(step);
00372 
00373       doMolly = simParams->mollyOn && doFullElectrostatics;
00374       // BEGIN LA
00375       doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00376       // END LA
00377 
00378       maxForceUsed = Results::normal;
00379       if ( doNonbonded ) maxForceUsed = Results::nbond;
00380       if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00381       if ( accelMDOn && (accelMDdihe || accelMDdual))  maxForceUsed = Results::amdf;
00382 
00383       // Migrate Atoms on stepsPerCycle
00384       doEnergy = ! ( step % energyFrequency );
00385       if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00386       if ( adaptTempOn ) doEnergy=1; 
00387 
00388       // Multigrator 
00389       if (simParams->multigratorOn) {
00390         doVirial = (!(step % energyFrequency) || ((simParams->outputPressure > 0) && !(step % simParams->outputPressure)) 
00391           || !(step % simParams->multigratorPressureFreq));
00392         doKineticEnergy = (!(step % energyFrequency) || !(step % simParams->multigratorTemperatureFreq));
00393         doMomenta = (simParams->outputMomenta > 0) && !(step % simParams->outputMomenta);
00394       } else {
00395         doVirial = 1;
00396         doKineticEnergy = 1;
00397         doMomenta = 1;
00398       }
00399 
00400       runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00401 
00402       rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
00403      
00404       if ( staleForces || doTcl || doColvars ) {
00405         if ( doNonbonded ) saveForce(Results::nbond);
00406         if ( doFullElectrostatics ) saveForce(Results::slow);
00407       }
00408 
00409       // reassignment based on full-step velocities
00410       if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00411         reassignVelocities(timestep,step);
00412         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00413         rattle1(-timestep,0);
00414       }
00415 
00416       if ( ! commOnly ) {
00417         langevinVelocitiesBBK1(timestep);
00418         newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00419         langevinVelocitiesBBK2(timestep);
00420       }
00421 
00422       // add drag to each atom's positions
00423       if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00424       if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00425 
00426       rattle1(timestep,1);
00427       if (doTcl || doColvars)  // include constraint forces
00428         computeGlobal->saveTotalForces(patch);
00429 
00430       submitHalfstep(step);
00431       if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00432 
00433       if ( ! commOnly ) {
00434         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00435       }
00436 
00437         // rattle2(timestep,step);
00438 
00439         submitReductions(step);
00440         submitCollections(step);
00441        //Update adaptive tempering temperature
00442         adaptTempUpdate(step);
00443 
00444       // Multigrator temperature and pressure steps
00445       multigratorTemperature(step, 1);
00446       multigratorPressure(step, 1);
00447       multigratorPressure(step, 2);
00448       multigratorTemperature(step, 2);
00449       doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
00450 
00451 #if CYCLE_BARRIER
00452         cycleBarrier(!((step+1) % stepsPerCycle), step);
00453 #elif PME_BARRIER
00454         cycleBarrier(doFullElectrostatics, step);
00455 #elif  STEP_BARRIER
00456         cycleBarrier(1, step);
00457 #endif
00458 
00459          if(Node::Object()->specialTracing || simParams->statsOn){
00460                  int bstep = simParams->traceStartStep;
00461                  int estep = bstep + simParams->numTraceSteps;
00462                  if(step == bstep || step == estep){
00463                          traceBarrier(step);
00464                  }                       
00465          }
00466 
00467 #ifdef MEASURE_NAMD_WITH_PAPI    
00468          if(simParams->papiMeasure) {
00469                  int bstep = simParams->papiMeasureStartStep;
00470                  int estep = bstep + simParams->numPapiMeasureSteps;
00471                  if(step == bstep || step==estep) {
00472                          papiMeasureBarrier(step);
00473                  }
00474          }
00475 #endif
00476           
00477         if(traceIsOn()){
00478             traceUserEvent(eventEndOfTimeStep);
00479             sprintf(traceNote, "%s%d",tracePrefix,step); 
00480             traceUserSuppliedNote(traceNote);
00481         }
00482         rebalanceLoad(step);
00483 
00484 #if PME_BARRIER
00485         // a step before PME
00486         cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00487 #endif
00488 
00489 #if USE_HPM
00490         if(step == START_HPM_STEP)
00491           (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
00492 
00493         if(step == STOP_HPM_STEP)
00494           (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
00495 #endif
00496 
00497     }
00498 }

void Sequencer::langevinPiston ( int   )  [protected]

Definition at line 1327 of file Sequencer.C.

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

Referenced by integrate().

01328 {
01329   if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
01330   {
01331    FullAtom *a = patch->atom.begin();
01332    int numAtoms = patch->numAtoms;
01333    Tensor factor = broadcast->positionRescaleFactor.get(step);
01334    // JCP FIX THIS!!!
01335    Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
01336    patch->lattice.rescale(factor);
01337    Molecule *mol = Node::Object()->molecule;
01338    if ( simParams->useGroupPressure )
01339    {
01340     int hgs;
01341     for ( int i = 0; i < numAtoms; i += hgs ) {
01342       int j;
01343       hgs = a[i].hydrogenGroupSize;
01344       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
01345         for ( j = i; j < (i+hgs); ++j ) {
01346           a[j].position = patch->lattice.apply_transform(
01347                                 a[j].fixedPosition,a[j].transform);
01348         }
01349         continue;
01350       }
01351       BigReal m_cm = 0;
01352       Position x_cm(0,0,0);
01353       Velocity v_cm(0,0,0);
01354       for ( j = i; j < (i+hgs); ++j ) {
01355         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01356         m_cm += a[j].mass;
01357         x_cm += a[j].mass * a[j].position;
01358         v_cm += a[j].mass * a[j].velocity;
01359       }
01360       x_cm /= m_cm;
01361       Position new_x_cm = x_cm;
01362       patch->lattice.rescale(new_x_cm,factor);
01363       Position delta_x_cm = new_x_cm - x_cm;
01364       v_cm /= m_cm;
01365       Velocity delta_v_cm;
01366       delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
01367       delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
01368       delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
01369       for ( j = i; j < (i+hgs); ++j ) {
01370         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01371           a[j].position = patch->lattice.apply_transform(
01372                                 a[j].fixedPosition,a[j].transform);
01373           continue;
01374         }
01375         if ( mol->is_atom_exPressure(a[j].id) ) continue;
01376         a[j].position += delta_x_cm;
01377         a[j].velocity += delta_v_cm;
01378       }
01379     }
01380    }
01381    else
01382    {
01383     for ( int i = 0; i < numAtoms; ++i )
01384     {
01385       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01386         a[i].position = patch->lattice.apply_transform(
01387                                 a[i].fixedPosition,a[i].transform);
01388         continue;
01389       }
01390       if ( mol->is_atom_exPressure(a[i].id) ) continue;
01391       patch->lattice.rescale(a[i].position,factor);
01392       a[i].velocity.x *= velFactor.x;
01393       a[i].velocity.y *= velFactor.y;
01394       a[i].velocity.z *= velFactor.z;
01395     }
01396    }
01397   }
01398 }

void Sequencer::langevinVelocities ( BigReal   )  [protected]

Definition at line 1069 of file Sequencer.C.

References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, adaptTempT, HomePatch::atom, ResizeArray< Elem >::begin(), BOLTZMANN, colvarproxy_namd::dt(), Random::gaussian_vector(), SimParameters::langevin_useBAOAB, SimParameters::langevinOn, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, Node::molecule, Patch::numAtoms, Node::Object(), partition(), patch, random, simParams, and FullAtom::velocity.

Referenced by integrate().

01070 {
01071 // This routine is used for the BAOAB integrator,
01072 // Ornstein-Uhlenbeck exact solve for the O-part.
01073 // See B. Leimkuhler and C. Matthews, AMRX (2012)
01074 // Routine originally written by JPhillips, with fresh errors by CMatthews June2012
01075 
01076   if ( simParams->langevinOn && simParams->langevin_useBAOAB )
01077   {
01078     FullAtom *a = patch->atom.begin();
01079     int numAtoms = patch->numAtoms;
01080     Molecule *molecule = Node::Object()->molecule;
01081     BigReal dt = dt_fs * 0.001;  // convert to ps
01082     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
01083     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
01084     {
01085         kbT = BOLTZMANN*adaptTempT;
01086     }
01087 
01088     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01089     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01090 
01091     for ( int i = 0; i < numAtoms; ++i )
01092     {
01093       BigReal dt_gamma = dt * a[i].langevinParam;
01094       if ( ! dt_gamma ) continue;
01095 
01096       BigReal f1 = exp( -dt_gamma );
01097       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT * 
01098                          ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
01099       a[i].velocity *= f1;
01100       a[i].velocity += f2 * random->gaussian_vector();
01101     }
01102   }
01103 }

void Sequencer::langevinVelocitiesBBK1 ( BigReal   )  [protected]

Definition at line 1105 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::drudeOn, colvarproxy_namd::dt(), SimParameters::langevin_useBAOAB, SimParameters::langevinOn, FullAtom::langevinParam, Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, and FullAtom::velocity.

Referenced by integrate().

01106 {
01107   if ( simParams->langevinOn && !simParams->langevin_useBAOAB )
01108   {
01109     FullAtom *a = patch->atom.begin();
01110     int numAtoms = patch->numAtoms;
01111     Molecule *molecule = Node::Object()->molecule;
01112     BigReal dt = dt_fs * 0.001;  // convert to ps
01113     int i;
01114 
01115     if (simParams->drudeOn) {
01116       for (i = 0;  i < numAtoms;  i++) {
01117 
01118         if (i < numAtoms-1 &&
01119             a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01120           //printf("*** Found Drude particle %d\n", a[i+1].id);
01121           // i+1 is a Drude particle with parent i
01122 
01123           // convert from Cartesian coordinates to (COM,bond) coordinates
01124           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
01125           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
01126           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
01127           BigReal dt_gamma;
01128 
01129           // use Langevin damping factor i for v_com
01130           dt_gamma = dt * a[i].langevinParam;
01131           if (dt_gamma != 0.0) {
01132             v_com *= ( 1. - 0.5 * dt_gamma );
01133           }
01134 
01135           // use Langevin damping factor i+1 for v_bnd
01136           dt_gamma = dt * a[i+1].langevinParam;
01137           if (dt_gamma != 0.0) {
01138             v_bnd *= ( 1. - 0.5 * dt_gamma );
01139           }
01140 
01141           // convert back
01142           a[i].velocity = v_com - m * v_bnd;
01143           a[i+1].velocity = v_bnd + a[i].velocity;
01144 
01145           i++;  // +1 from loop, we've updated both particles
01146         }
01147         else {
01148           BigReal dt_gamma = dt * a[i].langevinParam;
01149           if ( ! dt_gamma ) continue;
01150 
01151           a[i].velocity *= ( 1. - 0.5 * dt_gamma );
01152         }
01153 
01154       } // end for
01155     } // end if drudeOn
01156     else {
01157 
01158       for ( i = 0; i < numAtoms; ++i )
01159       {
01160         BigReal dt_gamma = dt * a[i].langevinParam;
01161         if ( ! dt_gamma ) continue;
01162 
01163         a[i].velocity *= ( 1. - 0.5 * dt_gamma );
01164       }
01165 
01166     } // end else
01167 
01168   } // end if langevinOn
01169 }

void Sequencer::langevinVelocitiesBBK2 ( BigReal   )  [protected]

Definition at line 1172 of file Sequencer.C.

References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, adaptTempT, HomePatch::atom, ResizeArray< Elem >::begin(), BOLTZMANN, SimParameters::drudeOn, SimParameters::drudeTemp, colvarproxy_namd::dt(), Random::gaussian_vector(), SimParameters::langevin_useBAOAB, SimParameters::langevinOn, FullAtom::langevinParam, SimParameters::langevinTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, FullAtom::mass, Node::molecule, Patch::numAtoms, Node::Object(), partition(), patch, random, rattle1(), simParams, and FullAtom::velocity.

Referenced by integrate().

01173 {
01174   if ( simParams->langevinOn && !simParams->langevin_useBAOAB ) 
01175   {
01176     rattle1(dt_fs,1);  // conserve momentum if gammas differ
01177 
01178     FullAtom *a = patch->atom.begin();
01179     int numAtoms = patch->numAtoms;
01180     Molecule *molecule = Node::Object()->molecule;
01181     BigReal dt = dt_fs * 0.001;  // convert to ps
01182     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
01183     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
01184     {
01185         kbT = BOLTZMANN*adaptTempT;
01186     }
01187     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01188     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01189     int i;
01190 
01191     if (simParams->drudeOn) {
01192       BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp);  // drude bond Temp
01193 
01194       for (i = 0;  i < numAtoms;  i++) {
01195 
01196         if (i < numAtoms-1 &&
01197             a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01198           //printf("*** Found Drude particle %d\n", a[i+1].id);
01199           // i+1 is a Drude particle with parent i
01200 
01201           // convert from Cartesian coordinates to (COM,bond) coordinates
01202           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
01203           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
01204           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
01205           BigReal dt_gamma;
01206 
01207           // use Langevin damping factor i for v_com
01208           dt_gamma = dt * a[i].langevinParam;
01209           if (dt_gamma != 0.0) {
01210             BigReal mass = a[i].mass + a[i+1].mass;
01211             v_com += random->gaussian_vector() *
01212               sqrt( 2 * dt_gamma * kbT *
01213                   ( a[i].partition ? tempFactor : 1.0 ) / mass );
01214             v_com /= ( 1. + 0.5 * dt_gamma );
01215           }
01216 
01217           // use Langevin damping factor i+1 for v_bnd
01218           dt_gamma = dt * a[i+1].langevinParam;
01219           if (dt_gamma != 0.0) {
01220             BigReal mass = a[i+1].mass * (1. - m);
01221             v_bnd += random->gaussian_vector() *
01222               sqrt( 2 * dt_gamma * kbT_bnd *
01223                   ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
01224             v_bnd /= ( 1. + 0.5 * dt_gamma );
01225           }
01226 
01227           // convert back
01228           a[i].velocity = v_com - m * v_bnd;
01229           a[i+1].velocity = v_bnd + a[i].velocity;
01230 
01231           i++;  // +1 from loop, we've updated both particles
01232         }
01233         else { 
01234           BigReal dt_gamma = dt * a[i].langevinParam;
01235           if ( ! dt_gamma ) continue;
01236 
01237           a[i].velocity += random->gaussian_vector() *
01238             sqrt( 2 * dt_gamma * kbT *
01239                 ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
01240           a[i].velocity /= ( 1. + 0.5 * dt_gamma );
01241         }
01242 
01243       } // end for
01244     } // end if drudeOn
01245     else {
01246 
01247       for ( i = 0; i < numAtoms; ++i )
01248       {
01249         BigReal dt_gamma = dt * a[i].langevinParam;
01250         if ( ! dt_gamma ) continue;
01251 
01252         a[i].velocity += random->gaussian_vector() *
01253           sqrt( 2 * dt_gamma * kbT *
01254               ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass );
01255         a[i].velocity /= ( 1. + 0.5 * dt_gamma );
01256       }
01257 
01258     } // end else
01259 
01260   } // end if langevinOn
01261 }

void Sequencer::maximumMove ( BigReal   )  [protected]

Definition at line 1641 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::cutoff, colvarproxy_namd::dt(), Node::enableEarlyExit(), endi(), CompAtomExt::id, iERROR(), iout, SimParameters::maximumMove, Patch::numAtoms, Node::Object(), patch, Patch::patchID, PDBVELFACTOR, simParams, terminate(), and TIMEFACTOR.

Referenced by integrate().

01642 {
01643   FullAtom *a = patch->atom.begin();
01644   int numAtoms = patch->numAtoms;
01645   if ( simParams->maximumMove ) {
01646     const BigReal dt = timestep / TIMEFACTOR;
01647     const BigReal maxvel = simParams->maximumMove / dt;
01648     const BigReal maxvel2 = maxvel * maxvel;
01649     for ( int i=0; i<numAtoms; ++i ) {
01650       if ( a[i].velocity.length2() > maxvel2 ) {
01651         a[i].velocity *= ( maxvel / a[i].velocity.length() );
01652       }
01653     }
01654   } else {
01655     const BigReal dt = timestep / TIMEFACTOR;
01656     const BigReal maxvel = simParams->cutoff / dt;
01657     const BigReal maxvel2 = maxvel * maxvel;
01658     int killme = 0;
01659     for ( int i=0; i<numAtoms; ++i ) {
01660       killme = killme || ( a[i].velocity.length2() > maxvel2 );
01661     }
01662     if ( killme ) {
01663       killme = 0;
01664       for ( int i=0; i<numAtoms; ++i ) {
01665         if ( a[i].velocity.length2() > maxvel2 ) {
01666           ++killme;
01667           iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
01668             << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
01669             << ( PDBVELFACTOR * maxvel ) << ", atom "
01670             << i << " of " << numAtoms << " on patch "
01671             << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
01672         }
01673       }
01674       iout << iERROR << 
01675         "Atoms moving too fast; simulation has become unstable ("
01676         << killme << " atoms on patch " << patch->patchID
01677         << " pe " << CkMyPe() << ").\n" << endi;
01678       Node::Object()->enableEarlyExit();
01679       terminate();
01680     }
01681   }
01682 }

void Sequencer::minimizationQuenchVelocity ( void   )  [protected]

Definition at line 1684 of file Sequencer.C.

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

Referenced by integrate().

01685 {
01686   if ( simParams->minimizeOn ) {
01687     FullAtom *a = patch->atom.begin();
01688     int numAtoms = patch->numAtoms;
01689     for ( int i=0; i<numAtoms; ++i ) {
01690       a[i].velocity = 0.;
01691     }
01692   }
01693 }

void Sequencer::minimize (  )  [protected]

Definition at line 543 of file Sequencer.C.

References HomePatch::atom, Patch::atomMapper, ResizeArray< Elem >::begin(), broadcast, SimParameters::colvarsOn, ComputeMgr::computeGlobalObject, Node::computeMgr, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, Flags::doNonbonded, ResizeArray< Elem >::end(), SimParameters::firstTimestep, Patch::flags, SimParameters::fullElectFrequency, SimParameters::GBISOn, SimParameters::LCPOOn, SimParameters::lonepairs, Flags::maxForceMerged, Flags::maxForceUsed, ControllerBroadcasts::minimizeCoefficient, minimizeMoveDownhill(), SimParameters::mollyOn, SimParameters::N, Results::nbond, newMinimizeDirection(), newMinimizePosition(), Results::normal, Node::Object(), patch, quenchVelocities(), HomePatch::rattle1(), rebalanceLoad(), AtomMapper::registerIDsFullAtom(), runComputeObjects(), saveForce(), ComputeGlobal::saveTotalForces(), simParams, Results::slow, Flags::step, GlobalMaster::step, SimParameters::stepsPerCycle, submitCollections(), submitMinimizeReductions(), SimParameters::tclForcesOn, and TIMEFACTOR.

Referenced by algorithm().

00543                          {
00544   const int numberOfSteps = simParams->N;
00545   const int stepsPerCycle = simParams->stepsPerCycle;
00546   int &step = patch->flags.step;
00547   step = simParams->firstTimestep;
00548 
00549   int &maxForceUsed = patch->flags.maxForceUsed;
00550   int &maxForceMerged = patch->flags.maxForceMerged;
00551   maxForceUsed = Results::normal;
00552   maxForceMerged = Results::normal;
00553   int &doNonbonded = patch->flags.doNonbonded;
00554   doNonbonded = 1;
00555   maxForceUsed = Results::nbond;
00556   maxForceMerged = Results::nbond;
00557   const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00558   int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00559   doFullElectrostatics = dofull;
00560   if ( dofull ) {
00561     maxForceMerged = Results::slow;
00562     maxForceUsed = Results::slow;
00563   }
00564   int &doMolly = patch->flags.doMolly;
00565   doMolly = simParams->mollyOn && doFullElectrostatics;
00566   // BEGIN LA
00567   int &doLoweAndersen = patch->flags.doLoweAndersen;
00568   doLoweAndersen = 0;
00569   // END LA
00570 
00571   int &doGBIS = patch->flags.doGBIS;
00572   doGBIS = simParams->GBISOn;
00573 
00574     int &doLCPO = patch->flags.doLCPO;
00575     doLCPO = simParams->LCPOOn;
00576 
00577     int doTcl = simParams->tclForcesOn;
00578         int doColvars = simParams->colvarsOn;
00579     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00580 
00581   int &doEnergy = patch->flags.doEnergy;
00582   doEnergy = 1;
00583 
00584   patch->rattle1(0.,0,0);  // enforce rigid bond constraints on initial positions
00585 
00586   if (simParams->lonepairs) {
00587     patch->atomMapper->registerIDsFullAtom(
00588                 patch->atom.begin(),patch->atom.end());
00589   }
00590 
00591   runComputeObjects(1,step<numberOfSteps); // must migrate here!
00592 
00593   if ( doTcl || doColvars ) {
00594     if ( doNonbonded ) saveForce(Results::nbond);
00595     if ( doFullElectrostatics ) saveForce(Results::slow);
00596     computeGlobal->saveTotalForces(patch);
00597   }
00598 
00599   BigReal fmax2 = TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00600 
00601   submitMinimizeReductions(step,fmax2);
00602   rebalanceLoad(step);
00603 
00604   int downhill = 1;  // start out just fixing bad contacts
00605   int minSeq = 0;
00606   for ( ++step; step <= numberOfSteps; ++step ) {
00607    BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00608    if ( downhill ) {
00609     if ( c ) minimizeMoveDownhill(fmax2);
00610     else {
00611       downhill = 0;
00612       fmax2 *= 10000.;
00613     }
00614    }
00615    if ( ! downhill ) {
00616     if ( ! c ) {  // new direction
00617       c = broadcast->minimizeCoefficient.get(minSeq++);
00618       newMinimizeDirection(c);  // v = c * v + f
00619       c = broadcast->minimizeCoefficient.get(minSeq++);
00620     }  // same direction
00621     newMinimizePosition(c);  // x = x + c * v
00622    }
00623 
00624     runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00625     if ( doTcl || doColvars ) {
00626       if ( doNonbonded ) saveForce(Results::nbond);
00627       if ( doFullElectrostatics ) saveForce(Results::slow);
00628       computeGlobal->saveTotalForces(patch);
00629     }
00630     submitMinimizeReductions(step,fmax2);
00631     submitCollections(step, 1);  // write out zeros for velocities
00632     rebalanceLoad(step);
00633   }
00634   quenchVelocities();  // zero out bogus velocity
00635 }

void Sequencer::minimizeMoveDownhill ( BigReal  fmax2  )  [protected]

Definition at line 638 of file Sequencer.C.

References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), f, Patch::f, SimParameters::fixedAtomsOn, CompAtom::hydrogenGroupSize, if(), j, Results::normal, Patch::numAtoms, patch, HomePatch::rattle1(), and simParams.

Referenced by minimize().

00638                                                   {
00639 
00640   FullAtom *a = patch->atom.begin();
00641   Force *f1 = patch->f[Results::normal].begin();  // includes nbond and slow
00642   int numAtoms = patch->numAtoms;
00643 
00644   for ( int i = 0; i < numAtoms; ++i ) {
00645     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
00646     Force f = f1[i];
00647     if ( f.length2() > fmax2 ) {
00648       a[i].position += ( 0.1 * f.unit() );
00649       int hgs = a[i].hydrogenGroupSize;  // 0 if not parent
00650       for ( int j=1; j<hgs; ++j ) {
00651         a[++i].position += ( 0.1 * f.unit() );
00652       }
00653     }
00654   }
00655 
00656   patch->rattle1(0.,0,0);
00657 }

void Sequencer::multigratorPressure ( int  step,
int  callNumber 
) [protected]

Definition at line 843 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), broadcast, calcFixVirial(), Flags::doFullElectrostatics, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, SimpleBroadcastObject< T >::get(), Tensor::identity(), SubmitReduction::item(), j, Patch::lattice, Vector::length2(), HomePatch::marginViolations, FullAtom::mass, SimParameters::multigratorOn, SimParameters::multigratorPressureFreq, SimParameters::N, NAMD_bug(), Results::nbond, Results::normal, Patch::numAtoms, Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, partition(), patch, CompAtom::position, ControllerBroadcasts::positionRescaleFactor, ControllerBroadcasts::positionRescaleFactor2, reduction, REDUCTION_ATOM_CHECKSUM, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_EXT_FORCE_NBOND, REDUCTION_EXT_FORCE_NORMAL, REDUCTION_EXT_FORCE_SLOW, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_MARGIN_VIOLATIONS, REDUCTION_MOMENTUM_SQUARED, REDUCTION_VIRIAL_NBOND, REDUCTION_VIRIAL_NORMAL, REDUCTION_VIRIAL_SLOW, Lattice::rescale(), runComputeObjects(), scalePositionsVelocities(), simParams, Results::slow, SimParameters::stepsPerCycle, SubmitReduction::submit(), SimParameters::useGroupPressure, FullAtom::velocity, ControllerBroadcasts::velocityRescaleTensor, and ControllerBroadcasts::velocityRescaleTensor2.

Referenced by integrate().

00843                                                             {
00844 // Calculate new positions, momenta, and volume using positionRescaleFactor and 
00845 // velocityRescaleTensor values returned from Controller::multigratorPressureCalcScale()
00846   if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
00847     FullAtom *a = patch->atom.begin();
00848     int numAtoms = patch->numAtoms;
00849 
00850     // Receive scaling factors from Controller
00851     Tensor scaleTensor    = (callNumber == 1) ? broadcast->positionRescaleFactor.get(step) : broadcast->positionRescaleFactor2.get(step);
00852     Tensor velScaleTensor = (callNumber == 1) ? broadcast->velocityRescaleTensor.get(step) : broadcast->velocityRescaleTensor2.get(step);
00853     Tensor posScaleTensor = scaleTensor;
00854     posScaleTensor -= Tensor::identity();
00855     if (simParams->useGroupPressure) {
00856       velScaleTensor -= Tensor::identity();
00857     }
00858 
00859     // Scale volume
00860     patch->lattice.rescale(scaleTensor);
00861     // Scale positions and velocities
00862     scalePositionsVelocities(posScaleTensor, velScaleTensor);
00863 
00864     if (!patch->flags.doFullElectrostatics) NAMD_bug("Sequencer::multigratorPressure, doFullElectrostatics must be true");
00865 
00866     // Calculate new forces
00867     // NOTE: We should not need to migrate here since any migration should have happened in the
00868     // previous call to runComputeObjects inside the MD loop in Sequencer::integrate()
00869     const int numberOfSteps = simParams->N;
00870     const int stepsPerCycle = simParams->stepsPerCycle;
00871     runComputeObjects(0 , step<numberOfSteps, 1);
00872 
00873     reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
00874     reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
00875 
00876     // Virials etc.
00877     Tensor virialNormal;
00878     Tensor momentumSqrSum;
00879     BigReal kineticEnergy = 0;
00880     if ( simParams->pairInteractionOn ) {
00881       if ( simParams->pairInteractionSelf ) {
00882         for ( int i = 0; i < numAtoms; ++i ) {
00883           if ( a[i].partition != 1 ) continue;
00884           kineticEnergy += a[i].mass * a[i].velocity.length2();
00885           virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
00886         }
00887       }
00888     } else {
00889       for ( int i = 0; i < numAtoms; ++i ) {
00890         if (a[i].mass < 0.01) continue;
00891         kineticEnergy += a[i].mass * a[i].velocity.length2();
00892         virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
00893       }
00894     }
00895     if (!simParams->useGroupPressure) momentumSqrSum = virialNormal;
00896     kineticEnergy *= 0.5;
00897     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
00898     ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virialNormal);
00899 
00900     if ( simParams->fixedAtomsOn ) {
00901       Tensor fixVirialNormal;
00902       Tensor fixVirialNbond;
00903       Tensor fixVirialSlow;
00904       Vector fixForceNormal = 0;
00905       Vector fixForceNbond = 0;
00906       Vector fixForceSlow = 0;
00907 
00908       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
00909 
00910       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, fixVirialNormal);
00911       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, fixVirialNbond);
00912       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, fixVirialSlow);
00913       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NORMAL, fixForceNormal);
00914       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NBOND, fixForceNbond);
00915       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_SLOW, fixForceSlow);
00916     }
00917 
00918     // Internal virial and group momentum
00919     Tensor intVirialNormal;
00920     Tensor intVirialNormal2;
00921     Tensor intVirialNbond;
00922     Tensor intVirialSlow;
00923     int hgs;
00924     for ( int i = 0; i < numAtoms; i += hgs ) {
00925       hgs = a[i].hydrogenGroupSize;
00926       int j;
00927       BigReal m_cm = 0;
00928       Position x_cm(0,0,0);
00929       Velocity v_cm(0,0,0);
00930       for ( j = i; j < (i+hgs); ++j ) {
00931         m_cm += a[j].mass;
00932         x_cm += a[j].mass * a[j].position;
00933         v_cm += a[j].mass * a[j].velocity;
00934       }
00935       if (simParams->useGroupPressure) momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
00936       x_cm /= m_cm;
00937       v_cm /= m_cm;
00938       if (simParams->fixedAtomsOn) NAMD_bug("Sequencer::multigratorPressure, simParams->fixedAtomsOn not implemented yet");
00939       if ( simParams->pairInteractionOn ) {
00940         if ( simParams->pairInteractionSelf ) {
00941           NAMD_bug("Sequencer::multigratorPressure, this part needs to be implemented correctly");
00942           for ( j = i; j < (i+hgs); ++j ) {
00943             if ( a[j].partition != 1 ) continue;
00944             BigReal mass = a[j].mass;
00945             Vector v = a[j].velocity;
00946             Vector dv = v - v_cm;
00947             intVirialNormal2.outerAdd (mass, v, dv);
00948             Vector dx = a[j].position - x_cm;
00949             intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
00950             intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
00951             intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
00952           }
00953         }
00954       } else {
00955         for ( j = i; j < (i+hgs); ++j ) {
00956           BigReal mass = a[j].mass;
00957           Vector v = a[j].velocity;
00958           Vector dv = v - v_cm;
00959           intVirialNormal2.outerAdd(mass, v, dv);
00960           Vector dx = a[j].position - x_cm;
00961           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
00962           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
00963           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
00964         }
00965       }
00966     }
00967 
00968     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
00969     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal2);
00970     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NBOND, intVirialNbond);
00971     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_SLOW, intVirialSlow);
00972     ADD_TENSOR_OBJECT(reduction, REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
00973 
00974     reduction->submit();
00975   }
00976 }

void Sequencer::multigratorTemperature ( int  step,
int  callNumber 
) [protected]

Definition at line 1006 of file Sequencer.C.

References ADD_TENSOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), broadcast, calcKineticEnergy(), SubmitReduction::item(), j, FullAtom::mass, MULTIGRATOR_REDUCTION_KINETIC_ENERGY, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, SimParameters::multigratorOn, SimParameters::multigratorPressureFreq, multigratorReduction, SimParameters::multigratorTemperatureFreq, Patch::numAtoms, Tensor::outerAdd(), patch, CompAtom::position, scaleVelocities(), simParams, SubmitReduction::submit(), SimParameters::useGroupPressure, FullAtom::velocity, ControllerBroadcasts::velocityRescaleFactor, and ControllerBroadcasts::velocityRescaleFactor2.

Referenced by integrate().

01006                                                                {
01007   if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
01008     // Scale velocities
01009     BigReal velScale = (callNumber == 1) ? broadcast->velocityRescaleFactor.get(step) : broadcast->velocityRescaleFactor2.get(step);
01010     scaleVelocities(velScale);
01011     // Calculate new kineticEnergy
01012     BigReal kineticEnergy = calcKineticEnergy();
01013     multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY) += kineticEnergy;
01014     if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
01015       // If this is a pressure cycle, calculate new momentum squared sum
01016       FullAtom *a = patch->atom.begin();
01017       int numAtoms = patch->numAtoms;
01018       Tensor momentumSqrSum;
01019       if (simParams->useGroupPressure) {
01020         int hgs;
01021         for ( int i = 0; i < numAtoms; i += hgs ) {
01022           hgs = a[i].hydrogenGroupSize;
01023           int j;
01024           BigReal m_cm = 0;
01025           Position x_cm(0,0,0);
01026           Velocity v_cm(0,0,0);
01027           for ( j = i; j < (i+hgs); ++j ) {
01028             m_cm += a[j].mass;
01029             x_cm += a[j].mass * a[j].position;
01030             v_cm += a[j].mass * a[j].velocity;
01031           }
01032           momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
01033         }
01034       } else {
01035         for ( int i = 0; i < numAtoms; i++) {
01036           momentumSqrSum.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01037         }
01038       }
01039       ADD_TENSOR_OBJECT(multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
01040     }
01041     // Submit reductions (kineticEnergy and, if applicable, momentumSqrSum)
01042     multigratorReduction->submit();
01043 
01044   }
01045 }

void Sequencer::newMinimizeDirection ( BigReal   )  [protected]

Definition at line 660 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::drudeHardWallOn, Patch::f, SimParameters::fixedAtomsOn, j, Vector::length2(), SubmitReduction::max(), min_reduction, HomePatch::minimize_rattle2(), Results::normal, Patch::numAtoms, patch, simParams, SubmitReduction::submit(), TIMEFACTOR, and FullAtom::velocity.

Referenced by minimize().

00660                                               {
00661   FullAtom *a = patch->atom.begin();
00662   Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
00663   const bool fixedAtomsOn = simParams->fixedAtomsOn;
00664   const bool drudeHardWallOn = simParams->drudeHardWallOn;
00665   int numAtoms = patch->numAtoms;
00666   BigReal maxv2 = 0.;
00667 
00668   for ( int i = 0; i < numAtoms; ++i ) {
00669     a[i].velocity *= c;
00670     a[i].velocity += f1[i];
00671     if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00672       a[i].velocity = a[i-1].velocity;
00673     }
00674     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00675     BigReal v2 = a[i].velocity.length2();
00676     if ( v2 > maxv2 ) maxv2 = v2;
00677   }
00678 
00679   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial); }
00680 
00681   maxv2 = 0.;
00682   for ( int i = 0; i < numAtoms; ++i ) {
00683     if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00684       a[i].velocity = a[i-1].velocity;
00685       if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00686     }
00687     BigReal v2 = a[i].velocity.length2();
00688     if ( v2 > maxv2 ) maxv2 = v2;
00689   }
00690 
00691   min_reduction->max(0,maxv2);
00692   min_reduction->submit();
00693 
00694   // prevent hydrogens from being left behind
00695   BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00696   // int adjustCount = 0;
00697   int hgs;
00698   for ( int i = 0; i < numAtoms; i += hgs ) {
00699     hgs = a[i].hydrogenGroupSize;
00700     BigReal minChildVel = a[i].velocity.length2();
00701     if ( minChildVel < fmax2 ) continue;
00702     int adjustChildren = 1;
00703     for ( int j = i+1; j < (i+hgs); ++j ) {
00704       if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
00705     }
00706     if ( adjustChildren ) {
00707       // if ( hgs > 1 ) ++adjustCount;
00708       for ( int j = i+1; j < (i+hgs); ++j ) {
00709         if (a[i].mass < 0.01) continue;  // lone pair
00710         a[j].velocity = a[i].velocity;
00711       }
00712     }
00713   }
00714   // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
00715 
00716 }

void Sequencer::newMinimizePosition ( BigReal   )  [protected]

Definition at line 719 of file Sequencer.C.

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

Referenced by minimize().

00719                                              {
00720   FullAtom *a = patch->atom.begin();
00721   int numAtoms = patch->numAtoms;
00722 
00723   for ( int i = 0; i < numAtoms; ++i ) {
00724     a[i].position += c * a[i].velocity;
00725   }
00726 
00727   if ( simParams->drudeHardWallOn ) {
00728     for ( int i = 1; i < numAtoms; ++i ) {
00729       if ( (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00730         a[i].position -= a[i-1].position;
00731       }
00732     }
00733   }
00734 
00735   patch->rattle1(0.,0,0);
00736 
00737   if ( simParams->drudeHardWallOn ) {
00738     for ( int i = 1; i < numAtoms; ++i ) {
00739       if ( (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
00740         a[i].position += a[i-1].position;
00741       }
00742     }
00743   }
00744 }

void Sequencer::newtonianVelocities ( BigReal  ,
const   BigReal,
const   BigReal,
const   BigReal,
const   int,
const   int,
const   int 
) [protected]

Definition at line 1048 of file Sequencer.C.

References addForceToMomentum(), addForceToMomentum3(), Results::nbond, Results::normal, and Results::slow.

Referenced by integrate().

01054 {
01055   // Deterministic velocity update, account for multigrator
01056   if (staleForces || (doNonbonded && doFullElectrostatics)) {
01057     addForceToMomentum3(stepscale*timestep, Results::normal, 0,
01058                         stepscale*nbondstep, Results::nbond, staleForces,
01059                         stepscale*slowstep, Results::slow, staleForces);
01060   } else {
01061     addForceToMomentum(stepscale*timestep);
01062     if (staleForces || doNonbonded)
01063       addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces);
01064     if (staleForces || doFullElectrostatics)
01065       addForceToMomentum(stepscale*slowstep, Results::slow, staleForces);
01066   }
01067 }

void Sequencer::quenchVelocities (  )  [protected]

Definition at line 747 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, and patch.

Referenced by minimize().

00747                                  {
00748   FullAtom *a = patch->atom.begin();
00749   int numAtoms = patch->numAtoms;
00750 
00751   for ( int i = 0; i < numAtoms; ++i ) {
00752     a[i].velocity = 0;
00753   }
00754 }

void Sequencer::rattle1 ( BigReal  ,
int   
) [protected]

Definition at line 1612 of file Sequencer.C.

References ADD_TENSOR_OBJECT, Node::enableEarlyExit(), endi(), iERROR(), iout, Node::Object(), patch, pressureProfileReduction, HomePatch::rattle1(), reduction, REDUCTION_VIRIAL_NORMAL, RIGID_NONE, SimParameters::rigidBonds, simParams, and terminate().

Referenced by integrate(), and langevinVelocitiesBBK2().

01613 {
01614   if ( simParams->rigidBonds != RIGID_NONE ) {
01615     Tensor virial;
01616     Tensor *vp = ( pressure ? &virial : 0 );
01617     if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
01618       iout << iERROR << 
01619         "Constraint failure; simulation has become unstable.\n" << endi;
01620       Node::Object()->enableEarlyExit();
01621       terminate();
01622     }
01623     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01624   }
01625 }

void Sequencer::reassignVelocities ( BigReal  ,
int   
) [protected]

Definition at line 1475 of file Sequencer.C.

References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BOLTZMANN, SimParameters::fixedAtomsOn, Random::gaussian_vector(), if(), SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, FullAtom::mass, NAMD_bug(), Patch::numAtoms, partition(), patch, random, SimParameters::reassignFreq, SimParameters::reassignHold, SimParameters::reassignIncr, SimParameters::reassignTemp, and simParams.

Referenced by integrate().

01476 {
01477   const int reassignFreq = simParams->reassignFreq;
01478   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01479     FullAtom *a = patch->atom.begin();
01480     int numAtoms = patch->numAtoms;
01481     BigReal newTemp = simParams->reassignTemp;
01482     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01483     if ( simParams->reassignIncr > 0.0 ) {
01484       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01485         newTemp = simParams->reassignHold;
01486     } else {
01487       if ( newTemp < simParams->reassignHold )
01488         newTemp = simParams->reassignHold;
01489     }
01490     BigReal kbT = BOLTZMANN * newTemp;
01491 
01492     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01493     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01494 
01495     for ( int i = 0; i < numAtoms; ++i )
01496     {
01497       a[i].velocity = ( ( simParams->fixedAtomsOn && a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
01498         sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
01499           * random->gaussian_vector() );
01500     }
01501   } else {
01502     NAMD_bug("Sequencer::reassignVelocities called improperly!");
01503   }
01504 }

void Sequencer::rebalanceLoad ( int  timestep  )  [protected]

Definition at line 2432 of file Sequencer.C.

References LdbCoordinator::getNumStepsToRun(), Patch::getPatchID(), ldbSteps, LdbCoordinator::Object(), pairlistsAreValid, patch, LdbCoordinator::rebalance(), and HomePatch::submitLoadStats().

Referenced by integrate(), and minimize().

02432                                           {
02433   if ( ! ldbSteps ) {
02434     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02435   }
02436   if ( ! --ldbSteps ) {
02437     patch->submitLoadStats(timestep);
02438     ldbCoordinator->rebalance(this,patch->getPatchID());
02439     pairlistsAreValid = 0;
02440   }
02441 }

void Sequencer::reinitVelocities ( void   )  [protected]

Definition at line 1506 of file Sequencer.C.

References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BOLTZMANN, SimParameters::drudeOn, SimParameters::fixedAtomsOn, Random::gaussian_vector(), SimParameters::initialTemp, SimParameters::lesFactor, SimParameters::lesOn, SimParameters::lesReduceTemp, FullAtom::mass, Patch::numAtoms, partition(), patch, random, simParams, and FullAtom::velocity.

Referenced by algorithm().

01507 {
01508   FullAtom *a = patch->atom.begin();
01509   int numAtoms = patch->numAtoms;
01510   BigReal newTemp = simParams->initialTemp;
01511   BigReal kbT = BOLTZMANN * newTemp;
01512 
01513   int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01514   BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01515 
01516   for ( int i = 0; i < numAtoms; ++i )
01517   {
01518     a[i].velocity = ( ( (simParams->fixedAtomsOn && a[i].atomFixed) || a[i].mass <= 0.) ? Vector(0,0,0) :
01519       sqrt( kbT * ( a[i].partition ? tempFactor : 1.0 ) / a[i].mass )
01520         * random->gaussian_vector() );
01521     if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass >= 0.001 ) {
01522       a[i+1].velocity = a[i].velocity;  // zero is good enough
01523       ++i;
01524     }
01525   }
01526 }

void Sequencer::reloadCharges (  )  [protected]

Definition at line 1538 of file Sequencer.C.

References HomePatch::atom, Molecule::atomcharge(), ResizeArray< Elem >::begin(), CompAtomExt::id, Node::molecule, Patch::numAtoms, Node::Object(), and patch.

Referenced by algorithm().

01539 {
01540   FullAtom *a = patch->atom.begin();
01541   int numAtoms = patch->numAtoms;
01542   Molecule *molecule = Node::Object()->molecule;
01543   for ( int i = 0; i < numAtoms; ++i )
01544   {
01545     a[i].charge = molecule->atomcharge(a[i].id);
01546   }
01547 }

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

Definition at line 1418 of file Sequencer.C.

References SimParameters::accelMDdihe, SimParameters::accelMDdual, SimParameters::accelMDLastStep, SimParameters::accelMDOn, ControllerBroadcasts::accelMDRescaleFactor, Results::amdf, broadcast, Patch::f, SimpleBroadcastObject< T >::get(), if(), NAMD_die(), Results::nbond, Results::normal, Patch::numAtoms, patch, simParams, and Results::slow.

Referenced by integrate().

01419 {
01420    if (!simParams->accelMDOn) return;
01421    if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
01422 
01423    Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
01424    const BigReal factor_dihe = accelMDfactor[0];
01425    const BigReal factor_tot  = accelMDfactor[1];
01426    const int numAtoms = patch->numAtoms;
01427 
01428    if (simParams->accelMDdihe && factor_tot <1 )
01429        NAMD_die("accelMD broadcasting error!\n");
01430    if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
01431        NAMD_die("accelMD broadcasting error!\n");
01432 
01433    if (simParams->accelMDdihe && factor_dihe < 1) {
01434         for (int i = 0; i < numAtoms; ++i)
01435           if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01436               patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);         
01437    }
01438 
01439    if ( !simParams->accelMDdihe && factor_tot < 1) {
01440         for (int i = 0; i < numAtoms; ++i)
01441             patch->f[Results::normal][i] *= factor_tot;
01442         if (doNonbonded) {
01443             for (int i = 0; i < numAtoms; ++i)
01444                  patch->f[Results::nbond][i] *= factor_tot;
01445         }
01446         if (doFullElectrostatics) {
01447             for (int i = 0; i < numAtoms; ++i)
01448                  patch->f[Results::slow][i] *= factor_tot;
01449         }
01450    }
01451 
01452    if (simParams->accelMDdual && factor_dihe < 1) {
01453         for (int i = 0; i < numAtoms; ++i)
01454            if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01455                patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
01456    }
01457 
01458 }

void Sequencer::rescaleVelocities ( int   )  [protected]

Definition at line 1400 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), broadcast, Patch::numAtoms, patch, SimParameters::rescaleFreq, rescaleVelocities_numTemps, simParams, and ControllerBroadcasts::velocityRescaleFactor.

Referenced by integrate().

01401 {
01402   const int rescaleFreq = simParams->rescaleFreq;
01403   if ( rescaleFreq > 0 ) {
01404     FullAtom *a = patch->atom.begin();
01405     int numAtoms = patch->numAtoms;
01406     ++rescaleVelocities_numTemps;
01407     if ( rescaleVelocities_numTemps == rescaleFreq ) {
01408       BigReal factor = broadcast->velocityRescaleFactor.get(step);
01409       for ( int i = 0; i < numAtoms; ++i )
01410       {
01411         a[i].velocity *= factor;
01412       }
01413       rescaleVelocities_numTemps = 0;
01414     }
01415   }
01416 }

void Sequencer::rescaleVelocitiesByFactor ( BigReal   )  [protected]

Definition at line 1528 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, and patch.

Referenced by algorithm().

01529 {
01530   FullAtom *a = patch->atom.begin();
01531   int numAtoms = patch->numAtoms;
01532   for ( int i = 0; i < numAtoms; ++i )
01533   {
01534     a[i].velocity *= factor;
01535   }
01536 }

void Sequencer::run ( void   ) 

Definition at line 97 of file Sequencer.C.

References awaken(), DebugM, Patch::getPatchID(), patch, PATCH_PRIORITY, and SEQ_STK_SZ.

Referenced by HomePatch::runSequencer().

00098 {
00099     // create a Thread and invoke it
00100     DebugM(4, "::run() - this = " << this << "\n" );
00101     thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
00102     CthSetStrategyDefault(thread);
00103     priority = PATCH_PRIORITY(patch->getPatchID());
00104     awaken();
00105 }

void Sequencer::runComputeObjects ( int  migration = 1,
int  pairlists = 0,
int  pressureStep = 0 
) [protected]

Definition at line 2259 of file Sequencer.C.

References Flags::doFullElectrostatics, Patch::flags, SimParameters::fullElectFrequency, if(), pairlistsAge, pairlistsAreValid, SimParameters::pairlistsPerCycle, patch, simParams, and SimParameters::stepsPerCycle.

Referenced by integrate(), minimize(), and multigratorPressure().

02260 {
02261   if ( migration ) pairlistsAreValid = 0;
02262 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
02263   if ( pairlistsAreValid &&
02264        ( patch->flags.doFullElectrostatics || ! simParams->fullElectFrequency )
02265                          && ( pairlistsAge > (
02266 #else
02267   if ( pairlistsAreValid && ( pairlistsAge > (
02268 #endif
02269          (simParams->stepsPerCycle - 1) / simParams->pairlistsPerCycle ) ) ) {
02270     pairlistsAreValid = 0;
02271   }
02272   if ( ! simParams->usePairlists ) pairlists = 0;
02273   patch->flags.usePairlists = pairlists || pairlistsAreValid;
02274   patch->flags.savePairlists =
02275         pairlists && ! pairlistsAreValid;
02276 
02277   if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
02278 
02279   patch->positionsReady(migration);  // updates flags.sequence
02280 
02281   int seq = patch->flags.sequence;
02282   int basePriority = ( (seq & 0xffff) << 15 )
02283                      + PATCH_PRIORITY(patch->getPatchID());
02284   if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
02285     priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
02286     suspend(); // until all deposit boxes close
02287     patch->gbisComputeAfterP1();
02288     priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
02289     suspend();
02290     patch->gbisComputeAfterP2();
02291     priority = basePriority + COMPUTE_HOME_PRIORITY;
02292     suspend();
02293   } else {
02294     priority = basePriority + COMPUTE_HOME_PRIORITY;
02295     suspend(); // until all deposit boxes close
02296   }
02297 
02298   if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
02299     pairlistsAreValid = 1;
02300     pairlistsAge = 0;
02301   }
02302   // For multigrator, do not age pairlist during pressure step
02303   // NOTE: for non-multigrator pressureStep = 0 always
02304   if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
02305 
02306   if (simParams->lonepairs) {
02307     {
02308       Tensor virial;
02309       patch->redistrib_lonepair_forces(Results::normal, &virial);
02310       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02311     }
02312     if (patch->flags.doNonbonded) {
02313       Tensor virial;
02314       patch->redistrib_lonepair_forces(Results::nbond, &virial);
02315       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02316     }
02317     if (patch->flags.doFullElectrostatics) {
02318       Tensor virial;
02319       patch->redistrib_lonepair_forces(Results::slow, &virial);
02320       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02321     }
02322   } else if (simParams->watmodel == WAT_TIP4) {
02323     {
02324       Tensor virial;
02325       patch->redistrib_tip4p_forces(Results::normal, &virial);
02326       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02327     }
02328     if (patch->flags.doNonbonded) {
02329       Tensor virial;
02330       patch->redistrib_tip4p_forces(Results::nbond, &virial);
02331       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02332     }
02333     if (patch->flags.doFullElectrostatics) {
02334       Tensor virial;
02335       patch->redistrib_tip4p_forces(Results::slow, &virial);
02336       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02337     }
02338   } else if (simParams->watmodel == WAT_SWM4) {
02339     {
02340       Tensor virial;
02341       patch->redistrib_swm4_forces(Results::normal, &virial);
02342       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02343     }
02344     if (patch->flags.doNonbonded) {
02345       Tensor virial;
02346       patch->redistrib_swm4_forces(Results::nbond, &virial);
02347       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02348     }
02349     if (patch->flags.doFullElectrostatics) {
02350       Tensor virial;
02351       patch->redistrib_swm4_forces(Results::slow, &virial);
02352       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02353     }
02354   }
02355 
02356   if ( patch->flags.doMolly ) {
02357     Tensor virial;
02358     patch->mollyMollify(&virial);
02359     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
02360   }
02361 
02362 
02363   // BEGIN LA
02364   if (patch->flags.doLoweAndersen) {
02365       patch->loweAndersenFinish();
02366   }
02367   // END LA
02368 #ifdef NAMD_CUDA_XXX
02369   int numAtoms = patch->numAtoms;
02370   FullAtom *a = patch->atom.begin();
02371   for ( int i=0; i<numAtoms; ++i ) {
02372     CkPrintf("%d %g %g %g\n", a[i].id,
02373         patch->f[Results::normal][i].x +
02374         patch->f[Results::nbond][i].x +
02375         patch->f[Results::slow][i].x,
02376         patch->f[Results::normal][i].y + 
02377         patch->f[Results::nbond][i].y +
02378         patch->f[Results::slow][i].y,
02379         patch->f[Results::normal][i].z +
02380         patch->f[Results::nbond][i].z +
02381         patch->f[Results::slow][i].z);
02382     CkPrintf("%d %g %g %g\n", a[i].id,
02383         patch->f[Results::normal][i].x,
02384         patch->f[Results::nbond][i].x,
02385         patch->f[Results::slow][i].x);
02386     CkPrintf("%d %g %g %g\n", a[i].id,
02387         patch->f[Results::normal][i].y,
02388         patch->f[Results::nbond][i].y,
02389         patch->f[Results::slow][i].y);
02390     CkPrintf("%d %g %g %g\n", a[i].id,
02391         patch->f[Results::normal][i].z,
02392         patch->f[Results::nbond][i].z,
02393         patch->f[Results::slow][i].z);
02394   }
02395 #endif
02396 
02397 #if PRINT_FORCES
02398   int numAtoms = patch->numAtoms;
02399   FullAtom *a = patch->atom.begin();
02400   for ( int i=0; i<numAtoms; ++i ) {
02401     float fxNo = patch->f[Results::normal][i].x;
02402     float fxNb = patch->f[Results::nbond][i].x;
02403     float fxSl = patch->f[Results::slow][i].x;
02404     float fyNo = patch->f[Results::normal][i].y;
02405     float fyNb = patch->f[Results::nbond][i].y;
02406     float fySl = patch->f[Results::slow][i].y;
02407     float fzNo = patch->f[Results::normal][i].z;
02408     float fzNb = patch->f[Results::nbond][i].z;
02409     float fzSl = patch->f[Results::slow][i].z;
02410     float fx = fxNo+fxNb+fxSl;
02411     float fy = fyNo+fyNb+fySl;
02412     float fz = fzNo+fzNb+fzSl;
02413 
02414                 float f = sqrt(fx*fx+fy*fy+fz*fz);
02415     int id = patch->pExt[i].id;
02416     int seq = patch->flags.sequence;
02417     float x = patch->p[i].position.x;
02418     float y = patch->p[i].position.y;
02419     float z = patch->p[i].position.z;
02420     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
02421     CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
02422     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
02423 //fxNo,fyNo,fzNo,
02424 //fxNb,fyNb,fzNb,
02425 //fxSl,fySl,fzSl,
02426 fx,fy,fz
02427 );
02428         }
02429 #endif
02430 }

void Sequencer::saveForce ( const int  ftag = Results::normal  )  [protected]

Definition at line 1567 of file Sequencer.C.

References patch, and HomePatch::saveForce().

Referenced by integrate(), and minimize().

01568 {
01569   patch->saveForce(ftag);
01570 }

void Sequencer::scalePositionsVelocities ( const Tensor posScale,
const Tensor velScale 
) [protected]

Definition at line 806 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), SimParameters::fixedAtomsOn, j, Patch::lattice, FullAtom::mass, NAMD_bug(), Patch::numAtoms, Lattice::origin(), patch, CompAtom::position, simParams, SimParameters::useGroupPressure, and FullAtom::velocity.

Referenced by multigratorPressure().

00806                                                                                        {
00807   FullAtom *a = patch->atom.begin();
00808   int numAtoms = patch->numAtoms;
00809   Position origin = patch->lattice.origin();
00810   if ( simParams->fixedAtomsOn ) {
00811     NAMD_bug("Sequencer::scalePositionsVelocities, fixed atoms not implemented");
00812   }
00813   if ( simParams->useGroupPressure ) {
00814     int hgs;
00815     for ( int i = 0; i < numAtoms; i += hgs ) {
00816       hgs = a[i].hydrogenGroupSize;
00817       Position pos_cm(0.0, 0.0, 0.0);
00818       Velocity vel_cm(0.0, 0.0, 0.0);
00819       BigReal m_cm = 0.0;
00820       for (int j=0;j < hgs;++j) {
00821         m_cm += a[i+j].mass;
00822         pos_cm += a[i+j].mass*a[i+j].position;
00823         vel_cm += a[i+j].mass*a[i+j].velocity;
00824       }
00825       pos_cm /= m_cm;
00826       vel_cm /= m_cm;
00827       pos_cm -= origin;
00828       Position dpos = posScale*pos_cm;
00829       Velocity dvel = velScale*vel_cm;
00830       for (int j=0;j < hgs;++j) {
00831         a[i+j].position += dpos;
00832         a[i+j].velocity += dvel;
00833       }
00834     }
00835   } else {
00836     for ( int i = 0; i < numAtoms; i++) {
00837       a[i].position += posScale*(a[i].position-origin);
00838       a[i].velocity = velScale*a[i].velocity;
00839     }
00840   }
00841 }

void Sequencer::scaleVelocities ( const BigReal  velScale  )  [protected]

Definition at line 978 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), Patch::numAtoms, and patch.

Referenced by multigratorTemperature().

00978                                                       {
00979   FullAtom *a = patch->atom.begin();
00980   int numAtoms = patch->numAtoms;
00981   for ( int i = 0; i < numAtoms; i++) {
00982     a[i].velocity *= velScale;
00983   }
00984 }

void Sequencer::submitCollections ( int  step,
int  zeroVel = 0 
) [protected]

Definition at line 2245 of file Sequencer.C.

References HomePatch::atom, collection, Output::coordinateNeeded(), Patch::f, Patch::flags, Output::forceNeeded(), Patch::lattice, Flags::maxForceUsed, patch, Results::slow, CollectionMgr::submitForces(), CollectionMgr::submitPositions(), CollectionMgr::submitVelocities(), and Output::velocityNeeded().

Referenced by algorithm(), integrate(), and minimize().

02246 {
02247   int prec = Output::coordinateNeeded(step);
02248   if ( prec )
02249     collection->submitPositions(step,patch->atom,patch->lattice,prec);
02250   if ( Output::velocityNeeded(step) )
02251     collection->submitVelocities(step,zeroVel,patch->atom);
02252   if ( Output::forceNeeded(step) ) {
02253     int maxForceUsed = patch->flags.maxForceUsed;
02254     if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
02255     collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
02256   }
02257 }

void Sequencer::submitHalfstep ( int   )  [protected]

Definition at line 1695 of file Sequencer.C.

References ADD_TENSOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), Lattice::c(), doKineticEnergy, Flags::doVirial, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, Patch::lattice, GlobalMaster::lattice, Vector::length2(), FullAtom::mass, SimParameters::multigratorOn, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, CompAtom::partition, partition(), patch, pressureProfileReduction, SimParameters::pressureProfileSlabs, reduction, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_MOMENTUM_SQUARED, REDUCTION_VIRIAL_NORMAL, simParams, SimParameters::useGroupPressure, FullAtom::velocity, Vector::x, Vector::y, z, and Vector::z.

Referenced by integrate().

01696 {
01697   // velocity-dependent quantities *** ONLY ***
01698   // positions are not at half-step when called
01699   FullAtom *a = patch->atom.begin();
01700   int numAtoms = patch->numAtoms;
01701 
01702 #if CMK_BLUEGENEL
01703   CmiNetworkProgressAfter (0);
01704 #endif
01705 
01706   // For non-Multigrator doKineticEnergy = 1 always
01707   Tensor momentumSqrSum;
01708   if (doKineticEnergy || patch->flags.doVirial)
01709   {
01710     BigReal kineticEnergy = 0;
01711     Tensor virial;
01712     if ( simParams->pairInteractionOn ) {
01713       if ( simParams->pairInteractionSelf ) {
01714         for ( int i = 0; i < numAtoms; ++i ) {
01715           if ( a[i].partition != 1 ) continue;
01716           kineticEnergy += a[i].mass * a[i].velocity.length2();
01717           virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01718         }
01719       }
01720     } else {
01721       for ( int i = 0; i < numAtoms; ++i ) {
01722         if (a[i].mass < 0.01) continue;
01723         kineticEnergy += a[i].mass * a[i].velocity.length2();
01724         virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01725       }
01726     }
01727     
01728     if (simParams->multigratorOn && !simParams->useGroupPressure) {
01729       momentumSqrSum = virial;
01730     }
01731     kineticEnergy *= 0.5 * 0.5;
01732     reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
01733     virial *= 0.5;
01734     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01735 #ifdef ALTVIRIAL
01736     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
01737 #endif
01738   }
01739  
01740   if (pressureProfileReduction) {
01741     int nslabs = simParams->pressureProfileSlabs;
01742     const Lattice &lattice = patch->lattice;
01743     BigReal idz = nslabs/lattice.c().z;
01744     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
01745     int useGroupPressure = simParams->useGroupPressure;
01746 
01747     // Compute kinetic energy partition, possibly subtracting off
01748     // internal kinetic energy if group pressure is enabled.
01749     // Since the regular pressure is 1/2 mvv and the internal kinetic
01750     // term that is subtracted off for the group pressure is
01751     // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
01752     // 1/2 m * v * v_cm.  The factor of 1/2 is because submitHalfstep
01753     // gets called twice per timestep.
01754     int hgs;
01755     for (int i=0; i<numAtoms; i += hgs) {
01756       int j, ppoffset;
01757       hgs = a[i].hydrogenGroupSize;
01758       int partition = a[i].partition;
01759 
01760       BigReal m_cm = 0;
01761       Velocity v_cm(0,0,0);
01762       for (j=i; j< i+hgs; ++j) {
01763         m_cm += a[j].mass;
01764         v_cm += a[j].mass * a[j].velocity;
01765       }
01766       v_cm /= m_cm;
01767       for (j=i; j < i+hgs; ++j) {
01768         BigReal mass = a[j].mass;
01769         if (! (useGroupPressure && j != i)) {
01770           BigReal z = a[j].position.z;
01771           int slab = (int)floor((z-zmin)*idz);
01772           if (slab < 0) slab += nslabs;
01773           else if (slab >= nslabs) slab -= nslabs;
01774           ppoffset = 3*(slab + partition*nslabs);
01775         }
01776         BigReal wxx, wyy, wzz;
01777         if (useGroupPressure) {
01778           wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
01779           wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
01780           wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
01781         } else {
01782           wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
01783           wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
01784           wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
01785         }
01786         pressureProfileReduction->item(ppoffset  ) += wxx;
01787         pressureProfileReduction->item(ppoffset+1) += wyy;
01788         pressureProfileReduction->item(ppoffset+2) += wzz;
01789       }
01790     }
01791   } 
01792 
01793   // For non-Multigrator doKineticEnergy = 1 always
01794   if (doKineticEnergy || patch->flags.doVirial)
01795   {
01796     BigReal intKineticEnergy = 0;
01797     Tensor intVirialNormal;
01798 
01799     int hgs;
01800     for ( int i = 0; i < numAtoms; i += hgs ) {
01801 
01802 #if CMK_BLUEGENEL
01803       CmiNetworkProgress ();
01804 #endif
01805 
01806       hgs = a[i].hydrogenGroupSize;
01807       int j;
01808       BigReal m_cm = 0;
01809       Velocity v_cm(0,0,0);
01810       for ( j = i; j < (i+hgs); ++j ) {
01811         m_cm += a[j].mass;
01812         v_cm += a[j].mass * a[j].velocity;
01813       }
01814       if (simParams->multigratorOn && simParams->useGroupPressure) {
01815         momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
01816       }
01817       v_cm /= m_cm;
01818       if ( simParams->pairInteractionOn ) {
01819         if ( simParams->pairInteractionSelf ) {
01820           for ( j = i; j < (i+hgs); ++j ) {
01821             if ( a[j].partition != 1 ) continue;
01822             BigReal mass = a[j].mass;
01823             Vector v = a[j].velocity;
01824             Vector dv = v - v_cm;
01825             intKineticEnergy += mass * (v * dv);
01826             intVirialNormal.outerAdd (mass, v, dv);
01827           }
01828         }
01829       } else {
01830         for ( j = i; j < (i+hgs); ++j ) {
01831           BigReal mass = a[j].mass;
01832           Vector v = a[j].velocity;
01833           Vector dv = v - v_cm;
01834           intKineticEnergy += mass * (v * dv);
01835           intVirialNormal.outerAdd(mass, v, dv);
01836         }
01837       }
01838     }
01839 
01840     intKineticEnergy *= 0.5 * 0.5;
01841     reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
01842     intVirialNormal *= 0.5;
01843     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
01844     if ( simParams->multigratorOn) {
01845       momentumSqrSum *= 0.5;
01846       ADD_TENSOR_OBJECT(reduction,REDUCTION_MOMENTUM_SQUARED,momentumSqrSum);
01847     }
01848   }
01849 
01850 }

void Sequencer::submitMinimizeReductions ( int  ,
BigReal  fmax2 
) [protected]

Definition at line 2107 of file Sequencer.C.

References HomePatch::atom, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), SimParameters::drudeBondLen, SimParameters::drudeHardWallOn, f, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, if(), SubmitReduction::item(), Vector::length2(), HomePatch::minimize_rattle2(), Results::nbond, Results::normal, Patch::numAtoms, patch, Patch::pExt, CompAtom::position, SimParameters::printBadContacts, reduction, REDUCTION_ATOM_CHECKSUM, Flags::sequence, simParams, Results::slow, TIMEFACTOR, and FullAtom::velocity.

Referenced by minimize().

02108 {
02109   FullAtom *a = patch->atom.begin();
02110   Force *f1 = patch->f[Results::normal].begin();
02111   Force *f2 = patch->f[Results::nbond].begin();
02112   Force *f3 = patch->f[Results::slow].begin();
02113   const bool fixedAtomsOn = simParams->fixedAtomsOn;
02114   const bool drudeHardWallOn = simParams->drudeHardWallOn;
02115   const double drudeBondLen = simParams->drudeBondLen;
02116   const double drudeBondLen2 = drudeBondLen * drudeBondLen;
02117   const double drudeStep = 0.1/(TIMEFACTOR*TIMEFACTOR);
02118   const double drudeMove = 0.01;
02119   const double drudeStep2 = drudeStep * drudeStep;
02120   const double drudeMove2 = drudeMove * drudeMove;
02121   int numAtoms = patch->numAtoms;
02122 
02123   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
02124 
02125   for ( int i = 0; i < numAtoms; ++i ) {
02126     f1[i] += f2[i] + f3[i];  // add all forces
02127     if ( drudeHardWallOn && i && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) { // drude particle
02128       if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
02129         if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
02130           a[i].position += drudeMove * f1[i].unit();
02131         } else {
02132           a[i].position += drudeStep * f1[i];
02133         }
02134         if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
02135           a[i].position = a[i-1].position + drudeBondLen * (a[i].position - a[i-1].position).unit();
02136         }
02137       }
02138       Vector netf = f1[i-1] + f1[i];
02139       if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
02140       f1[i-1] = netf;
02141       f1[i] = 0.;
02142     }
02143     if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
02144   }
02145 
02146   f2 = f3 = 0;  // included in f1
02147 
02148   BigReal maxv2 = 0.;
02149 
02150   for ( int i = 0; i < numAtoms; ++i ) {
02151     BigReal v2 = a[i].velocity.length2();
02152     if ( v2 > 0. ) {
02153       if ( v2 > maxv2 ) maxv2 = v2;
02154     } else {
02155       v2 = f1[i].length2();
02156       if ( v2 > maxv2 ) maxv2 = v2;
02157     }
02158   }
02159 
02160   if ( fmax2 > 10. * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR )
02161   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial, true /* forces */); }
02162 
02163   BigReal fdotf = 0;
02164   BigReal fdotv = 0;
02165   BigReal vdotv = 0;
02166   int numHuge = 0;
02167   for ( int i = 0; i < numAtoms; ++i ) {
02168     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
02169     if ( drudeHardWallOn && (a[i].mass > 0.01) && ((a[i].mass < 1.0)) ) continue; // drude particle
02170     Force f = f1[i];
02171     BigReal ff = f * f;
02172     if ( ff > fmax2 ) {
02173       if (simParams->printBadContacts) {
02174         CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
02175       }
02176       ++numHuge;
02177       // pad scaling so minimizeMoveDownhill() doesn't miss them
02178       BigReal fmult = 1.01 * sqrt(fmax2/ff);
02179       f *= fmult;  ff = f * f;
02180       f1[i] *= fmult;
02181     }
02182     fdotf += ff;
02183     fdotv += f * a[i].velocity;
02184     vdotv += a[i].velocity * a[i].velocity;
02185   }
02186 
02187   reduction->item(REDUCTION_MIN_F_DOT_F) += fdotf;
02188   reduction->item(REDUCTION_MIN_F_DOT_V) += fdotv;
02189   reduction->item(REDUCTION_MIN_V_DOT_V) += vdotv;
02190   reduction->item(REDUCTION_MIN_HUGE_COUNT) += numHuge;
02191 
02192   {
02193     Tensor intVirialNormal;
02194     Tensor intVirialNbond;
02195     Tensor intVirialSlow;
02196 
02197     int hgs;
02198     for ( int i = 0; i < numAtoms; i += hgs ) {
02199       hgs = a[i].hydrogenGroupSize;
02200       int j;
02201       BigReal m_cm = 0;
02202       Position x_cm(0,0,0);
02203       for ( j = i; j < (i+hgs); ++j ) {
02204         m_cm += a[j].mass;
02205         x_cm += a[j].mass * a[j].position;
02206       }
02207       x_cm /= m_cm;
02208       for ( j = i; j < (i+hgs); ++j ) {
02209         BigReal mass = a[j].mass;
02210         // net force treated as zero for fixed atoms
02211         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
02212         Vector dx = a[j].position - x_cm;
02213         intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02214         intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02215         intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02216       }
02217     }
02218 
02219     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02220     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
02221     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
02222   }
02223 
02224   if ( simParams->fixedAtomsOn ) {
02225     Tensor fixVirialNormal;
02226     Tensor fixVirialNbond;
02227     Tensor fixVirialSlow;
02228     Vector fixForceNormal = 0;
02229     Vector fixForceNbond = 0;
02230     Vector fixForceSlow = 0;
02231 
02232     calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
02233 
02234     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
02235     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
02236     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
02237     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
02238     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
02239     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
02240   }
02241 
02242   reduction->submit();
02243 }

void Sequencer::submitMomentum ( int  step  )  [protected]

Definition at line 756 of file Sequencer.C.

References ADD_VECTOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), SubmitReduction::item(), FullAtom::mass, Patch::numAtoms, patch, reduction, REDUCTION_HALFSTEP_MOMENTUM, REDUCTION_MOMENTUM_MASS, simParams, and SimParameters::zeroMomentumAlt.

Referenced by integrate().

00756                                        {
00757 
00758   FullAtom *a = patch->atom.begin();
00759   const int numAtoms = patch->numAtoms;
00760 
00761   Vector momentum = 0;
00762   BigReal mass = 0;
00763 if ( simParams->zeroMomentumAlt ) {
00764   for ( int i = 0; i < numAtoms; ++i ) {
00765     momentum += a[i].mass * a[i].velocity;
00766     mass += 1.;
00767   }
00768 } else {
00769   for ( int i = 0; i < numAtoms; ++i ) {
00770     momentum += a[i].mass * a[i].velocity;
00771     mass += a[i].mass;
00772   }
00773 }
00774 
00775   ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
00776   reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
00777 }

void Sequencer::submitReductions ( int   )  [protected]

Definition at line 1872 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, HomePatch::atom, ResizeArray< Elem >::begin(), cross(), doKineticEnergy, doMomenta, Flags::doVirial, SimParameters::drudeOn, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, Patch::lattice, Vector::length2(), HomePatch::marginViolations, FullAtom::mass, Results::nbond, Results::normal, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, partition(), patch, CompAtom::position, reduction, REDUCTION_ANGULAR_MOMENTUM, REDUCTION_ATOM_CHECKSUM, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_INT_VIRIAL_NBOND, REDUCTION_INT_VIRIAL_NORMAL, REDUCTION_INT_VIRIAL_SLOW, REDUCTION_MARGIN_VIOLATIONS, REDUCTION_MOMENTUM, simParams, Results::slow, and FullAtom::velocity.

Referenced by integrate().

01873 {
01874   FullAtom *a = patch->atom.begin();
01875   int numAtoms = patch->numAtoms;
01876 
01877 #if CMK_BLUEGENEL
01878   CmiNetworkProgressAfter(0);
01879 #endif
01880 
01881   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01882   reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
01883 
01884   // For non-Multigrator doKineticEnergy = 1 always
01885   if (doKineticEnergy || doMomenta || patch->flags.doVirial)
01886   {
01887     BigReal kineticEnergy = 0;
01888     Vector momentum = 0;
01889     Vector angularMomentum = 0;
01890     Vector o = patch->lattice.origin();
01891     int i;
01892     if ( simParams->pairInteractionOn ) {
01893       if ( simParams->pairInteractionSelf ) {
01894         for (i = 0; i < numAtoms; ++i ) {
01895           if ( a[i].partition != 1 ) continue;
01896           kineticEnergy += a[i].mass * a[i].velocity.length2();
01897           momentum += a[i].mass * a[i].velocity;
01898           angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01899         }
01900       }
01901     } else {
01902       for (i = 0; i < numAtoms; ++i ) {
01903         kineticEnergy += a[i].mass * a[i].velocity.length2();
01904         momentum += a[i].mass * a[i].velocity;
01905         angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
01906       }
01907       if (simParams->drudeOn) {
01908         BigReal drudeComKE = 0.;
01909         BigReal drudeBondKE = 0.;
01910 
01911         for (i = 0;  i < numAtoms;  i++) {
01912           if (i < numAtoms-1 &&
01913               a[i+1].mass < 1.0 && a[i+1].mass >= 0.001) {
01914             // i+1 is a Drude particle with parent i
01915 
01916             // convert from Cartesian coordinates to (COM,bond) coordinates
01917             BigReal m_com = (a[i].mass + a[i+1].mass);  // mass of COM
01918             BigReal m = a[i+1].mass / m_com;  // mass ratio
01919             BigReal m_bond = a[i+1].mass * (1. - m);  // mass of bond
01920             Vector v_bond = a[i+1].velocity - a[i].velocity;  // vel of bond
01921             Vector v_com = a[i].velocity + m * v_bond;  // vel of COM
01922 
01923             drudeComKE += m_com * v_com.length2();
01924             drudeBondKE += m_bond * v_bond.length2();
01925 
01926             i++;  // +1 from loop, we've updated both particles
01927           }
01928           else {
01929             drudeComKE += a[i].mass * a[i].velocity.length2();
01930           }
01931         } // end for
01932 
01933         drudeComKE *= 0.5;
01934         drudeBondKE *= 0.5;
01935         reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY)
01936           += drudeComKE;
01937         reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY)
01938           += drudeBondKE;
01939       } // end drudeOn
01940 
01941     } // end else
01942 
01943     kineticEnergy *= 0.5;
01944     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
01945     ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
01946     ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);  
01947   }
01948 
01949 #ifdef ALTVIRIAL
01950   // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
01951   {
01952     Tensor altVirial;
01953     for ( int i = 0; i < numAtoms; ++i ) {
01954       altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
01955     }
01956     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
01957   }
01958   {
01959     Tensor altVirial;
01960     for ( int i = 0; i < numAtoms; ++i ) {
01961       altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
01962     }
01963     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
01964   }
01965   {
01966     Tensor altVirial;
01967     for ( int i = 0; i < numAtoms; ++i ) {
01968       altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
01969     }
01970     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
01971   }
01972 #endif
01973 
01974   // For non-Multigrator doKineticEnergy = 1 always
01975   if (doKineticEnergy || patch->flags.doVirial)
01976   {
01977     BigReal intKineticEnergy = 0;
01978     Tensor intVirialNormal;
01979     Tensor intVirialNbond;
01980     Tensor intVirialSlow;
01981 
01982     int hgs;
01983     for ( int i = 0; i < numAtoms; i += hgs ) {
01984 #if CMK_BLUEGENEL
01985       CmiNetworkProgress();
01986 #endif
01987       hgs = a[i].hydrogenGroupSize;
01988       int j;
01989       BigReal m_cm = 0;
01990       Position x_cm(0,0,0);
01991       Velocity v_cm(0,0,0);
01992       for ( j = i; j < (i+hgs); ++j ) {
01993         m_cm += a[j].mass;
01994         x_cm += a[j].mass * a[j].position;
01995         v_cm += a[j].mass * a[j].velocity;
01996       }
01997       x_cm /= m_cm;
01998       v_cm /= m_cm;
01999       int fixedAtomsOn = simParams->fixedAtomsOn;
02000       if ( simParams->pairInteractionOn ) {
02001         int pairInteractionSelf = simParams->pairInteractionSelf;
02002         for ( j = i; j < (i+hgs); ++j ) {
02003           if ( a[j].partition != 1 &&
02004                ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
02005           // net force treated as zero for fixed atoms
02006           if ( fixedAtomsOn && a[j].atomFixed ) continue;
02007           BigReal mass = a[j].mass;
02008           Vector v = a[j].velocity;
02009           Vector dv = v - v_cm;
02010           intKineticEnergy += mass * (v * dv);
02011           Vector dx = a[j].position - x_cm;
02012           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02013           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02014           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02015         }
02016       } else {
02017         for ( j = i; j < (i+hgs); ++j ) {
02018           // net force treated as zero for fixed atoms
02019           if ( fixedAtomsOn && a[j].atomFixed ) continue;
02020           BigReal mass = a[j].mass;
02021           Vector v = a[j].velocity;
02022           Vector dv = v - v_cm;
02023           intKineticEnergy += mass * (v * dv);
02024           Vector dx = a[j].position - x_cm;
02025           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02026           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02027           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02028         }
02029       }
02030     }
02031 
02032     intKineticEnergy *= 0.5;
02033     reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
02034     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02035     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
02036     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
02037   }
02038 
02039   if (pressureProfileReduction && simParams->useGroupPressure) {
02040     // subtract off internal virial term, calculated as for intVirial.
02041     int nslabs = simParams->pressureProfileSlabs;
02042     const Lattice &lattice = patch->lattice;
02043     BigReal idz = nslabs/lattice.c().z;
02044     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
02045     int useGroupPressure = simParams->useGroupPressure;
02046 
02047     int hgs;
02048     for (int i=0; i<numAtoms; i += hgs) {
02049       int j;
02050       hgs = a[i].hydrogenGroupSize;
02051       BigReal m_cm = 0;
02052       Position x_cm(0,0,0);
02053       for (j=i; j< i+hgs; ++j) {
02054         m_cm += a[j].mass;
02055         x_cm += a[j].mass * a[j].position;
02056       }
02057       x_cm /= m_cm;
02058       
02059       BigReal z = a[i].position.z;
02060       int slab = (int)floor((z-zmin)*idz);
02061       if (slab < 0) slab += nslabs;
02062       else if (slab >= nslabs) slab -= nslabs;
02063       int partition = a[i].partition;
02064       int ppoffset = 3*(slab + nslabs*partition);
02065       for (j=i; j < i+hgs; ++j) {
02066         BigReal mass = a[j].mass;
02067         Vector dx = a[j].position - x_cm;
02068         const Vector &fnormal = patch->f[Results::normal][j];
02069         const Vector &fnbond  = patch->f[Results::nbond][j];
02070         const Vector &fslow   = patch->f[Results::slow][j];
02071         BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
02072         BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
02073         BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
02074         pressureProfileReduction->item(ppoffset  ) -= wxx;
02075         pressureProfileReduction->item(ppoffset+1) -= wyy;
02076         pressureProfileReduction->item(ppoffset+2) -= wzz;
02077       }
02078     }
02079   }
02080 
02081   // For non-Multigrator doVirial = 1 always
02082   if (patch->flags.doVirial)
02083   {
02084     if ( simParams->fixedAtomsOn ) {
02085       Tensor fixVirialNormal;
02086       Tensor fixVirialNbond;
02087       Tensor fixVirialSlow;
02088       Vector fixForceNormal = 0;
02089       Vector fixForceNbond = 0;
02090       Vector fixForceSlow = 0;
02091 
02092       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
02093 
02094       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
02095       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
02096       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
02097       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
02098       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
02099       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
02100     }
02101   }
02102 
02103   reduction->submit();
02104   if (pressureProfileReduction) pressureProfileReduction->submit();
02105 }

void Sequencer::suspend ( void   ) 

Definition at line 107 of file Sequencer.C.

References HomePatch::ldObjHandle, LdbCoordinator::Object(), patch, LdbCoordinator::pauseWork(), and LdbCoordinator::startWork().

Referenced by HomePatch::doAtomMigration(), and LdbCoordinator::rebalance().

00108 {
00109     LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
00110     CthSuspend();
00111     LdbCoordinator::Object()->startWork(patch->ldObjHandle);
00112 }

void Sequencer::tcoupleVelocities ( BigReal  ,
int   
) [protected]

Definition at line 1549 of file Sequencer.C.

References HomePatch::atom, ResizeArray< Elem >::begin(), broadcast, colvarproxy_namd::dt(), Node::molecule, Patch::numAtoms, Node::Object(), patch, simParams, ControllerBroadcasts::tcoupleCoefficient, SimParameters::tCoupleOn, and FullAtom::velocity.

Referenced by integrate().

01550 {
01551   if ( simParams->tCoupleOn )
01552   {
01553     FullAtom *a = patch->atom.begin();
01554     int numAtoms = patch->numAtoms;
01555     BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
01556     Molecule *molecule = Node::Object()->molecule;
01557     BigReal dt = dt_fs * 0.001;  // convert to ps
01558     coefficient *= dt;
01559     for ( int i = 0; i < numAtoms; ++i )
01560     {
01561       BigReal f1 = exp( coefficient * a[i].langevinParam );
01562       a[i].velocity *= f1;
01563     }
01564   }
01565 }

void Sequencer::terminate ( void   )  [protected]

Definition at line 2460 of file Sequencer.C.

References HomePatch::ldObjHandle, LdbCoordinator::Object(), patch, and LdbCoordinator::pauseWork().

Referenced by algorithm(), hardWallDrude(), maximumMove(), and rattle1().

02460                           {
02461   LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
02462   CthFree(thread);
02463   CthSuspend();
02464 }

void Sequencer::traceBarrier ( int   )  [protected]

Definition at line 2450 of file Sequencer.C.

References broadcast, SimpleBroadcastObject< T >::get(), and ControllerBroadcasts::traceBarrier.

Referenced by integrate().

02450                                     {
02451         broadcast->traceBarrier.get(step);
02452 }


Friends And Related Function Documentation

friend class HomePatch [friend]

Definition at line 24 of file Sequencer.h.


Member Data Documentation

BigReal Sequencer::adaptTempT [protected]

Definition at line 81 of file Sequencer.h.

Referenced by adaptTempUpdate(), integrate(), langevinVelocities(), and langevinVelocitiesBBK2().

int Sequencer::berendsenPressure_count [protected]

Definition at line 92 of file Sequencer.h.

Referenced by algorithm(), berendsenPressure(), HomePatch::recvCheckpointLoad(), and Sequencer().

ControllerBroadcasts* Sequencer::broadcast [protected]

Definition at line 126 of file Sequencer.h.

Referenced by adaptTempUpdate(), algorithm(), berendsenPressure(), correctMomentum(), cycleBarrier(), langevinPiston(), minimize(), multigratorPressure(), multigratorTemperature(), rescaleaccelMD(), rescaleVelocities(), Sequencer(), tcoupleVelocities(), traceBarrier(), and ~Sequencer().

int Sequencer::checkpoint_berendsenPressure_count [protected]

Definition at line 93 of file Sequencer.h.

Referenced by algorithm().

CollectionMgr* const Sequencer::collection [protected]

Definition at line 125 of file Sequencer.h.

Referenced by submitCollections().

int Sequencer::doKineticEnergy [protected]

Definition at line 108 of file Sequencer.h.

Referenced by integrate(), submitHalfstep(), and submitReductions().

int Sequencer::doMomenta [protected]

Definition at line 109 of file Sequencer.h.

Referenced by integrate(), and submitReductions().

int Sequencer::ldbSteps [protected]

Definition at line 128 of file Sequencer.h.

Referenced by rebalanceLoad().

SubmitReduction* Sequencer::min_reduction [protected]

Definition at line 39 of file Sequencer.h.

Referenced by newMinimizeDirection(), Sequencer(), and ~Sequencer().

SubmitReduction* Sequencer::multigratorReduction [protected]

Definition at line 107 of file Sequencer.h.

Referenced by multigratorTemperature(), Sequencer(), and ~Sequencer().

int Sequencer::pairlistsAge [protected]

Definition at line 43 of file Sequencer.h.

Referenced by runComputeObjects().

int Sequencer::pairlistsAreValid [protected]

Definition at line 42 of file Sequencer.h.

Referenced by algorithm(), rebalanceLoad(), and runComputeObjects().

HomePatch* const Sequencer::patch [protected]

Definition at line 121 of file Sequencer.h.

Referenced by addForceToMomentum(), addForceToMomentum3(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), algorithm(), berendsenPressure(), calcFixVirial(), calcKineticEnergy(), correctMomentum(), hardWallDrude(), integrate(), langevinPiston(), langevinVelocities(), langevinVelocitiesBBK1(), langevinVelocitiesBBK2(), maximumMove(), minimizationQuenchVelocity(), minimize(), minimizeMoveDownhill(), multigratorPressure(), multigratorTemperature(), newMinimizeDirection(), newMinimizePosition(), quenchVelocities(), rattle1(), reassignVelocities(), rebalanceLoad(), reinitVelocities(), reloadCharges(), rescaleaccelMD(), rescaleVelocities(), rescaleVelocitiesByFactor(), run(), runComputeObjects(), saveForce(), scalePositionsVelocities(), scaleVelocities(), Sequencer(), submitCollections(), submitHalfstep(), submitMinimizeReductions(), submitMomentum(), submitReductions(), suspend(), tcoupleVelocities(), and terminate().

SubmitReduction* Sequencer::pressureProfileReduction [protected]

Definition at line 123 of file Sequencer.h.

Referenced by hardWallDrude(), rattle1(), Sequencer(), submitHalfstep(), and ~Sequencer().

Random* Sequencer::random [protected]

Definition at line 119 of file Sequencer.h.

Referenced by langevinVelocities(), langevinVelocitiesBBK2(), reassignVelocities(), reinitVelocities(), Sequencer(), and ~Sequencer().

SubmitReduction* Sequencer::reduction [protected]

Definition at line 122 of file Sequencer.h.

Referenced by hardWallDrude(), multigratorPressure(), rattle1(), Sequencer(), submitHalfstep(), submitMinimizeReductions(), submitMomentum(), submitReductions(), and ~Sequencer().

int Sequencer::rescaleVelocities_numTemps [protected]

Definition at line 86 of file Sequencer.h.

Referenced by rescaleVelocities(), and Sequencer().

SimParameters* const Sequencer::simParams [protected]

Definition at line 120 of file Sequencer.h.

Referenced by adaptTempUpdate(), addMovDragToPosition(), addRotDragToPosition(), algorithm(), berendsenPressure(), calcFixVirial(), calcKineticEnergy(), correctMomentum(), hardWallDrude(), integrate(), langevinPiston(), langevinVelocities(), langevinVelocitiesBBK1(), langevinVelocitiesBBK2(), maximumMove(), minimizationQuenchVelocity(), minimize(), minimizeMoveDownhill(), multigratorPressure(), multigratorTemperature(), newMinimizeDirection(), newMinimizePosition(), rattle1(), reassignVelocities(), reinitVelocities(), rescaleaccelMD(), rescaleVelocities(), runComputeObjects(), scalePositionsVelocities(), Sequencer(), submitHalfstep(), submitMinimizeReductions(), submitMomentum(), submitReductions(), and tcoupleVelocities().

int Sequencer::slowFreq [protected]

Definition at line 95 of file Sequencer.h.

Referenced by integrate(), and langevinPiston().


The documentation for this class was generated from the following files:
Generated on Wed Sep 20 01:17:20 2017 for NAMD by  doxygen 1.4.7