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 rescaleSoluteCharges (BigReal)
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 stochRescaleVelocities (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 stochRescale_count
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 66 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, rescaleSoluteCharges(), rescaleVelocities_numTemps, simParams, SimParameters::soluteScalingFactorCharge, SimParameters::soluteScalingOn, Random::split(), stochRescale_count, and ReductionMgr::willSubmit().

00066                                  :
00067         simParams(Node::Object()->simParameters),
00068         patch(p),
00069         collection(CollectionMgr::Object()),
00070         ldbSteps(0)
00071 {
00072     broadcast = new ControllerBroadcasts(& patch->ldObjHandle);
00073     reduction = ReductionMgr::Object()->willSubmit(
00074                   simParams->accelMDOn ? REDUCTIONS_AMD : REDUCTIONS_BASIC );
00075     min_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MINIMIZER,1);
00076     if (simParams->pressureProfileOn) {
00077       int ntypes = simParams->pressureProfileAtomTypes;
00078       int nslabs = simParams->pressureProfileSlabs;
00079       pressureProfileReduction = 
00080         ReductionMgr::Object()->willSubmit(
00081                 REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
00082     } else {
00083       pressureProfileReduction = NULL;
00084     }
00085     if (simParams->multigratorOn) {
00086       multigratorReduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
00087     } else {
00088       multigratorReduction = NULL;
00089     }
00090     ldbCoordinator = (LdbCoordinator::Object());
00091     random = new Random(simParams->randomSeed);
00092     random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
00093 
00094     // Is soluteScaling enabled?
00095     if (simParams->soluteScalingOn) {
00096       // If so, we must "manually" perform charge scaling on startup because
00097       // Sequencer will not get a scripting task for initial charge scaling.
00098       // Subsequent rescalings will take place through a scripting task.
00099       rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
00100     }
00101 
00102     rescaleVelocities_numTemps = 0;
00103     stochRescale_count = 0;
00104     berendsenPressure_count = 0;
00105 //    patch->write_tip4_props();
00106 }

Sequencer::~Sequencer ( void   )  [virtual]

Definition at line 108 of file Sequencer.C.

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

00109 {
00110     delete broadcast;
00111     delete reduction;
00112     delete min_reduction;
00113     if (pressureProfileReduction) delete pressureProfileReduction;
00114     delete random;
00115     if (multigratorReduction) delete multigratorReduction;
00116 }


Member Function Documentation

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

Definition at line 1740 of file Sequencer.C.

References ControllerBroadcasts::adaptTemperature, SimParameters::adaptTempFreq, SimParameters::adaptTempLastStep, SimParameters::adaptTempOn, adaptTempT, broadcast, SimParameters::firstTimestep, SimpleBroadcastObject< T >::get(), SimParameters::langevinOn, SimParameters::langevinTemp, and simParams.

Referenced by integrate().

01741 {
01742    //check if adaptive tempering is enabled and in the right timestep range
01743    if (!simParams->adaptTempOn) return;
01744    if ( (step < simParams->adaptTempFirstStep ) || 
01745      ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) {
01746         if (simParams->langevinOn) // restore langevin temperature
01747             adaptTempT = simParams->langevinTemp;
01748         return;
01749    }
01750    // Get Updated Temperature
01751    if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
01752      // Blocking receive for the updated adaptive tempering temperature.
01753      adaptTempT = broadcast->adaptTemperature.get(step);
01754 }

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

Definition at line 1904 of file Sequencer.C.

References HomePatch::addForceToMomentum(), ResizeArray< Elem >::begin(), ResizeArray< Elem >::const_begin(), colvarproxy_namd::dt(), Patch::f, Patch::flags, NAMD_EVENT_RANGE_2, Patch::numAtoms, patch, and TIMEFACTOR.

Referenced by newtonianVelocities().

01906       {
01907   NAMD_EVENT_RANGE_2(patch->flags.event_on,
01908       NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
01909 #if CMK_BLUEGENEL
01910   CmiNetworkProgressAfter (0);
01911 #endif
01912   const BigReal dt = timestep / TIMEFACTOR;
01913   FullAtom *atom_arr  = patch->atom.begin();
01914   ForceList *f_use = (useSaved ? patch->f_saved : patch->f);
01915   const Force *force_arr = f_use[ftag].const_begin();
01916   patch->addForceToMomentum(atom_arr, force_arr, dt, patch->numAtoms);
01917 }

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 1919 of file Sequencer.C.

References HomePatch::addForceToMomentum3(), ResizeArray< Elem >::begin(), ResizeArray< Elem >::const_begin(), Patch::f, Patch::flags, NAMD_EVENT_RANGE_2, Patch::numAtoms, patch, and TIMEFACTOR.

Referenced by newtonianVelocities().

01923       {
01924   NAMD_EVENT_RANGE_2(patch->flags.event_on,
01925       NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
01926 #if CMK_BLUEGENEL
01927   CmiNetworkProgressAfter (0);
01928 #endif
01929   const BigReal dt1 = timestep1 / TIMEFACTOR;
01930   const BigReal dt2 = timestep2 / TIMEFACTOR;
01931   const BigReal dt3 = timestep3 / TIMEFACTOR;
01932   ForceList *f_use1 = (useSaved1 ? patch->f_saved : patch->f);
01933   ForceList *f_use2 = (useSaved2 ? patch->f_saved : patch->f);
01934   ForceList *f_use3 = (useSaved3 ? patch->f_saved : patch->f);
01935   FullAtom *atom_arr  = patch->atom.begin();
01936   const Force *force_arr1 = f_use1[ftag1].const_begin();
01937   const Force *force_arr2 = f_use2[ftag2].const_begin();
01938   const Force *force_arr3 = f_use3[ftag3].const_begin();
01939   patch->addForceToMomentum3 (atom_arr, force_arr1, force_arr2, force_arr3,
01940       dt1, dt2, dt3, patch->numAtoms);
01941 }

void Sequencer::addMovDragToPosition ( BigReal  timestep  )  [protected]

Definition at line 718 of file Sequencer.C.

References 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().

00718                                                      {
00719   FullAtom *atom = patch->atom.begin();
00720   int numAtoms = patch->numAtoms;
00721   Molecule *molecule = Node::Object()->molecule;   // need its methods
00722   const BigReal movDragGlobVel = simParams->movDragGlobVel;
00723   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00724   Vector movDragVel, dragIncrement;
00725   for ( int i = 0; i < numAtoms; ++i )
00726   {
00727     // skip if fixed atom or zero drag attribute
00728     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00729          || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
00730     molecule->get_movdrag_params(movDragVel, atom[i].id);
00731     dragIncrement = movDragGlobVel * movDragVel * dt;
00732     atom[i].position += dragIncrement;
00733   }
00734 }

void Sequencer::addRotDragToPosition ( BigReal  timestep  )  [protected]

Definition at line 737 of file Sequencer.C.

References 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().

00737                                                      {
00738   FullAtom *atom = patch->atom.begin();
00739   int numAtoms = patch->numAtoms;
00740   Molecule *molecule = Node::Object()->molecule;   // need its methods
00741   const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
00742   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
00743   BigReal rotDragVel, dAngle;
00744   Vector atomRadius;
00745   Vector rotDragAxis, rotDragPivot, dragIncrement;
00746   for ( int i = 0; i < numAtoms; ++i )
00747   {
00748     // skip if fixed atom or zero drag attribute
00749     if ( (simParams->fixedAtomsOn && atom[i].atomFixed) 
00750          || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
00751     molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
00752     dAngle = rotDragGlobVel * rotDragVel * dt;
00753     rotDragAxis /= rotDragAxis.length();
00754     atomRadius = atom[i].position - rotDragPivot;
00755     dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
00756     atom[i].position += dragIncrement;
00757   }
00758 }

void Sequencer::addVelocityToPosition ( BigReal  timestep  )  [protected]

Definition at line 1943 of file Sequencer.C.

References HomePatch::addVelocityToPosition(), ResizeArray< Elem >::begin(), colvarproxy_namd::dt(), Patch::flags, NAMD_EVENT_RANGE_2, Patch::numAtoms, patch, and TIMEFACTOR.

Referenced by integrate().

01944 {
01945   NAMD_EVENT_RANGE_2(patch->flags.event_on,
01946       NamdProfileEvent::ADD_VELOCITY_TO_POSITION);
01947 #if CMK_BLUEGENEL
01948   CmiNetworkProgressAfter (0);
01949 #endif
01950   const BigReal dt = timestep / TIMEFACTOR;
01951   FullAtom *atom_arr  = patch->atom.begin();
01952   patch->addVelocityToPosition(atom_arr, dt, patch->numAtoms);
01953 }

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

Definition at line 146 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(), rescaleSoluteCharges(), 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_RESCALESOLUTECHARGES, SCRIPT_RESCALEVELS, SCRIPT_REVERT, SCRIPT_RUN, SimParameters::scriptArg1, ControllerBroadcasts::scriptBarrier, simParams, SimParameters::soluteScalingFactorCharge, submitCollections(), and terminate().

00147 {
00148   int scriptTask;
00149   int scriptSeq = 0;
00150   // Blocking receive for the script barrier.
00151   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
00152     switch ( scriptTask ) {
00153       case SCRIPT_OUTPUT:
00154         submitCollections(FILE_OUTPUT);
00155         break;
00156       case SCRIPT_FORCEOUTPUT:
00157         submitCollections(FORCE_OUTPUT);
00158         break;
00159       case SCRIPT_MEASURE:
00160         submitCollections(EVAL_MEASURE);
00161         break;
00162       case SCRIPT_REINITVELS:
00163         reinitVelocities();
00164         break;
00165       case SCRIPT_RESCALEVELS:
00166         rescaleVelocitiesByFactor(simParams->scriptArg1);
00167         break;
00168       case SCRIPT_RESCALESOLUTECHARGES:
00169         rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
00170         break;
00171       case SCRIPT_RELOADCHARGES:
00172         reloadCharges();
00173         break;
00174       case SCRIPT_CHECKPOINT:
00175         patch->checkpoint();
00176         checkpoint_berendsenPressure_count = berendsenPressure_count;
00177         break;
00178       case SCRIPT_REVERT:
00179         patch->revert();
00180         berendsenPressure_count = checkpoint_berendsenPressure_count;
00181         pairlistsAreValid = 0;
00182         break;
00183       case SCRIPT_CHECKPOINT_STORE:
00184       case SCRIPT_CHECKPOINT_LOAD:
00185       case SCRIPT_CHECKPOINT_SWAP:
00186       case SCRIPT_CHECKPOINT_FREE:
00187         patch->exchangeCheckpoint(scriptTask,berendsenPressure_count);
00188         break;
00189       case SCRIPT_ATOMSENDRECV:
00190       case SCRIPT_ATOMSEND:
00191       case SCRIPT_ATOMRECV:
00192         patch->exchangeAtoms(scriptTask);
00193         break;
00194       case SCRIPT_MINIMIZE:
00195         minimize();
00196         break;
00197       case SCRIPT_RUN:
00198       case SCRIPT_CONTINUE:
00199   //
00200   // DJH: Call a cleaned up version of integrate().
00201   //
00202   // We could test for simulation options and call a more basic version
00203   // of integrate() where we can avoid performing most tests.
00204   //
00205         integrate(scriptTask);
00206         break;
00207       default:
00208         NAMD_bug("Unknown task in Sequencer::algorithm");
00209     }
00210   }
00211   submitCollections(END_OF_RUN);
00212   terminate();
00213 }

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  step  )  [protected]

Definition at line 1534 of file Sequencer.C.

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

Referenced by integrate().

01535 {
01536   if ( simParams->berendsenPressureOn ) {
01537   berendsenPressure_count += 1;
01538   const int freq = simParams->berendsenPressureFreq;
01539   if ( ! (berendsenPressure_count % freq ) ) {
01540    berendsenPressure_count = 0;
01541    FullAtom *a = patch->atom.begin();
01542    int numAtoms = patch->numAtoms;
01543    // Blocking receive for the updated lattice scaling factor.
01544    Tensor factor = broadcast->positionRescaleFactor.get(step);
01545    patch->lattice.rescale(factor);
01546    if ( simParams->useGroupPressure )
01547    {
01548     int hgs;
01549     for ( int i = 0; i < numAtoms; i += hgs ) {
01550       int j;
01551       hgs = a[i].hydrogenGroupSize;
01552       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
01553         for ( j = i; j < (i+hgs); ++j ) {
01554           a[j].position = patch->lattice.apply_transform(
01555                                 a[j].fixedPosition,a[j].transform);
01556         }
01557         continue;
01558       }
01559       BigReal m_cm = 0;
01560       Position x_cm(0,0,0);
01561       for ( j = i; j < (i+hgs); ++j ) {
01562         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01563         m_cm += a[j].mass;
01564         x_cm += a[j].mass * a[j].position;
01565       }
01566       x_cm /= m_cm;
01567       Position new_x_cm = x_cm;
01568       patch->lattice.rescale(new_x_cm,factor);
01569       Position delta_x_cm = new_x_cm - x_cm;
01570       for ( j = i; j < (i+hgs); ++j ) {
01571         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01572           a[j].position = patch->lattice.apply_transform(
01573                                 a[j].fixedPosition,a[j].transform);
01574           continue;
01575         }
01576         a[j].position += delta_x_cm;
01577       }
01578     }
01579    }
01580    else
01581    {
01582     for ( int i = 0; i < numAtoms; ++i )
01583     {
01584       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01585         a[i].position = patch->lattice.apply_transform(
01586                                 a[i].fixedPosition,a[i].transform);
01587         continue;
01588       }
01589       patch->lattice.rescale(a[i].position,factor);
01590     }
01591    }
01592   }
01593   } else {
01594     berendsenPressure_count = 0;
01595   }
01596 }

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

Definition at line 2242 of file Sequencer.C.

References CompAtomExt::atomFixed, 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(), submitMinimizeReductions(), and submitReductions().

02243                                                                        {
02244 
02245   FullAtom *a = patch->atom.begin();
02246   int numAtoms = patch->numAtoms;
02247 
02248   for ( int j = 0; j < numAtoms; j++ ) {
02249     if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
02250       Vector dx = a[j].fixedPosition;
02251       // all negative because fixed atoms cancels these forces
02252       fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
02253       fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
02254       fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
02255       fixForceNormal -= patch->f[Results::normal][j];
02256       fixForceNbond -= patch->f[Results::nbond][j];
02257       fixForceSlow -= patch->f[Results::slow][j];
02258     }
02259   }
02260 }

BigReal Sequencer::calcKineticEnergy (  )  [protected]

Definition at line 1227 of file Sequencer.C.

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

Referenced by multigratorTemperature().

