NAMD
Public Member Functions | Protected Member Functions | Protected Attributes | Friends | List of all members
Sequencer Class Reference

#include <Sequencer.h>

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

66  :
67  simParams(Node::Object()->simParameters),
68  patch(p),
70  ldbSteps(0)
71 {
77  int ntypes = simParams->pressureProfileAtomTypes;
78  int nslabs = simParams->pressureProfileSlabs;
81  REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
82  } else {
84  }
85  if (simParams->multigratorOn) {
87  } else {
88  multigratorReduction = NULL;
89  }
90  ldbCoordinator = (LdbCoordinator::Object());
93 
94  // Is soluteScaling enabled?
96  // If so, we must "manually" perform charge scaling on startup because
97  // Sequencer will not get a scripting task for initial charge scaling.
98  // Subsequent rescalings will take place through a scripting task.
100  }
101 
103  stochRescale_count = 0;
105 // patch->write_tip4_props();
106 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
SubmitReduction * multigratorReduction
Definition: Sequencer.h:119
static CollectionMgr * Object()
Definition: CollectionMgr.h:30
BigReal soluteScalingFactorCharge
SubmitReduction * pressureProfileReduction
Definition: Sequencer.h:135
SubmitReduction * reduction
Definition: Sequencer.h:134
SubmitReduction * min_reduction
Definition: Sequencer.h:39
static PatchMap * Object()
Definition: PatchMap.h:27
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
int pressureProfileSlabs
LDObjHandle ldObjHandle
Definition: HomePatch.h:479
void split(int iStream, int numStreams)
Definition: Random.h:77
Definition: Random.h:37
int rescaleVelocities_numTemps
Definition: Sequencer.h:87
int ldbSteps
Definition: Sequencer.h:140
static LdbCoordinator * Object()
int berendsenPressure_count
Definition: Sequencer.h:104
Random * random
Definition: Sequencer.h:131
PatchID getPatchID()
Definition: Patch.h:114
unsigned int randomSeed
int pressureProfileAtomTypes
int numPatches(void) const
Definition: PatchMap.h:59
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
CollectionMgr *const collection
Definition: Sequencer.h:137
void rescaleSoluteCharges(BigReal)
Definition: Sequencer.C:1833
Bool pressureProfileOn
int stochRescale_count
Definition: Sequencer.h:100
SimParameters *const simParams
Definition: Sequencer.h:132
Sequencer::~Sequencer ( void  )
virtual

Definition at line 108 of file Sequencer.C.

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

109 {
110  delete broadcast;
111  delete reduction;
112  delete min_reduction;
114  delete random;
116 }
SubmitReduction * multigratorReduction
Definition: Sequencer.h:119
SubmitReduction * pressureProfileReduction
Definition: Sequencer.h:135
SubmitReduction * reduction
Definition: Sequencer.h:134
SubmitReduction * min_reduction
Definition: Sequencer.h:39
Random * random
Definition: Sequencer.h:131
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138

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

1741 {
1742  //check if adaptive tempering is enabled and in the right timestep range
1743  if (!simParams->adaptTempOn) return;
1744  if ( (step < simParams->adaptTempFirstStep ) ||
1746  if (simParams->langevinOn) // restore langevin temperature
1748  return;
1749  }
1750  // Get Updated Temperature
1751  if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
1752  // Blocking receive for the updated adaptive tempering temperature.
1754 }
SimpleBroadcastObject< BigReal > adaptTemperature
Definition: Broadcasts.h:86
BigReal langevinTemp
BigReal adaptTempT
Definition: Sequencer.h:82
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
SimParameters *const simParams
Definition: Sequencer.h:132
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< T >::begin(), ResizeArray< T >::const_begin(), colvarproxy_namd::dt(), Patch::f, Patch::flags, NAMD_EVENT_RANGE_2, Patch::numAtoms, patch, and TIMEFACTOR.

Referenced by newtonianVelocities().

1906  {
1907  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1908  NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
1909 #if CMK_BLUEGENEL
1910  CmiNetworkProgressAfter (0);
1911 #endif
1912  const BigReal dt = timestep / TIMEFACTOR;
1913  FullAtom *atom_arr = patch->atom.begin();
1914  ForceList *f_use = (useSaved ? patch->f_saved : patch->f);
1915  const Force *force_arr = f_use[ftag].const_begin();
1916  patch->addForceToMomentum(atom_arr, force_arr, dt, patch->numAtoms);
1917 }
HomePatch *const patch
Definition: Sequencer.h:133
Definition: Vector.h:64
void addForceToMomentum(FullAtom *__restrict atom_arr, const Force *__restrict force_arr, const BigReal dt, int num_atoms)
Definition: HomePatch.C:1933
Flags flags
Definition: Patch.h:127
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
const_iterator const_begin(void) const
Definition: ResizeArray.h:39
#define TIMEFACTOR
Definition: common.h:48
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
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< T >::begin(), ResizeArray< T >::const_begin(), Patch::f, Patch::flags, NAMD_EVENT_RANGE_2, Patch::numAtoms, patch, and TIMEFACTOR.

Referenced by newtonianVelocities().

1923  {
1924  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1925  NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
1926 #if CMK_BLUEGENEL
1927  CmiNetworkProgressAfter (0);
1928 #endif
1929  const BigReal dt1 = timestep1 / TIMEFACTOR;
1930  const BigReal dt2 = timestep2 / TIMEFACTOR;
1931  const BigReal dt3 = timestep3 / TIMEFACTOR;
1932  ForceList *f_use1 = (useSaved1 ? patch->f_saved : patch->f);
1933  ForceList *f_use2 = (useSaved2 ? patch->f_saved : patch->f);
1934  ForceList *f_use3 = (useSaved3 ? patch->f_saved : patch->f);
1935  FullAtom *atom_arr = patch->atom.begin();
1936  const Force *force_arr1 = f_use1[ftag1].const_begin();
1937  const Force *force_arr2 = f_use2[ftag2].const_begin();
1938  const Force *force_arr3 = f_use3[ftag3].const_begin();
1939  patch->addForceToMomentum3 (atom_arr, force_arr1, force_arr2, force_arr3,
1940  dt1, dt2, dt3, patch->numAtoms);
1941 }
HomePatch *const patch
Definition: Sequencer.h:133
Definition: Vector.h:64
Flags flags
Definition: Patch.h:127
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
const_iterator const_begin(void) const
Definition: ResizeArray.h:39
#define TIMEFACTOR
Definition: common.h:48
void addForceToMomentum3(FullAtom *__restrict atom_arr, const Force *__restrict force_arr1, const Force *__restrict force_arr2, const Force *__restrict force_arr3, const BigReal dt1, const BigReal dt2, const BigReal dt3, int num_atoms)
Definition: HomePatch.C:1962
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::addMovDragToPosition ( BigReal  timestep)
protected

Definition at line 718 of file Sequencer.C.

References CompAtomExt::atomFixed, ResizeArray< T >::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().

718  {
719  FullAtom *atom = patch->atom.begin();
720  int numAtoms = patch->numAtoms;
721  Molecule *molecule = Node::Object()->molecule; // need its methods
722  const BigReal movDragGlobVel = simParams->movDragGlobVel;
723  const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
724  Vector movDragVel, dragIncrement;
725  for ( int i = 0; i < numAtoms; ++i )
726  {
727  // skip if fixed atom or zero drag attribute
728  if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
729  || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
730  molecule->get_movdrag_params(movDragVel, atom[i].id);
731  dragIncrement = movDragGlobVel * movDragVel * dt;
732  atom[i].position += dragIncrement;
733  }
734 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
Bool is_atom_movdragged(int atomnum) const
Definition: Molecule.h:1214
unsigned int atomFixed
Definition: NamdTypes.h:96
Definition: Vector.h:64
Position position
Definition: NamdTypes.h:53
int numAtoms
Definition: Patch.h:144
void get_movdrag_params(Vector &v, int atomnum) const
Definition: Molecule.h:1321
BigReal movDragGlobVel
#define TIMEFACTOR
Definition: common.h:48
SimParameters *const simParams
Definition: Sequencer.h:132
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::addRotDragToPosition ( BigReal  timestep)
protected

Definition at line 737 of file Sequencer.C.

References CompAtomExt::atomFixed, ResizeArray< T >::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().

737  {
738  FullAtom *atom = patch->atom.begin();
739  int numAtoms = patch->numAtoms;
740  Molecule *molecule = Node::Object()->molecule; // need its methods
741  const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
742  const BigReal dt = timestep / TIMEFACTOR; // MUST be as in the integrator!
743  BigReal rotDragVel, dAngle;
744  Vector atomRadius;
745  Vector rotDragAxis, rotDragPivot, dragIncrement;
746  for ( int i = 0; i < numAtoms; ++i )
747  {
748  // skip if fixed atom or zero drag attribute
749  if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
750  || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
751  molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
752  dAngle = rotDragGlobVel * rotDragVel * dt;
753  rotDragAxis /= rotDragAxis.length();
754  atomRadius = atom[i].position - rotDragPivot;
755  dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
756  atom[i].position += dragIncrement;
757  }
758 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
unsigned int atomFixed
Definition: NamdTypes.h:96
Definition: Vector.h:64
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
Position position
Definition: NamdTypes.h:53
BigReal rotDragGlobVel
BigReal length(void) const
Definition: Vector.h:169
Bool is_atom_rotdragged(int atomnum) const
Definition: Molecule.h:1230
int numAtoms
Definition: Patch.h:144
void get_rotdrag_params(BigReal &v, Vector &a, Vector &p, int atomnum) const
Definition: Molecule.h:1327
#define TIMEFACTOR
Definition: common.h:48
SimParameters *const simParams
Definition: Sequencer.h:132
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::addVelocityToPosition ( BigReal  timestep)
protected

Definition at line 1943 of file Sequencer.C.

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

Referenced by integrate().

1944 {
1945  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1946  NamdProfileEvent::ADD_VELOCITY_TO_POSITION);
1947 #if CMK_BLUEGENEL
1948  CmiNetworkProgressAfter (0);
1949 #endif
1950  const BigReal dt = timestep / TIMEFACTOR;
1951  FullAtom *atom_arr = patch->atom.begin();
1952  patch->addVelocityToPosition(atom_arr, dt, patch->numAtoms);
1953 }
HomePatch *const patch
Definition: Sequencer.h:133
void addVelocityToPosition(FullAtom *__restrict atom_arr, const BigReal dt, int num_atoms)
Definition: HomePatch.C:2001
Flags flags
Definition: Patch.h:127
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
#define TIMEFACTOR
Definition: common.h:48
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::algorithm ( void  )
protectedvirtual

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

147 {
148  int scriptTask;
149  int scriptSeq = 0;
150  // Blocking receive for the script barrier.
151  while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
152  switch ( scriptTask ) {
153  case SCRIPT_OUTPUT:
155  break;
156  case SCRIPT_FORCEOUTPUT:
158  break;
159  case SCRIPT_MEASURE:
161  break;
162  case SCRIPT_REINITVELS:
164  break;
165  case SCRIPT_RESCALEVELS:
167  break;
170  break;
172  reloadCharges();
173  break;
174  case SCRIPT_CHECKPOINT:
175  patch->checkpoint();
177  break;
178  case SCRIPT_REVERT:
179  patch->revert();
181  pairlistsAreValid = 0;
182  break;
188  break;
189  case SCRIPT_ATOMSENDRECV:
190  case SCRIPT_ATOMSEND:
191  case SCRIPT_ATOMRECV:
192  patch->exchangeAtoms(scriptTask);
193  break;
194  case SCRIPT_MINIMIZE:
195  minimize();
196  break;
197  case SCRIPT_RUN:
198  case SCRIPT_CONTINUE:
199  //
200  // DJH: Call a cleaned up version of integrate().
201  //
202  // We could test for simulation options and call a more basic version
203  // of integrate() where we can avoid performing most tests.
204  //
205  integrate(scriptTask);
206  break;
207  default:
208  NAMD_bug("Unknown task in Sequencer::algorithm");
209  }
210  }
212  terminate();
213 }
HomePatch *const patch
Definition: Sequencer.h:133
void terminate(void)
Definition: Sequencer.C:2909
BigReal soluteScalingFactorCharge
void integrate(int)
Definition: Sequencer.C:218
#define FILE_OUTPUT
Definition: Output.h:25
#define EVAL_MEASURE
Definition: Output.h:27
void exchangeCheckpoint(int scriptTask, int &bpc)
Definition: HomePatch.C:3478
void revert(void)
Definition: HomePatch.C:3460
void submitCollections(int step, int zeroVel=0)
Definition: Sequencer.C:2643
void NAMD_bug(const char *err_msg)
Definition: common.C:123
SimpleBroadcastObject< int > scriptBarrier
Definition: Broadcasts.h:83
BigReal scriptArg1
#define END_OF_RUN
Definition: Output.h:26
int berendsenPressure_count
Definition: Sequencer.h:104
void checkpoint(void)
Definition: HomePatch.C:3450
void reinitVelocities(void)
Definition: Sequencer.C:1788
int checkpoint_berendsenPressure_count
Definition: Sequencer.h:105
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
void rescaleSoluteCharges(BigReal)
Definition: Sequencer.C:1833
int pairlistsAreValid
Definition: Sequencer.h:42
SimParameters *const simParams
Definition: Sequencer.h:132
void rescaleVelocitiesByFactor(BigReal)
Definition: Sequencer.C:1811
void reloadCharges()
Definition: Sequencer.C:1821
#define FORCE_OUTPUT
Definition: Output.h:28
void minimize()
Definition: Sequencer.C:760
void exchangeAtoms(int scriptTask)
Definition: HomePatch.C:3574
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().

29  {
30  CthAwakenPrio(thread, CK_QUEUEING_IFIFO, PRIORITY_SIZE, &priority);
31  }
#define PRIORITY_SIZE
Definition: Priorities.h:13
void Sequencer::berendsenPressure ( int  step)
protected

Definition at line 1534 of file Sequencer.C.

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

Referenced by integrate().

1535 {
1536  if ( simParams->berendsenPressureOn ) {
1538  const int freq = simParams->berendsenPressureFreq;
1539  if ( ! (berendsenPressure_count % freq ) ) {
1541  FullAtom *a = patch->atom.begin();
1542  int numAtoms = patch->numAtoms;
1543  // Blocking receive for the updated lattice scaling factor.
1544  Tensor factor = broadcast->positionRescaleFactor.get(step);
1545  patch->lattice.rescale(factor);
1546  if ( simParams->useGroupPressure )
1547  {
1548  int hgs;
1549  for ( int i = 0; i < numAtoms; i += hgs ) {
1550  int j;
1551  hgs = a[i].hydrogenGroupSize;
1552  if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
1553  for ( j = i; j < (i+hgs); ++j ) {
1555  a[j].fixedPosition,a[j].transform);
1556  }
1557  continue;
1558  }
1559  BigReal m_cm = 0;
1560  Position x_cm(0,0,0);
1561  for ( j = i; j < (i+hgs); ++j ) {
1562  if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
1563  m_cm += a[j].mass;
1564  x_cm += a[j].mass * a[j].position;
1565  }
1566  x_cm /= m_cm;
1567  Position new_x_cm = x_cm;
1568  patch->lattice.rescale(new_x_cm,factor);
1569  Position delta_x_cm = new_x_cm - x_cm;
1570  for ( j = i; j < (i+hgs); ++j ) {
1571  if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
1573  a[j].fixedPosition,a[j].transform);
1574  continue;
1575  }
1576  a[j].position += delta_x_cm;
1577  }
1578  }
1579  }
1580  else
1581  {
1582  for ( int i = 0; i < numAtoms; ++i )
1583  {
1584  if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
1586  a[i].fixedPosition,a[i].transform);
1587  continue;
1588  }
1589  patch->lattice.rescale(a[i].position,factor);
1590  }
1591  }
1592  }
1593  } else {
1595  }
1596 }
HomePatch *const patch
Definition: Sequencer.h:133
Bool berendsenPressureOn
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
unsigned int atomFixed
Definition: NamdTypes.h:96
Lattice & lattice
Definition: Patch.h:126
Definition: Vector.h:64
Position position
Definition: NamdTypes.h:53
int berendsenPressureFreq
unsigned int groupFixed
Definition: NamdTypes.h:97
void rescale(Tensor factor)
Definition: Lattice.h:60
int numAtoms
Definition: Patch.h:144
int berendsenPressure_count
Definition: Sequencer.h:104
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
Position apply_transform(Position data, const Transform &t) const
Definition: Lattice.h:132
Definition: Tensor.h:15
Mass mass
Definition: NamdTypes.h:108
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
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< T >::begin(), Patch::f, SimParameters::fixedAtomsOn, FullAtom::fixedPosition, Results::nbond, Results::normal, Patch::numAtoms, Tensor::outerAdd(), patch, simParams, and Results::slow.

Referenced by multigratorPressure(), submitMinimizeReductions(), and submitReductions().

2243  {
2244 
2245  FullAtom *a = patch->atom.begin();
2246  int numAtoms = patch->numAtoms;
2247 
2248  for ( int j = 0; j < numAtoms; j++ ) {
2249  if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
2250  Vector dx = a[j].fixedPosition;
2251  // all negative because fixed atoms cancels these forces
2252  fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
2253  fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
2254  fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
2255  fixForceNormal -= patch->f[Results::normal][j];
2256  fixForceNbond -= patch->f[Results::nbond][j];
2257  fixForceSlow -= patch->f[Results::slow][j];
2258  }
2259  }
2260 }
HomePatch *const patch
Definition: Sequencer.h:133
unsigned int atomFixed
Definition: NamdTypes.h:96
Position fixedPosition
Definition: NamdTypes.h:102
Definition: Vector.h:64
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
Definition: Tensor.h:255
int numAtoms
Definition: Patch.h:144
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
SimParameters *const simParams
Definition: Sequencer.h:132
iterator begin(void)
Definition: ResizeArray.h:36
BigReal Sequencer::calcKineticEnergy ( )
protected

Definition at line 1227 of file Sequencer.C.

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

Referenced by multigratorTemperature().

1227  {
1228  FullAtom *a = patch->atom.begin();
1229  int numAtoms = patch->numAtoms;
1230  BigReal kineticEnergy = 0.0;
1231  if ( simParams->pairInteractionOn ) {
1232  if ( simParams->pairInteractionSelf ) {
1233  for (int i = 0; i < numAtoms; ++i ) {
1234  if ( a[i].partition != 1 ) continue;
1235  kineticEnergy += a[i].mass * a[i].velocity.length2();
1236  }
1237  }
1238  } else {
1239  for (int i = 0; i < numAtoms; ++i ) {
1240  kineticEnergy += a[i].mass * a[i].velocity.length2();
1241  }
1242  }
1243  kineticEnergy *= 0.5;
1244  return kineticEnergy;
1245 }
HomePatch *const patch
Definition: Sequencer.h:133
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Velocity velocity
Definition: NamdTypes.h:101
Bool pairInteractionOn
int numAtoms
Definition: Patch.h:144
BigReal length2(void) const
Definition: Vector.h:173
Mass mass
Definition: NamdTypes.h:108
Bool pairInteractionSelf
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::correctMomentum ( int  step,
BigReal  drifttime 
)
protected

Definition at line 1016 of file Sequencer.C.

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

Referenced by integrate().

1016  {
1017 
1018  //
1019  // DJH: This test should be done in SimParameters.
1020  //
1021  if ( simParams->fixedAtomsOn )
1022  NAMD_die("Cannot zero momentum when fixed atoms are present.");
1023 
1024  // Blocking receive for the momentum correction vector.
1025  const Vector dv = broadcast->momentumCorrection.get(step);
1026 
1027  const Vector dx = dv * ( drifttime / TIMEFACTOR );
1028 
1029  FullAtom *a = patch->atom.begin();
1030  const int numAtoms = patch->numAtoms;
1031 
1032 if ( simParams->zeroMomentumAlt ) {
1033  for ( int i = 0; i < numAtoms; ++i ) {
1034  a[i].velocity += dv * a[i].recipMass;
1035  a[i].position += dx * a[i].recipMass;
1036  }
1037 } else {
1038  for ( int i = 0; i < numAtoms; ++i ) {
1039  a[i].velocity += dv;
1040  a[i].position += dx;
1041  }
1042 }
1043 
1044 }
HomePatch *const patch
Definition: Sequencer.h:133
SimpleBroadcastObject< Vector > momentumCorrection
Definition: Broadcasts.h:79
Definition: Vector.h:64
int numAtoms
Definition: Patch.h:144
void NAMD_die(const char *err_msg)
Definition: common.C:83
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
#define TIMEFACTOR
Definition: common.h:48
SimParameters *const simParams
Definition: Sequencer.h:132
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::cycleBarrier ( int  doBarrier,
int  step 
)
protected

Definition at line 2889 of file Sequencer.C.

References broadcast.

Referenced by integrate().

2889  {
2890 #if USE_BARRIER
2891  if (doBarrier)
2892  // Blocking receive for the cycle barrier.
2893  broadcast->cycleBarrier.get(step);
2894 #endif
2895 }
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
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().