01227                                      {
01228   FullAtom *a = patch->atom.begin();
01229   int numAtoms = patch->numAtoms;
01230   BigReal kineticEnergy = 0.0;
01231   if ( simParams->pairInteractionOn ) {
01232     if ( simParams->pairInteractionSelf ) {
01233       for (int i = 0; i < numAtoms; ++i ) {
01234         if ( a[i].partition != 1 ) continue;
01235         kineticEnergy += a[i].mass * a[i].velocity.length2();
01236       }
01237     }
01238   } else {
01239     for (int i = 0; i < numAtoms; ++i ) {
01240       kineticEnergy += a[i].mass * a[i].velocity.length2();
01241     }
01242   }
01243   kineticEnergy *= 0.5;
01244   return kineticEnergy;
01245 }

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

Definition at line 1016 of file Sequencer.C.

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

Referenced by integrate().

01016                                                            {
01017 
01018   //
01019   // DJH: This test should be done in SimParameters.
01020   //
01021   if ( simParams->fixedAtomsOn )
01022     NAMD_die("Cannot zero momentum when fixed atoms are present.");
01023 
01024   // Blocking receive for the momentum correction vector.
01025   const Vector dv = broadcast->momentumCorrection.get(step);
01026 
01027   const Vector dx = dv * ( drifttime / TIMEFACTOR );
01028 
01029   FullAtom *a = patch->atom.begin();
01030   const int numAtoms = patch->numAtoms;
01031 
01032 if ( simParams->zeroMomentumAlt ) {
01033   for ( int i = 0; i < numAtoms; ++i ) {
01034     a[i].velocity += dv * a[i].recipMass;
01035     a[i].position += dx * a[i].recipMass;
01036   }
01037 } else {
01038   for ( int i = 0; i < numAtoms; ++i ) {
01039     a[i].velocity += dv;
01040     a[i].position += dx;
01041   }
01042 }
01043 
01044 }

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

Definition at line 2878 of file Sequencer.C.

References broadcast.

Referenced by integrate().

02878                                                     {
02879 #if USE_BARRIER
02880         if (doBarrier)
02881           // Blocking receive for the cycle barrier.
02882           broadcast->cycleBarrier.get(step);
02883 #endif
02884 }

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

Definition at line 1955 of file Sequencer.C.

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

Referenced by integrate().

01956 {
01957   if ( simParams->drudeHardWallOn ) {
01958     Tensor virial;
01959     Tensor *vp = ( pressure ? &virial : 0 );
01960     if ( patch->hardWallDrude(dt, vp, pressureProfileReduction) ) {
01961       iout << iERROR << "Constraint failure in HardWallDrude(); "
01962         << "simulation may become unstable.\n" << endi;
01963       Node::Object()->enableEarlyExit();
01964       terminate();
01965     }
01966     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01967   }
01968 }

void Sequencer::integrate ( int  scriptTask  )  [protected]

Definition at line 218 of file Sequencer.C.

References SimParameters::accelMDdihe, SimParameters::accelMDdual, SimParameters::accelMDOn, SimParameters::adaptTempOn, adaptTempT, adaptTempUpdate(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), Results::amdf, Patch::atomMapper, ResizeArray< Elem >::begin(), berendsenPressure(), SimParameters::colvarsOn, SimParameters::commOnly, ComputeMgr::computeGlobalObject, Node::computeMgr, correctMomentum(), cycleBarrier(), D_MSG, 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, Patch::getPatchID(), 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, NAMD_EVENT_START, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NAMD_PROFILE_START, NAMD_PROFILE_STOP, NamdProfileEventStr, 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, SPECIAL_PATCH_ID, SimParameters::statsOn, Flags::step, GlobalMaster::step, SimParameters::stepsPerCycle, stochRescaleVelocities(), submitCollections(), submitHalfstep(), submitMomentum(), submitReductions(), SimParameters::tclForcesOn, tcoupleVelocities(), TIMER_DONE, TIMER_INIT_WIDTH, TIMER_REPORT, TIMER_START, TIMER_STOP, traceBarrier(), SimParameters::traceStartStep, and SimParameters::zeroMomentum.

Referenced by algorithm().

00218                                         {
00219     char traceNote[24];
00220     char tracePrefix[20];
00221     sprintf(tracePrefix, "p:%d,s:",patch->patchID);
00222 //    patch->write_tip4_props();
00223 
00224     //
00225     // DJH: Copy all data into SOA (structure of arrays)
00226     // from AOS (array of structures) data structure.
00227     //
00228     //patch->copy_all_to_SOA();
00229 
00230 #ifdef TIMER_COLLECTION
00231     TimerSet& t = patch->timerSet;
00232 #endif
00233     TIMER_INIT_WIDTH(t, KICK, simParams->timerBinWidth);
00234     TIMER_INIT_WIDTH(t, MAXMOVE, simParams->timerBinWidth);
00235     TIMER_INIT_WIDTH(t, DRIFT, simParams->timerBinWidth);
00236     TIMER_INIT_WIDTH(t, PISTON, simParams->timerBinWidth);
00237     TIMER_INIT_WIDTH(t, SUBMITHALF, simParams->timerBinWidth);
00238     TIMER_INIT_WIDTH(t, VELBBK1, simParams->timerBinWidth);
00239     TIMER_INIT_WIDTH(t, VELBBK2, simParams->timerBinWidth);
00240     TIMER_INIT_WIDTH(t, RATTLE1, simParams->timerBinWidth);
00241     TIMER_INIT_WIDTH(t, SUBMITFULL, simParams->timerBinWidth);
00242     TIMER_INIT_WIDTH(t, SUBMITCOLLECT, simParams->timerBinWidth);
00243 
00244     int &step = patch->flags.step;
00245     step = simParams->firstTimestep;
00246 
00247     // drag switches
00248     const Bool rotDragOn = simParams->rotDragOn;
00249     const Bool movDragOn = simParams->movDragOn;
00250 
00251     const int commOnly = simParams->commOnly;
00252 
00253     int &maxForceUsed = patch->flags.maxForceUsed;
00254     int &maxForceMerged = patch->flags.maxForceMerged;
00255     maxForceUsed = Results::normal;
00256     maxForceMerged = Results::normal;
00257 
00258     const int numberOfSteps = simParams->N;
00259     const int stepsPerCycle = simParams->stepsPerCycle;
00260     const BigReal timestep = simParams->dt;
00261 
00262     // what MTS method?
00263     const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
00264 
00265     const int nonbondedFrequency = simParams->nonbondedFrequency;
00266     slowFreq = nonbondedFrequency;
00267     const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
00268     int &doNonbonded = patch->flags.doNonbonded;
00269     doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
00270     if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
00271     if ( doNonbonded ) maxForceUsed = Results::nbond;
00272 
00273     // Do we do full electrostatics?
00274     const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00275     const int fullElectFrequency = simParams->fullElectFrequency;
00276     if ( dofull ) slowFreq = fullElectFrequency;
00277     const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
00278     int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00279     doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
00280     if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
00281                                         maxForceMerged = Results::slow;
00282     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00283 
00284 //#ifndef UPPER_BOUND
00285     const Bool accelMDOn = simParams->accelMDOn;
00286     const Bool accelMDdihe = simParams->accelMDdihe;
00287     const Bool accelMDdual = simParams->accelMDdual;
00288     if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
00289 
00290     // Is adaptive tempering on?
00291     const Bool adaptTempOn = simParams->adaptTempOn;
00292     adaptTempT = simParams->initialTemp;
00293     if (simParams->langevinOn)
00294         adaptTempT = simParams->langevinTemp;
00295     else if (simParams->rescaleFreq > 0)
00296         adaptTempT = simParams->rescaleTemp;
00297         
00298 
00299     int &doMolly = patch->flags.doMolly;
00300     doMolly = simParams->mollyOn && doFullElectrostatics;
00301     // BEGIN LA
00302     int &doLoweAndersen = patch->flags.doLoweAndersen;
00303     doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00304     // END LA
00305 
00306     int &doGBIS = patch->flags.doGBIS;
00307     doGBIS = simParams->GBISOn;
00308 
00309     int &doLCPO = patch->flags.doLCPO;
00310     doLCPO = simParams->LCPOOn;
00311 
00312     int zeroMomentum = simParams->zeroMomentum;
00313     
00314     // Do we need to return forces to TCL script or Colvar module?
00315     int doTcl = simParams->tclForcesOn;
00316         int doColvars = simParams->colvarsOn;
00317 //#endif
00318     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00319 
00320     // Bother to calculate energies?
00321     int &doEnergy = patch->flags.doEnergy;
00322     int energyFrequency = simParams->outputEnergies;
00323 #ifndef UPPER_BOUND
00324     const int reassignFreq = simParams->reassignFreq;
00325 #endif
00326 
00327     int &doVirial = patch->flags.doVirial;
00328     doVirial = 1;
00329 
00330   if ( scriptTask == SCRIPT_RUN ) {
00331 
00332 #ifndef UPPER_BOUND
00333 //    printf("Doing initial rattle\n");
00334 #ifndef UPPER_BOUND
00335 D_MSG("rattle1()");
00336     TIMER_START(t, RATTLE1);
00337     rattle1(0.,0);  // enforce rigid bond constraints on initial positions
00338     TIMER_STOP(t, RATTLE1);
00339 #endif
00340 
00341     if (simParams->lonepairs) {
00342       patch->atomMapper->registerIDsFullAtom(
00343                 patch->atom.begin(),patch->atom.end());
00344     }
00345 
00346     if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00347        reassignVelocities(timestep,step);
00348     }
00349 #endif
00350 
00351     doEnergy = ! ( step % energyFrequency );
00352 #ifndef UPPER_BOUND
00353     if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00354     //Update energy every timestep for adaptive tempering
00355     if ( adaptTempOn ) doEnergy=1;
00356 #endif
00357 D_MSG("runComputeObjects()");
00358     runComputeObjects(1,step<numberOfSteps); // must migrate here!
00359 #ifndef UPPER_BOUND
00360     rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD 
00361     adaptTempUpdate(step); // update adaptive tempering temperature
00362 #endif
00363 
00364 #ifndef UPPER_BOUND
00365     if ( staleForces || doTcl || doColvars ) {
00366       if ( doNonbonded ) saveForce(Results::nbond);
00367       if ( doFullElectrostatics ) saveForce(Results::slow);
00368     }
00369     if ( ! commOnly ) {
00370 D_MSG("newtonianVelocities()");
00371       TIMER_START(t, KICK);
00372       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
00373       TIMER_STOP(t, KICK);
00374     }
00375     minimizationQuenchVelocity();
00376 #ifndef UPPER_BOUND
00377 D_MSG("rattle1()");
00378     TIMER_START(t, RATTLE1);
00379     rattle1(-timestep,0);
00380     TIMER_STOP(t, RATTLE1);
00381 #endif
00382 D_MSG("submitHalfstep()");
00383     TIMER_START(t, SUBMITHALF);
00384     submitHalfstep(step);
00385     TIMER_STOP(t, SUBMITHALF);
00386     if ( ! commOnly ) {
00387 D_MSG("newtonianVelocities()");
00388       TIMER_START(t, KICK);
00389       newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
00390       TIMER_STOP(t, KICK);
00391     }
00392 D_MSG("rattle1()");
00393     TIMER_START(t, RATTLE1);
00394     rattle1(timestep,1);
00395     TIMER_STOP(t, RATTLE1);
00396     if (doTcl || doColvars)  // include constraint forces
00397       computeGlobal->saveTotalForces(patch);
00398 D_MSG("submitHalfstep()");
00399     TIMER_START(t, SUBMITHALF);
00400     submitHalfstep(step);
00401     TIMER_STOP(t, SUBMITHALF);
00402     if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00403     if ( ! commOnly ) {
00404 D_MSG("newtonianVelocities()");
00405       TIMER_START(t, KICK);
00406       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
00407       TIMER_STOP(t, KICK);
00408     }
00409 #endif
00410 D_MSG("submitReductions()");
00411     TIMER_START(t, SUBMITFULL);
00412     submitReductions(step);
00413     TIMER_STOP(t, SUBMITFULL);
00414 #ifndef UPPER_BOUND
00415     if(0){ // if(traceIsOn()){
00416         traceUserEvent(eventEndOfTimeStep);
00417         sprintf(traceNote, "%s%d",tracePrefix,step); 
00418         traceUserSuppliedNote(traceNote);
00419     }
00420 #endif
00421     rebalanceLoad(step);
00422 
00423   } // scriptTask == SCRIPT_RUN
00424 
00425 #ifndef UPPER_BOUND
00426   bool doMultigratorRattle = false;
00427 #endif
00428 
00429   //
00430   // DJH: There are a lot of mod operations below and elsewhere to
00431   // test step number against the frequency of something happening.
00432   // Mod and integer division are expensive!
00433   // Might be better to replace with counters and test equality.
00434   //
00435 #if 0
00436     for(int i = 0; i < NamdProfileEvent::EventsCount; i++)
00437         CkPrintf("-------------- [%d] %s -------------\n", i, NamdProfileEventStr[i]);
00438 #endif
00439 
00440 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
00441   int& eon = patch->flags.event_on;
00442   int epid = (simParams->beginEventPatchID <= patch->getPatchID()
00443       && patch->getPatchID() <= simParams->endEventPatchID);
00444   int beginStep = simParams->beginEventStep;
00445   int endStep = simParams->endEventStep;
00446   bool controlProfiling = patch->getPatchID() == 0;
00447 #endif
00448  
00449     for ( ++step; step <= numberOfSteps; ++step )
00450     {
00451 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)
00452       eon = epid && (beginStep < step && step <= endStep);
00453 
00454       if (controlProfiling && step == beginStep) {
00455         NAMD_PROFILE_START();
00456       }
00457       if (controlProfiling && step == endStep) {
00458         NAMD_PROFILE_STOP();
00459       }
00460       char buf[32];
00461       sprintf(buf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::INTEGRATE_1], patch->getPatchID());
00462       NAMD_EVENT_START_EX(eon, NamdProfileEvent::INTEGRATE_1, buf);
00463 #endif
00464 #ifndef UPPER_BOUND
00465 
00466       rescaleVelocities(step);
00467       tcoupleVelocities(timestep,step);
00468       stochRescaleVelocities(timestep,step);
00469       berendsenPressure(step);
00470 
00471       if ( ! commOnly ) {
00472         TIMER_START(t, KICK);
00473         newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics); 
00474         TIMER_STOP(t, KICK);
00475       }
00476 
00477       // We do RATTLE here if multigrator thermostat was applied in the previous step
00478       if (doMultigratorRattle) rattle1(timestep, 1);
00479 
00480       /* reassignment based on half-step velocities
00481          if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00482          addVelocityToPosition(0.5*timestep);
00483          reassignVelocities(timestep,step);
00484          addVelocityToPosition(0.5*timestep);
00485          rattle1(0.,0);
00486          rattle1(-timestep,0);
00487          addVelocityToPosition(-1.0*timestep);
00488          rattle1(timestep,0);
00489          } */
00490 
00491       TIMER_START(t, MAXMOVE);
00492       maximumMove(timestep);
00493       TIMER_STOP(t, MAXMOVE);
00494 
00495       NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_1);  // integrate 1
00496 
00497       if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) {
00498         if ( ! commOnly ) {
00499           TIMER_START(t, DRIFT);
00500           addVelocityToPosition(0.5*timestep);
00501           TIMER_STOP(t, DRIFT);
00502         }
00503         // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
00504         langevinVelocities(timestep);
00505 
00506         // There is a blocking receive inside of langevinPiston()
00507         // that might suspend the current thread of execution,
00508         // so split profiling around this conditional block.
00509         langevinPiston(step);
00510 
00511         if ( ! commOnly ) {
00512           TIMER_START(t, DRIFT);
00513           addVelocityToPosition(0.5*timestep);
00514           TIMER_STOP(t, DRIFT);
00515         }
00516       } else {
00517         // If Langevin is not used, take full time step directly instread of two half steps      
00518         if ( ! commOnly ) {
00519           TIMER_START(t, DRIFT);
00520           addVelocityToPosition(timestep); 
00521           TIMER_STOP(t, DRIFT);
00522         }
00523       }
00524 
00525       NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_2);
00526 
00527       // impose hard wall potential for Drude bond length
00528       hardWallDrude(timestep, 1);
00529 
00530       minimizationQuenchVelocity();
00531 #endif // UPPER_BOUND
00532 
00533       doNonbonded = !(step%nonbondedFrequency);
00534       doFullElectrostatics = (dofull && !(step%fullElectFrequency));
00535 
00536 #ifndef UPPER_BOUND 
00537       if ( zeroMomentum && doFullElectrostatics ) {
00538         // There is a blocking receive inside of correctMomentum().
00539         correctMomentum(step,slowstep);
00540       }
00541 
00542       // There are NO sends in submitHalfstep() just local summation 
00543       // into the Reduction struct.
00544       TIMER_START(t, SUBMITHALF);
00545       submitHalfstep(step);
00546       TIMER_STOP(t, SUBMITHALF);
00547 
00548       doMolly = simParams->mollyOn && doFullElectrostatics;
00549       // BEGIN LA
00550       doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
00551       // END LA
00552 
00553       maxForceUsed = Results::normal;
00554       if ( doNonbonded ) maxForceUsed = Results::nbond;
00555       if ( doFullElectrostatics ) maxForceUsed = Results::slow;
00556       if ( accelMDOn && (accelMDdihe || accelMDdual))  maxForceUsed = Results::amdf;
00557 
00558       // Migrate Atoms on stepsPerCycle
00559       doEnergy = ! ( step % energyFrequency );
00560       if ( accelMDOn && !accelMDdihe ) doEnergy=1;
00561       if ( adaptTempOn ) doEnergy=1; 
00562 
00563       // Multigrator 
00564       if (simParams->multigratorOn) {
00565         doVirial = (!(step % energyFrequency) || ((simParams->outputPressure > 0) && !(step % simParams->outputPressure)) 
00566           || !(step % simParams->multigratorPressureFreq));
00567         doKineticEnergy = (!(step % energyFrequency) || !(step % simParams->multigratorTemperatureFreq));
00568         doMomenta = (simParams->outputMomenta > 0) && !(step % simParams->outputMomenta);
00569       } else {
00570         doVirial = 1;
00571         doKineticEnergy = 1;
00572         doMomenta = 1;
00573       }
00574 #endif
00575       NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_2);  // integrate 2
00576 
00577       // The current thread of execution will suspend in runComputeObjects().
00578       runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00579 
00580       NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_3);
00581 
00582 #ifndef UPPER_BOUND
00583       rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
00584      
00585       if ( staleForces || doTcl || doColvars ) {
00586         if ( doNonbonded ) saveForce(Results::nbond);
00587         if ( doFullElectrostatics ) saveForce(Results::slow);
00588       }
00589 
00590       // reassignment based on full-step velocities
00591       if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
00592         reassignVelocities(timestep,step);
00593         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00594         rattle1(-timestep,0);
00595       }
00596 
00597       if ( ! commOnly ) {
00598         TIMER_START(t, VELBBK1);
00599         langevinVelocitiesBBK1(timestep);
00600         TIMER_STOP(t, VELBBK1);
00601         TIMER_START(t, KICK);
00602         newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00603         TIMER_STOP(t, KICK);
00604         TIMER_START(t, VELBBK2);
00605         langevinVelocitiesBBK2(timestep);
00606         TIMER_STOP(t, VELBBK2);
00607       }
00608 
00609       // add drag to each atom's positions
00610       if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
00611       if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
00612 
00613       TIMER_START(t, RATTLE1);
00614       rattle1(timestep,1);
00615       TIMER_STOP(t, RATTLE1);
00616       if (doTcl || doColvars)  // include constraint forces
00617         computeGlobal->saveTotalForces(patch);
00618 
00619       TIMER_START(t, SUBMITHALF);
00620       submitHalfstep(step);
00621       TIMER_STOP(t, SUBMITHALF);
00622       if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
00623 
00624       if ( ! commOnly ) {
00625         TIMER_START(t, KICK);
00626         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
00627         TIMER_STOP(t, KICK);
00628       }
00629 
00630         // rattle2(timestep,step);
00631 #endif
00632 
00633         TIMER_START(t, SUBMITFULL);
00634         submitReductions(step);
00635         TIMER_STOP(t, SUBMITFULL);
00636         TIMER_START(t, SUBMITCOLLECT);
00637         submitCollections(step);
00638         TIMER_STOP(t, SUBMITCOLLECT);
00639 #ifndef UPPER_BOUND
00640        //Update adaptive tempering temperature
00641         adaptTempUpdate(step);
00642 
00643       // Multigrator temperature and pressure steps
00644       multigratorTemperature(step, 1);
00645       multigratorPressure(step, 1);
00646       multigratorPressure(step, 2);
00647       multigratorTemperature(step, 2);
00648       doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
00649 
00650       NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_3); // integrate 3
00651 #endif
00652 
00653 #if CYCLE_BARRIER
00654         cycleBarrier(!((step+1) % stepsPerCycle), step);
00655 #elif PME_BARRIER
00656         cycleBarrier(doFullElectrostatics, step);
00657 #elif  STEP_BARRIER
00658         cycleBarrier(1, step);
00659 #endif
00660 
00661 #ifndef UPPER_BOUND
00662          if(Node::Object()->specialTracing || simParams->statsOn){
00663                  int bstep = simParams->traceStartStep;
00664                  int estep = bstep + simParams->numTraceSteps;
00665                  if(step == bstep || step == estep){
00666                          traceBarrier(step);
00667                  }                       
00668          }
00669 
00670 #ifdef MEASURE_NAMD_WITH_PAPI    
00671          if(simParams->papiMeasure) {
00672                  int bstep = simParams->papiMeasureStartStep;
00673                  int estep = bstep + simParams->numPapiMeasureSteps;
00674                  if(step == bstep || step==estep) {
00675                          papiMeasureBarrier(step);
00676                  }
00677          }
00678 #endif
00679           
00680         if(0){ // if(traceIsOn()){
00681             traceUserEvent(eventEndOfTimeStep);
00682             sprintf(traceNote, "%s%d",tracePrefix,step); 
00683             traceUserSuppliedNote(traceNote);
00684         }
00685 #endif // UPPER_BOUND
00686         rebalanceLoad(step);
00687 
00688 #if PME_BARRIER
00689         // a step before PME
00690         cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
00691 #endif
00692 
00693 #if USE_HPM
00694         if(step == START_HPM_STEP)
00695           (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
00696 
00697         if(step == STOP_HPM_STEP)
00698           (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
00699 #endif
00700 
00701     }
00702 
00703   TIMER_DONE(t);
00704 #ifdef TIMER_COLLECTION
00705   if (patch->patchID == SPECIAL_PATCH_ID) {
00706     printf("Timer collection reporting in microseconds for "
00707         "Patch %d\n", patch->patchID);
00708     TIMER_REPORT(t);
00709   }
00710 #endif // TIMER_COLLECTION
00711     //
00712     // DJH: Copy updates of SOA back into AOS.
00713     //
00714     //patch->copy_updates_to_AOS();
00715 }

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