1956 {
1957  if ( simParams->drudeHardWallOn ) {
1958  Tensor virial;
1959  Tensor *vp = ( pressure ? &virial : 0 );
1960  if ( patch->hardWallDrude(dt, vp, pressureProfileReduction) ) {
1961  iout << iERROR << "Constraint failure in HardWallDrude(); "
1962  << "simulation may become unstable.\n" << endi;
1964  terminate();
1965  }
1966  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
1967  }
1968 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
void terminate(void)
Definition: Sequencer.C:2909
SubmitReduction * pressureProfileReduction
Definition: Sequencer.h:135
SubmitReduction * reduction
Definition: Sequencer.h:134
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
#define iout
Definition: InfoStream.h:87
int hardWallDrude(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2024
Definition: Tensor.h:15
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void enableEarlyExit(void)
Definition: Node.C:1365
SimParameters *const simParams
Definition: Sequencer.h:132
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< T >::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< T >::end(), NamdProfileEvent::EventsCount, 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, SimParameters::singleTopology, 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().

218  {
219  char traceNote[24];
220  char tracePrefix[20];
221  sprintf(tracePrefix, "p:%d,s:",patch->patchID);
222 // patch->write_tip4_props();
223 
224  //
225  // DJH: Copy all data into SOA (structure of arrays)
226  // from AOS (array of structures) data structure.
227  //
228  //patch->copy_all_to_SOA();
229 
230 #ifdef TIMER_COLLECTION
231  TimerSet& t = patch->timerSet;
232 #endif
233  TIMER_INIT_WIDTH(t, KICK, simParams->timerBinWidth);
234  TIMER_INIT_WIDTH(t, MAXMOVE, simParams->timerBinWidth);
235  TIMER_INIT_WIDTH(t, DRIFT, simParams->timerBinWidth);
236  TIMER_INIT_WIDTH(t, PISTON, simParams->timerBinWidth);
237  TIMER_INIT_WIDTH(t, SUBMITHALF, simParams->timerBinWidth);
238  TIMER_INIT_WIDTH(t, VELBBK1, simParams->timerBinWidth);
239  TIMER_INIT_WIDTH(t, VELBBK2, simParams->timerBinWidth);
240  TIMER_INIT_WIDTH(t, RATTLE1, simParams->timerBinWidth);
241  TIMER_INIT_WIDTH(t, SUBMITFULL, simParams->timerBinWidth);
242  TIMER_INIT_WIDTH(t, SUBMITCOLLECT, simParams->timerBinWidth);
243 
244  int &step = patch->flags.step;
245  step = simParams->firstTimestep;
246 
247  // drag switches
248  const Bool rotDragOn = simParams->rotDragOn;
249  const Bool movDragOn = simParams->movDragOn;
250 
251  const int commOnly = simParams->commOnly;
252 
253  int &maxForceUsed = patch->flags.maxForceUsed;
254  int &maxForceMerged = patch->flags.maxForceMerged;
255  maxForceUsed = Results::normal;
256  maxForceMerged = Results::normal;
257 
258  const int numberOfSteps = simParams->N;
259  const int stepsPerCycle = simParams->stepsPerCycle;
260  const BigReal timestep = simParams->dt;
261 
262  // what MTS method?
263  const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
264 
265  const int nonbondedFrequency = simParams->nonbondedFrequency;
266  slowFreq = nonbondedFrequency;
267  const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
268  int &doNonbonded = patch->flags.doNonbonded;
269  doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
270  if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
271  if ( doNonbonded ) maxForceUsed = Results::nbond;
272 
273  // Do we do full electrostatics?
274  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
275  const int fullElectFrequency = simParams->fullElectFrequency;
276  if ( dofull ) slowFreq = fullElectFrequency;
277  const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
278  int &doFullElectrostatics = patch->flags.doFullElectrostatics;
279  doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
280  if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
281  maxForceMerged = Results::slow;
282  if ( doFullElectrostatics ) maxForceUsed = Results::slow;
283 
284 //#ifndef UPPER_BOUND
285  const Bool accelMDOn = simParams->accelMDOn;
286  const Bool accelMDdihe = simParams->accelMDdihe;
287  const Bool accelMDdual = simParams->accelMDdual;
288  if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
289 
290  // Is adaptive tempering on?
291  const Bool adaptTempOn = simParams->adaptTempOn;
293  if (simParams->langevinOn)
295  else if (simParams->rescaleFreq > 0)
297 
298 
299  int &doMolly = patch->flags.doMolly;
300  doMolly = simParams->mollyOn && doFullElectrostatics;
301  // BEGIN LA
302  int &doLoweAndersen = patch->flags.doLoweAndersen;
303  doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
304  // END LA
305 
306  int &doGBIS = patch->flags.doGBIS;
307  doGBIS = simParams->GBISOn;
308 
309  int &doLCPO = patch->flags.doLCPO;
310  doLCPO = simParams->LCPOOn;
311 
312  int zeroMomentum = simParams->zeroMomentum;
313 
314  // Do we need to return forces to TCL script or Colvar module?
315  int doTcl = simParams->tclForcesOn;
316  int doColvars = simParams->colvarsOn;
317 //#endif
319 
320  // Bother to calculate energies?
321  int &doEnergy = patch->flags.doEnergy;
322  int energyFrequency = simParams->outputEnergies;
323 #ifndef UPPER_BOUND
324  const int reassignFreq = simParams->reassignFreq;
325 #endif
326 
327  int &doVirial = patch->flags.doVirial;
328  doVirial = 1;
329 
330  if ( scriptTask == SCRIPT_RUN ) {
331 
332 #ifndef UPPER_BOUND
333 // printf("Doing initial rattle\n");
334 #ifndef UPPER_BOUND
335 D_MSG("rattle1()");
336  TIMER_START(t, RATTLE1);
337  rattle1(0.,0); // enforce rigid bond constraints on initial positions
338  TIMER_STOP(t, RATTLE1);
339 #endif
340 
343  patch->atom.begin(),patch->atom.end());
344  }
345 
346  if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
347  reassignVelocities(timestep,step);
348  }
349 #endif
350 
351  doEnergy = ! ( step % energyFrequency );
352 #ifndef UPPER_BOUND
353  if ( accelMDOn && !accelMDdihe ) doEnergy=1;
354  //Update energy every timestep for adaptive tempering
355  if ( adaptTempOn ) doEnergy=1;
356 #endif
357 D_MSG("runComputeObjects()");
358  runComputeObjects(1,step<numberOfSteps); // must migrate here!
359 #ifndef UPPER_BOUND
360  rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
361  adaptTempUpdate(step); // update adaptive tempering temperature
362 #endif
363 
364 #ifndef UPPER_BOUND
365  if ( staleForces || doTcl || doColvars ) {
366  if ( doNonbonded ) saveForce(Results::nbond);
367  if ( doFullElectrostatics ) saveForce(Results::slow);
368  }
369  if ( ! commOnly ) {
370 D_MSG("newtonianVelocities()");
371  TIMER_START(t, KICK);
372  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
373  TIMER_STOP(t, KICK);
374  }
376 #ifndef UPPER_BOUND
377 D_MSG("rattle1()");
378  TIMER_START(t, RATTLE1);
379  rattle1(-timestep,0);
380  TIMER_STOP(t, RATTLE1);
381 #endif
382 D_MSG("submitHalfstep()");
383  TIMER_START(t, SUBMITHALF);
384  submitHalfstep(step);
385  TIMER_STOP(t, SUBMITHALF);
386  if ( ! commOnly ) {
387 D_MSG("newtonianVelocities()");
388  TIMER_START(t, KICK);
389  newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
390  TIMER_STOP(t, KICK);
391  }
392 D_MSG("rattle1()");
393  TIMER_START(t, RATTLE1);
394  rattle1(timestep,1);
395  TIMER_STOP(t, RATTLE1);
396  if (doTcl || doColvars) // include constraint forces
397  computeGlobal->saveTotalForces(patch);
398 D_MSG("submitHalfstep()");
399  TIMER_START(t, SUBMITHALF);
400  submitHalfstep(step);
401  TIMER_STOP(t, SUBMITHALF);
402  if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
403  if ( ! commOnly ) {
404 D_MSG("newtonianVelocities()");
405  TIMER_START(t, KICK);
406  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
407  TIMER_STOP(t, KICK);
408  }
409 #endif
410 D_MSG("submitReductions()");
411  TIMER_START(t, SUBMITFULL);
412  submitReductions(step);
413  TIMER_STOP(t, SUBMITFULL);
414 #ifndef UPPER_BOUND
415  if(0){ // if(traceIsOn()){
416  traceUserEvent(eventEndOfTimeStep);
417  sprintf(traceNote, "%s%d",tracePrefix,step);
418  traceUserSuppliedNote(traceNote);
419  }
420 #endif
421  rebalanceLoad(step);
422 
423  } // scriptTask == SCRIPT_RUN
424 
425 #ifndef UPPER_BOUND
426  bool doMultigratorRattle = false;
427 #endif
428 
429  //
430  // DJH: There are a lot of mod operations below and elsewhere to
431  // test step number against the frequency of something happening.
432  // Mod and integer division are expensive!
433  // Might be better to replace with counters and test equality.
434  //
435 #if 0
436  for(int i = 0; i < NamdProfileEvent::EventsCount; i++)
437  CkPrintf("-------------- [%d] %s -------------\n", i, NamdProfileEventStr[i]);
438 #endif
439 
440 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
441  int& eon = patch->flags.event_on;
442  int epid = (simParams->beginEventPatchID <= patch->getPatchID()
443  && patch->getPatchID() <= simParams->endEventPatchID);
444  int beginStep = simParams->beginEventStep;
445  int endStep = simParams->endEventStep;
446  bool controlProfiling = patch->getPatchID() == 0;
447 #endif
448 
449  for ( ++step; step <= numberOfSteps; ++step )
450  {
451 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
452  eon = epid && (beginStep < step && step <= endStep);
453 
454  if (controlProfiling && step == beginStep) {
456  }
457  if (controlProfiling && step == endStep) {
459  }
460  char buf[32];
461  sprintf(buf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::INTEGRATE_1], patch->getPatchID());
462  NAMD_EVENT_START_EX(eon, NamdProfileEvent::INTEGRATE_1, buf);
463 #endif
464 #ifndef UPPER_BOUND
465 
466  rescaleVelocities(step);
467  tcoupleVelocities(timestep,step);
468  stochRescaleVelocities(timestep,step);
469  berendsenPressure(step);
470 
471  if ( ! commOnly ) {
472  TIMER_START(t, KICK);
473  newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
474  TIMER_STOP(t, KICK);
475  }
476 
477  // We do RATTLE here if multigrator thermostat was applied in the previous step
478  if (doMultigratorRattle) rattle1(timestep, 1);
479 
480  /* reassignment based on half-step velocities
481  if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
482  addVelocityToPosition(0.5*timestep);
483  reassignVelocities(timestep,step);
484  addVelocityToPosition(0.5*timestep);
485  rattle1(0.,0);
486  rattle1(-timestep,0);
487  addVelocityToPosition(-1.0*timestep);
488  rattle1(timestep,0);
489  } */
490 
491  TIMER_START(t, MAXMOVE);
492  maximumMove(timestep);
493  TIMER_STOP(t, MAXMOVE);
494 
495  NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_1); // integrate 1
496 
498  if ( ! commOnly ) {
499  TIMER_START(t, DRIFT);
500  addVelocityToPosition(0.5*timestep);
501  TIMER_STOP(t, DRIFT);
502  }
503  // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
504  langevinVelocities(timestep);
505 
506  // There is a blocking receive inside of langevinPiston()
507  // that might suspend the current thread of execution,
508  // so split profiling around this conditional block.
509  langevinPiston(step);
510 
511  if ( ! commOnly ) {
512  TIMER_START(t, DRIFT);
513  addVelocityToPosition(0.5*timestep);
514  TIMER_STOP(t, DRIFT);
515  }
516  } else {
517  // If Langevin is not used, take full time step directly instread of two half steps
518  if ( ! commOnly ) {
519  TIMER_START(t, DRIFT);
520  addVelocityToPosition(timestep);
521  TIMER_STOP(t, DRIFT);
522  }
523  }
524 
525  NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_2);
526 
527  // impose hard wall potential for Drude bond length
528  hardWallDrude(timestep, 1);
529 
531 #endif // UPPER_BOUND
532 
533  doNonbonded = !(step%nonbondedFrequency);
534  doFullElectrostatics = (dofull && !(step%fullElectFrequency));
535 
536 #ifndef UPPER_BOUND
537  if ( zeroMomentum && doFullElectrostatics ) {
538  // There is a blocking receive inside of correctMomentum().
539  correctMomentum(step,slowstep);
540  }
541 
542  // There are NO sends in submitHalfstep() just local summation
543  // into the Reduction struct.
544  TIMER_START(t, SUBMITHALF);
545  submitHalfstep(step);
546  TIMER_STOP(t, SUBMITHALF);
547 
548  doMolly = simParams->mollyOn && doFullElectrostatics;
549  // BEGIN LA
550  doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
551  // END LA
552 
553  maxForceUsed = Results::normal;
554  if ( doNonbonded ) maxForceUsed = Results::nbond;
555  if ( doFullElectrostatics ) maxForceUsed = Results::slow;
556  if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
557 
558  // Migrate Atoms on stepsPerCycle
559  doEnergy = ! ( step % energyFrequency );
560  if ( accelMDOn && !accelMDdihe ) doEnergy=1;
561  if ( adaptTempOn ) doEnergy=1;
562 
563  // Multigrator
564  if (simParams->multigratorOn) {
565  doVirial = (!(step % energyFrequency) || ((simParams->outputPressure > 0) && !(step % simParams->outputPressure))
566  || !(step % simParams->multigratorPressureFreq));
567  doKineticEnergy = (!(step % energyFrequency) || !(step % simParams->multigratorTemperatureFreq));
569  } else {
570  doVirial = 1;
571  doKineticEnergy = 1;
572  doMomenta = 1;
573  }
574 #endif
575  NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_2); // integrate 2
576 
577  // The current thread of execution will suspend in runComputeObjects().
578  runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
579 
580  NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_3);
581 
582 #ifndef UPPER_BOUND
583  rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
584 
585  if ( staleForces || doTcl || doColvars ) {
586  if ( doNonbonded ) saveForce(Results::nbond);
587  if ( doFullElectrostatics ) saveForce(Results::slow);
588  }
589 
590  // reassignment based on full-step velocities
591  if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
592  reassignVelocities(timestep,step);
593  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
594  rattle1(-timestep,0);
595  }
596 
597  if ( ! commOnly ) {
598  TIMER_START(t, VELBBK1);
599  langevinVelocitiesBBK1(timestep);
600  TIMER_STOP(t, VELBBK1);
601  TIMER_START(t, KICK);
602  newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
603  TIMER_STOP(t, KICK);
604  TIMER_START(t, VELBBK2);
605  langevinVelocitiesBBK2(timestep);
606  TIMER_STOP(t, VELBBK2);
607  }
608 
609  // add drag to each atom's positions
610  if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
611  if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
612 
613  TIMER_START(t, RATTLE1);
614  rattle1(timestep,1);
615  TIMER_STOP(t, RATTLE1);
616  if (doTcl || doColvars) // include constraint forces
617  computeGlobal->saveTotalForces(patch);
618 
619  TIMER_START(t, SUBMITHALF);
620  submitHalfstep(step);
621  TIMER_STOP(t, SUBMITHALF);
622  if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
623 
624  if ( ! commOnly ) {
625  TIMER_START(t, KICK);
626  newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
627  TIMER_STOP(t, KICK);
628  }
629 
630  // rattle2(timestep,step);
631 #endif
632 
633  TIMER_START(t, SUBMITFULL);
634  submitReductions(step);
635  TIMER_STOP(t, SUBMITFULL);
636  TIMER_START(t, SUBMITCOLLECT);
637  submitCollections(step);
638  TIMER_STOP(t, SUBMITCOLLECT);
639 #ifndef UPPER_BOUND
640  //Update adaptive tempering temperature
641  adaptTempUpdate(step);
642 
643  // Multigrator temperature and pressure steps
644  multigratorTemperature(step, 1);
645  multigratorPressure(step, 1);
646  multigratorPressure(step, 2);
647  multigratorTemperature(step, 2);
648  doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
649 
650  NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_3); // integrate 3
651 #endif
652 
653 #if CYCLE_BARRIER
654  cycleBarrier(!((step+1) % stepsPerCycle), step);
655 #elif PME_BARRIER
656  cycleBarrier(doFullElectrostatics, step);
657 #elif STEP_BARRIER
658  cycleBarrier(1, step);
659 #endif
660 
661 #ifndef UPPER_BOUND
662  if(Node::Object()->specialTracing || simParams->statsOn){
663  int bstep = simParams->traceStartStep;
664  int estep = bstep + simParams->numTraceSteps;
665  if(step == bstep || step == estep){
666  traceBarrier(step);
667  }
668  }
669 
670 #ifdef MEASURE_NAMD_WITH_PAPI
671  if(simParams->papiMeasure) {
672  int bstep = simParams->papiMeasureStartStep;
673  int estep = bstep + simParams->numPapiMeasureSteps;
674  if(step == bstep || step==estep) {
675  papiMeasureBarrier(step);
676  }
677  }
678 #endif
679 
680  if(0){ // if(traceIsOn()){
681  traceUserEvent(eventEndOfTimeStep);
682  sprintf(traceNote, "%s%d",tracePrefix,step);
683  traceUserSuppliedNote(traceNote);
684  }
685 #endif // UPPER_BOUND
686  rebalanceLoad(step);
687 
688 #if PME_BARRIER
689  // a step before PME
690  cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
691 #endif
692 
693 #if USE_HPM
694  if(step == START_HPM_STEP)
695  (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
696 
697  if(step == STOP_HPM_STEP)
698  (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
699 #endif
700 
701  }
702 
703  TIMER_DONE(t);
704 #ifdef TIMER_COLLECTION
705  if (patch->patchID == SPECIAL_PATCH_ID) {
706  printf("Timer collection reporting in microseconds for "
707  "Patch %d\n", patch->patchID);
708  TIMER_REPORT(t);
709  }
710 #endif // TIMER_COLLECTION
711  //
712  // DJH: Copy updates of SOA back into AOS.
713  //
714  //patch->copy_updates_to_AOS();
715 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
void rescaleVelocities(int)
Definition: Sequencer.C:1678
int doKineticEnergy
Definition: Sequencer.h:120
void minimizationQuenchVelocity(void)
Definition: Sequencer.C:2072
void tcoupleVelocities(BigReal, int)
Definition: Sequencer.C:1848
void addMovDragToPosition(BigReal)
Definition: Sequencer.C:718
#define NAMD_EVENT_STOP(eon, id)
void addVelocityToPosition(BigReal)
Definition: Sequencer.C:1943
int eventEndOfTimeStep
Definition: Node.C:286
void maximumMove(BigReal)
Definition: Sequencer.C:2027
void cycleBarrier(int, int)
Definition: Sequencer.C:2889
void addRotDragToPosition(BigReal)
Definition: Sequencer.C:737
void saveForce(const int ftag=Results::normal)
Definition: Sequencer.C:1893
void registerIDsFullAtom(const FullAtom *begin, const FullAtom *end)
Definition: AtomMap.C:50
void langevinVelocitiesBBK2(BigReal)
Definition: Sequencer.C:1427
int slowFreq
Definition: Sequencer.h:107
void langevinVelocitiesBBK1(BigReal)
Definition: Sequencer.C:1354
char const *const NamdProfileEventStr[]
void rattle1(BigReal, int)
Definition: Sequencer.C:1970
void saveTotalForces(HomePatch *)
void rebalanceLoad(int timestep)
Definition: Sequencer.C:2878
void submitHalfstep(int)
Definition: Sequencer.C:2083
int doLoweAndersen
Definition: PatchTypes.h:26
void langevinPiston(int)
Definition: Sequencer.C:1598
AtomMapper * atomMapper
Definition: Patch.h:152
Flags flags
Definition: Patch.h:127
void submitCollections(int step, int zeroVel=0)
Definition: Sequencer.C:2643
void adaptTempUpdate(int)
Definition: Sequencer.C:1740
Bool langevin_useBAOAB
#define TIMER_START(T, TYPE)
Definition: HomePatch.h:265
#define NAIVE
Definition: SimParameters.h:49
#define NAMD_PROFILE_START()
#define NAMD_EVENT_START(eon, id)
BigReal rescaleTemp
#define TIMER_REPORT(T)
Definition: HomePatch.h:268
void multigratorPressure(int step, int callNumber)
Definition: Sequencer.C:1084
int doEnergy
Definition: PatchTypes.h:20
void berendsenPressure(int)
Definition: Sequencer.C:1534
void submitMomentum(int step)
Definition: Sequencer.C:993
int doFullElectrostatics
Definition: PatchTypes.h:23
iterator end(void)
Definition: ResizeArray.h:37
void runComputeObjects(int migration=1, int pairlists=0, int pressureStep=0)
Definition: Sequencer.C:2668
void rescaleaccelMD(int, int, int)
Definition: Sequencer.C:1697
int Bool
Definition: common.h:133
BigReal langevinTemp
MTSChoices MTSAlgorithm
BigReal adaptTempT
Definition: Sequencer.h:82
int maxForceUsed
Definition: PatchTypes.h:31
#define D_MSG(t)
Definition: Debug.h:149
void traceBarrier(int)
Definition: Sequencer.C:2897
int doNonbonded
Definition: PatchTypes.h:22
#define TIMER_INIT_WIDTH(T, TYPE, WIDTH)
Definition: HomePatch.h:264
void reassignVelocities(BigReal, int)
Definition: Sequencer.C:1756
ComputeGlobal * computeGlobalObject
Definition: ComputeMgr.h:97
void langevinVelocities(BigReal)
Definition: Sequencer.C:1317
void hardWallDrude(BigReal, int)
Definition: Sequencer.C:1955
#define TIMER_STOP(T, TYPE)
Definition: HomePatch.h:266
PatchID getPatchID()
Definition: Patch.h:114
void multigratorTemperature(int step, int callNumber)
Definition: Sequencer.C:1247
BigReal initialTemp
#define NAMD_EVENT_START_EX(eon, id, str)
void stochRescaleVelocities(BigReal, int)
Definition: Sequencer.C:1872
const PatchID patchID
Definition: Patch.h:143
#define NAMD_PROFILE_STOP()
int doVirial
Definition: PatchTypes.h:21
int doLCPO
Definition: PatchTypes.h:29
void newtonianVelocities(BigReal, const BigReal, const BigReal, const BigReal, const int, const int, const int)
Definition: Sequencer.C:1293
int doMomenta
Definition: Sequencer.h:121
#define TIMER_DONE(T)
Definition: HomePatch.h:267
int multigratorPressureFreq
#define SPECIAL_PATCH_ID
Definition: Sequencer.C:64
void correctMomentum(int step, BigReal drifttime)
Definition: Sequencer.C:1016
int doGBIS
Definition: PatchTypes.h:28
ComputeMgr * computeMgr
Definition: Node.h:169
int maxForceMerged
Definition: PatchTypes.h:32
void submitReductions(int)
Definition: Sequencer.C:2262
SimParameters *const simParams
Definition: Sequencer.h:132
int doMolly
Definition: PatchTypes.h:24
int multigratorTemperatureFreq
double BigReal
Definition: common.h:114
int step
Definition: PatchTypes.h:16
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::langevinPiston ( int  step)
protected

Definition at line 1598 of file Sequencer.C.

References Lattice::apply_transform(), CompAtomExt::atomFixed, ResizeArray< T >::begin(), broadcast, SimParameters::fixedAtomsOn, for(), SimpleBroadcastObject< T >::get(), CompAtomExt::groupFixed, CompAtom::hydrogenGroupSize, CompAtomExt::id, Molecule::is_atom_exPressure(), 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().

1599 {
1600  if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
1601  {
1602  //
1603  // DJH: Loops below simplify if we lift out special cases of fixed atoms
1604  // and pressure excluded atoms and make them their own branch.
1605  //
1606  FullAtom *a = patch->atom.begin();
1607  int numAtoms = patch->numAtoms;
1608  // Blocking receive for the updated lattice scaling factor.
1609  Tensor factor = broadcast->positionRescaleFactor.get(step);
1610  TIMER_START(patch->timerSet, PISTON);
1611  // JCP FIX THIS!!!
1612  Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
1613  patch->lattice.rescale(factor);
1614  Molecule *mol = Node::Object()->molecule;
1615  if ( simParams->useGroupPressure )
1616  {
1617  int hgs;
1618  for ( int i = 0; i < numAtoms; i += hgs ) {
1619  int j;
1620  hgs = a[i].hydrogenGroupSize;
1621  if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
1622  for ( j = i; j < (i+hgs); ++j ) {
1624  a[j].fixedPosition,a[j].transform);
1625  }
1626  continue;
1627  }
1628  BigReal m_cm = 0;
1629  Position x_cm(0,0,0);
1630  Velocity v_cm(0,0,0);
1631  for ( j = i; j < (i+hgs); ++j ) {
1632  if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
1633  m_cm += a[j].mass;
1634  x_cm += a[j].mass * a[j].position;
1635  v_cm += a[j].mass * a[j].velocity;
1636  }
1637  x_cm /= m_cm;
1638  Position new_x_cm = x_cm;
1639  patch->lattice.rescale(new_x_cm,factor);
1640  Position delta_x_cm = new_x_cm - x_cm;
1641  v_cm /= m_cm;
1642  Velocity delta_v_cm;
1643  delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
1644  delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
1645  delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
1646  for ( j = i; j < (i+hgs); ++j ) {
1647  if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
1649  a[j].fixedPosition,a[j].transform);
1650  continue;
1651  }
1652  if ( mol->is_atom_exPressure(a[j].id) ) continue;
1653  a[j].position += delta_x_cm;
1654  a[j].velocity += delta_v_cm;
1655  }
1656  }
1657  }
1658  else
1659  {
1660  for ( int i = 0; i < numAtoms; ++i )
1661  {
1662  if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
1664  a[i].fixedPosition,a[i].transform);
1665  continue;
1666  }
1667  if ( mol->is_atom_exPressure(a[i].id) ) continue;
1668  patch->lattice.rescale(a[i].position,factor);
1669  a[i].velocity.x *= velFactor.x;
1670  a[i].velocity.y *= velFactor.y;
1671  a[i].velocity.z *= velFactor.z;
1672  }
1673  }
1674  TIMER_STOP(patch->timerSet, PISTON);
1675  }
1676 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
unsigned int atomFixed
Definition: NamdTypes.h:96
Lattice & lattice
Definition: Patch.h:126
Definition: Vector.h:64
int slowFreq
Definition: Sequencer.h:107
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
Velocity velocity
Definition: NamdTypes.h:101
unsigned int groupFixed
Definition: NamdTypes.h:97
void rescale(Tensor factor)
Definition: Lattice.h:60
#define TIMER_START(T, TYPE)
Definition: HomePatch.h:265
int numAtoms
Definition: Patch.h:144
BigReal x
Definition: Vector.h:66
BigReal xx
Definition: Tensor.h:17
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
BigReal zz
Definition: Tensor.h:19
#define TIMER_STOP(T, TYPE)
Definition: HomePatch.h:266
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
Position apply_transform(Position data, const Transform &t) const
Definition: Lattice.h:132
Definition: Tensor.h:15
BigReal y
Definition: Vector.h:66
Mass mass
Definition: NamdTypes.h:108
BigReal yy
Definition: Tensor.h:18
Bool is_atom_exPressure(int atomnum) const
Definition: Molecule.h:1447
SimParameters *const simParams
Definition: Sequencer.h:132
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
for(int i=0;i< n1;++i)
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::langevinVelocities ( BigReal  dt_fs)
protected

Definition at line 1317 of file Sequencer.C.

References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, adaptTempT, ResizeArray< T >::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().

1318 {
1319 // This routine is used for the BAOAB integrator,
1320 // Ornstein-Uhlenbeck exact solve for the O-part.
1321 // See B. Leimkuhler and C. Matthews, AMRX (2012)
1322 // Routine originally written by JPhillips, with fresh errors by CMatthews June2012
1323 
1325  {
1326  FullAtom *a = patch->atom.begin();
1327  int numAtoms = patch->numAtoms;
1328  Molecule *molecule = Node::Object()->molecule;
1329  BigReal dt = dt_fs * 0.001; // convert to ps
1332  {
1333  kbT = BOLTZMANN*adaptTempT;
1334  }
1335 
1336  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1337  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1338 
1339  for ( int i = 0; i < numAtoms; ++i )
1340  {
1341  BigReal dt_gamma = dt * a[i].langevinParam;
1342  if ( ! dt_gamma ) continue;
1343 
1344  BigReal f1 = exp( -dt_gamma );
1345  BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
1346  ( a[i].partition ? tempFactor : 1.0 ) *
1347  a[i].recipMass );
1348  a[i].velocity *= f1;
1349  a[i].velocity += f2 * random->gaussian_vector();
1350  }
1351  }
1352 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
Vector gaussian_vector(void)
Definition: Random.h:208
Real langevinParam
Definition: NamdTypes.h:110
#define BOLTZMANN
Definition: common.h:47
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Velocity velocity
Definition: NamdTypes.h:101
Bool langevin_useBAOAB
Bool adaptTempLangevin
BigReal langevinTemp
int numAtoms
Definition: Patch.h:144
BigReal adaptTempT
Definition: Sequencer.h:82
Random * random
Definition: Sequencer.h:131
SimParameters *const simParams
Definition: Sequencer.h:132
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::langevinVelocitiesBBK1 ( BigReal  dt_fs)
protected

Definition at line 1354 of file Sequencer.C.

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

1355 {
1356  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1357  NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1);
1359  {
1360  FullAtom *a = patch->atom.begin();
1361  int numAtoms = patch->numAtoms;
1362  Molecule *molecule = Node::Object()->molecule;
1363  BigReal dt = dt_fs * 0.001; // convert to ps
1364  int i;
1365 
1366  if (simParams->drudeOn) {
1367  for (i = 0; i < numAtoms; i++) {
1368 
1369  if (i < numAtoms-1 &&
1370  a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
1371  //printf("*** Found Drude particle %d\n", a[i+1].id);
1372  // i+1 is a Drude particle with parent i
1373 
1374  // convert from Cartesian coordinates to (COM,bond) coordinates
1375  BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
1376  Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
1377  Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
1378  BigReal dt_gamma;
1379 
1380  // use Langevin damping factor i for v_com
1381  dt_gamma = dt * a[i].langevinParam;
1382  if (dt_gamma != 0.0) {
1383  v_com *= ( 1. - 0.5 * dt_gamma );
1384  }
1385 
1386  // use Langevin damping factor i+1 for v_bnd
1387  dt_gamma = dt * a[i+1].langevinParam;
1388  if (dt_gamma != 0.0) {
1389  v_bnd *= ( 1. - 0.5 * dt_gamma );
1390  }
1391 
1392  // convert back
1393  a[i].velocity = v_com - m * v_bnd;
1394  a[i+1].velocity = v_bnd + a[i].velocity;
1395 
1396  i++; // +1 from loop, we've updated both particles
1397  }
1398  else {
1399  BigReal dt_gamma = dt * a[i].langevinParam;
1400  if ( ! dt_gamma ) continue;
1401 
1402  a[i].velocity *= ( 1. - 0.5 * dt_gamma );
1403  }
1404 
1405  } // end for
1406  } // end if drudeOn
1407  else {
1408 
1409  //
1410  // DJH: The conditional inside loop prevents vectorization and doesn't
1411  // avoid much work since addition and multiplication are cheap.
1412  //
1413  for ( i = 0; i < numAtoms; ++i )
1414  {
1415  BigReal dt_gamma = dt * a[i].langevinParam;
1416  if ( ! dt_gamma ) continue;
1417 
1418  a[i].velocity *= ( 1. - 0.5 * dt_gamma );
1419  }
1420 
1421  } // end else
1422 
1423  } // end if langevinOn
1424 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
Real langevinParam
Definition: NamdTypes.h:110
Definition: Vector.h:64
Velocity velocity
Definition: NamdTypes.h:101
Flags flags
Definition: Patch.h:127
Bool langevin_useBAOAB
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
Mass mass
Definition: NamdTypes.h:108
SimParameters *const simParams
Definition: Sequencer.h:132
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::langevinVelocitiesBBK2 ( BigReal  dt_fs)
protected

Definition at line 1427 of file Sequencer.C.

References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, adaptTempT, ResizeArray< T >::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().

1428 {
1429  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1430  NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2);
1432  {
1433  //
1434  // DJH: This call is expensive. Avoid calling when gammas don't differ.
1435  // Set flag in SimParameters and make this call conditional.
1436  //
1437  TIMER_START(patch->timerSet, RATTLE1);
1438  rattle1(dt_fs,1); // conserve momentum if gammas differ
1439  TIMER_STOP(patch->timerSet, RATTLE1);
1440 
1441  FullAtom *a = patch->atom.begin();
1442  int numAtoms = patch->numAtoms;
1443  Molecule *molecule = Node::Object()->molecule;
1444  BigReal dt = dt_fs * 0.001; // convert to ps
1447  {
1448  kbT = BOLTZMANN*adaptTempT;
1449  }
1450  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1451  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1452  int i;
1453 
1454  if (simParams->drudeOn) {
1455  BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp); // drude bond Temp
1456 
1457  for (i = 0; i < numAtoms; i++) {
1458 
1459  if (i < numAtoms-1 &&
1460  a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
1461  //printf("*** Found Drude particle %d\n", a[i+1].id);
1462  // i+1 is a Drude particle with parent i
1463 
1464  // convert from Cartesian coordinates to (COM,bond) coordinates
1465  BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass); // mass ratio
1466  Vector v_bnd = a[i+1].velocity - a[i].velocity; // vel of bond
1467  Vector v_com = a[i].velocity + m * v_bnd; // vel of COM
1468  BigReal dt_gamma;
1469 
1470  // use Langevin damping factor i for v_com
1471  dt_gamma = dt * a[i].langevinParam;
1472  if (dt_gamma != 0.0) {
1473  BigReal mass = a[i].mass + a[i+1].mass;
1474  v_com += random->gaussian_vector() *
1475  sqrt( 2 * dt_gamma * kbT *
1476  ( a[i].partition ? tempFactor : 1.0 ) / mass );
1477  v_com /= ( 1. + 0.5 * dt_gamma );
1478  }
1479 
1480  // use Langevin damping factor i+1 for v_bnd
1481  dt_gamma = dt * a[i+1].langevinParam;
1482  if (dt_gamma != 0.0) {
1483  BigReal mass = a[i+1].mass * (1. - m);
1484  v_bnd += random->gaussian_vector() *
1485  sqrt( 2 * dt_gamma * kbT_bnd *
1486  ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
1487  v_bnd /= ( 1. + 0.5 * dt_gamma );
1488  }
1489 
1490  // convert back
1491  a[i].velocity = v_com - m * v_bnd;
1492  a[i+1].velocity = v_bnd + a[i].velocity;
1493 
1494  i++; // +1 from loop, we've updated both particles
1495  }
1496  else {
1497  BigReal dt_gamma = dt * a[i].langevinParam;
1498  if ( ! dt_gamma ) continue;
1499 
1500  a[i].velocity += random->gaussian_vector() *
1501  sqrt( 2 * dt_gamma * kbT *
1502  ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
1503  a[i].velocity /= ( 1. + 0.5 * dt_gamma );
1504  }
1505 
1506  } // end for
1507  } // end if drudeOn
1508  else {
1509 
1510  //
1511  // DJH: For case using same gamma (the Langevin parameter),
1512  // no partitions (e.g. FEP), and no adaptive tempering (adaptTempMD),
1513  // we can precompute constants. Then by lifting the RNG from the
1514  // loop (filling up an array of random numbers), we can vectorize
1515  // loop and simplify arithmetic to just addition and multiplication.
1516  //
1517  for ( i = 0; i < numAtoms; ++i )
1518  {
1519  BigReal dt_gamma = dt * a[i].langevinParam;
1520  if ( ! dt_gamma ) continue;
1521 
1522  a[i].velocity += random->gaussian_vector() *
1523  sqrt( 2 * dt_gamma * kbT *
1524  ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
1525  a[i].velocity /= ( 1. + 0.5 * dt_gamma );
1526  }
1527 
1528  } // end else
1529 
1530  } // end if langevinOn
1531 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
Vector gaussian_vector(void)
Definition: Random.h:208
Real langevinParam
Definition: NamdTypes.h:110
#define BOLTZMANN
Definition: common.h:47
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Definition: Vector.h:64
void rattle1(BigReal, int)
Definition: Sequencer.C:1970
Velocity velocity
Definition: NamdTypes.h:101
Flags flags
Definition: Patch.h:127
Bool langevin_useBAOAB
#define TIMER_START(T, TYPE)
Definition: HomePatch.h:265
Bool adaptTempLangevin
BigReal langevinTemp
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
BigReal adaptTempT
Definition: Sequencer.h:82
Random * random
Definition: Sequencer.h:131
#define TIMER_STOP(T, TYPE)
Definition: HomePatch.h:266
Mass mass
Definition: NamdTypes.h:108
SimParameters *const simParams
Definition: Sequencer.h:132
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
BigReal drudeTemp
void Sequencer::maximumMove ( BigReal  timestep)
protected

Definition at line 2027 of file Sequencer.C.

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

2028 {
2029  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::MAXIMUM_MOVE);
2030 
2031  FullAtom *a = patch->atom.begin();
2032  int numAtoms = patch->numAtoms;
2033  if ( simParams->maximumMove ) {
2034  const BigReal dt = timestep / TIMEFACTOR;
2035  const BigReal maxvel = simParams->maximumMove / dt;
2036  const BigReal maxvel2 = maxvel * maxvel;
2037  for ( int i=0; i<numAtoms; ++i ) {
2038  if ( a[i].velocity.length2() > maxvel2 ) {
2039  a[i].velocity *= ( maxvel / a[i].velocity.length() );
2040  }
2041  }
2042  } else {
2043  const BigReal dt = timestep / TIMEFACTOR;
2044  const BigReal maxvel = simParams->cutoff / dt;
2045  const BigReal maxvel2 = maxvel * maxvel;
2046  int killme = 0;
2047  for ( int i=0; i<numAtoms; ++i ) {
2048  killme = killme || ( a[i].velocity.length2() > maxvel2 );
2049  }
2050  if ( killme ) {
2051  killme = 0;
2052  for ( int i=0; i<numAtoms; ++i ) {
2053  if ( a[i].velocity.length2() > maxvel2 ) {
2054  ++killme;
2055  iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
2056  << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
2057  << ( PDBVELFACTOR * maxvel ) << ", atom "
2058  << i << " of " << numAtoms << " on patch "
2059  << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
2060  }
2061  }
2062  iout << iERROR <<
2063  "Atoms moving too fast; simulation has become unstable ("
2064  << killme << " atoms on patch " << patch->patchID
2065  << " pe " << CkMyPe() << ").\n" << endi;
2067  terminate();
2068  }
2069  }
2070 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
void terminate(void)
Definition: Sequencer.C:2909
#define iout
Definition: InfoStream.h:87
Velocity velocity
Definition: NamdTypes.h:101
BigReal length(void) const
Definition: Vector.h:169
Flags flags
Definition: Patch.h:127
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
BigReal length2(void) const
Definition: Vector.h:173
const PatchID patchID
Definition: Patch.h:143
#define PDBVELFACTOR
Definition: common.h:50
#define TIMEFACTOR
Definition: common.h:48
BigReal maximumMove
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
void enableEarlyExit(void)
Definition: Node.C:1365
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::minimizationQuenchVelocity ( void  )
protected

Definition at line 2072 of file Sequencer.C.

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

Referenced by integrate().

2073 {
2074  if ( simParams->minimizeOn ) {
2075  FullAtom *a = patch->atom.begin();
2076  int numAtoms = patch->numAtoms;
2077  for ( int i=0; i<numAtoms; ++i ) {
2078  a[i].velocity = 0.;
2079  }
2080  }
2081 }
HomePatch *const patch
Definition: Sequencer.h:133
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
SimParameters *const simParams
Definition: Sequencer.h:132
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::minimize ( )
protected

Definition at line 760 of file Sequencer.C.

References Patch::atomMapper, ResizeArray< T >::begin(), broadcast, SimParameters::colvarsOn, ComputeMgr::computeGlobalObject, Node::computeMgr, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, Flags::doNonbonded, ResizeArray< T >::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, SimParameters::singleTopology, Results::slow, Flags::step, GlobalMaster::step, SimParameters::stepsPerCycle, submitCollections(), submitMinimizeReductions(), SimParameters::tclForcesOn, and TIMEFACTOR.

Referenced by algorithm().

760  {
761  //
762  // DJH: Copy all data into SOA (structure of arrays)
763  // from AOS (array of structures) data structure.
764  //
765  //patch->copy_all_to_SOA();
766 
767  const int numberOfSteps = simParams->N;
768  const int stepsPerCycle = simParams->stepsPerCycle;
769  int &step = patch->flags.step;
770  step = simParams->firstTimestep;
771 
772  int &maxForceUsed = patch->flags.maxForceUsed;
773  int &maxForceMerged = patch->flags.maxForceMerged;
774  maxForceUsed = Results::normal;
775  maxForceMerged = Results::normal;
776  int &doNonbonded = patch->flags.doNonbonded;
777  doNonbonded = 1;
778  maxForceUsed = Results::nbond;
779  maxForceMerged = Results::nbond;
780  const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
781  int &doFullElectrostatics = patch->flags.doFullElectrostatics;
782  doFullElectrostatics = dofull;
783  if ( dofull ) {
784  maxForceMerged = Results::slow;
785  maxForceUsed = Results::slow;
786  }
787  int &doMolly = patch->flags.doMolly;
788  doMolly = simParams->mollyOn && doFullElectrostatics;
789  // BEGIN LA
790  int &doLoweAndersen = patch->flags.doLoweAndersen;
791  doLoweAndersen = 0;
792  // END LA
793 
794  int &doGBIS = patch->flags.doGBIS;
795  doGBIS = simParams->GBISOn;
796 
797  int &doLCPO = patch->flags.doLCPO;
798  doLCPO = simParams->LCPOOn;
799 
800  int doTcl = simParams->tclForcesOn;
801  int doColvars = simParams->colvarsOn;
803 
804  int &doEnergy = patch->flags.doEnergy;
805  doEnergy = 1;
806 
807  patch->rattle1(0.,0,0); // enforce rigid bond constraints on initial positions
808 
811  patch->atom.begin(),patch->atom.end());
812  }
813 
814  runComputeObjects(1,step<numberOfSteps); // must migrate here!
815 
816  if ( doTcl || doColvars ) {
817  if ( doNonbonded ) saveForce(Results::nbond);
818  if ( doFullElectrostatics ) saveForce(Results::slow);
819  computeGlobal->saveTotalForces(patch);
820  }
821 
823 
824  submitMinimizeReductions(step,fmax2);
825  rebalanceLoad(step);
826 
827  int downhill = 1; // start out just fixing bad contacts
828  int minSeq = 0;
829  for ( ++step; step <= numberOfSteps; ++step ) {
830 
831  // Blocking receive for the minimization coefficient.
832  BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
833 
834  if ( downhill ) {
835  if ( c ) minimizeMoveDownhill(fmax2);
836  else {
837  downhill = 0;
838  fmax2 *= 10000.;
839  }
840  }
841  if ( ! downhill ) {
842  if ( ! c ) { // new direction
843 
844  // Blocking receive for the minimization coefficient.
845  c = broadcast->minimizeCoefficient.get(minSeq++);
846 
847  newMinimizeDirection(c); // v = c * v + f
848 
849  // Blocking receive for the minimization coefficient.
850  c = broadcast->minimizeCoefficient.get(minSeq++);
851 
852  } // same direction
853  newMinimizePosition(c); // x = x + c * v
854  }
855 
856  runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
857  if ( doTcl || doColvars ) {
858  if ( doNonbonded ) saveForce(Results::nbond);
859  if ( doFullElectrostatics ) saveForce(Results::slow);
860  computeGlobal->saveTotalForces(patch);
861  }
862  submitMinimizeReductions(step,fmax2);
863  submitCollections(step, 1); // write out zeros for velocities
864  rebalanceLoad(step);
865  }
866  quenchVelocities(); // zero out bogus velocity
867 
868  //
869  // DJH: Copy updates of SOA back into AOS.
870  //
871  //patch->copy_updates_to_AOS();
872 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
void saveForce(const int ftag=Results::normal)
Definition: Sequencer.C:1893
void registerIDsFullAtom(const FullAtom *begin, const FullAtom *end)
Definition: AtomMap.C:50
void newMinimizeDirection(BigReal)
Definition: Sequencer.C:897
void newMinimizePosition(BigReal)
Definition: Sequencer.C:956
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2382
void saveTotalForces(HomePatch *)
void rebalanceLoad(int timestep)
Definition: Sequencer.C:2878
int doLoweAndersen
Definition: PatchTypes.h:26
void minimizeMoveDownhill(BigReal fmax2)
Definition: Sequencer.C:875
AtomMapper * atomMapper
Definition: Patch.h:152
Flags flags
Definition: Patch.h:127
void submitCollections(int step, int zeroVel=0)
Definition: Sequencer.C:2643
int doEnergy
Definition: PatchTypes.h:20
int doFullElectrostatics
Definition: PatchTypes.h:23
iterator end(void)
Definition: ResizeArray.h:37
void runComputeObjects(int migration=1, int pairlists=0, int pressureStep=0)
Definition: Sequencer.C:2668
int maxForceUsed
Definition: PatchTypes.h:31
int doNonbonded
Definition: PatchTypes.h:22
SimpleBroadcastObject< BigReal > minimizeCoefficient
Definition: Broadcasts.h:78
ComputeGlobal * computeGlobalObject
Definition: ComputeMgr.h:97
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
int doLCPO
Definition: PatchTypes.h:29
void submitMinimizeReductions(int, BigReal fmax2)
Definition: Sequencer.C:2505
#define TIMEFACTOR
Definition: common.h:48
int doGBIS
Definition: PatchTypes.h:28
ComputeMgr * computeMgr
Definition: Node.h:169
int maxForceMerged
Definition: PatchTypes.h:32
void quenchVelocities()
Definition: Sequencer.C:984
SimParameters *const simParams
Definition: Sequencer.h:132
int doMolly
Definition: PatchTypes.h:24
double BigReal
Definition: common.h:114
int step
Definition: PatchTypes.h:16
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::minimizeMoveDownhill ( BigReal  fmax2)
protected