Definition at line 1598 of file Sequencer.C.

References Lattice::apply_transform(), CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), broadcast, SimParameters::fixedAtomsOn, for(), SimpleBroadcastObject< T >::get(), CompAtomExt::groupFixed, CompAtom::hydrogenGroupSize, CompAtomExt::id, 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, TIMER_START, TIMER_STOP, SimParameters::useGroupPressure, FullAtom::velocity, Vector::x, Tensor::xx, Vector::y, Tensor::yy, Vector::z, and Tensor::zz.

Referenced by integrate().

01599 {
01600   if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
01601   {
01602     //
01603     // DJH: Loops below simplify if we lift out special cases of fixed atoms
01604     // and pressure excluded atoms and make them their own branch.
01605     //
01606    FullAtom *a = patch->atom.begin();
01607    int numAtoms = patch->numAtoms;
01608    // Blocking receive for the updated lattice scaling factor.
01609    Tensor factor = broadcast->positionRescaleFactor.get(step);
01610    TIMER_START(patch->timerSet, PISTON);
01611    // JCP FIX THIS!!!
01612    Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
01613    patch->lattice.rescale(factor);
01614    Molecule *mol = Node::Object()->molecule;
01615    if ( simParams->useGroupPressure )
01616    {
01617     int hgs;
01618     for ( int i = 0; i < numAtoms; i += hgs ) {
01619       int j;
01620       hgs = a[i].hydrogenGroupSize;
01621       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
01622         for ( j = i; j < (i+hgs); ++j ) {
01623           a[j].position = patch->lattice.apply_transform(
01624                                 a[j].fixedPosition,a[j].transform);
01625         }
01626         continue;
01627       }
01628       BigReal m_cm = 0;
01629       Position x_cm(0,0,0);
01630       Velocity v_cm(0,0,0);
01631       for ( j = i; j < (i+hgs); ++j ) {
01632         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
01633         m_cm += a[j].mass;
01634         x_cm += a[j].mass * a[j].position;
01635         v_cm += a[j].mass * a[j].velocity;
01636       }
01637       x_cm /= m_cm;
01638       Position new_x_cm = x_cm;
01639       patch->lattice.rescale(new_x_cm,factor);
01640       Position delta_x_cm = new_x_cm - x_cm;
01641       v_cm /= m_cm;
01642       Velocity delta_v_cm;
01643       delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
01644       delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
01645       delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
01646       for ( j = i; j < (i+hgs); ++j ) {
01647         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
01648           a[j].position = patch->lattice.apply_transform(
01649                                 a[j].fixedPosition,a[j].transform);
01650           continue;
01651         }
01652         if ( mol->is_atom_exPressure(a[j].id) ) continue;
01653         a[j].position += delta_x_cm;
01654         a[j].velocity += delta_v_cm;
01655       }
01656     }
01657    }
01658    else
01659    {
01660     for ( int i = 0; i < numAtoms; ++i )
01661     {
01662       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
01663         a[i].position = patch->lattice.apply_transform(
01664                                 a[i].fixedPosition,a[i].transform);
01665         continue;
01666       }
01667       if ( mol->is_atom_exPressure(a[i].id) ) continue;
01668       patch->lattice.rescale(a[i].position,factor);
01669       a[i].velocity.x *= velFactor.x;
01670       a[i].velocity.y *= velFactor.y;
01671       a[i].velocity.z *= velFactor.z;
01672     }
01673    }
01674    TIMER_STOP(patch->timerSet, PISTON);
01675   }
01676 }

void Sequencer::langevinVelocities ( BigReal  dt_fs  )  [protected]

Definition at line 1317 of file Sequencer.C.

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

Referenced by integrate().

01318 {
01319 // This routine is used for the BAOAB integrator,
01320 // Ornstein-Uhlenbeck exact solve for the O-part.
01321 // See B. Leimkuhler and C. Matthews, AMRX (2012)
01322 // Routine originally written by JPhillips, with fresh errors by CMatthews June2012
01323 
01324   if ( simParams->langevinOn && simParams->langevin_useBAOAB )
01325   {
01326     FullAtom *a = patch->atom.begin();
01327     int numAtoms = patch->numAtoms;
01328     Molecule *molecule = Node::Object()->molecule;
01329     BigReal dt = dt_fs * 0.001;  // convert to ps
01330     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
01331     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
01332     {
01333         kbT = BOLTZMANN*adaptTempT;
01334     }
01335 
01336     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01337     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01338 
01339     for ( int i = 0; i < numAtoms; ++i )
01340     {
01341       BigReal dt_gamma = dt * a[i].langevinParam;
01342       if ( ! dt_gamma ) continue;
01343 
01344       BigReal f1 = exp( -dt_gamma );
01345       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT * 
01346                          ( a[i].partition ? tempFactor : 1.0 ) * 
01347                          a[i].recipMass );
01348       a[i].velocity *= f1;
01349       a[i].velocity += f2 * random->gaussian_vector();
01350     }
01351   }
01352 }

void Sequencer::langevinVelocitiesBBK1 ( BigReal  dt_fs  )  [protected]

Definition at line 1354 of file Sequencer.C.

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

Referenced by integrate().

01355 {
01356   NAMD_EVENT_RANGE_2(patch->flags.event_on,
01357       NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1);
01358   if ( simParams->langevinOn && !simParams->langevin_useBAOAB )
01359   {
01360     FullAtom *a = patch->atom.begin();
01361     int numAtoms = patch->numAtoms;
01362     Molecule *molecule = Node::Object()->molecule;
01363     BigReal dt = dt_fs * 0.001;  // convert to ps
01364     int i;
01365 
01366     if (simParams->drudeOn) {
01367       for (i = 0;  i < numAtoms;  i++) {
01368 
01369         if (i < numAtoms-1 &&
01370             a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
01371           //printf("*** Found Drude particle %d\n", a[i+1].id);
01372           // i+1 is a Drude particle with parent i
01373 
01374           // convert from Cartesian coordinates to (COM,bond) coordinates
01375           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
01376           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
01377           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
01378           BigReal dt_gamma;
01379 
01380           // use Langevin damping factor i for v_com
01381           dt_gamma = dt * a[i].langevinParam;
01382           if (dt_gamma != 0.0) {
01383             v_com *= ( 1. - 0.5 * dt_gamma );
01384           }
01385 
01386           // use Langevin damping factor i+1 for v_bnd
01387           dt_gamma = dt * a[i+1].langevinParam;
01388           if (dt_gamma != 0.0) {
01389             v_bnd *= ( 1. - 0.5 * dt_gamma );
01390           }
01391 
01392           // convert back
01393           a[i].velocity = v_com - m * v_bnd;
01394           a[i+1].velocity = v_bnd + a[i].velocity;
01395 
01396           i++;  // +1 from loop, we've updated both particles
01397         }
01398         else {
01399           BigReal dt_gamma = dt * a[i].langevinParam;
01400           if ( ! dt_gamma ) continue;
01401 
01402           a[i].velocity *= ( 1. - 0.5 * dt_gamma );
01403         }
01404 
01405       } // end for
01406     } // end if drudeOn
01407     else {
01408 
01409       //
01410       // DJH: The conditional inside loop prevents vectorization and doesn't
01411       // avoid much work since addition and multiplication are cheap.
01412       //
01413       for ( i = 0; i < numAtoms; ++i )
01414       {
01415         BigReal dt_gamma = dt * a[i].langevinParam;
01416         if ( ! dt_gamma ) continue;
01417 
01418         a[i].velocity *= ( 1. - 0.5 * dt_gamma );
01419       }
01420 
01421     } // end else
01422 
01423   } // end if langevinOn
01424 }