Definition at line 875 of file Sequencer.C.

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

Referenced by minimize().

875  {
876 
877  FullAtom *a = patch->atom.begin();
878  Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
879  int numAtoms = patch->numAtoms;
880 
881  for ( int i = 0; i < numAtoms; ++i ) {
882  if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
883  Force f = f1[i];
884  if ( f.length2() > fmax2 ) {
885  a[i].position += ( 0.1 * f.unit() );
886  int hgs = a[i].hydrogenGroupSize; // 0 if not parent
887  for ( int j=1; j<hgs; ++j ) {
888  a[++i].position += ( 0.1 * f.unit() );
889  }
890  }
891  }
892 
893  patch->rattle1(0.,0,0);
894 }
HomePatch *const patch
Definition: Sequencer.h:133
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
unsigned int atomFixed
Definition: NamdTypes.h:96
Definition: Vector.h:64
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2382
Position position
Definition: NamdTypes.h:53
int numAtoms
Definition: Patch.h:144
BigReal length2(void) const
Definition: Vector.h:173
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
SimParameters *const simParams
Definition: Sequencer.h:132
Vector unit(void) const
Definition: Vector.h:182
iterator begin(void)
Definition: ResizeArray.h:36
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< T >::begin(), broadcast, calcFixVirial(), Flags::doFullElectrostatics, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, SimpleBroadcastObject< T >::get(), CompAtom::hydrogenGroupSize, Tensor::identity(), SubmitReduction::item(), 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().

1084  {
1085 // Calculate new positions, momenta, and volume using positionRescaleFactor and
1086 // velocityRescaleTensor values returned from Controller::multigratorPressureCalcScale()
1088  FullAtom *a = patch->atom.begin();
1089  int numAtoms = patch->numAtoms;
1090 
1091  // Blocking receive (get) scaling factors from Controller
1092  Tensor scaleTensor = (callNumber == 1) ? broadcast->positionRescaleFactor.get(step) : broadcast->positionRescaleFactor2.get(step);
1093  Tensor velScaleTensor = (callNumber == 1) ? broadcast->velocityRescaleTensor.get(step) : broadcast->velocityRescaleTensor2.get(step);
1094  Tensor posScaleTensor = scaleTensor;
1095  posScaleTensor -= Tensor::identity();
1096  if (simParams->useGroupPressure) {
1097  velScaleTensor -= Tensor::identity();
1098  }
1099 
1100  // Scale volume
1101  patch->lattice.rescale(scaleTensor);
1102  // Scale positions and velocities
1103  scalePositionsVelocities(posScaleTensor, velScaleTensor);
1104 
1105  if (!patch->flags.doFullElectrostatics) NAMD_bug("Sequencer::multigratorPressure, doFullElectrostatics must be true");
1106 
1107  // Calculate new forces
1108  // NOTE: We should not need to migrate here since any migration should have happened in the
1109  // previous call to runComputeObjects inside the MD loop in Sequencer::integrate()
1110  const int numberOfSteps = simParams->N;
1111  const int stepsPerCycle = simParams->stepsPerCycle;
1112  runComputeObjects(0 , step<numberOfSteps, 1);
1113 
1114  reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
1116 
1117  // Virials etc.
1118  Tensor virialNormal;
1119  Tensor momentumSqrSum;
1120  BigReal kineticEnergy = 0;
1121  if ( simParams->pairInteractionOn ) {
1122  if ( simParams->pairInteractionSelf ) {
1123  for ( int i = 0; i < numAtoms; ++i ) {
1124  if ( a[i].partition != 1 ) continue;
1125  kineticEnergy += a[i].mass * a[i].velocity.length2();
1126  virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1127  }
1128  }
1129  } else {
1130  for ( int i = 0; i < numAtoms; ++i ) {
1131  if (a[i].mass < 0.01) continue;
1132  kineticEnergy += a[i].mass * a[i].velocity.length2();
1133  virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1134  }
1135  }
1136  if (!simParams->useGroupPressure) momentumSqrSum = virialNormal;
1137  kineticEnergy *= 0.5;
1139  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virialNormal);
1140 
1141  if ( simParams->fixedAtomsOn ) {
1142  Tensor fixVirialNormal;
1143  Tensor fixVirialNbond;
1144  Tensor fixVirialSlow;
1145  Vector fixForceNormal = 0;
1146  Vector fixForceNbond = 0;
1147  Vector fixForceSlow = 0;
1148 
1149  calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
1150 
1151  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, fixVirialNormal);
1152  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, fixVirialNbond);
1153  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, fixVirialSlow);
1154  ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NORMAL, fixForceNormal);
1155  ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NBOND, fixForceNbond);
1156  ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_SLOW, fixForceSlow);
1157  }
1158 
1159  // Internal virial and group momentum
1160  Tensor intVirialNormal;
1161  Tensor intVirialNormal2;
1162  Tensor intVirialNbond;
1163  Tensor intVirialSlow;
1164  int hgs;
1165  for ( int i = 0; i < numAtoms; i += hgs ) {
1166  hgs = a[i].hydrogenGroupSize;
1167  int j;
1168  BigReal m_cm = 0;
1169  Position x_cm(0,0,0);
1170  Velocity v_cm(0,0,0);
1171  for ( j = i; j < (i+hgs); ++j ) {
1172  m_cm += a[j].mass;
1173  x_cm += a[j].mass * a[j].position;
1174  v_cm += a[j].mass * a[j].velocity;
1175  }
1176  if (simParams->useGroupPressure) momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
1177  x_cm /= m_cm;
1178  v_cm /= m_cm;
1179  if (simParams->fixedAtomsOn) NAMD_bug("Sequencer::multigratorPressure, simParams->fixedAtomsOn not implemented yet");
1180  if ( simParams->pairInteractionOn ) {
1181  if ( simParams->pairInteractionSelf ) {
1182  NAMD_bug("Sequencer::multigratorPressure, this part needs to be implemented correctly");
1183  for ( j = i; j < (i+hgs); ++j ) {
1184  if ( a[j].partition != 1 ) continue;
1185  BigReal mass = a[j].mass;
1186  Vector v = a[j].velocity;
1187  Vector dv = v - v_cm;
1188  intVirialNormal2.outerAdd (mass, v, dv);
1189  Vector dx = a[j].position - x_cm;
1190  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
1191  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
1192  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
1193  }
1194  }
1195  } else {
1196  for ( j = i; j < (i+hgs); ++j ) {
1197  BigReal mass = a[j].mass;
1198  Vector v = a[j].velocity;
1199  Vector dv = v - v_cm;
1200  intVirialNormal2.outerAdd(mass, v, dv);
1201  Vector dx = a[j].position - x_cm;
1202  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
1203  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
1204  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
1205  }
1206  }
1207  }
1208 
1209  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
1210  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal2);
1211  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NBOND, intVirialNbond);
1212  ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_SLOW, intVirialSlow);
1213  ADD_TENSOR_OBJECT(reduction, REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
1214 
1215  reduction->submit();
1216  }
1217 }
HomePatch *const patch
Definition: Sequencer.h:133
SubmitReduction * reduction
Definition: Sequencer.h:134
int marginViolations
Definition: HomePatch.h:333
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
Lattice & lattice
Definition: Patch.h:126
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
BigReal & item(int i)
Definition: ReductionMgr.h:312
Position position
Definition: NamdTypes.h:53
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
Definition: Tensor.h:255
Velocity velocity
Definition: NamdTypes.h:101
Bool pairInteractionOn
Flags flags
Definition: Patch.h:127
void rescale(Tensor factor)
Definition: Lattice.h:60
void calcFixVirial(Tensor &fixVirialNormal, Tensor &fixVirialNbond, Tensor &fixVirialSlow, Vector &fixForceNormal, Vector &fixForceNbond, Vector &fixForceSlow)
Definition: Sequencer.C:2242
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int doFullElectrostatics
Definition: PatchTypes.h:23
void runComputeObjects(int migration=1, int pairlists=0, int pressureStep=0)
Definition: Sequencer.C:2668
SimpleBroadcastObject< Tensor > velocityRescaleTensor2
Definition: Broadcasts.h:72
int numAtoms
Definition: Patch.h:144
void scalePositionsVelocities(const Tensor &posScale, const Tensor &velScale)
Definition: Sequencer.C:1047
static Tensor identity(BigReal v1=1.0)
Definition: Tensor.h:31
SimpleBroadcastObject< Tensor > positionRescaleFactor
Definition: Broadcasts.h:69
BigReal length2(void) const
Definition: Vector.h:173
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
Definition: Tensor.h:15
Mass mass
Definition: NamdTypes.h:108
#define ADD_VECTOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:27
Bool pairInteractionSelf
int multigratorPressureFreq
void submit(void)
Definition: ReductionMgr.h:323
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
SimpleBroadcastObject< Tensor > positionRescaleFactor2
Definition: Broadcasts.h:74
SimParameters *const simParams
Definition: Sequencer.h:132
SimpleBroadcastObject< Tensor > velocityRescaleTensor
Definition: Broadcasts.h:71
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::multigratorTemperature ( int  step,
int  callNumber 
)
protected

Definition at line 1247 of file Sequencer.C.

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

1247  {
1249  // Blocking receive (get) velocity scaling factor.
1250  BigReal velScale = (callNumber == 1) ? broadcast->velocityRescaleFactor.get(step) : broadcast->velocityRescaleFactor2.get(step);
1251  scaleVelocities(velScale);
1252  // Calculate new kineticEnergy
1253  BigReal kineticEnergy = calcKineticEnergy();
1255  if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
1256  // If this is a pressure cycle, calculate new momentum squared sum
1257  FullAtom *a = patch->atom.begin();
1258  int numAtoms = patch->numAtoms;
1259  Tensor momentumSqrSum;
1260  if (simParams->useGroupPressure) {
1261  int hgs;
1262  for ( int i = 0; i < numAtoms; i += hgs ) {
1263  hgs = a[i].hydrogenGroupSize;
1264  int j;
1265  BigReal m_cm = 0;
1266  Position x_cm(0,0,0);
1267  Velocity v_cm(0,0,0);
1268  for ( j = i; j < (i+hgs); ++j ) {
1269  m_cm += a[j].mass;
1270  x_cm += a[j].mass * a[j].position;
1271  v_cm += a[j].mass * a[j].velocity;
1272  }
1273  momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
1274  }
1275  } else {
1276  for ( int i = 0; i < numAtoms; i++) {
1277  momentumSqrSum.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
1278  }
1279  }
1280  ADD_TENSOR_OBJECT(multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
1281  }
1282  // Submit reductions (kineticEnergy and, if applicable, momentumSqrSum)
1284 
1285  }
1286 }
HomePatch *const patch
Definition: Sequencer.h:133
SubmitReduction * multigratorReduction
Definition: Sequencer.h:119
void scaleVelocities(const BigReal velScale)
Definition: Sequencer.C:1219
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
BigReal & item(int i)
Definition: ReductionMgr.h:312
Position position
Definition: NamdTypes.h:53
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
Definition: Tensor.h:255
Velocity velocity
Definition: NamdTypes.h:101
BigReal calcKineticEnergy()
Definition: Sequencer.C:1227
int numAtoms
Definition: Patch.h:144
SimpleBroadcastObject< BigReal > velocityRescaleFactor2
Definition: Broadcasts.h:73
SimpleBroadcastObject< BigReal > velocityRescaleFactor
Definition: Broadcasts.h:68
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
Definition: Tensor.h:15
Mass mass
Definition: NamdTypes.h:108
int multigratorPressureFreq
void submit(void)
Definition: ReductionMgr.h:323
SimParameters *const simParams
Definition: Sequencer.h:132
int multigratorTemperatureFreq
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::newMinimizeDirection ( BigReal  c)
protected

Definition at line 897 of file Sequencer.C.

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

Referenced by minimize().

897  {
898  FullAtom *a = patch->atom.begin();
899  Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
900  const bool fixedAtomsOn = simParams->fixedAtomsOn;
901  const bool drudeHardWallOn = simParams->drudeHardWallOn;
902  int numAtoms = patch->numAtoms;
903  BigReal maxv2 = 0.;
904 
905  for ( int i = 0; i < numAtoms; ++i ) {
906  a[i].velocity *= c;
907  a[i].velocity += f1[i];
908  if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
909  a[i].velocity = a[i-1].velocity;
910  }
911  if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
912  BigReal v2 = a[i].velocity.length2();
913  if ( v2 > maxv2 ) maxv2 = v2;
914  }
915 
916  { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial); }
917 
918  maxv2 = 0.;
919  for ( int i = 0; i < numAtoms; ++i ) {
920  if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
921  a[i].velocity = a[i-1].velocity;
922  }
923  if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
924  BigReal v2 = a[i].velocity.length2();
925  if ( v2 > maxv2 ) maxv2 = v2;
926  }
927 
928  min_reduction->max(0,maxv2);
930 
931  // prevent hydrogens from being left behind
932  BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
933  // int adjustCount = 0;
934  int hgs;
935  for ( int i = 0; i < numAtoms; i += hgs ) {
936  hgs = a[i].hydrogenGroupSize;
937  BigReal minChildVel = a[i].velocity.length2();
938  if ( minChildVel < fmax2 ) continue;
939  int adjustChildren = 1;
940  for ( int j = i+1; j < (i+hgs); ++j ) {
941  if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
942  }
943  if ( adjustChildren ) {
944  // if ( hgs > 1 ) ++adjustCount;
945  for ( int j = i+1; j < (i+hgs); ++j ) {
946  if (a[i].mass < 0.01) continue; // lone pair
947  a[j].velocity = a[i].velocity;
948  }
949  }
950  }
951  // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
952 
953 }
HomePatch *const patch
Definition: Sequencer.h:133
void max(int i, BigReal v)
Definition: ReductionMgr.h:315
void minimize_rattle2(const BigReal, Tensor *virial, bool forces=false)
Definition: HomePatch.C:2963
SubmitReduction * min_reduction
Definition: Sequencer.h:39
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
Definition: Vector.h:64
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
BigReal length2(void) const
Definition: Vector.h:173
Definition: Tensor.h:15
#define TIMEFACTOR
Definition: common.h:48
void submit(void)
Definition: ReductionMgr.h:323
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::newMinimizePosition ( BigReal  c)
protected

Definition at line 956 of file Sequencer.C.

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

Referenced by minimize().

956  {
957  FullAtom *a = patch->atom.begin();
958  int numAtoms = patch->numAtoms;
959 
960  for ( int i = 0; i < numAtoms; ++i ) {
961  a[i].position += c * a[i].velocity;
962  }
963 
964  if ( simParams->drudeHardWallOn ) {
965  for ( int i = 1; i < numAtoms; ++i ) {
966  if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
967  a[i].position -= a[i-1].position;
968  }
969  }
970  }
971 
972  patch->rattle1(0.,0,0);
973 
974  if ( simParams->drudeHardWallOn ) {
975  for ( int i = 1; i < numAtoms; ++i ) {
976  if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
977  a[i].position += a[i-1].position;
978  }
979  }
980  }
981 }
HomePatch *const patch
Definition: Sequencer.h:133
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2382
Position position
Definition: NamdTypes.h:53
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
SimParameters *const simParams
Definition: Sequencer.h:132
iterator begin(void)
Definition: ResizeArray.h:36
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().

1299 {
1300  NAMD_EVENT_RANGE_2(patch->flags.event_on,
1301  NamdProfileEvent::NEWTONIAN_VELOCITIES);
1302 
1303  // Deterministic velocity update, account for multigrator
1304  if (staleForces || (doNonbonded && doFullElectrostatics)) {
1305  addForceToMomentum3(stepscale*timestep, Results::normal, 0,
1306  stepscale*nbondstep, Results::nbond, staleForces,
1307  stepscale*slowstep, Results::slow, staleForces);
1308  } else {
1309  addForceToMomentum(stepscale*timestep);
1310  if (staleForces || doNonbonded)
1311  addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces);
1312  if (staleForces || doFullElectrostatics)
1313  addForceToMomentum(stepscale*slowstep, Results::slow, staleForces);
1314  }
1315 }
HomePatch *const patch
Definition: Sequencer.h:133
void addForceToMomentum(BigReal, const int ftag=Results::normal, const int useSaved=0)
Definition: Sequencer.C:1904
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)
Definition: Sequencer.C:1919
Flags flags
Definition: Patch.h:127
#define NAMD_EVENT_RANGE_2(eon, id)
void Sequencer::quenchVelocities ( )
protected

Definition at line 984 of file Sequencer.C.

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

Referenced by minimize().

984  {
985  FullAtom *a = patch->atom.begin();
986  int numAtoms = patch->numAtoms;
987 
988  for ( int i = 0; i < numAtoms; ++i ) {
989  a[i].velocity = 0;
990  }
991 }
HomePatch *const patch
Definition: Sequencer.h:133
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::rattle1 ( BigReal  dt,
int  pressure 
)
protected

Definition at line 1970 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ResizeArray< T >::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().

1971 {
1972  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::RATTLE1);
1973  if ( simParams->rigidBonds != RIGID_NONE ) {
1974  Tensor virial;
1975  Tensor *vp = ( pressure ? &virial : 0 );
1976  if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
1977  iout << iERROR <<
1978  "Constraint failure; simulation has become unstable.\n" << endi;
1980  terminate();
1981  }
1982 #if 0
1983  printf("virial = %g %g %g %g %g %g %g %g %g\n",
1984  virial.xx, virial.xy, virial.xz,
1985  virial.yx, virial.yy, virial.yz,
1986  virial.zx, virial.zy, virial.zz);
1987 #endif
1988  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
1989 #if 0
1990  {
1991  const FullAtom *a = patch->atom.const_begin();
1992  for (int n=0; n < patch->numAtoms; n++) {
1993  printf("pos[%d] = %g %g %g\n", n,
1994  a[n].position.x, a[n].position.y, a[n].position.z);
1995  }
1996  for (int n=0; n < patch->numAtoms; n++) {
1997  printf("vel[%d] = %g %g %g\n", n,
1998  a[n].velocity.x, a[n].velocity.y, a[n].velocity.z);
1999  }
2000  if (pressure) {
2001  for (int n=0; n < patch->numAtoms; n++) {
2002  printf("force[%d] = %g %g %g\n", n,
2003  patch->f[Results::normal][n].x,
2004  patch->f[Results::normal][n].y,
2005  patch->f[Results::normal][n].z);
2006  }
2007  }
2008  }
2009 #endif
2010  }
2011 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
BigReal zy
Definition: Tensor.h:19
void terminate(void)
Definition: Sequencer.C:2909
SubmitReduction * pressureProfileReduction
Definition: Sequencer.h:135
SubmitReduction * reduction
Definition: Sequencer.h:134
BigReal xz
Definition: Tensor.h:17
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
int rattle1(const BigReal, Tensor *virial, SubmitReduction *)
Definition: HomePatch.C:2382
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
BigReal yz
Definition: Tensor.h:18
#define iout
Definition: InfoStream.h:87
Velocity velocity
Definition: NamdTypes.h:101
Flags flags
Definition: Patch.h:127
BigReal yx
Definition: Tensor.h:18
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
BigReal xx
Definition: Tensor.h:17
BigReal zz
Definition: Tensor.h:19
Definition: Tensor.h:15
BigReal xy
Definition: Tensor.h:17
BigReal y
Definition: Vector.h:66
BigReal yy
Definition: Tensor.h:18
const_iterator const_begin(void) const
Definition: ResizeArray.h:39
std::ostream & iERROR(std::ostream &s)
Definition: InfoStream.C:109
infostream & endi(infostream &s)
Definition: InfoStream.C:38
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
void enableEarlyExit(void)
Definition: Node.C:1365
#define RIGID_NONE
Definition: SimParameters.h:77
SimParameters *const simParams
Definition: Sequencer.h:132
BigReal zx
Definition: Tensor.h:19
void Sequencer::reassignVelocities ( BigReal  timestep,
int  step 
)
protected

Definition at line 1756 of file Sequencer.C.

References CompAtomExt::atomFixed, ResizeArray< T >::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().

1757 {
1758  const int reassignFreq = simParams->reassignFreq;
1759  if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
1760  FullAtom *a = patch->atom.begin();
1761  int numAtoms = patch->numAtoms;
1762  BigReal newTemp = simParams->reassignTemp;
1763  newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
1764  if ( simParams->reassignIncr > 0.0 ) {
1765  if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
1766  newTemp = simParams->reassignHold;
1767  } else {
1768  if ( newTemp < simParams->reassignHold )
1769  newTemp = simParams->reassignHold;
1770  }
1771  BigReal kbT = BOLTZMANN * newTemp;
1772 
1773  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1774  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1775 
1776  for ( int i = 0; i < numAtoms; ++i )
1777  {
1778  a[i].velocity = ( ( simParams->fixedAtomsOn &&
1779  a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
1780  sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) *
1781  random->gaussian_vector() );
1782  }
1783  } else {
1784  NAMD_bug("Sequencer::reassignVelocities called improperly!");
1785  }
1786 }
HomePatch *const patch
Definition: Sequencer.h:133
#define BOLTZMANN
Definition: common.h:47
unsigned int atomFixed
Definition: NamdTypes.h:96
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Definition: Vector.h:64
BigReal reassignTemp
if(ComputeNonbondedUtil::goMethod==2)
Velocity velocity
Definition: NamdTypes.h:101
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int numAtoms
Definition: Patch.h:144
BigReal reassignIncr
Random * random
Definition: Sequencer.h:131
Mass mass
Definition: NamdTypes.h:108
BigReal reassignHold
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::rebalanceLoad ( int  timestep)
protected

Definition at line 2878 of file Sequencer.C.

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

Referenced by integrate(), and minimize().

2878  {
2879  if ( ! ldbSteps ) {
2881  }
2882  if ( ! --ldbSteps ) {
2883  patch->submitLoadStats(timestep);
2884  ldbCoordinator->rebalance(this,patch->getPatchID());
2885  pairlistsAreValid = 0;
2886  }
2887 }
HomePatch *const patch
Definition: Sequencer.h:133
void submitLoadStats(int timestep)
Definition: HomePatch.C:3620
void rebalance(Sequencer *seq, PatchID id)
int ldbSteps
Definition: Sequencer.h:140
static LdbCoordinator * Object()
PatchID getPatchID()
Definition: Patch.h:114
int getNumStepsToRun(void)
int pairlistsAreValid
Definition: Sequencer.h:42
void Sequencer::reinitVelocities ( void  )
protected

Definition at line 1788 of file Sequencer.C.

References CompAtomExt::atomFixed, ResizeArray< T >::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().

1789 {
1790  FullAtom *a = patch->atom.begin();
1791  int numAtoms = patch->numAtoms;
1792  BigReal newTemp = simParams->initialTemp;
1793  BigReal kbT = BOLTZMANN * newTemp;
1794 
1795  int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
1796  BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
1797 
1798  for ( int i = 0; i < numAtoms; ++i )
1799  {
1800  a[i].velocity = ( ( (simParams->fixedAtomsOn && a[i].atomFixed) ||
1801  a[i].mass <= 0.) ? Vector(0,0,0) :
1802  sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) *
1803  random->gaussian_vector() );
1804  if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass > 0.05 ) {
1805  a[i+1].velocity = a[i].velocity; // zero is good enough
1806  ++i;
1807  }
1808  }
1809 }
HomePatch *const patch
Definition: Sequencer.h:133
#define BOLTZMANN
Definition: common.h:47
unsigned int atomFixed
Definition: NamdTypes.h:96
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Definition: Vector.h:64
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
Random * random
Definition: Sequencer.h:131
BigReal initialTemp
Mass mass
Definition: NamdTypes.h:108
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::reloadCharges ( )
protected

Definition at line 1821 of file Sequencer.C.

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

Referenced by algorithm().

1822 {
1823  FullAtom *a = patch->atom.begin();
1824  int numAtoms = patch->numAtoms;
1825  Molecule *molecule = Node::Object()->molecule;
1826  for ( int i = 0; i < numAtoms; ++i )
1827  {
1828  a[i].charge = molecule->atomcharge(a[i].id);
1829  }
1830 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
Charge charge
Definition: NamdTypes.h:54
int numAtoms
Definition: Patch.h:144
Real atomcharge(int anum) const
Definition: Molecule.h:1048
Molecule * molecule
Definition: Node.h:176
iterator begin(void)
Definition: ResizeArray.h:36
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().

1698 {
1699  if (!simParams->accelMDOn) return;
1700  if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
1701 
1702  // Blocking receive for the Accelerated MD scaling factors.
1703  Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
1704  const BigReal factor_dihe = accelMDfactor[0];
1705  const BigReal factor_tot = accelMDfactor[1];
1706  const int numAtoms = patch->numAtoms;
1707 
1708  if (simParams->accelMDdihe && factor_tot <1 )
1709  NAMD_die("accelMD broadcasting error!\n");
1710  if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
1711  NAMD_die("accelMD broadcasting error!\n");
1712 
1713  if (simParams->accelMDdihe && factor_dihe < 1) {
1714  for (int i = 0; i < numAtoms; ++i)
1715  if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
1716  patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);
1717  }
1718 
1719  if ( !simParams->accelMDdihe && factor_tot < 1) {
1720  for (int i = 0; i < numAtoms; ++i)
1721  patch->f[Results::normal][i] *= factor_tot;
1722  if (doNonbonded) {
1723  for (int i = 0; i < numAtoms; ++i)
1724  patch->f[Results::nbond][i] *= factor_tot;
1725  }
1726  if (doFullElectrostatics) {
1727  for (int i = 0; i < numAtoms; ++i)
1728  patch->f[Results::slow][i] *= factor_tot;
1729  }
1730  }
1731 
1732  if (simParams->accelMDdual && factor_dihe < 1) {
1733  for (int i = 0; i < numAtoms; ++i)
1734  if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
1735  patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
1736  }
1737 
1738 }
HomePatch *const patch
Definition: Sequencer.h:133
BigReal accelMDLastStep
Definition: Vector.h:64
if(ComputeNonbondedUtil::goMethod==2)
int numAtoms
Definition: Patch.h:144
void NAMD_die(const char *err_msg)
Definition: common.C:83
SimpleBroadcastObject< Vector > accelMDRescaleFactor
Definition: Broadcasts.h:85
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
void Sequencer::rescaleSoluteCharges ( BigReal  factor)
protected

Definition at line 1833 of file Sequencer.C.

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

Referenced by algorithm(), and Sequencer().