void Sequencer::langevinVelocitiesBBK2 ( BigReal  dt_fs  )  [protected]

Definition at line 1427 of file Sequencer.C.

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

Referenced by integrate().

01428 {
01429   NAMD_EVENT_RANGE_2(patch->flags.event_on,
01430       NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2);
01431   if ( simParams->langevinOn && !simParams->langevin_useBAOAB ) 
01432   {
01433     //
01434     // DJH: This call is expensive. Avoid calling when gammas don't differ.
01435     // Set flag in SimParameters and make this call conditional.
01436     // 
01437     TIMER_START(patch->timerSet, RATTLE1);
01438     rattle1(dt_fs,1);  // conserve momentum if gammas differ
01439     TIMER_STOP(patch->timerSet, RATTLE1);
01440 
01441     FullAtom *a = patch->atom.begin();
01442     int numAtoms = patch->numAtoms;
01443     Molecule *molecule = Node::Object()->molecule;
01444     BigReal dt = dt_fs * 0.001;  // convert to ps
01445     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
01446     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
01447     {
01448         kbT = BOLTZMANN*adaptTempT;
01449     }
01450     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01451     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01452     int i;
01453 
01454     if (simParams->drudeOn) {
01455       BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp);  // drude bond Temp
01456 
01457       for (i = 0;  i < numAtoms;  i++) {
01458 
01459         if (i < numAtoms-1 &&
01460             a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
01461           //printf("*** Found Drude particle %d\n", a[i+1].id);
01462           // i+1 is a Drude particle with parent i
01463 
01464           // convert from Cartesian coordinates to (COM,bond) coordinates
01465           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
01466           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
01467           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
01468           BigReal dt_gamma;
01469 
01470           // use Langevin damping factor i for v_com
01471           dt_gamma = dt * a[i].langevinParam;
01472           if (dt_gamma != 0.0) {
01473             BigReal mass = a[i].mass + a[i+1].mass;
01474             v_com += random->gaussian_vector() *
01475               sqrt( 2 * dt_gamma * kbT *
01476                   ( a[i].partition ? tempFactor : 1.0 ) / mass );
01477             v_com /= ( 1. + 0.5 * dt_gamma );
01478           }
01479 
01480           // use Langevin damping factor i+1 for v_bnd
01481           dt_gamma = dt * a[i+1].langevinParam;
01482           if (dt_gamma != 0.0) {
01483             BigReal mass = a[i+1].mass * (1. - m);
01484             v_bnd += random->gaussian_vector() *
01485               sqrt( 2 * dt_gamma * kbT_bnd *
01486                   ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
01487             v_bnd /= ( 1. + 0.5 * dt_gamma );
01488           }
01489 
01490           // convert back
01491           a[i].velocity = v_com - m * v_bnd;
01492           a[i+1].velocity = v_bnd + a[i].velocity;
01493 
01494           i++;  // +1 from loop, we've updated both particles
01495         }
01496         else { 
01497           BigReal dt_gamma = dt * a[i].langevinParam;
01498           if ( ! dt_gamma ) continue;
01499 
01500           a[i].velocity += random->gaussian_vector() *
01501             sqrt( 2 * dt_gamma * kbT *
01502                 ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
01503           a[i].velocity /= ( 1. + 0.5 * dt_gamma );
01504         }
01505 
01506       } // end for
01507     } // end if drudeOn
01508     else {
01509 
01510       //
01511       // DJH: For case using same gamma (the Langevin parameter),
01512       // no partitions (e.g. FEP), and no adaptive tempering (adaptTempMD),
01513       // we can precompute constants. Then by lifting the RNG from the
01514       // loop (filling up an array of random numbers), we can vectorize
01515       // loop and simplify arithmetic to just addition and multiplication.
01516       //
01517       for ( i = 0; i < numAtoms; ++i )
01518       {
01519         BigReal dt_gamma = dt * a[i].langevinParam;
01520         if ( ! dt_gamma ) continue;
01521 
01522         a[i].velocity += random->gaussian_vector() *
01523           sqrt( 2 * dt_gamma * kbT *
01524               ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
01525         a[i].velocity /= ( 1. + 0.5 * dt_gamma );
01526       }
01527 
01528     } // end else
01529 
01530   } // end if langevinOn
01531 }

void Sequencer::maximumMove ( BigReal  timestep  )  [protected]

Definition at line 2027 of file Sequencer.C.

References ResizeArray< Elem >::begin(), SimParameters::cutoff, colvarproxy_namd::dt(), Node::enableEarlyExit(), endi(), Patch::flags, CompAtomExt::id, iERROR(), iout, Vector::length(), Vector::length2(), SimParameters::maximumMove, NAMD_EVENT_RANGE_2, Patch::numAtoms, Node::Object(), patch, Patch::patchID, PDBVELFACTOR, simParams, terminate(), TIMEFACTOR, and FullAtom::velocity.

Referenced by integrate().

02028 {
02029   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::MAXIMUM_MOVE);
02030 
02031   FullAtom *a = patch->atom.begin();
02032   int numAtoms = patch->numAtoms;
02033   if ( simParams->maximumMove ) {
02034     const BigReal dt = timestep / TIMEFACTOR;
02035     const BigReal maxvel = simParams->maximumMove / dt;
02036     const BigReal maxvel2 = maxvel * maxvel;
02037     for ( int i=0; i<numAtoms; ++i ) {
02038       if ( a[i].velocity.length2() > maxvel2 ) {
02039         a[i].velocity *= ( maxvel / a[i].velocity.length() );
02040       }
02041     }
02042   } else {
02043     const BigReal dt = timestep / TIMEFACTOR;
02044     const BigReal maxvel = simParams->cutoff / dt;
02045     const BigReal maxvel2 = maxvel * maxvel;
02046     int killme = 0;
02047     for ( int i=0; i<numAtoms; ++i ) {
02048       killme = killme || ( a[i].velocity.length2() > maxvel2 );
02049     }
02050     if ( killme ) {
02051       killme = 0;
02052       for ( int i=0; i<numAtoms; ++i ) {
02053         if ( a[i].velocity.length2() > maxvel2 ) {
02054           ++killme;
02055           iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
02056             << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
02057             << ( PDBVELFACTOR * maxvel ) << ", atom "
02058             << i << " of " << numAtoms << " on patch "
02059             << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
02060         }
02061       }
02062       iout << iERROR << 
02063         "Atoms moving too fast; simulation has become unstable ("
02064         << killme << " atoms on patch " << patch->patchID
02065         << " pe " << CkMyPe() << ").\n" << endi;
02066       Node::Object()->enableEarlyExit();
02067       terminate();
02068     }
02069   }
02070 }

void Sequencer::minimizationQuenchVelocity ( void   )  [protected]

Definition at line 2072 of file Sequencer.C.

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

Referenced by integrate().

02073 {
02074   if ( simParams->minimizeOn ) {
02075     FullAtom *a = patch->atom.begin();
02076     int numAtoms = patch->numAtoms;
02077     for ( int i=0; i<numAtoms; ++i ) {
02078       a[i].velocity = 0.;
02079     }
02080   }
02081 }

void Sequencer::minimize (  )  [protected]

Definition at line 760 of file Sequencer.C.

References 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, SimpleBroadcastObject< T >::get(), 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().

00760                          {
00761     //
00762     // DJH: Copy all data into SOA (structure of arrays)
00763     // from AOS (array of structures) data structure.
00764     //
00765     //patch->copy_all_to_SOA();
00766 
00767   const int numberOfSteps = simParams->N;
00768   const int stepsPerCycle = simParams->stepsPerCycle;
00769   int &step = patch->flags.step;
00770   step = simParams->firstTimestep;
00771 
00772   int &maxForceUsed = patch->flags.maxForceUsed;
00773   int &maxForceMerged = patch->flags.maxForceMerged;
00774   maxForceUsed = Results::normal;
00775   maxForceMerged = Results::normal;
00776   int &doNonbonded = patch->flags.doNonbonded;
00777   doNonbonded = 1;
00778   maxForceUsed = Results::nbond;
00779   maxForceMerged = Results::nbond;
00780   const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
00781   int &doFullElectrostatics = patch->flags.doFullElectrostatics;
00782   doFullElectrostatics = dofull;
00783   if ( dofull ) {
00784     maxForceMerged = Results::slow;
00785     maxForceUsed = Results::slow;
00786   }
00787   int &doMolly = patch->flags.doMolly;
00788   doMolly = simParams->mollyOn && doFullElectrostatics;
00789   // BEGIN LA
00790   int &doLoweAndersen = patch->flags.doLoweAndersen;
00791   doLoweAndersen = 0;
00792   // END LA
00793 
00794   int &doGBIS = patch->flags.doGBIS;
00795   doGBIS = simParams->GBISOn;
00796 
00797     int &doLCPO = patch->flags.doLCPO;
00798     doLCPO = simParams->LCPOOn;
00799 
00800     int doTcl = simParams->tclForcesOn;
00801         int doColvars = simParams->colvarsOn;
00802     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
00803 
00804   int &doEnergy = patch->flags.doEnergy;
00805   doEnergy = 1;
00806 
00807   patch->rattle1(0.,0,0);  // enforce rigid bond constraints on initial positions
00808 
00809   if (simParams->lonepairs) {
00810     patch->atomMapper->registerIDsFullAtom(
00811                 patch->atom.begin(),patch->atom.end());
00812   }
00813 
00814   runComputeObjects(1,step<numberOfSteps); // must migrate here!
00815 
00816   if ( doTcl || doColvars ) {
00817     if ( doNonbonded ) saveForce(Results::nbond);
00818     if ( doFullElectrostatics ) saveForce(Results::slow);
00819     computeGlobal->saveTotalForces(patch);
00820   }
00821 
00822   BigReal fmax2 = TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00823 
00824   submitMinimizeReductions(step,fmax2);
00825   rebalanceLoad(step);
00826 
00827   int downhill = 1;  // start out just fixing bad contacts
00828   int minSeq = 0;
00829   for ( ++step; step <= numberOfSteps; ++step ) {
00830 
00831    // Blocking receive for the minimization coefficient.
00832    BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
00833 
00834    if ( downhill ) {
00835     if ( c ) minimizeMoveDownhill(fmax2);
00836     else {
00837       downhill = 0;
00838       fmax2 *= 10000.;
00839     }
00840    }
00841    if ( ! downhill ) {
00842     if ( ! c ) {  // new direction
00843 
00844       // Blocking receive for the minimization coefficient.
00845       c = broadcast->minimizeCoefficient.get(minSeq++);
00846 
00847       newMinimizeDirection(c);  // v = c * v + f
00848 
00849       // Blocking receive for the minimization coefficient.
00850       c = broadcast->minimizeCoefficient.get(minSeq++);
00851 
00852     }  // same direction
00853     newMinimizePosition(c);  // x = x + c * v
00854    }
00855 
00856     runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
00857     if ( doTcl || doColvars ) {
00858       if ( doNonbonded ) saveForce(Results::nbond);
00859       if ( doFullElectrostatics ) saveForce(Results::slow);
00860       computeGlobal->saveTotalForces(patch);
00861     }
00862     submitMinimizeReductions(step,fmax2);
00863     submitCollections(step, 1);  // write out zeros for velocities
00864     rebalanceLoad(step);
00865   }
00866   quenchVelocities();  // zero out bogus velocity
00867 
00868     //
00869     // DJH: Copy updates of SOA back into AOS.
00870     //
00871     //patch->copy_updates_to_AOS();
00872 }

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

Definition at line 875 of file Sequencer.C.

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

Referenced by minimize().

00875                                                   {
00876 
00877   FullAtom *a = patch->atom.begin();
00878   Force *f1 = patch->f[Results::normal].begin();  // includes nbond and slow
00879   int numAtoms = patch->numAtoms;
00880 
00881   for ( int i = 0; i < numAtoms; ++i ) {
00882     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
00883     Force f = f1[i];
00884     if ( f.length2() > fmax2 ) {
00885       a[i].position += ( 0.1 * f.unit() );
00886       int hgs = a[i].hydrogenGroupSize;  // 0 if not parent
00887       for ( int j=1; j<hgs; ++j ) {
00888         a[++i].position += ( 0.1 * f.unit() );
00889       }
00890     }
00891   }
00892 
00893   patch->rattle1(0.,0,0);
00894 }

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

(stepstepsPerCycle)

Definition at line 1084 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, ResizeArray< Elem >::begin(), broadcast, calcFixVirial(), Flags::doFullElectrostatics, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, SimpleBroadcastObject< T >::get(), CompAtom::hydrogenGroupSize, 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_MARGIN_VIOLATIONS, Lattice::rescale(), runComputeObjects(), scalePositionsVelocities(), simParams, Results::slow, SimParameters::stepsPerCycle, SubmitReduction::submit(), SimParameters::useGroupPressure, FullAtom::velocity, ControllerBroadcasts::velocityRescaleTensor, and ControllerBroadcasts::velocityRescaleTensor2.

Referenced by integrate().

01084                                                             {
01085 // Calculate new positions, momenta, and volume using positionRescaleFactor and 
01086 // velocityRescaleTensor values returned from Controller::multigratorPressureCalcScale()
01087   if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
01088     FullAtom *a = patch->atom.begin();
01089     int numAtoms = patch->numAtoms;
01090 
01091     // Blocking receive (get) scaling factors from Controller
01092     Tensor scaleTensor    = (callNumber == 1) ? broadcast->positionRescaleFactor.get(step) : broadcast->positionRescaleFactor2.get(step);
01093     Tensor velScaleTensor = (callNumber == 1) ? broadcast->velocityRescaleTensor.get(step) : broadcast->velocityRescaleTensor2.get(step);
01094     Tensor posScaleTensor = scaleTensor;
01095     posScaleTensor -= Tensor::identity();
01096     if (simParams->useGroupPressure) {
01097       velScaleTensor -= Tensor::identity();
01098     }
01099 
01100     // Scale volume
01101     patch->lattice.rescale(scaleTensor);
01102     // Scale positions and velocities
01103     scalePositionsVelocities(posScaleTensor, velScaleTensor);
01104 
01105     if (!patch->flags.doFullElectrostatics) NAMD_bug("Sequencer::multigratorPressure, doFullElectrostatics must be true");
01106 
01107     // Calculate new forces
01108     // NOTE: We should not need to migrate here since any migration should have happened in the
01109     // previous call to runComputeObjects inside the MD loop in Sequencer::integrate()
01110     const int numberOfSteps = simParams->N;
01111     const int stepsPerCycle = simParams->stepsPerCycle;
01112     runComputeObjects(0 , step<numberOfSteps, 1);
01113 
01114     reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
01115     reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
01116 
01117     // Virials etc.
01118     Tensor virialNormal;
01119     Tensor momentumSqrSum;
01120     BigReal kineticEnergy = 0;
01121     if ( simParams->pairInteractionOn ) {
01122       if ( simParams->pairInteractionSelf ) {
01123         for ( int i = 0; i < numAtoms; ++i ) {
01124           if ( a[i].partition != 1 ) continue;
01125           kineticEnergy += a[i].mass * a[i].velocity.length2();
01126           virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01127         }
01128       }
01129     } else {
01130       for ( int i = 0; i < numAtoms; ++i ) {
01131         if (a[i].mass < 0.01) continue;
01132         kineticEnergy += a[i].mass * a[i].velocity.length2();
01133         virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01134       }
01135     }
01136     if (!simParams->useGroupPressure) momentumSqrSum = virialNormal;
01137     kineticEnergy *= 0.5;
01138     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
01139     ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virialNormal);
01140 
01141     if ( simParams->fixedAtomsOn ) {
01142       Tensor fixVirialNormal;
01143       Tensor fixVirialNbond;
01144       Tensor fixVirialSlow;
01145       Vector fixForceNormal = 0;
01146       Vector fixForceNbond = 0;
01147       Vector fixForceSlow = 0;
01148 
01149       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
01150 
01151       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, fixVirialNormal);
01152       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, fixVirialNbond);
01153       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, fixVirialSlow);
01154       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NORMAL, fixForceNormal);
01155       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NBOND, fixForceNbond);
01156       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_SLOW, fixForceSlow);
01157     }
01158 
01159     // Internal virial and group momentum
01160     Tensor intVirialNormal;
01161     Tensor intVirialNormal2;
01162     Tensor intVirialNbond;
01163     Tensor intVirialSlow;
01164     int hgs;
01165     for ( int i = 0; i < numAtoms; i += hgs ) {
01166       hgs = a[i].hydrogenGroupSize;
01167       int j;
01168       BigReal m_cm = 0;
01169       Position x_cm(0,0,0);
01170       Velocity v_cm(0,0,0);
01171       for ( j = i; j < (i+hgs); ++j ) {
01172         m_cm += a[j].mass;
01173         x_cm += a[j].mass * a[j].position;
01174         v_cm += a[j].mass * a[j].velocity;
01175       }
01176       if (simParams->useGroupPressure) momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
01177       x_cm /= m_cm;
01178       v_cm /= m_cm;
01179       if (simParams->fixedAtomsOn) NAMD_bug("Sequencer::multigratorPressure, simParams->fixedAtomsOn not implemented yet");
01180       if ( simParams->pairInteractionOn ) {
01181         if ( simParams->pairInteractionSelf ) {
01182           NAMD_bug("Sequencer::multigratorPressure, this part needs to be implemented correctly");
01183           for ( j = i; j < (i+hgs); ++j ) {
01184             if ( a[j].partition != 1 ) continue;
01185             BigReal mass = a[j].mass;
01186             Vector v = a[j].velocity;
01187             Vector dv = v - v_cm;
01188             intVirialNormal2.outerAdd (mass, v, dv);
01189             Vector dx = a[j].position - x_cm;
01190             intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01191             intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01192             intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01193           }
01194         }
01195       } else {
01196         for ( j = i; j < (i+hgs); ++j ) {
01197           BigReal mass = a[j].mass;
01198           Vector v = a[j].velocity;
01199           Vector dv = v - v_cm;
01200           intVirialNormal2.outerAdd(mass, v, dv);
01201           Vector dx = a[j].position - x_cm;
01202           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
01203           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
01204           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
01205         }
01206       }
01207     }
01208 
01209     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
01210     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal2);
01211     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NBOND, intVirialNbond);
01212     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_SLOW, intVirialSlow);
01213     ADD_TENSOR_OBJECT(reduction, REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
01214 
01215     reduction->submit();
01216   }
01217 }

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

Definition at line 1247 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ResizeArray< Elem >::begin(), broadcast, calcKineticEnergy(), SimpleBroadcastObject< T >::get(), CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, FullAtom::mass, MULTIGRATOR_REDUCTION_KINETIC_ENERGY, 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().

01247                                                                {
01248   if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
01249     // Blocking receive (get) velocity scaling factor.
01250     BigReal velScale = (callNumber == 1) ? broadcast->velocityRescaleFactor.get(step) : broadcast->velocityRescaleFactor2.get(step);
01251     scaleVelocities(velScale);
01252     // Calculate new kineticEnergy
01253     BigReal kineticEnergy = calcKineticEnergy();
01254     multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY) += kineticEnergy;
01255     if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
01256       // If this is a pressure cycle, calculate new momentum squared sum
01257       FullAtom *a = patch->atom.begin();
01258       int numAtoms = patch->numAtoms;
01259       Tensor momentumSqrSum;
01260       if (simParams->useGroupPressure) {
01261         int hgs;
01262         for ( int i = 0; i < numAtoms; i += hgs ) {
01263           hgs = a[i].hydrogenGroupSize;
01264           int j;
01265           BigReal m_cm = 0;
01266           Position x_cm(0,0,0);
01267           Velocity v_cm(0,0,0);
01268           for ( j = i; j < (i+hgs); ++j ) {
01269             m_cm += a[j].mass;
01270             x_cm += a[j].mass * a[j].position;
01271             v_cm += a[j].mass * a[j].velocity;
01272           }
01273           momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
01274         }
01275       } else {
01276         for ( int i = 0; i < numAtoms; i++) {
01277           momentumSqrSum.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
01278         }
01279       }
01280       ADD_TENSOR_OBJECT(multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
01281     }
01282     // Submit reductions (kineticEnergy and, if applicable, momentumSqrSum)
01283     multigratorReduction->submit();
01284 
01285   }
01286 }

void Sequencer::newMinimizeDirection ( BigReal  c  )  [protected]

Definition at line 897 of file Sequencer.C.

References ResizeArray< Elem >::begin(), SimParameters::drudeHardWallOn, Patch::f, SimParameters::fixedAtomsOn, CompAtom::hydrogenGroupSize, 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().

00897                                               {
00898   FullAtom *a = patch->atom.begin();
00899   Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
00900   const bool fixedAtomsOn = simParams->fixedAtomsOn;
00901   const bool drudeHardWallOn = simParams->drudeHardWallOn;
00902   int numAtoms = patch->numAtoms;
00903   BigReal maxv2 = 0.;
00904 
00905   for ( int i = 0; i < numAtoms; ++i ) {
00906     a[i].velocity *= c;
00907     a[i].velocity += f1[i];
00908     if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
00909       a[i].velocity = a[i-1].velocity;
00910     }
00911     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00912     BigReal v2 = a[i].velocity.length2();
00913     if ( v2 > maxv2 ) maxv2 = v2;
00914   }
00915 
00916   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial); }
00917 
00918   maxv2 = 0.;
00919   for ( int i = 0; i < numAtoms; ++i ) {
00920     if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
00921       a[i].velocity = a[i-1].velocity;
00922     }
00923     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
00924     BigReal v2 = a[i].velocity.length2();
00925     if ( v2 > maxv2 ) maxv2 = v2;
00926   }
00927 
00928   min_reduction->max(0,maxv2);
00929   min_reduction->submit();
00930 
00931   // prevent hydrogens from being left behind
00932   BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
00933   // int adjustCount = 0;
00934   int hgs;
00935   for ( int i = 0; i < numAtoms; i += hgs ) {
00936     hgs = a[i].hydrogenGroupSize;
00937     BigReal minChildVel = a[i].velocity.length2();
00938     if ( minChildVel < fmax2 ) continue;
00939     int adjustChildren = 1;
00940     for ( int j = i+1; j < (i+hgs); ++j ) {
00941       if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
00942     }
00943     if ( adjustChildren ) {
00944       // if ( hgs > 1 ) ++adjustCount;
00945       for ( int j = i+1; j < (i+hgs); ++j ) {
00946         if (a[i].mass < 0.01) continue;  // lone pair
00947         a[j].velocity = a[i].velocity;
00948       }
00949     }
00950   }
00951   // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
00952 
00953 }

void Sequencer::newMinimizePosition ( BigReal  c  )  [protected]

Definition at line 956 of file Sequencer.C.

References ResizeArray< Elem >::begin(), SimParameters::drudeHardWallOn, Patch::numAtoms, patch, CompAtom::position, HomePatch::rattle1(), simParams, and FullAtom::velocity.

Referenced by minimize().

00956                                              {
00957   FullAtom *a = patch->atom.begin();
00958   int numAtoms = patch->numAtoms;
00959 
00960   for ( int i = 0; i < numAtoms; ++i ) {
00961     a[i].position += c * a[i].velocity;
00962   }
00963 
00964   if ( simParams->drudeHardWallOn ) {
00965     for ( int i = 1; i < numAtoms; ++i ) {
00966       if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
00967         a[i].position -= a[i-1].position;
00968       }
00969     }
00970   }
00971 
00972   patch->rattle1(0.,0,0);
00973 
00974   if ( simParams->drudeHardWallOn ) {
00975     for ( int i = 1; i < numAtoms; ++i ) {
00976       if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
00977         a[i].position += a[i-1].position;
00978       }
00979     }
00980   }
00981 }

void Sequencer::newtonianVelocities ( BigReal  stepscale,
const BigReal  timestep,
const BigReal  nbondstep,
const BigReal  slowstep,
const int  staleForces,
const int  doNonbonded,
const int  doFullElectrostatics 
) [protected]

Definition at line 1293 of file Sequencer.C.

References addForceToMomentum(), addForceToMomentum3(), Patch::flags, NAMD_EVENT_RANGE_2, Results::nbond, Results::normal, patch, and Results::slow.

Referenced by integrate().

01299 {
01300   NAMD_EVENT_RANGE_2(patch->flags.event_on,
01301       NamdProfileEvent::NEWTONIAN_VELOCITIES);
01302 
01303   // Deterministic velocity update, account for multigrator
01304   if (staleForces || (doNonbonded && doFullElectrostatics)) {
01305     addForceToMomentum3(stepscale*timestep, Results::normal, 0,
01306                         stepscale*nbondstep, Results::nbond, staleForces,
01307                         stepscale*slowstep, Results::slow, staleForces);
01308   } else {
01309     addForceToMomentum(stepscale*timestep);
01310     if (staleForces || doNonbonded)
01311       addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces);
01312     if (staleForces || doFullElectrostatics)
01313       addForceToMomentum(stepscale*slowstep, Results::slow, staleForces);
01314   }
01315 }

void Sequencer::quenchVelocities (  )  [protected]

Definition at line 984 of file Sequencer.C.

References ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity.

Referenced by minimize().

00984                                  {
00985   FullAtom *a = patch->atom.begin();
00986   int numAtoms = patch->numAtoms;
00987 
00988   for ( int i = 0; i < numAtoms; ++i ) {
00989     a[i].velocity = 0;
00990   }
00991 }

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

Definition at line 1970 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ResizeArray< Elem >::const_begin(), Node::enableEarlyExit(), endi(), Patch::f, Patch::flags, iERROR(), iout, NAMD_EVENT_RANGE_2, Results::normal, Patch::numAtoms, Node::Object(), patch, CompAtom::position, pressureProfileReduction, HomePatch::rattle1(), reduction, RIGID_NONE, SimParameters::rigidBonds, simParams, terminate(), FullAtom::velocity, Tensor::xx, Tensor::xy, Tensor::xz, Vector::y, Tensor::yx, Tensor::yy, Tensor::yz, Vector::z, Tensor::zx, Tensor::zy, and Tensor::zz.

Referenced by integrate(), and langevinVelocitiesBBK2().

01971 {
01972   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::RATTLE1);
01973   if ( simParams->rigidBonds != RIGID_NONE ) {
01974     Tensor virial;
01975     Tensor *vp = ( pressure ? &virial : 0 );
01976     if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
01977       iout << iERROR << 
01978         "Constraint failure; simulation has become unstable.\n" << endi;
01979       Node::Object()->enableEarlyExit();
01980       terminate();
01981     }
01982 #if 0
01983     printf("virial = %g %g %g  %g %g %g  %g %g %g\n",
01984         virial.xx, virial.xy, virial.xz,
01985         virial.yx, virial.yy, virial.yz,
01986         virial.zx, virial.zy, virial.zz);
01987 #endif
01988     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
01989 #if 0
01990     {
01991       const FullAtom *a = patch->atom.const_begin();
01992       for (int n=0;  n < patch->numAtoms;  n++) {
01993         printf("pos[%d] =  %g %g %g\n", n,
01994             a[n].position.x, a[n].position.y, a[n].position.z);
01995       }
01996       for (int n=0;  n < patch->numAtoms;  n++) {
01997         printf("vel[%d] =  %g %g %g\n", n,
01998             a[n].velocity.x, a[n].velocity.y, a[n].velocity.z);
01999       }
02000       if (pressure) {
02001         for (int n=0;  n < patch->numAtoms;  n++) {
02002           printf("force[%d] =  %g %g %g\n", n,
02003               patch->f[Results::normal][n].x,
02004               patch->f[Results::normal][n].y,
02005               patch->f[Results::normal][n].z);
02006         }
02007       }
02008     }
02009 #endif
02010   }
02011 }

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

Definition at line 1756 of file Sequencer.C.

References 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, simParams, and FullAtom::velocity.

Referenced by integrate().