1834 {
1835  FullAtom *a = patch->atom.begin();
1836  int numAtoms = patch->numAtoms;
1837  Molecule *molecule = Node::Object()->molecule;
1838  BigReal sqrt_factor = sqrt(factor);
1839  // apply scaling to the original charge (stored in molecule)
1840  // of just the marked solute atoms
1841  for ( int i = 0; i < numAtoms; ++i ) {
1842  if (molecule->get_ss_type(a[i].id)) {
1843  a[i].charge = sqrt_factor * molecule->atomcharge(a[i].id);
1844  }
1845  }
1846 }
static Node * Object()
Definition: Node.h:86
HomePatch *const patch
Definition: Sequencer.h:133
Charge charge
Definition: NamdTypes.h:54
unsigned char get_ss_type(int anum) const
Definition: Molecule.h:1351
int numAtoms
Definition: Patch.h:144
Real atomcharge(int anum) const
Definition: Molecule.h:1048
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::rescaleVelocities ( int  step)
protected

Definition at line 1678 of file Sequencer.C.

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

Referenced by integrate().

1679 {
1680  const int rescaleFreq = simParams->rescaleFreq;
1681  if ( rescaleFreq > 0 ) {
1682  FullAtom *a = patch->atom.begin();
1683  int numAtoms = patch->numAtoms;
1685  if ( rescaleVelocities_numTemps == rescaleFreq ) {
1686  // Blocking receive for the velcity scaling factor.
1687  BigReal factor = broadcast->velocityRescaleFactor.get(step);
1688  for ( int i = 0; i < numAtoms; ++i )
1689  {
1690  a[i].velocity *= factor;
1691  }
1693  }
1694  }
1695 }
HomePatch *const patch
Definition: Sequencer.h:133
Velocity velocity
Definition: NamdTypes.h:101
int rescaleVelocities_numTemps
Definition: Sequencer.h:87
int numAtoms
Definition: Patch.h:144
SimpleBroadcastObject< BigReal > velocityRescaleFactor
Definition: Broadcasts.h:68
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::rescaleVelocitiesByFactor ( BigReal  factor)
protected

Definition at line 1811 of file Sequencer.C.

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

Referenced by algorithm().

1812 {
1813  FullAtom *a = patch->atom.begin();
1814  int numAtoms = patch->numAtoms;
1815  for ( int i = 0; i < numAtoms; ++i )
1816  {
1817  a[i].velocity *= factor;
1818  }
1819 }
HomePatch *const patch
Definition: Sequencer.h:133
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
iterator begin(void)
Definition: ResizeArray.h:36
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().

127 {
128  // create a Thread and invoke it
129  DebugM(4, "::run() - this = " << this << "\n" );
130  thread = CthCreate((CthVoidFn)&(threadRun),(void*)(this),SEQ_STK_SZ);
131  CthSetStrategyDefault(thread);
132  priority = PATCH_PRIORITY(patch->getPatchID());
133  awaken();
134 }
HomePatch *const patch
Definition: Sequencer.h:133
#define DebugM(x, y)
Definition: Debug.h:59
#define SEQ_STK_SZ
Definition: Thread.h:11
void awaken(void)
Definition: Sequencer.h:29
PatchID getPatchID()
Definition: Patch.h:114
#define PATCH_PRIORITY(PID)
Definition: Priorities.h:25
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().

2669 {
2670  if ( migration ) pairlistsAreValid = 0;
2671 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
2672  if ( pairlistsAreValid &&
2674  && ( pairlistsAge > (
2675 #else
2676  if ( pairlistsAreValid && ( pairlistsAge > (
2677 #endif
2679  pairlistsAreValid = 0;
2680  }
2681  if ( ! simParams->usePairlists ) pairlists = 0;
2682  patch->flags.usePairlists = pairlists || pairlistsAreValid;
2684  pairlists && ! pairlistsAreValid;
2685 
2686  if ( simParams->singleTopology ) patch->reposition_all_alchpairs();
2687  if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
2688 
2689  //
2690  // DJH: Copy updates of SOA back into AOS.
2691  // The positionsReady() routine starts force computation and atom migration.
2692  //
2693  // We could reduce amount of copying here by checking migration status
2694  // and copying velocities only when migrating. Some types of simulation
2695  // always require velocities, such as Lowe-Anderson.
2696  //
2697  //patch->copy_updates_to_AOS();
2698 
2699  patch->positionsReady(migration); // updates flags.sequence
2700 
2701  int seq = patch->flags.sequence;
2702  int basePriority = ( (seq & 0xffff) << 15 )
2704  if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
2705  priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
2706  suspend(); // until all deposit boxes close
2708  priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
2709  suspend();
2711  priority = basePriority + COMPUTE_HOME_PRIORITY;
2712  suspend();
2713  } else {
2714  priority = basePriority + COMPUTE_HOME_PRIORITY;
2715  suspend(); // until all deposit boxes close
2716  }
2717 
2718  //
2719  // DJH: Copy all data into SOA from AOS.
2720  //
2721  // We need everything copied after atom migration.
2722  // When doing force computation without atom migration,
2723  // all data except forces will already be up-to-date in SOA
2724  // (except maybe for some special types of simulation).
2725  //
2726  //patch->copy_all_to_SOA();
2727 
2728  //
2729  // DJH: Copy forces to SOA.
2730  // Force available after suspend() has returned.
2731  //
2732  //patch->copy_forces_to_SOA();
2733 
2735  pairlistsAreValid = 1;
2736  pairlistsAge = 0;
2737  }
2738  // For multigrator, do not age pairlist during pressure step
2739  // NOTE: for non-multigrator pressureStep = 0 always
2740  if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
2741 
2742  if (simParams->lonepairs) {
2743  {
2744  Tensor virial;
2745  patch->redistrib_lonepair_forces(Results::normal, &virial);
2746  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
2747  }
2748  if (patch->flags.doNonbonded) {
2749  Tensor virial;
2750  patch->redistrib_lonepair_forces(Results::nbond, &virial);
2751  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
2752  }
2754  Tensor virial;
2755  patch->redistrib_lonepair_forces(Results::slow, &virial);
2756  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
2757  }
2758  } else if (simParams->watmodel == WAT_TIP4) {
2759  {
2760  Tensor virial;
2761  patch->redistrib_tip4p_forces(Results::normal, &virial);
2762  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
2763  }
2764  if (patch->flags.doNonbonded) {
2765  Tensor virial;
2766  patch->redistrib_tip4p_forces(Results::nbond, &virial);
2767  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
2768  }
2770  Tensor virial;
2771  patch->redistrib_tip4p_forces(Results::slow, &virial);
2772  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
2773  }
2774  } else if (simParams->watmodel == WAT_SWM4) {
2775  {
2776  Tensor virial;
2777  patch->redistrib_swm4_forces(Results::normal, &virial);
2778  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
2779  }
2780  if (patch->flags.doNonbonded) {
2781  Tensor virial;
2782  patch->redistrib_swm4_forces(Results::nbond, &virial);
2783  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
2784  }
2786  Tensor virial;
2787  patch->redistrib_swm4_forces(Results::slow, &virial);
2788  ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
2789  }
2790  }
2791 
2792  if (simParams->singleTopology) {
2793  patch->redistrib_alchpair_forces(Results::normal);
2794  if (patch->flags.doNonbonded) {
2795  patch->redistrib_alchpair_forces(Results::nbond);
2796  }
2798  patch->redistrib_alchpair_forces(Results::slow);
2799  }
2800  }
2801 
2802  if ( patch->flags.doMolly ) {
2803  Tensor virial;
2804  patch->mollyMollify(&virial);
2805  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
2806  }
2807 
2808 
2809  // BEGIN LA
2810  if (patch->flags.doLoweAndersen) {
2812  }
2813  // END LA
2814 #ifdef NAMD_CUDA_XXX
2815  int numAtoms = patch->numAtoms;
2816  FullAtom *a = patch->atom.begin();
2817  for ( int i=0; i<numAtoms; ++i ) {
2818  CkPrintf("%d %g %g %g\n", a[i].id,
2819  patch->f[Results::normal][i].x +
2820  patch->f[Results::nbond][i].x +
2821  patch->f[Results::slow][i].x,
2822  patch->f[Results::normal][i].y +
2823  patch->f[Results::nbond][i].y +
2824  patch->f[Results::slow][i].y,
2825  patch->f[Results::normal][i].z +
2826  patch->f[Results::nbond][i].z +
2827  patch->f[Results::slow][i].z);
2828  CkPrintf("%d %g %g %g\n", a[i].id,
2829  patch->f[Results::normal][i].x,
2830  patch->f[Results::nbond][i].x,
2831  patch->f[Results::slow][i].x);
2832  CkPrintf("%d %g %g %g\n", a[i].id,
2833  patch->f[Results::normal][i].y,
2834  patch->f[Results::nbond][i].y,
2835  patch->f[Results::slow][i].y);
2836  CkPrintf("%d %g %g %g\n", a[i].id,
2837  patch->f[Results::normal][i].z,
2838  patch->f[Results::nbond][i].z,
2839  patch->f[Results::slow][i].z);
2840  }
2841 #endif
2842 
2843 #if PRINT_FORCES
2844  int numAtoms = patch->numAtoms;
2845  FullAtom *a = patch->atom.begin();
2846  for ( int i=0; i<numAtoms; ++i ) {
2847  float fxNo = patch->f[Results::normal][i].x;
2848  float fxNb = patch->f[Results::nbond][i].x;
2849  float fxSl = patch->f[Results::slow][i].x;
2850  float fyNo = patch->f[Results::normal][i].y;
2851  float fyNb = patch->f[Results::nbond][i].y;
2852  float fySl = patch->f[Results::slow][i].y;
2853  float fzNo = patch->f[Results::normal][i].z;
2854  float fzNb = patch->f[Results::nbond][i].z;
2855  float fzSl = patch->f[Results::slow][i].z;
2856  float fx = fxNo+fxNb+fxSl;
2857  float fy = fyNo+fyNb+fySl;
2858  float fz = fzNo+fzNb+fzSl;
2859 
2860  float f = sqrt(fx*fx+fy*fy+fz*fz);
2861  int id = patch->pExt[i].id;
2862  int seq = patch->flags.sequence;
2863  float x = patch->p[i].position.x;
2864  float y = patch->p[i].position.y;
2865  float z = patch->p[i].position.z;
2866  //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
2867  CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
2868  //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
2869 //fxNo,fyNo,fzNo,
2870 //fxNb,fyNb,fzNb,
2871 //fxSl,fySl,fzSl,
2872 fx,fy,fz
2873 );
2874  }
2875 #endif
2876 }
HomePatch *const patch
Definition: Sequencer.h:133
#define GB1_COMPUTE_HOME_PRIORITY
Definition: Priorities.h:56
SubmitReduction * reduction
Definition: Sequencer.h:134
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
int savePairlists
Definition: PatchTypes.h:39
void gbisComputeAfterP2()
Definition: HomePatch.C:3175
int usePairlists
Definition: PatchTypes.h:38
if(ComputeNonbondedUtil::goMethod==2)
int doLoweAndersen
Definition: PatchTypes.h:26
void positionsReady(int doMigration=0)
Definition: HomePatch.C:925
Flags flags
Definition: Patch.h:127
int pairlistsAge
Definition: Sequencer.h:43
#define COMPUTE_HOME_PRIORITY
Definition: Priorities.h:76
int doFullElectrostatics
Definition: PatchTypes.h:23
#define WAT_SWM4
Definition: common.h:192
void mollyMollify(Tensor *virial)
Definition: HomePatch.C:3387
gridSize z
CompAtomList p
Definition: Patch.h:146
int numAtoms
Definition: Patch.h:144
int sequence
Definition: PatchTypes.h:18
void gbisComputeAfterP1()
Definition: HomePatch.C:3147
int doNonbonded
Definition: PatchTypes.h:22
PatchID getPatchID()
Definition: Patch.h:114
void suspend(void)
Definition: Sequencer.C:136
#define GB2_COMPUTE_HOME_PRIORITY
Definition: Priorities.h:64
Definition: Tensor.h:15
gridSize y
#define WAT_TIP4
Definition: common.h:191
int pairlistsAreValid
Definition: Sequencer.h:42
int doGBIS
Definition: PatchTypes.h:28
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
gridSize x
void loweAndersenFinish()
Definition: HomePatch.C:3113
SimParameters *const simParams
Definition: Sequencer.h:132
CompAtomExtList pExt
Definition: Patch.h:174
int doMolly
Definition: PatchTypes.h:24
#define PATCH_PRIORITY(PID)
Definition: Priorities.h:25
iterator begin(void)
Definition: ResizeArray.h:36
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().

1894 {
1895  patch->saveForce(ftag);
1896 }
HomePatch *const patch
Definition: Sequencer.h:133
void saveForce(const int ftag=Results::normal)
Definition: HomePatch.C:1336
void Sequencer::scalePositionsVelocities ( const Tensor posScale,
const Tensor velScale 
)
protected

Definition at line 1047 of file Sequencer.C.

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

Referenced by multigratorPressure().

1047  {
1048  FullAtom *a = patch->atom.begin();
1049  int numAtoms = patch->numAtoms;
1050  Position origin = patch->lattice.origin();
1051  if ( simParams->fixedAtomsOn ) {
1052  NAMD_bug("Sequencer::scalePositionsVelocities, fixed atoms not implemented");
1053  }
1054  if ( simParams->useGroupPressure ) {
1055  int hgs;
1056  for ( int i = 0; i < numAtoms; i += hgs ) {
1057  hgs = a[i].hydrogenGroupSize;
1058  Position pos_cm(0.0, 0.0, 0.0);
1059  Velocity vel_cm(0.0, 0.0, 0.0);
1060  BigReal m_cm = 0.0;
1061  for (int j=0;j < hgs;++j) {
1062  m_cm += a[i+j].mass;
1063  pos_cm += a[i+j].mass*a[i+j].position;
1064  vel_cm += a[i+j].mass*a[i+j].velocity;
1065  }
1066  pos_cm /= m_cm;
1067  vel_cm /= m_cm;
1068  pos_cm -= origin;
1069  Position dpos = posScale*pos_cm;
1070  Velocity dvel = velScale*vel_cm;
1071  for (int j=0;j < hgs;++j) {
1072  a[i+j].position += dpos;
1073  a[i+j].velocity += dvel;
1074  }
1075  }
1076  } else {
1077  for ( int i = 0; i < numAtoms; i++) {
1078  a[i].position += posScale*(a[i].position-origin);
1079  a[i].velocity = velScale*a[i].velocity;
1080  }
1081  }
1082 }
HomePatch *const patch
Definition: Sequencer.h:133
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
Lattice & lattice
Definition: Patch.h:126
Definition: Vector.h:64
Position position
Definition: NamdTypes.h:53
Velocity velocity
Definition: NamdTypes.h:101
Vector origin() const
Definition: Lattice.h:262
void NAMD_bug(const char *err_msg)
Definition: common.C:123
int numAtoms
Definition: Patch.h:144
Mass mass
Definition: NamdTypes.h:108
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
void Sequencer::scaleVelocities ( const BigReal  velScale)
protected

Definition at line 1219 of file Sequencer.C.

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

Referenced by multigratorTemperature().

1219  {
1220  FullAtom *a = patch->atom.begin();
1221  int numAtoms = patch->numAtoms;
1222  for ( int i = 0; i < numAtoms; i++) {
1223  a[i].velocity *= velScale;
1224  }
1225 }
HomePatch *const patch
Definition: Sequencer.h:133
Velocity velocity
Definition: NamdTypes.h:101
int numAtoms
Definition: Patch.h:144
iterator begin(void)
Definition: ResizeArray.h:36
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_fsThe integration increment, in fs (not currently used)
stepThe current timestep

Definition at line 1872 of file Sequencer.C.

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

Referenced by integrate().