01757 {
01758   const int reassignFreq = simParams->reassignFreq;
01759   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
01760     FullAtom *a = patch->atom.begin();
01761     int numAtoms = patch->numAtoms;
01762     BigReal newTemp = simParams->reassignTemp;
01763     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
01764     if ( simParams->reassignIncr > 0.0 ) {
01765       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
01766         newTemp = simParams->reassignHold;
01767     } else {
01768       if ( newTemp < simParams->reassignHold )
01769         newTemp = simParams->reassignHold;
01770     }
01771     BigReal kbT = BOLTZMANN * newTemp;
01772 
01773     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01774     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01775 
01776     for ( int i = 0; i < numAtoms; ++i )
01777     {
01778       a[i].velocity = ( ( simParams->fixedAtomsOn && 
01779             a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
01780           sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) * 
01781           random->gaussian_vector() );
01782     }
01783   } else {
01784     NAMD_bug("Sequencer::reassignVelocities called improperly!");
01785   }
01786 }

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

Definition at line 2867 of file Sequencer.C.

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

Referenced by integrate(), and minimize().

02867                                           {
02868   if ( ! ldbSteps ) {
02869     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
02870   }
02871   if ( ! --ldbSteps ) {
02872     patch->submitLoadStats(timestep);
02873     ldbCoordinator->rebalance(this,patch->getPatchID());
02874     pairlistsAreValid = 0;
02875   }
02876 }

void Sequencer::reinitVelocities ( void   )  [protected]

Definition at line 1788 of file Sequencer.C.

References 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().

01789 {
01790   FullAtom *a = patch->atom.begin();
01791   int numAtoms = patch->numAtoms;
01792   BigReal newTemp = simParams->initialTemp;
01793   BigReal kbT = BOLTZMANN * newTemp;
01794 
01795   int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
01796   BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
01797 
01798   for ( int i = 0; i < numAtoms; ++i )
01799   {
01800     a[i].velocity = ( ( (simParams->fixedAtomsOn && a[i].atomFixed) || 
01801           a[i].mass <= 0.) ? Vector(0,0,0) :
01802         sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) * 
01803         random->gaussian_vector() );
01804     if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass > 0.05 ) {
01805       a[i+1].velocity = a[i].velocity;  // zero is good enough
01806       ++i;
01807     }
01808   }
01809 }

void Sequencer::reloadCharges (  )  [protected]

Definition at line 1821 of file Sequencer.C.

References Molecule::atomcharge(), ResizeArray< Elem >::begin(), CompAtom::charge, Node::molecule, Patch::numAtoms, Node::Object(), and patch.

Referenced by algorithm().

01822 {
01823   FullAtom *a = patch->atom.begin();
01824   int numAtoms = patch->numAtoms;
01825   Molecule *molecule = Node::Object()->molecule;
01826   for ( int i = 0; i < numAtoms; ++i )
01827   {
01828     a[i].charge = molecule->atomcharge(a[i].id);
01829   }
01830 }

void Sequencer::rescaleaccelMD ( int  step,
int  doNonbonded,
int  doFullElectrostatics 
) [protected]

Definition at line 1697 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().

01698 {
01699    if (!simParams->accelMDOn) return;
01700    if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
01701 
01702    // Blocking receive for the Accelerated MD scaling factors.
01703    Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
01704    const BigReal factor_dihe = accelMDfactor[0];
01705    const BigReal factor_tot  = accelMDfactor[1];
01706    const int numAtoms = patch->numAtoms;
01707 
01708    if (simParams->accelMDdihe && factor_tot <1 )
01709        NAMD_die("accelMD broadcasting error!\n");
01710    if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
01711        NAMD_die("accelMD broadcasting error!\n");
01712 
01713    if (simParams->accelMDdihe && factor_dihe < 1) {
01714         for (int i = 0; i < numAtoms; ++i)
01715           if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01716               patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);         
01717    }
01718 
01719    if ( !simParams->accelMDdihe && factor_tot < 1) {
01720         for (int i = 0; i < numAtoms; ++i)
01721             patch->f[Results::normal][i] *= factor_tot;
01722         if (doNonbonded) {
01723             for (int i = 0; i < numAtoms; ++i)
01724                  patch->f[Results::nbond][i] *= factor_tot;
01725         }
01726         if (doFullElectrostatics) {
01727             for (int i = 0; i < numAtoms; ++i)
01728                  patch->f[Results::slow][i] *= factor_tot;
01729         }
01730    }
01731 
01732    if (simParams->accelMDdual && factor_dihe < 1) {
01733         for (int i = 0; i < numAtoms; ++i)
01734            if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
01735                patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
01736    }
01737 
01738 }

void Sequencer::rescaleSoluteCharges ( BigReal  factor  )  [protected]

Definition at line 1833 of file Sequencer.C.

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

Referenced by algorithm(), and Sequencer().

01834 {
01835   FullAtom *a = patch->atom.begin();
01836   int numAtoms = patch->numAtoms;
01837   Molecule *molecule = Node::Object()->molecule;
01838   BigReal sqrt_factor = sqrt(factor);
01839   // apply scaling to the original charge (stored in molecule)
01840   // of just the marked solute atoms
01841   for ( int i = 0; i < numAtoms; ++i ) {
01842     if (molecule->get_ss_type(a[i].id)) {
01843       a[i].charge = sqrt_factor * molecule->atomcharge(a[i].id);
01844     }
01845   }
01846 }

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

Definition at line 1678 of file Sequencer.C.

References ResizeArray< Elem >::begin(), broadcast, SimpleBroadcastObject< T >::get(), Patch::numAtoms, patch, SimParameters::rescaleFreq, rescaleVelocities_numTemps, simParams, FullAtom::velocity, and ControllerBroadcasts::velocityRescaleFactor.

Referenced by integrate().

01679 {
01680   const int rescaleFreq = simParams->rescaleFreq;
01681   if ( rescaleFreq > 0 ) {
01682     FullAtom *a = patch->atom.begin();
01683     int numAtoms = patch->numAtoms;
01684     ++rescaleVelocities_numTemps;
01685     if ( rescaleVelocities_numTemps == rescaleFreq ) {
01686       // Blocking receive for the velcity scaling factor.
01687       BigReal factor = broadcast->velocityRescaleFactor.get(step);
01688       for ( int i = 0; i < numAtoms; ++i )
01689       {
01690         a[i].velocity *= factor;
01691       }
01692       rescaleVelocities_numTemps = 0;
01693     }
01694   }
01695 }

void Sequencer::rescaleVelocitiesByFactor ( BigReal  factor  )  [protected]

Definition at line 1811 of file Sequencer.C.

References ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity.

Referenced by algorithm().

01812 {
01813   FullAtom *a = patch->atom.begin();
01814   int numAtoms = patch->numAtoms;
01815   for ( int i = 0; i < numAtoms; ++i )
01816   {
01817     a[i].velocity *= factor;
01818   }
01819 }

void Sequencer::run ( void   ) 

Definition at line 126 of file Sequencer.C.

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

Referenced by HomePatch::runSequencer().

00127 {
00128     // create a Thread and invoke it
00129     DebugM(4, "::run() - this = " << this << "\n" );
00130     thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
00131     CthSetStrategyDefault(thread);
00132     priority = PATCH_PRIORITY(patch->getPatchID());
00133     awaken();
00134 }

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

Definition at line 2668 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().

02669 {
02670   if ( migration ) pairlistsAreValid = 0;
02671 #if defined(NAMD_CUDA) || defined(NAMD_MIC)
02672   if ( pairlistsAreValid &&
02673        ( patch->flags.doFullElectrostatics || ! simParams->fullElectFrequency )
02674                          && ( pairlistsAge > (
02675 #else
02676   if ( pairlistsAreValid && ( pairlistsAge > (
02677 #endif
02678          (simParams->stepsPerCycle - 1) / simParams->pairlistsPerCycle ) ) ) {
02679     pairlistsAreValid = 0;
02680   }
02681   if ( ! simParams->usePairlists ) pairlists = 0;
02682   patch->flags.usePairlists = pairlists || pairlistsAreValid;
02683   patch->flags.savePairlists =
02684         pairlists && ! pairlistsAreValid;
02685 
02686   if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
02687 
02688   //
02689   // DJH: Copy updates of SOA back into AOS.
02690   // The positionsReady() routine starts force computation and atom migration.
02691   //
02692   // We could reduce amount of copying here by checking migration status
02693   // and copying velocities only when migrating. Some types of simulation
02694   // always require velocities, such as Lowe-Anderson.
02695   //
02696   //patch->copy_updates_to_AOS();
02697 
02698   patch->positionsReady(migration);  // updates flags.sequence
02699 
02700   int seq = patch->flags.sequence;
02701   int basePriority = ( (seq & 0xffff) << 15 )
02702                      + PATCH_PRIORITY(patch->getPatchID());
02703   if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
02704     priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
02705     suspend(); // until all deposit boxes close
02706     patch->gbisComputeAfterP1();
02707     priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
02708     suspend();
02709     patch->gbisComputeAfterP2();
02710     priority = basePriority + COMPUTE_HOME_PRIORITY;
02711     suspend();
02712   } else {
02713     priority = basePriority + COMPUTE_HOME_PRIORITY;
02714     suspend(); // until all deposit boxes close
02715   }
02716 
02717   //
02718   // DJH: Copy all data into SOA from AOS.
02719   //
02720   // We need everything copied after atom migration.
02721   // When doing force computation without atom migration,
02722   // all data except forces will already be up-to-date in SOA
02723   // (except maybe for some special types of simulation).
02724   //
02725   //patch->copy_all_to_SOA();
02726 
02727   //
02728   // DJH: Copy forces to SOA.
02729   // Force available after suspend() has returned.
02730   //
02731   //patch->copy_forces_to_SOA();
02732 
02733   if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
02734     pairlistsAreValid = 1;
02735     pairlistsAge = 0;
02736   }
02737   // For multigrator, do not age pairlist during pressure step
02738   // NOTE: for non-multigrator pressureStep = 0 always
02739   if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
02740 
02741   if (simParams->lonepairs) {
02742     {
02743       Tensor virial;
02744       patch->redistrib_lonepair_forces(Results::normal, &virial);
02745       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02746     }
02747     if (patch->flags.doNonbonded) {
02748       Tensor virial;
02749       patch->redistrib_lonepair_forces(Results::nbond, &virial);
02750       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02751     }
02752     if (patch->flags.doFullElectrostatics) {
02753       Tensor virial;
02754       patch->redistrib_lonepair_forces(Results::slow, &virial);
02755       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02756     }
02757   } else if (simParams->watmodel == WAT_TIP4) {
02758     {
02759       Tensor virial;
02760       patch->redistrib_tip4p_forces(Results::normal, &virial);
02761       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02762     }
02763     if (patch->flags.doNonbonded) {
02764       Tensor virial;
02765       patch->redistrib_tip4p_forces(Results::nbond, &virial);
02766       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02767     }
02768     if (patch->flags.doFullElectrostatics) {
02769       Tensor virial;
02770       patch->redistrib_tip4p_forces(Results::slow, &virial);
02771       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02772     }
02773   } else if (simParams->watmodel == WAT_SWM4) {
02774     {
02775       Tensor virial;
02776       patch->redistrib_swm4_forces(Results::normal, &virial);
02777       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
02778     }
02779     if (patch->flags.doNonbonded) {
02780       Tensor virial;
02781       patch->redistrib_swm4_forces(Results::nbond, &virial);
02782       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
02783     }
02784     if (patch->flags.doFullElectrostatics) {
02785       Tensor virial;
02786       patch->redistrib_swm4_forces(Results::slow, &virial);
02787       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
02788     }
02789   }
02790 
02791   if ( patch->flags.doMolly ) {
02792     Tensor virial;
02793     patch->mollyMollify(&virial);
02794     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
02795   }
02796 
02797 
02798   // BEGIN LA
02799   if (patch->flags.doLoweAndersen) {
02800       patch->loweAndersenFinish();
02801   }
02802   // END LA
02803 #ifdef NAMD_CUDA_XXX
02804   int numAtoms = patch->numAtoms;
02805   FullAtom *a = patch->atom.begin();
02806   for ( int i=0; i<numAtoms; ++i ) {
02807     CkPrintf("%d %g %g %g\n", a[i].id,
02808         patch->f[Results::normal][i].x +
02809         patch->f[Results::nbond][i].x +
02810         patch->f[Results::slow][i].x,
02811         patch->f[Results::normal][i].y + 
02812         patch->f[Results::nbond][i].y +
02813         patch->f[Results::slow][i].y,
02814         patch->f[Results::normal][i].z +
02815         patch->f[Results::nbond][i].z +
02816         patch->f[Results::slow][i].z);
02817     CkPrintf("%d %g %g %g\n", a[i].id,
02818         patch->f[Results::normal][i].x,
02819         patch->f[Results::nbond][i].x,
02820         patch->f[Results::slow][i].x);
02821     CkPrintf("%d %g %g %g\n", a[i].id,
02822         patch->f[Results::normal][i].y,
02823         patch->f[Results::nbond][i].y,
02824         patch->f[Results::slow][i].y);
02825     CkPrintf("%d %g %g %g\n", a[i].id,
02826         patch->f[Results::normal][i].z,
02827         patch->f[Results::nbond][i].z,
02828         patch->f[Results::slow][i].z);
02829   }
02830 #endif
02831 
02832 #if PRINT_FORCES
02833   int numAtoms = patch->numAtoms;
02834   FullAtom *a = patch->atom.begin();
02835   for ( int i=0; i<numAtoms; ++i ) {
02836     float fxNo = patch->f[Results::normal][i].x;
02837     float fxNb = patch->f[Results::nbond][i].x;
02838     float fxSl = patch->f[Results::slow][i].x;
02839     float fyNo = patch->f[Results::normal][i].y;
02840     float fyNb = patch->f[Results::nbond][i].y;
02841     float fySl = patch->f[Results::slow][i].y;
02842     float fzNo = patch->f[Results::normal][i].z;
02843     float fzNb = patch->f[Results::nbond][i].z;
02844     float fzSl = patch->f[Results::slow][i].z;
02845     float fx = fxNo+fxNb+fxSl;
02846     float fy = fyNo+fyNb+fySl;
02847     float fz = fzNo+fzNb+fzSl;
02848 
02849                 float f = sqrt(fx*fx+fy*fy+fz*fz);
02850     int id = patch->pExt[i].id;
02851     int seq = patch->flags.sequence;
02852     float x = patch->p[i].position.x;
02853     float y = patch->p[i].position.y;
02854     float z = patch->p[i].position.z;
02855     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
02856     CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
02857     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
02858 //fxNo,fyNo,fzNo,
02859 //fxNb,fyNb,fzNb,
02860 //fxSl,fySl,fzSl,
02861 fx,fy,fz
02862 );
02863         }
02864 #endif
02865 }

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

Definition at line 1893 of file Sequencer.C.

References patch, and HomePatch::saveForce().

Referenced by integrate(), and minimize().

01894 {
01895   patch->saveForce(ftag);
01896 }

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

Definition at line 1047 of file Sequencer.C.

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

Referenced by multigratorPressure().

01047                                                                                        {
01048   FullAtom *a = patch->atom.begin();
01049   int numAtoms = patch->numAtoms;
01050   Position origin = patch->lattice.origin();
01051   if ( simParams->fixedAtomsOn ) {
01052     NAMD_bug("Sequencer::scalePositionsVelocities, fixed atoms not implemented");
01053   }
01054   if ( simParams->useGroupPressure ) {
01055     int hgs;
01056     for ( int i = 0; i < numAtoms; i += hgs ) {
01057       hgs = a[i].hydrogenGroupSize;
01058       Position pos_cm(0.0, 0.0, 0.0);
01059       Velocity vel_cm(0.0, 0.0, 0.0);
01060       BigReal m_cm = 0.0;
01061       for (int j=0;j < hgs;++j) {
01062         m_cm += a[i+j].mass;
01063         pos_cm += a[i+j].mass*a[i+j].position;
01064         vel_cm += a[i+j].mass*a[i+j].velocity;
01065       }
01066       pos_cm /= m_cm;
01067       vel_cm /= m_cm;
01068       pos_cm -= origin;
01069       Position dpos = posScale*pos_cm;
01070       Velocity dvel = velScale*vel_cm;
01071       for (int j=0;j < hgs;++j) {
01072         a[i+j].position += dpos;
01073         a[i+j].velocity += dvel;
01074       }
01075     }
01076   } else {
01077     for ( int i = 0; i < numAtoms; i++) {
01078       a[i].position += posScale*(a[i].position-origin);
01079       a[i].velocity = velScale*a[i].velocity;
01080     }
01081   }
01082 }

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

Definition at line 1219 of file Sequencer.C.

References ResizeArray< Elem >::begin(), Patch::numAtoms, patch, and FullAtom::velocity.

Referenced by multigratorTemperature().

01219                                                       {
01220   FullAtom *a = patch->atom.begin();
01221   int numAtoms = patch->numAtoms;
01222   for ( int i = 0; i < numAtoms; i++) {
01223     a[i].velocity *= velScale;
01224   }
01225 }

void Sequencer::stochRescaleVelocities ( BigReal  dt_fs,
int  step 
) [protected]

When doing stochastic velocity rescaling, every stochRescaleFreq steps we receive the globally computed rescaling coefficient and apply it to the velocities of all the atoms in our patch.

Rescale velocities with the scale factor sent from the Controller.

Parameters:
dt_fs The integration increment, in fs (not currently used)
step The current timestep

Definition at line 1872 of file Sequencer.C.

References ResizeArray< Elem >::begin(), broadcast, DebugM, SimpleBroadcastObject< T >::get(), Patch::numAtoms, patch, simParams, stochRescale_count, ControllerBroadcasts::stochRescaleCoefficient, SimParameters::stochRescaleFreq, SimParameters::stochRescaleOn, and FullAtom::velocity.

Referenced by integrate().

01873 {
01874   if ( simParams->stochRescaleOn )
01875   {
01876     FullAtom *a = patch->atom.begin();
01877     int numAtoms = patch->numAtoms;
01878     ++stochRescale_count;
01879     if ( stochRescale_count == simParams->stochRescaleFreq )
01880     {
01881       // Blocking receive for the temperature coupling coefficient.
01882       BigReal coefficient = broadcast->stochRescaleCoefficient.get(step);
01883       DebugM(4, "stochastically rescaling velocities at step " << step << " by " << coefficient << "\n");
01884       for ( int i = 0; i < numAtoms; ++i )
01885       {
01886         a[i].velocity *= coefficient;
01887       }
01888       stochRescale_count = 0;
01889     }
01890   }
01891 }

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

Definition at line 2643 of file Sequencer.C.

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

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

02644 {
02645   //
02646   // DJH: Copy updates of SOA back into AOS.
02647   // Do we need to update everything or is it safe to just update
02648   // positions and velocities separately, as needed?
02649   //
02650   //patch->copy_updates_to_AOS();
02651 
02652   NAMD_EVENT_RANGE_2(patch->flags.event_on,
02653       NamdProfileEvent::SUBMIT_COLLECTIONS);
02654   int prec = Output::coordinateNeeded(step);
02655   if ( prec ) {
02656     collection->submitPositions(step,patch->atom,patch->lattice,prec);
02657   }
02658   if ( Output::velocityNeeded(step) ) {
02659     collection->submitVelocities(step,zeroVel,patch->atom);
02660   }
02661   if ( Output::forceNeeded(step) ) {
02662     int maxForceUsed = patch->flags.maxForceUsed;
02663     if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
02664     collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
02665   }
02666 }

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

Definition at line 2083 of file Sequencer.C.

References ADD_TENSOR_OBJECT, 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, NAMD_EVENT_RANGE_2, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, CompAtom::partition, partition(), patch, CompAtom::position, pressureProfileReduction, SimParameters::pressureProfileSlabs, reduction, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, simParams, SimParameters::useGroupPressure, FullAtom::velocity, Vector::x, Vector::y, and Vector::z.

Referenced by integrate().

02084 {
02085   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::SUBMIT_HALFSTEP);
02086 
02087   // velocity-dependent quantities *** ONLY ***
02088   // positions are not at half-step when called
02089   FullAtom *a = patch->atom.begin();
02090   int numAtoms = patch->numAtoms;
02091 
02092 #if CMK_BLUEGENEL
02093   CmiNetworkProgressAfter (0);
02094 #endif
02095 
02096   // For non-Multigrator doKineticEnergy = 1 always
02097   Tensor momentumSqrSum;
02098   if (doKineticEnergy || patch->flags.doVirial)
02099   {
02100     BigReal kineticEnergy = 0;
02101     Tensor virial;
02102     if ( simParams->pairInteractionOn ) {
02103       if ( simParams->pairInteractionSelf ) {
02104         for ( int i = 0; i < numAtoms; ++i ) {
02105           if ( a[i].partition != 1 ) continue;
02106           kineticEnergy += a[i].mass * a[i].velocity.length2();
02107           virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
02108         }
02109       }
02110     } else {
02111       for ( int i = 0; i < numAtoms; ++i ) {
02112         if (a[i].mass < 0.01) continue;
02113         kineticEnergy += a[i].mass * a[i].velocity.length2();
02114         virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
02115       }
02116     }
02117     
02118     if (simParams->multigratorOn && !simParams->useGroupPressure) {
02119       momentumSqrSum = virial;
02120     }
02121     kineticEnergy *= 0.5 * 0.5;
02122     reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
02123     virial *= 0.5;
02124     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
02125 #ifdef ALTVIRIAL
02126     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
02127 #endif
02128   }
02129  
02130   if (pressureProfileReduction) {
02131     int nslabs = simParams->pressureProfileSlabs;
02132     const Lattice &lattice = patch->lattice;
02133     BigReal idz = nslabs/lattice.c().z;
02134     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
02135     int useGroupPressure = simParams->useGroupPressure;
02136 
02137     // Compute kinetic energy partition, possibly subtracting off
02138     // internal kinetic energy if group pressure is enabled.
02139     // Since the regular pressure is 1/2 mvv and the internal kinetic
02140     // term that is subtracted off for the group pressure is
02141     // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
02142     // 1/2 m * v * v_cm.  The factor of 1/2 is because submitHalfstep
02143     // gets called twice per timestep.
02144     int hgs;
02145     for (int i=0; i<numAtoms; i += hgs) {
02146       int j, ppoffset;
02147       hgs = a[i].hydrogenGroupSize;
02148       int partition = a[i].partition;
02149 
02150       BigReal m_cm = 0;
02151       Velocity v_cm(0,0,0);
02152       for (j=i; j< i+hgs; ++j) {
02153         m_cm += a[j].mass;
02154         v_cm += a[j].mass * a[j].velocity;
02155       }
02156       v_cm /= m_cm;
02157       for (j=i; j < i+hgs; ++j) {
02158         BigReal mass = a[j].mass;
02159         if (! (useGroupPressure && j != i)) {
02160           BigReal z = a[j].position.z;
02161           int slab = (int)floor((z-zmin)*idz);
02162           if (slab < 0) slab += nslabs;
02163           else if (slab >= nslabs) slab -= nslabs;
02164           ppoffset = 3*(slab + partition*nslabs);
02165         }
02166         BigReal wxx, wyy, wzz;
02167         if (useGroupPressure) {
02168           wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
02169           wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
02170           wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
02171         } else {
02172           wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
02173           wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
02174           wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
02175         }
02176         pressureProfileReduction->item(ppoffset  ) += wxx;
02177         pressureProfileReduction->item(ppoffset+1) += wyy;
02178         pressureProfileReduction->item(ppoffset+2) += wzz;
02179       }
02180     }
02181   } 
02182 
02183   // For non-Multigrator doKineticEnergy = 1 always
02184   if (doKineticEnergy || patch->flags.doVirial)
02185   {
02186     BigReal intKineticEnergy = 0;
02187     Tensor intVirialNormal;
02188 
02189     int hgs;
02190     for ( int i = 0; i < numAtoms; i += hgs ) {
02191 
02192 #if CMK_BLUEGENEL
02193       CmiNetworkProgress ();
02194 #endif
02195 
02196       hgs = a[i].hydrogenGroupSize;
02197       int j;
02198       BigReal m_cm = 0;
02199       Velocity v_cm(0,0,0);
02200       for ( j = i; j < (i+hgs); ++j ) {
02201         m_cm += a[j].mass;
02202         v_cm += a[j].mass * a[j].velocity;
02203       }
02204       if (simParams->multigratorOn && simParams->useGroupPressure) {
02205         momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
02206       }
02207       v_cm /= m_cm;
02208       if ( simParams->pairInteractionOn ) {
02209         if ( simParams->pairInteractionSelf ) {
02210           for ( j = i; j < (i+hgs); ++j ) {
02211             if ( a[j].partition != 1 ) continue;
02212             BigReal mass = a[j].mass;
02213             Vector v = a[j].velocity;
02214             Vector dv = v - v_cm;
02215             intKineticEnergy += mass * (v * dv);
02216             intVirialNormal.outerAdd (mass, v, dv);
02217           }
02218         }
02219       } else {
02220         for ( j = i; j < (i+hgs); ++j ) {
02221           BigReal mass = a[j].mass;
02222           Vector v = a[j].velocity;
02223           Vector dv = v - v_cm;
02224           intKineticEnergy += mass * (v * dv);
02225           intVirialNormal.outerAdd(mass, v, dv);
02226         }
02227       }
02228     }
02229 
02230     intKineticEnergy *= 0.5 * 0.5;
02231     reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
02232     intVirialNormal *= 0.5;
02233     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02234     if ( simParams->multigratorOn) {
02235       momentumSqrSum *= 0.5;
02236       ADD_TENSOR_OBJECT(reduction,REDUCTION_MOMENTUM_SQUARED,momentumSqrSum);
02237     }
02238   }
02239 
02240 }

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

Definition at line 2505 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), calcFixVirial(), SimParameters::drudeBondLen, SimParameters::drudeHardWallOn, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, Vector::length2(), FullAtom::mass, HomePatch::minimize_rattle2(), Results::nbond, Results::normal, Patch::numAtoms, Tensor::outerAdd(), patch, Patch::pExt, CompAtom::position, SimParameters::printBadContacts, reduction, REDUCTION_ATOM_CHECKSUM, REDUCTION_MIN_F_DOT_F, REDUCTION_MIN_F_DOT_V, REDUCTION_MIN_HUGE_COUNT, REDUCTION_MIN_V_DOT_V, Flags::sequence, simParams, Results::slow, SubmitReduction::submit(), TIMEFACTOR, Vector::unit(), and FullAtom::velocity.

Referenced by minimize().

02506 {
02507   FullAtom *a = patch->atom.begin();
02508   Force *f1 = patch->f[Results::normal].begin();
02509   Force *f2 = patch->f[Results::nbond].begin();
02510   Force *f3 = patch->f[Results::slow].begin();
02511   const bool fixedAtomsOn = simParams->fixedAtomsOn;
02512   const bool drudeHardWallOn = simParams->drudeHardWallOn;
02513   const double drudeBondLen = simParams->drudeBondLen;
02514   const double drudeBondLen2 = drudeBondLen * drudeBondLen;
02515   const double drudeStep = 0.1/(TIMEFACTOR*TIMEFACTOR);
02516   const double drudeMove = 0.01;
02517   const double drudeStep2 = drudeStep * drudeStep;
02518   const double drudeMove2 = drudeMove * drudeMove;
02519   int numAtoms = patch->numAtoms;
02520 
02521   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
02522 
02523   for ( int i = 0; i < numAtoms; ++i ) {
02524     f1[i] += f2[i] + f3[i];  // add all forces
02525     if ( drudeHardWallOn && i && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) { // drude particle
02526       if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
02527         if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
02528           a[i].position += drudeMove * f1[i].unit();
02529         } else {
02530           a[i].position += drudeStep * f1[i];
02531         }
02532         if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
02533           a[i].position = a[i-1].position + drudeBondLen * (a[i].position - a[i-1].position).unit();
02534         }
02535       }
02536       Vector netf = f1[i-1] + f1[i];
02537       if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
02538       f1[i-1] = netf;
02539       f1[i] = 0.;
02540     }
02541     if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
02542   }
02543 
02544   f2 = f3 = 0;  // included in f1
02545 
02546   BigReal maxv2 = 0.;
02547 
02548   for ( int i = 0; i < numAtoms; ++i ) {
02549     BigReal v2 = a[i].velocity.length2();
02550     if ( v2 > 0. ) {
02551       if ( v2 > maxv2 ) maxv2 = v2;
02552     } else {
02553       v2 = f1[i].length2();
02554       if ( v2 > maxv2 ) maxv2 = v2;
02555     }
02556   }
02557 
02558   if ( fmax2 > 10. * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR )
02559   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial, true /* forces */); }
02560 
02561   BigReal fdotf = 0;
02562   BigReal fdotv = 0;
02563   BigReal vdotv = 0;
02564   int numHuge = 0;
02565   for ( int i = 0; i < numAtoms; ++i ) {
02566     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
02567     if ( drudeHardWallOn && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) continue; // drude particle
02568     Force f = f1[i];
02569     BigReal ff = f * f;
02570     if ( ff > fmax2 ) {
02571       if (simParams->printBadContacts) {
02572         CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
02573       }
02574       ++numHuge;
02575       // pad scaling so minimizeMoveDownhill() doesn't miss them
02576       BigReal fmult = 1.01 * sqrt(fmax2/ff);
02577       f *= fmult;  ff = f * f;
02578       f1[i] *= fmult;
02579     }
02580     fdotf += ff;
02581     fdotv += f * a[i].velocity;
02582     vdotv += a[i].velocity * a[i].velocity;
02583   }
02584 
02585   reduction->item(REDUCTION_MIN_F_DOT_F) += fdotf;
02586   reduction->item(REDUCTION_MIN_F_DOT_V) += fdotv;
02587   reduction->item(REDUCTION_MIN_V_DOT_V) += vdotv;
02588   reduction->item(REDUCTION_MIN_HUGE_COUNT) += numHuge;
02589 
02590   {
02591     Tensor intVirialNormal;
02592     Tensor intVirialNbond;
02593     Tensor intVirialSlow;
02594 
02595     int hgs;
02596     for ( int i = 0; i < numAtoms; i += hgs ) {
02597       hgs = a[i].hydrogenGroupSize;
02598       int j;
02599       BigReal m_cm = 0;
02600       Position x_cm(0,0,0);
02601       for ( j = i; j < (i+hgs); ++j ) {
02602         m_cm += a[j].mass;
02603         x_cm += a[j].mass * a[j].position;
02604       }
02605       x_cm /= m_cm;
02606       for ( j = i; j < (i+hgs); ++j ) {
02607         BigReal mass = a[j].mass;
02608         // net force treated as zero for fixed atoms
02609         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
02610         Vector dx = a[j].position - x_cm;
02611         intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02612         intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02613         intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02614       }
02615     }
02616 
02617     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02618     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
02619     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
02620   }
02621 
02622   if ( simParams->fixedAtomsOn ) {
02623     Tensor fixVirialNormal;
02624     Tensor fixVirialNbond;
02625     Tensor fixVirialSlow;
02626     Vector fixForceNormal = 0;
02627     Vector fixForceNbond = 0;
02628     Vector fixForceSlow = 0;
02629 
02630     calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
02631 
02632     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
02633     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
02634     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
02635     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
02636     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
02637     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
02638   }
02639 
02640   reduction->submit();
02641 }

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