1873 {
1874  if ( simParams->stochRescaleOn )
1875  {
1876  FullAtom *a = patch->atom.begin();
1877  int numAtoms = patch->numAtoms;
1880  {
1881  // Blocking receive for the temperature coupling coefficient.
1882  BigReal coefficient = broadcast->stochRescaleCoefficient.get(step);
1883  DebugM(4, "stochastically rescaling velocities at step " << step << " by " << coefficient << "\n");
1884  for ( int i = 0; i < numAtoms; ++i )
1885  {
1886  a[i].velocity *= coefficient;
1887  }
1888  stochRescale_count = 0;
1889  }
1890  }
1891 }
HomePatch *const patch
Definition: Sequencer.h:133
#define DebugM(x, y)
Definition: Debug.h:59
Velocity velocity
Definition: NamdTypes.h:101
SimpleBroadcastObject< BigReal > stochRescaleCoefficient
Definition: Broadcasts.h:77
int numAtoms
Definition: Patch.h:144
ControllerBroadcasts * broadcast
Definition: Sequencer.h:138
int stochRescale_count
Definition: Sequencer.h:100
SimParameters *const simParams
Definition: Sequencer.h:132
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
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().

2644 {
2645  //
2646  // DJH: Copy updates of SOA back into AOS.
2647  // Do we need to update everything or is it safe to just update
2648  // positions and velocities separately, as needed?
2649  //
2650  //patch->copy_updates_to_AOS();
2651 
2652  NAMD_EVENT_RANGE_2(patch->flags.event_on,
2653  NamdProfileEvent::SUBMIT_COLLECTIONS);
2654  int prec = Output::coordinateNeeded(step);
2655  if ( prec ) {
2656  collection->submitPositions(step,patch->atom,patch->lattice,prec);
2657  }
2658  if ( Output::velocityNeeded(step) ) {
2659  collection->submitVelocities(step,zeroVel,patch->atom);
2660  }
2661  if ( Output::forceNeeded(step) ) {
2662  int maxForceUsed = patch->flags.maxForceUsed;
2663  if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
2664  collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
2665  }
2666 }
HomePatch *const patch
Definition: Sequencer.h:133
static int velocityNeeded(int)
Definition: Output.C:365
static int coordinateNeeded(int)
Definition: Output.C:198
Lattice & lattice
Definition: Patch.h:126
void submitVelocities(int seq, int zero, FullAtomList &a)
static int forceNeeded(int)
Definition: Output.C:460
Flags flags
Definition: Patch.h:127
#define NAMD_EVENT_RANGE_2(eon, id)
int maxForceUsed
Definition: PatchTypes.h:31
void submitPositions(int seq, FullAtomList &a, Lattice l, int prec)
void submitForces(int seq, FullAtomList &a, int maxForceUsed, ForceList *f)
CollectionMgr *const collection
Definition: Sequencer.h:137
ForceList f[Results::maxNumForces]
Definition: Patch.h:207
void Sequencer::submitHalfstep ( int  step)
protected

Definition at line 2083 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ResizeArray< T >::begin(), Lattice::c(), doKineticEnergy, Flags::doVirial, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), Patch::lattice, GlobalMaster::lattice, Vector::length2(), FullAtom::mass, SimParameters::multigratorOn, NAMD_EVENT_RANGE_2, Patch::numAtoms, Lattice::origin(), Tensor::outerAdd(), SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, partition(), CompAtom::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, Vector::z, and z.

Referenced by integrate().

2084 {
2085  NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::SUBMIT_HALFSTEP);
2086 
2087  // velocity-dependent quantities *** ONLY ***
2088  // positions are not at half-step when called
2089  FullAtom *a = patch->atom.begin();
2090  int numAtoms = patch->numAtoms;
2091 
2092 #if CMK_BLUEGENEL
2093  CmiNetworkProgressAfter (0);
2094 #endif
2095 
2096  // For non-Multigrator doKineticEnergy = 1 always
2097  Tensor momentumSqrSum;
2099  {
2100  BigReal kineticEnergy = 0;
2101  Tensor virial;
2102  if ( simParams->pairInteractionOn ) {
2103  if ( simParams->pairInteractionSelf ) {
2104  for ( int i = 0; i < numAtoms; ++i ) {
2105  if ( a[i].partition != 1 ) continue;
2106  kineticEnergy += a[i].mass * a[i].velocity.length2();
2107  virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
2108  }
2109  }
2110  } else {
2111  for ( int i = 0; i < numAtoms; ++i ) {
2112  if (a[i].mass < 0.01) continue;
2113  kineticEnergy += a[i].mass * a[i].velocity.length2();
2114  virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
2115  }
2116  }
2117 
2119  momentumSqrSum = virial;
2120  }
2121  kineticEnergy *= 0.5 * 0.5;
2123  virial *= 0.5;
2124  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
2125 #ifdef ALTVIRIAL
2126  ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
2127 #endif
2128  }
2129 
2131  int nslabs = simParams->pressureProfileSlabs;
2132  const Lattice &lattice = patch->lattice;
2133  BigReal idz = nslabs/lattice.c().z;
2134  BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
2135  int useGroupPressure = simParams->useGroupPressure;
2136 
2137  // Compute kinetic energy partition, possibly subtracting off
2138  // internal kinetic energy if group pressure is enabled.
2139  // Since the regular pressure is 1/2 mvv and the internal kinetic
2140  // term that is subtracted off for the group pressure is
2141  // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
2142  // 1/2 m * v * v_cm. The factor of 1/2 is because submitHalfstep
2143  // gets called twice per timestep.
2144  int hgs;
2145  for (int i=0; i<numAtoms; i += hgs) {
2146  int j, ppoffset;
2147  hgs = a[i].hydrogenGroupSize;
2148  int partition = a[i].partition;
2149 
2150  BigReal m_cm = 0;
2151  Velocity v_cm(0,0,0);
2152  for (j=i; j< i+hgs; ++j) {
2153  m_cm += a[j].mass;
2154  v_cm += a[j].mass * a[j].velocity;
2155  }
2156  v_cm /= m_cm;
2157  for (j=i; j < i+hgs; ++j) {
2158  BigReal mass = a[j].mass;
2159  if (! (useGroupPressure && j != i)) {
2160  BigReal z = a[j].position.z;
2161  int slab = (int)floor((z-zmin)*idz);
2162  if (slab < 0) slab += nslabs;
2163  else if (slab >= nslabs) slab -= nslabs;
2164  ppoffset = 3*(slab + partition*nslabs);
2165  }
2166  BigReal wxx, wyy, wzz;
2167  if (useGroupPressure) {
2168  wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
2169  wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
2170  wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
2171  } else {
2172  wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
2173  wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
2174  wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
2175  }
2176  pressureProfileReduction->item(ppoffset ) += wxx;
2177  pressureProfileReduction->item(ppoffset+1) += wyy;
2178  pressureProfileReduction->item(ppoffset+2) += wzz;
2179  }
2180  }
2181  }
2182 
2183  // For non-Multigrator doKineticEnergy = 1 always
2185  {
2186  BigReal intKineticEnergy = 0;
2187  Tensor intVirialNormal;
2188 
2189  int hgs;
2190  for ( int i = 0; i < numAtoms; i += hgs ) {
2191 
2192 #if CMK_BLUEGENEL
2193  CmiNetworkProgress ();
2194 #endif
2195 
2196  hgs = a[i].hydrogenGroupSize;
2197  int j;
2198  BigReal m_cm = 0;
2199  Velocity v_cm(0,0,0);
2200  for ( j = i; j < (i+hgs); ++j ) {
2201  m_cm += a[j].mass;
2202  v_cm += a[j].mass * a[j].velocity;
2203  }
2205  momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
2206  }
2207  v_cm /= m_cm;
2208  if ( simParams->pairInteractionOn ) {
2209  if ( simParams->pairInteractionSelf ) {
2210  for ( j = i; j < (i+hgs); ++j ) {
2211  if ( a[j].partition != 1 ) continue;
2212  BigReal mass = a[j].mass;
2213  Vector v = a[j].velocity;
2214  Vector dv = v - v_cm;
2215  intKineticEnergy += mass * (v * dv);
2216  intVirialNormal.outerAdd (mass, v, dv);
2217  }
2218  }
2219  } else {
2220  for ( j = i; j < (i+hgs); ++j ) {
2221  BigReal mass = a[j].mass;
2222  Vector v = a[j].velocity;
2223  Vector dv = v - v_cm;
2224  intKineticEnergy += mass * (v * dv);
2225  intVirialNormal.outerAdd(mass, v, dv);
2226  }
2227  }
2228  }
2229 
2230  intKineticEnergy *= 0.5 * 0.5;
2232  intVirialNormal *= 0.5;
2233  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
2234  if ( simParams->multigratorOn) {
2235  momentumSqrSum *= 0.5;
2236  ADD_TENSOR_OBJECT(reduction,REDUCTION_MOMENTUM_SQUARED,momentumSqrSum);
2237  }
2238  }
2239 
2240 }
HomePatch *const patch
Definition: Sequencer.h:133
int doKineticEnergy
Definition: Sequencer.h:120
unsigned char partition
Definition: NamdTypes.h:56
SubmitReduction * pressureProfileReduction
Definition: Sequencer.h:135
SubmitReduction * reduction
Definition: Sequencer.h:134
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
Lattice & lattice
Definition: Patch.h:126
static void partition(int *order, const FullAtom *atoms, int begin, int end)
Definition: SortAtoms.C:42
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
BigReal & item(int i)
Definition: ReductionMgr.h:312
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
Definition: Tensor.h:255
Velocity velocity
Definition: NamdTypes.h:101
int pressureProfileSlabs
Vector origin() const
Definition: Lattice.h:262
Bool pairInteractionOn
Flags flags
Definition: Patch.h:127
gridSize z
int numAtoms
Definition: Patch.h:144
#define NAMD_EVENT_RANGE_2(eon, id)
BigReal x
Definition: Vector.h:66
BigReal length2(void) const
Definition: Vector.h:173
Definition: Tensor.h:15
int doVirial
Definition: PatchTypes.h:21
BigReal y
Definition: Vector.h:66
Mass mass
Definition: NamdTypes.h:108
Bool pairInteractionSelf
SimParameters *const simParams
Definition: Sequencer.h:132
Vector c() const
Definition: Lattice.h:254
double BigReal
Definition: common.h:114
iterator begin(void)
Definition: ResizeArray.h:36
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< T >::begin(), calcFixVirial(), SimParameters::drudeBondLen, SimParameters::drudeHardWallOn, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), 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().

2506 {
2507  FullAtom *a = patch->atom.begin();
2508  Force *f1 = patch->f[Results::normal].begin();
2509  Force *f2 = patch->f[Results::nbond].begin();
2510  Force *f3 = patch->f[Results::slow].begin();
2511  const bool fixedAtomsOn = simParams->fixedAtomsOn;
2512  const bool drudeHardWallOn = simParams->drudeHardWallOn;
2513  const double drudeBondLen = simParams->drudeBondLen;
2514  const double drudeBondLen2 = drudeBondLen * drudeBondLen;
2515  const double drudeStep = 0.1/(TIMEFACTOR*TIMEFACTOR);
2516  const double drudeMove = 0.01;
2517  const double drudeStep2 = drudeStep * drudeStep;
2518  const double drudeMove2 = drudeMove * drudeMove;
2519  int numAtoms = patch->numAtoms;
2520 
2521  reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
2522 
2523  for ( int i = 0; i < numAtoms; ++i ) {
2524  f1[i] += f2[i] + f3[i]; // add all forces
2525  if ( drudeHardWallOn && i && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) { // drude particle
2526  if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
2527  if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
2528  a[i].position += drudeMove * f1[i].unit();
2529  } else {
2530  a[i].position += drudeStep * f1[i];
2531  }
2532  if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
2533  a[i].position = a[i-1].position + drudeBondLen * (a[i].position - a[i-1].position).unit();
2534  }
2535  }
2536  Vector netf = f1[i-1] + f1[i];
2537  if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
2538  f1[i-1] = netf;
2539  f1[i] = 0.;
2540  }
2541  if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
2542  }
2543 
2544  f2 = f3 = 0; // included in f1
2545 
2546  BigReal maxv2 = 0.;
2547 
2548  for ( int i = 0; i < numAtoms; ++i ) {
2549  BigReal v2 = a[i].velocity.length2();
2550  if ( v2 > 0. ) {
2551  if ( v2 > maxv2 ) maxv2 = v2;
2552  } else {
2553  v2 = f1[i].length2();
2554  if ( v2 > maxv2 ) maxv2 = v2;
2555  }
2556  }
2557 
2558  if ( fmax2 > 10. * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR )
2559  { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial, true /* forces */); }
2560 
2561  BigReal fdotf = 0;
2562  BigReal fdotv = 0;
2563  BigReal vdotv = 0;
2564  int numHuge = 0;
2565  for ( int i = 0; i < numAtoms; ++i ) {
2566  if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
2567  if ( drudeHardWallOn && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) continue; // drude particle
2568  Force f = f1[i];
2569  BigReal ff = f * f;
2570  if ( ff > fmax2 ) {
2571  if (simParams->printBadContacts) {
2572  CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
2573  }
2574  ++numHuge;
2575  // pad scaling so minimizeMoveDownhill() doesn't miss them
2576  BigReal fmult = 1.01 * sqrt(fmax2/ff);
2577  f *= fmult; ff = f * f;
2578  f1[i] *= fmult;
2579  }
2580  fdotf += ff;
2581  fdotv += f * a[i].velocity;
2582  vdotv += a[i].velocity * a[i].velocity;
2583  }
2584 
2589 
2590  {
2591  Tensor intVirialNormal;
2592  Tensor intVirialNbond;
2593  Tensor intVirialSlow;
2594 
2595  int hgs;
2596  for ( int i = 0; i < numAtoms; i += hgs ) {
2597  hgs = a[i].hydrogenGroupSize;
2598  int j;
2599  BigReal m_cm = 0;
2600  Position x_cm(0,0,0);
2601  for ( j = i; j < (i+hgs); ++j ) {
2602  m_cm += a[j].mass;
2603  x_cm += a[j].mass * a[j].position;
2604  }
2605  x_cm /= m_cm;
2606  for ( j = i; j < (i+hgs); ++j ) {
2607  BigReal mass = a[j].mass;
2608  // net force treated as zero for fixed atoms
2609  if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
2610  Vector dx = a[j].position - x_cm;
2611  intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
2612  intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
2613  intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
2614  }
2615  }
2616 
2617  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
2618  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
2619  ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
2620  }
2621 
2622  if ( simParams->fixedAtomsOn ) {
2623  Tensor fixVirialNormal;
2624  Tensor fixVirialNbond;
2625  Tensor fixVirialSlow;
2626  Vector fixForceNormal = 0;
2627  Vector fixForceNbond = 0;
2628  Vector fixForceSlow = 0;
2629 
2630  calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
2631 
2632  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
2633  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
2634  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
2635  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
2636  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
2637  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
2638  }
2639 
2640  reduction->submit();
2641 }
HomePatch *const patch
Definition: Sequencer.h:133
void minimize_rattle2(const BigReal, Tensor *virial, bool forces=false)
Definition: HomePatch.C:2963
SubmitReduction * reduction
Definition: Sequencer.h:134
unsigned int hydrogenGroupSize
Definition: NamdTypes.h:58
unsigned int atomFixed
Definition: NamdTypes.h:96
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
BigReal & item(int i)
Definition: ReductionMgr.h:312
Position position
Definition: NamdTypes.h:53
void outerAdd(BigReal scale, const Vector &v1, const Vector &v2)
Definition: Tensor.h:255
Velocity velocity
Definition: NamdTypes.h:101
Flags flags
Definition: Patch.h:127
void calcFixVirial(Tensor &fixVirialNormal, Tensor &fixVirialNbond, Tensor &fixVirialSlow, Vector &fixForceNormal, Vector &fixForceNbond, Vector &fixForceSlow)
Definition: Sequencer.C:2242
BigReal drudeBondLen
int numAtoms
Definition: Patch.h:144
int sequence
Definition: PatchTypes.h:18
BigReal length2(void) const
Definition: Vector.h:173