Definition at line 993 of file Sequencer.C.

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

Referenced by integrate().

00993                                        {
00994 
00995   FullAtom *a = patch->atom.begin();
00996   const int numAtoms = patch->numAtoms;
00997 
00998   Vector momentum = 0;
00999   BigReal mass = 0;
01000 if ( simParams->zeroMomentumAlt ) {
01001   for ( int i = 0; i < numAtoms; ++i ) {
01002     momentum += a[i].mass * a[i].velocity;
01003     mass += 1.;
01004   }
01005 } else {
01006   for ( int i = 0; i < numAtoms; ++i ) {
01007     momentum += a[i].mass * a[i].velocity;
01008     mass += a[i].mass;
01009   }
01010 }
01011 
01012   ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
01013   reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
01014 }

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

Definition at line 2262 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, ResizeArray< Elem >::begin(), Lattice::c(), calcFixVirial(), cross(), doKineticEnergy, doMomenta, Flags::doVirial, SimParameters::drudeOn, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), j, GlobalMaster::lattice, Patch::lattice, Vector::length2(), HomePatch::marginViolations, FullAtom::mass, NAMD_EVENT_RANGE_2, Results::nbond, Results::normal, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, CompAtom::partition, partition(), patch, CompAtom::position, pressureProfileReduction, SimParameters::pressureProfileSlabs, reduction, REDUCTION_ATOM_CHECKSUM, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY, REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_MARGIN_VIOLATIONS, simParams, Results::slow, SubmitReduction::submit(), SimParameters::useGroupPressure, FullAtom::velocity, Vector::x, Vector::y, and Vector::z.

Referenced by integrate().

02263 {
02264 #ifndef UPPER_BOUND
02265   NAMD_EVENT_RANGE_2(patch->flags.event_on,
02266       NamdProfileEvent::SUBMIT_REDUCTIONS);
02267   FullAtom *a = patch->atom.begin();
02268 #endif
02269   int numAtoms = patch->numAtoms;
02270 
02271 #if CMK_BLUEGENEL
02272   CmiNetworkProgressAfter(0);
02273 #endif
02274 
02275   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
02276   reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
02277 
02278 #ifndef UPPER_BOUND
02279   // For non-Multigrator doKineticEnergy = 1 always
02280   if (doKineticEnergy || doMomenta || patch->flags.doVirial)
02281   {
02282     BigReal kineticEnergy = 0;
02283     Vector momentum = 0;
02284     Vector angularMomentum = 0;
02285     Vector o = patch->lattice.origin();
02286     int i;
02287     if ( simParams->pairInteractionOn ) {
02288       if ( simParams->pairInteractionSelf ) {
02289         for (i = 0; i < numAtoms; ++i ) {
02290           if ( a[i].partition != 1 ) continue;
02291           kineticEnergy += a[i].mass * a[i].velocity.length2();
02292           momentum += a[i].mass * a[i].velocity;
02293           angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
02294         }
02295       }
02296     } else {
02297       for (i = 0; i < numAtoms; ++i ) {
02298         kineticEnergy += a[i].mass * a[i].velocity.length2();
02299         momentum += a[i].mass * a[i].velocity;
02300         angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
02301       }
02302       if (simParams->drudeOn) {
02303         BigReal drudeComKE = 0.;
02304         BigReal drudeBondKE = 0.;
02305 
02306         for (i = 0;  i < numAtoms;  i++) {
02307           if (i < numAtoms-1 &&
02308               a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
02309             // i+1 is a Drude particle with parent i
02310 
02311             // convert from Cartesian coordinates to (COM,bond) coordinates
02312             BigReal m_com = (a[i].mass + a[i+1].mass);  // mass of COM
02313             BigReal m = a[i+1].mass / m_com;  // mass ratio
02314             BigReal m_bond = a[i+1].mass * (1. - m);  // mass of bond
02315             Vector v_bond = a[i+1].velocity - a[i].velocity;  // vel of bond
02316             Vector v_com = a[i].velocity + m * v_bond;  // vel of COM
02317 
02318             drudeComKE += m_com * v_com.length2();
02319             drudeBondKE += m_bond * v_bond.length2();
02320 
02321             i++;  // +1 from loop, we've updated both particles
02322           }
02323           else {
02324             drudeComKE += a[i].mass * a[i].velocity.length2();
02325           }
02326         } // end for
02327 
02328         drudeComKE *= 0.5;
02329         drudeBondKE *= 0.5;
02330         reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY)
02331           += drudeComKE;
02332         reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY)
02333           += drudeBondKE;
02334       } // end drudeOn
02335 
02336     } // end else
02337 
02338     kineticEnergy *= 0.5;
02339     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
02340     ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
02341     ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);  
02342   }
02343 
02344 #ifdef ALTVIRIAL
02345   // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
02346   {
02347     Tensor altVirial;
02348     for ( int i = 0; i < numAtoms; ++i ) {
02349       altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
02350     }
02351     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
02352   }
02353   {
02354     Tensor altVirial;
02355     for ( int i = 0; i < numAtoms; ++i ) {
02356       altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
02357     }
02358     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
02359   }
02360   {
02361     Tensor altVirial;
02362     for ( int i = 0; i < numAtoms; ++i ) {
02363       altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
02364     }
02365     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
02366   }
02367 #endif
02368 
02369   // For non-Multigrator doKineticEnergy = 1 always
02370   if (doKineticEnergy || patch->flags.doVirial)
02371   {
02372     BigReal intKineticEnergy = 0;
02373     Tensor intVirialNormal;
02374     Tensor intVirialNbond;
02375     Tensor intVirialSlow;
02376 
02377     int hgs;
02378     for ( int i = 0; i < numAtoms; i += hgs ) {
02379 #if CMK_BLUEGENEL
02380       CmiNetworkProgress();
02381 #endif
02382       hgs = a[i].hydrogenGroupSize;
02383       int j;
02384       BigReal m_cm = 0;
02385       Position x_cm(0,0,0);
02386       Velocity v_cm(0,0,0);
02387       for ( j = i; j < (i+hgs); ++j ) {
02388         m_cm += a[j].mass;
02389         x_cm += a[j].mass * a[j].position;
02390         v_cm += a[j].mass * a[j].velocity;
02391       }
02392       x_cm /= m_cm;
02393       v_cm /= m_cm;
02394       int fixedAtomsOn = simParams->fixedAtomsOn;
02395       if ( simParams->pairInteractionOn ) {
02396         int pairInteractionSelf = simParams->pairInteractionSelf;
02397         for ( j = i; j < (i+hgs); ++j ) {
02398           if ( a[j].partition != 1 &&
02399                ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
02400           // net force treated as zero for fixed atoms
02401           if ( fixedAtomsOn && a[j].atomFixed ) continue;
02402           BigReal mass = a[j].mass;
02403           Vector v = a[j].velocity;
02404           Vector dv = v - v_cm;
02405           intKineticEnergy += mass * (v * dv);
02406           Vector dx = a[j].position - x_cm;
02407           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02408           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02409           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02410         }
02411       } else {
02412         for ( j = i; j < (i+hgs); ++j ) {
02413           // net force treated as zero for fixed atoms
02414           if ( fixedAtomsOn && a[j].atomFixed ) continue;
02415           BigReal mass = a[j].mass;
02416           Vector v = a[j].velocity;
02417           Vector dv = v - v_cm;
02418           intKineticEnergy += mass * (v * dv);
02419           Vector dx = a[j].position - x_cm;
02420           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
02421           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
02422           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
02423         }
02424       }
02425     }
02426 
02427     intKineticEnergy *= 0.5;
02428     reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
02429     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
02430     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
02431     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
02432   }
02433 
02434   if (pressureProfileReduction && simParams->useGroupPressure) {
02435     // subtract off internal virial term, calculated as for intVirial.
02436     int nslabs = simParams->pressureProfileSlabs;
02437     const Lattice &lattice = patch->lattice;
02438     BigReal idz = nslabs/lattice.c().z;
02439     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
02440     int useGroupPressure = simParams->useGroupPressure;
02441 
02442     int hgs;
02443     for (int i=0; i<numAtoms; i += hgs) {
02444       int j;
02445       hgs = a[i].hydrogenGroupSize;
02446       BigReal m_cm = 0;
02447       Position x_cm(0,0,0);
02448       for (j=i; j< i+hgs; ++j) {
02449         m_cm += a[j].mass;
02450         x_cm += a[j].mass * a[j].position;
02451       }
02452       x_cm /= m_cm;
02453       
02454       BigReal z = a[i].position.z;
02455       int slab = (int)floor((z-zmin)*idz);
02456       if (slab < 0) slab += nslabs;
02457       else if (slab >= nslabs) slab -= nslabs;
02458       int partition = a[i].partition;
02459       int ppoffset = 3*(slab + nslabs*partition);
02460       for (j=i; j < i+hgs; ++j) {
02461         BigReal mass = a[j].mass;
02462         Vector dx = a[j].position - x_cm;
02463         const Vector &fnormal = patch->f[Results::normal][j];
02464         const Vector &fnbond  = patch->f[Results::nbond][j];
02465         const Vector &fslow   = patch->f[Results::slow][j];
02466         BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
02467         BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
02468         BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
02469         pressureProfileReduction->item(ppoffset  ) -= wxx;
02470         pressureProfileReduction->item(ppoffset+1) -= wyy;
02471         pressureProfileReduction->item(ppoffset+2) -= wzz;
02472       }
02473     }
02474   }
02475 
02476   // For non-Multigrator doVirial = 1 always
02477   if (patch->flags.doVirial)
02478   {
02479     if ( simParams->fixedAtomsOn ) {
02480       Tensor fixVirialNormal;
02481       Tensor fixVirialNbond;
02482       Tensor fixVirialSlow;
02483       Vector fixForceNormal = 0;
02484       Vector fixForceNbond = 0;
02485       Vector fixForceSlow = 0;
02486 
02487       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
02488 
02489       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
02490       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
02491       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
02492       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
02493       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
02494       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
02495     }
02496   }
02497 #endif // UPPER_BOUND
02498 
02499   reduction->submit();
02500 #ifndef UPPER_BOUND
02501   if (pressureProfileReduction) pressureProfileReduction->submit();
02502 #endif
02503 }

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

Definition at line 1848 of file Sequencer.C.

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

Referenced by integrate().

01849 {
01850   if ( simParams->tCoupleOn )
01851   {
01852     FullAtom *a = patch->atom.begin();
01853     int numAtoms = patch->numAtoms;
01854     // Blocking receive for the temperature coupling coefficient.
01855     BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
01856     Molecule *molecule = Node::Object()->molecule;
01857     BigReal dt = dt_fs * 0.001;  // convert to ps
01858     coefficient *= dt;
01859     for ( int i = 0; i < numAtoms; ++i )
01860     {
01861       BigReal f1 = exp( coefficient * a[i].langevinParam );
01862       a[i].velocity *= f1;
01863     }
01864   }
01865 }

void Sequencer::terminate ( void   )  [protected]

Definition at line 2898 of file Sequencer.C.

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

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

02898                           {
02899   LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
02900   CthFree(thread);
02901   CthSuspend();
02902 }

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

Definition at line 2886 of file Sequencer.C.

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

Referenced by integrate().

02886                                     {
02887         // Blocking receive for the trace barrier.
02888         broadcast->traceBarrier.get(step);
02889 }


Friends And Related Function Documentation

friend class HomePatch [friend]

Definition at line 24 of file Sequencer.h.


Member Data Documentation

Definition at line 82 of file Sequencer.h.

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

Definition at line 105 of file Sequencer.h.

Referenced by algorithm().

Definition at line 137 of file Sequencer.h.

Referenced by submitCollections().

int Sequencer::doKineticEnergy [protected]

Definition at line 120 of file Sequencer.h.

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

int Sequencer::doMomenta [protected]

Definition at line 121 of file Sequencer.h.

Referenced by integrate(), and submitReductions().

int Sequencer::ldbSteps [protected]

Definition at line 140 of file Sequencer.h.

Referenced by rebalanceLoad().

Definition at line 39 of file Sequencer.h.

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

Definition at line 119 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().

Definition at line 42 of file Sequencer.h.

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

HomePatch* const Sequencer::patch [protected]
Random* Sequencer::random [protected]

Definition at line 87 of file Sequencer.h.

Referenced by rescaleVelocities(), and Sequencer().

SimParameters* const Sequencer::simParams [protected]
int Sequencer::slowFreq [protected]

Definition at line 107 of file Sequencer.h.

Referenced by integrate(), and langevinPiston().

Count time steps until next stochastic velocity rescaling.

Definition at line 101 of file Sequencer.h.

Referenced by Sequencer(), and stochRescaleVelocities().


The documentation for this class was generated from the following files:

Generated on 20 Oct 2019 for NAMD by  doxygen 1.6.1