#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_CUDA_SOA (int scriptTask)

void	initialize_integrate_CUDA_SOA (int scriptTask, int step, BigReal timestep, int numberOfSteps, int nbondstep, int slowstep, int maxForceUsed)

void	suspendULTs ()

void	wakeULTs ()

void	runComputeObjectsCUDA (int doMigration, int doGlobal, int pairlists, int nstep, int startup)

void	constructDevicePatchMap ()

void	updateDevicePatchMap (int startup)

void	printDevicePatchMap ()

void	clearDevicePatchMap ()

void	updateDeviceData (const int startup, const int maxForceUsed, const int doGlobal)

void	doMigrationGPU (const int startup, const int doGlobal, const int updatePatchMap)

void	monteCarloPressureControl (const int step, const int doMigration, const int doEnergy, const int doVirial, const int maxForceNumber, const int doGlobal)

void	integrate_SOA (int)

void	rattle1_SOA (BigReal, int)

void	addForceToMomentum_SOA (const double scaling, double dt_normal, double dt_nbond, double dt_slow, int maxForceNumber)

void	addVelocityToPosition_SOA (const double dt)

void	submitHalfstep_SOA ()

void	submitReductions_SOA ()

void	submitCollections_SOA (int step, int zeroVel=0)

void	maximumMove_SOA (const double dt, const double maxvel2)

void	langevinVelocitiesBBK1_SOA (BigReal timestep)

void	langevinVelocitiesBBK2_SOA (BigReal timestep)

void	berendsenPressure_SOA (int step)

void	langevinPiston_SOA (int step)

void	stochRescaleVelocities_SOA (int step)

void	runComputeObjects_SOA (int migration, int pairlists, int step)

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 (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
SubmitReduction *	min_reduction

int	pairlistsAreValid

int	pairlistsAge

int	pairlistsAgeLimit

BigReal	adaptTempT

int	rescaleVelocities_numTemps

int	stochRescale_count

int	berendsenPressure_count

int	checkpoint_berendsenPressure_count

int	slowFreq

SubmitReduction *	multigratorReduction

int	doKineticEnergy

int	doMomenta

Random *	random

SimParameters *const	simParams

HomePatch *const	patch

SubmitReduction *	reduction

SubmitReduction *	pressureProfileReduction

CollectionMgr *const	collection

ControllerBroadcasts *	broadcast

int	ldbSteps

bool	masterThread

Friends
class	HomePatch

class	SequencerCUDA

Detailed Description

Definition at line 43 of file Sequencer.h.

Constructor & Destructor Documentation

◆ Sequencer()

Sequencer::Sequencer ( HomePatch * p )

Definition at line 171 of file Sequencer.C.

References SimParameters::accelMDOn, berendsenPressure_count, broadcast, constructDevicePatchMap(), SimParameters::CUDASOAintegrateMode, deviceCUDA, DeviceCUDA::getDeviceID(), Patch::getPatchID(), HomePatch::ldObjHandle, masterThread, min_reduction, MULTIGRATOR_REDUCTION_MAX_RESERVED, SimParameters::multigratorOn, multigratorReduction, PatchMap::numPatches(), PatchMap::Object(), LdbCoordinator::Object(), ReductionMgr::Object(), pairlistsAgeLimit, SimParameters::pairlistsPerCycle, 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, SequencerCUDA, simParams, SimParameters::soluteScalingFactorCharge, SimParameters::soluteScalingOn, Random::split(), SimParameters::stepsPerCycle, stochRescale_count, and ReductionMgr::willSubmit().

                                  :
         simParams(Node::Object()->simParameters),
         patch(p),
         collection(CollectionMgr::Object()),
         ldbSteps(0),
         pairlistsAreValid(0),
         pairlistsAge(0),
         pairlistsAgeLimit(0)
 {
     pairlistsAgeLimit =
       (simParams->stepsPerCycle - 1) / simParams->pairlistsPerCycle;
     broadcast = new ControllerBroadcasts(& patch->ldObjHandle);
     reduction = ReductionMgr::Object()->willSubmit(
                   simParams->accelMDOn ? REDUCTIONS_AMD : REDUCTIONS_BASIC );
     min_reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MINIMIZER,1);
     if (simParams->pressureProfileOn) {
       int ntypes = simParams->pressureProfileAtomTypes;
       int nslabs = simParams->pressureProfileSlabs;
       pressureProfileReduction =
         ReductionMgr::Object()->willSubmit(
                 REDUCTIONS_PPROF_INTERNAL, 3*nslabs*ntypes);
     } else {
       pressureProfileReduction = NULL;
     }
     if (simParams->multigratorOn) {
       multigratorReduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_MULTIGRATOR,MULTIGRATOR_REDUCTION_MAX_RESERVED);
     } else {
       multigratorReduction = NULL;
     }
     ldbCoordinator = (LdbCoordinator::Object());
     random = new Random(simParams->randomSeed);
     random->split(patch->getPatchID()+1,PatchMap::Object()->numPatches()+1);
 
     // Is soluteScaling enabled?
     if (simParams->soluteScalingOn) {
       // If so, we must "manually" perform charge scaling on startup because
       // Sequencer will not get a scripting task for initial charge scaling.
       // Subsequent rescalings will take place through a scripting task.
       rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
     }
 
     rescaleVelocities_numTemps = 0;
     stochRescale_count = 0;
     berendsenPressure_count = 0;
     masterThread = true;
 //    patch->write_tip4_props();
 #if (defined(NAMD_CUDA) || defined(NAMD_HIP)) && defined(SEQUENCER_SOA) && defined(NODEGROUP_FORCE_REGISTER)
     if ( simParams->CUDASOAintegrateMode ){
 #if 0
       CUDASequencer = new SequencerCUDA(deviceCUDA->getDeviceID(),
                                         simParams);
 #else
       CUDASequencer = SequencerCUDA::InstanceInit(deviceCUDA->getDeviceID(),
                                         simParams);
        
       constructDevicePatchMap();
 #endif
       CUDASequencer->patchData->reduction->zero();
     }
 #endif
 }

◆ ~Sequencer()

Sequencer::~Sequencer ( void )

virtual

Definition at line 233 of file Sequencer.C.

References broadcast, clearDevicePatchMap(), SimParameters::CUDASOAintegrateMode, min_reduction, multigratorReduction, pressureProfileReduction, random, reduction, and simParams.

 {
     delete broadcast;
     delete reduction;
     delete min_reduction;
     if (pressureProfileReduction) delete pressureProfileReduction;
     delete random;
     if (multigratorReduction) delete multigratorReduction;
 #if (defined(NAMD_CUDA) || defined(NAMD_HIP)) && defined(SEQUENCER_SOA) && defined(NODEGROUP_FORCE_REGISTER)
     if ( simParams->CUDASOAintegrateMode ){
       delete CUDASequencer;
       clearDevicePatchMap();
     }
 #endif
 }

Member Function Documentation

◆ adaptTempUpdate()

void Sequencer::adaptTempUpdate ( int step )

protected

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

 {
    //check if adaptive tempering is enabled and in the right timestep range
    if (!simParams->adaptTempOn) return;
    if ( (step < simParams->adaptTempFirstStep ) ||
      ( simParams->adaptTempLastStep > 0 && step > simParams->adaptTempLastStep )) {
         if (simParams->langevinOn) // restore langevin temperature
             adaptTempT = simParams->langevinTemp;
         return;
    }
    // Get Updated Temperature
    if ( !(step % simParams->adaptTempFreq ) && (step > simParams->firstTimestep ))
      // Blocking receive for the updated adaptive tempering temperature.
      adaptTempT = broadcast->adaptTemperature.get(step);
 }

◆ addForceToMomentum()

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

protected

Definition at line 5641 of file Sequencer.C.

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

Referenced by newtonianVelocities().

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

◆ addForceToMomentum3()

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

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

Referenced by newtonianVelocities().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::ADD_FORCE_TO_MOMENTUM);
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter (0);
 #endif
   const BigReal dt1 = timestep1 / TIMEFACTOR;
   const BigReal dt2 = timestep2 / TIMEFACTOR;
   const BigReal dt3 = timestep3 / TIMEFACTOR;
   ForceList *f_use1 = (useSaved1 ? patch->f_saved : patch->f);
   ForceList *f_use2 = (useSaved2 ? patch->f_saved : patch->f);
   ForceList *f_use3 = (useSaved3 ? patch->f_saved : patch->f);
   FullAtom *atom_arr  = patch->atom.begin();
   const Force *force_arr1 = f_use1[ftag1].const_begin();
   const Force *force_arr2 = f_use2[ftag2].const_begin();
   const Force *force_arr3 = f_use3[ftag3].const_begin();
   patch->addForceToMomentum3 (atom_arr, force_arr1, force_arr2, force_arr3,
       dt1, dt2, dt3, patch->numAtoms);
 }

◆ addForceToMomentum_SOA()

void Sequencer::addForceToMomentum_SOA	(	const double	scaling,
		double	dt_normal,
		double	dt_nbond,
		double	dt_slow,
		int	maxForceNumber
	)

protected

Definition at line 2707 of file Sequencer.C.

References PatchDataSOA::f_nbond_x, PatchDataSOA::f_nbond_y, PatchDataSOA::f_nbond_z, PatchDataSOA::f_normal_x, PatchDataSOA::f_normal_y, PatchDataSOA::f_normal_z, PatchDataSOA::f_slow_x, PatchDataSOA::f_slow_y, PatchDataSOA::f_slow_z, Patch::flags, NAMD_EVENT_RANGE_2, Results::nbond, Results::normal, PatchDataSOA::numAtoms, patch, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, PatchDataSOA::recipMass, Results::slow, PatchDataSOA::vel_x, PatchDataSOA::vel_y, and PatchDataSOA::vel_z.

Referenced by integrate_SOA().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::ADD_FORCE_TO_MOMENTUM_SOA);
 
 #ifdef SOA_SIMPLIFY_PARAMS
   const double * __restrict recipMass = patch->patchDataSOA.recipMass;
   // force    Results::normal = 0
   const double * __restrict f_normal_x = patch->patchDataSOA.f_normal_x;
   const double * __restrict f_normal_y = patch->patchDataSOA.f_normal_y;
   const double * __restrict f_normal_z = patch->patchDataSOA.f_normal_z;
   // force    Results::nbond  = 1
   const double * __restrict f_nbond_x = patch->patchDataSOA.f_nbond_x;
   const double * __restrict f_nbond_y = patch->patchDataSOA.f_nbond_y;
   const double * __restrict f_nbond_z = patch->patchDataSOA.f_nbond_z;
   // force    Results::slow   = 2
   const double * __restrict f_slow_x = patch->patchDataSOA.f_slow_x;
   const double * __restrict f_slow_y = patch->patchDataSOA.f_slow_y;
   const double * __restrict f_slow_z = patch->patchDataSOA.f_slow_z;
   double       * __restrict vel_x = patch->patchDataSOA.vel_x;
   double       * __restrict vel_y = patch->patchDataSOA.vel_y;
   double       * __restrict vel_z = patch->patchDataSOA.vel_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
   //
   // We could combine each case into a single loop with breaks,
   // with all faster forces also summed, like addForceToMomentum3().
   //
   // Things to consider:
   // - Do we always use acceleration (f/m) instead of just plain force?
   //   Then we could instead buffer accel_slow, accel_nbond, etc.
   // - We will always need one multiply, since each dt includes
   //   also a scaling factor.
   //
   
 #if 0
     if(this->patch->getPatchID() == 538){
       // fprintf(stderr, "Old Positions %lf %lf %lf\n", patch->patchDataSOA.pos_x[43], patch->patchDataSOA.pos_y[43], patch->patchDataSOA.pos_z[43]);
       // fprintf(stderr, "Old Velocities %lf %lf %lf\n", vel_x[43], vel_y[43], vel_z[ 43]);
       // fprintf(stderr, "Adding forces %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", 
       //   f_slow_x[43],   f_slow_y[43],   f_slow_z[43], 
       //   f_nbond_x[43],  f_nbond_y[43],  f_nbond_z[43], 
       //   f_normal_x[43], f_normal_y[43], f_normal_z[43]);
       fprintf(stderr, "Old Positions %lf %lf %lf\n", patch->patchDataSOA.pos_x[0], patch->patchDataSOA.pos_y[0], patch->patchDataSOA.pos_z[0]);
       fprintf(stderr, "Old Velocities %lf %lf %lf\n", vel_x[0], vel_y[0], vel_z[ 0]);
       fprintf(stderr, "Adding forces %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", 
         f_slow_x[43],   f_slow_y[43],   f_slow_z[43], 
         f_nbond_x[43],  f_nbond_y[43],  f_nbond_z[43], 
         f_normal_x[43], f_normal_y[43], f_normal_z[43]);
     }
 #endif
   switch (maxForceNumber) {
     case Results::slow:
       dt_slow *= scaling;
       for (int i=0;  i < numAtoms;  i++) {
         vel_x[i] += f_slow_x[i] * recipMass[i] * dt_slow;
         vel_y[i] += f_slow_y[i] * recipMass[i] * dt_slow;
         vel_z[i] += f_slow_z[i] * recipMass[i] * dt_slow;
       }
       // fall through because we will always have the "faster" forces
     case Results::nbond:
       dt_nbond *= scaling;
       for (int i=0;  i < numAtoms;  i++) {
         vel_x[i] += f_nbond_x[i] * recipMass[i] * dt_nbond;
         vel_y[i] += f_nbond_y[i] * recipMass[i] * dt_nbond;
         vel_z[i] += f_nbond_z[i] * recipMass[i] * dt_nbond;
       }
       // fall through because we will always have the "faster" forces
     case Results::normal:
       dt_normal *= scaling;
       for (int i=0;  i < numAtoms;  i++) {
         vel_x[i] += f_normal_x[i] * recipMass[i] * dt_normal;
         vel_y[i] += f_normal_y[i] * recipMass[i] * dt_normal;
         vel_z[i] += f_normal_z[i] * recipMass[i] * dt_normal;
       }
   }
 }

◆ addMovDragToPosition()

void Sequencer::addMovDragToPosition ( BigReal timestep )

protected

Definition at line 4429 of file Sequencer.C.

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

                                                      {
   FullAtom *atom = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   Molecule *molecule = Node::Object()->molecule;   // need its methods
   const BigReal movDragGlobVel = simParams->movDragGlobVel;
   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
   Vector movDragVel, dragIncrement;
   for ( int i = 0; i < numAtoms; ++i )
   {
     // skip if fixed atom or zero drag attribute
     if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
          || !(molecule->is_atom_movdragged(atom[i].id)) ) continue;
     molecule->get_movdrag_params(movDragVel, atom[i].id);
     dragIncrement = movDragGlobVel * movDragVel * dt;
     atom[i].position += dragIncrement;
   }
 }

◆ addRotDragToPosition()

void Sequencer::addRotDragToPosition ( BigReal timestep )

protected

Definition at line 4448 of file Sequencer.C.

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

                                                      {
   FullAtom *atom = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   Molecule *molecule = Node::Object()->molecule;   // need its methods
   const BigReal rotDragGlobVel = simParams->rotDragGlobVel;
   const BigReal dt = timestep / TIMEFACTOR;   // MUST be as in the integrator!
   BigReal rotDragVel, dAngle;
   Vector atomRadius;
   Vector rotDragAxis, rotDragPivot, dragIncrement;
   for ( int i = 0; i < numAtoms; ++i )
   {
     // skip if fixed atom or zero drag attribute
     if ( (simParams->fixedAtomsOn && atom[i].atomFixed)
          || !(molecule->is_atom_rotdragged(atom[i].id)) ) continue;
     molecule->get_rotdrag_params(rotDragVel, rotDragAxis, rotDragPivot, atom[i].id);
     dAngle = rotDragGlobVel * rotDragVel * dt;
     rotDragAxis /= rotDragAxis.length();
     atomRadius = atom[i].position - rotDragPivot;
     dragIncrement = cross(rotDragAxis, atomRadius) * dAngle;
     atom[i].position += dragIncrement;
   }
 }

◆ addVelocityToPosition()

void Sequencer::addVelocityToPosition ( BigReal timestep )

protected

Definition at line 5680 of file Sequencer.C.

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

Referenced by integrate().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::ADD_VELOCITY_TO_POSITION);
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter (0);
 #endif
   const BigReal dt = timestep / TIMEFACTOR;
   FullAtom *atom_arr  = patch->atom.begin();
   patch->addVelocityToPosition(atom_arr, dt, patch->numAtoms);
 }

◆ addVelocityToPosition_SOA()

void Sequencer::addVelocityToPosition_SOA ( const double dt )

protected

Parameters

dt	scaled timestep

Definition at line 2810 of file Sequencer.C.

References Patch::flags, NAMD_EVENT_RANGE_2, PatchDataSOA::numAtoms, patch, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, PatchDataSOA::vel_x, PatchDataSOA::vel_y, and PatchDataSOA::vel_z.

Referenced by integrate_SOA().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::ADD_VELOCITY_TO_POSITION_SOA);
 #ifdef SOA_SIMPLIFY_PARAMS
   const double * __restrict vel_x = patch->patchDataSOA.vel_x;
   const double * __restrict vel_y = patch->patchDataSOA.vel_y;
   const double * __restrict vel_z = patch->patchDataSOA.vel_z;
   double *       __restrict pos_x = patch->patchDataSOA.pos_x;
   double *       __restrict pos_y = patch->patchDataSOA.pos_y;
   double *       __restrict pos_z = patch->patchDataSOA.pos_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
   for (int i=0;  i < numAtoms;  i++) {
     pos_x[i] += vel_x[i] * dt;
     pos_y[i] += vel_y[i] * dt;
     pos_z[i] += vel_z[i] * dt;
   }
 #if 0
     if(this->patch->getPatchID() == 538){
       fprintf(stderr, "New Positions %lf %lf %lf\n",  pos_x[43], pos_y[43], pos_z[43]);
       fprintf(stderr, "New Velocities %lf %lf %lf\n", vel_x[43], vel_y[43], vel_z[43]);
     }
 #endif
 
 }

◆ algorithm()

void Sequencer::algorithm ( void )

protectedvirtual

Definition at line 277 of file Sequencer.C.

References berendsenPressure_count, broadcast, HomePatch::checkpoint(), checkpoint_berendsenPressure_count, SimParameters::CUDASOAintegrate, END_OF_RUN, EVAL_MEASURE, HomePatch::exchangeAtoms(), HomePatch::exchangeCheckpoint(), FILE_OUTPUT, FORCE_OUTPUT, SimpleBroadcastObject< T >::get(), integrate(), integrate_CUDA_SOA(), integrate_SOA(), minimize(), NAMD_bug(), NAMD_die(), 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::SOAintegrateOn, SimParameters::soluteScalingFactorCharge, submitCollections(), and terminate().

 {
   int scriptTask;
   int scriptSeq = 0;
   // Blocking receive for the script barrier.
   while ( (scriptTask = broadcast->scriptBarrier.get(scriptSeq++)) != SCRIPT_END ) {
     switch ( scriptTask ) {
       case SCRIPT_OUTPUT:
         submitCollections(FILE_OUTPUT);
         break;
       case SCRIPT_FORCEOUTPUT:
         submitCollections(FORCE_OUTPUT);
         break;
       case SCRIPT_MEASURE:
         submitCollections(EVAL_MEASURE);
         break;
       case SCRIPT_REINITVELS:
         reinitVelocities();
         break;
       case SCRIPT_RESCALEVELS:
         rescaleVelocitiesByFactor(simParams->scriptArg1);
         break;
       case SCRIPT_RESCALESOLUTECHARGES:
         rescaleSoluteCharges(simParams->soluteScalingFactorCharge);
         break;
       case SCRIPT_RELOADCHARGES:
         reloadCharges();
         break;
       case SCRIPT_CHECKPOINT:
         patch->checkpoint();
         checkpoint_berendsenPressure_count = berendsenPressure_count;
         break;
       case SCRIPT_REVERT:
         patch->revert();
         berendsenPressure_count = checkpoint_berendsenPressure_count;
         pairlistsAreValid = 0;
         break;
       case SCRIPT_CHECKPOINT_STORE:
       case SCRIPT_CHECKPOINT_LOAD:
       case SCRIPT_CHECKPOINT_SWAP:
       case SCRIPT_CHECKPOINT_FREE:
         patch->exchangeCheckpoint(scriptTask,berendsenPressure_count);
         break;
       case SCRIPT_ATOMSENDRECV:
       case SCRIPT_ATOMSEND:
       case SCRIPT_ATOMRECV:
         patch->exchangeAtoms(scriptTask);
         break;
       case SCRIPT_MINIMIZE:
 #if 0
         if (simParams->CUDASOAintegrate){
           NAMD_die("Minimization is currently not supported on the GPU integrator\n");
         }
 #endif
               minimize();
               break;
       case SCRIPT_RUN:
       case SCRIPT_CONTINUE:
   //
   // DJH: Call a cleaned up version of integrate().
   //
   // We could test for simulation options and call a more basic version
   // of integrate() where we can avoid performing most tests.
   //
 #ifdef SEQUENCER_SOA
         if ( simParams->SOAintegrateOn ) {
 #ifdef NODEGROUP_FORCE_REGISTER
 
           if(simParams->CUDASOAintegrate) integrate_CUDA_SOA(scriptTask);
           else {
 #endif
             integrate_SOA(scriptTask);
 #ifdef NODEGROUP_FORCE_REGISTER
           }
 #endif
         }
         else
 #endif
         integrate(scriptTask);
         break;
       default:
         NAMD_bug("Unknown task in Sequencer::algorithm");
     }
   }
   submitCollections(END_OF_RUN);
   terminate();
 }

◆ awaken()

void Sequencer::awaken ( void )

inline

Definition at line 53 of file Sequencer.h.

References PRIORITY_SIZE.

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

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

◆ berendsenPressure()

void Sequencer::berendsenPressure ( int step )

protected

Definition at line 5276 of file Sequencer.C.

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

Referenced by integrate().

 {
   if ( simParams->berendsenPressureOn ) {
   berendsenPressure_count += 1;
   const int freq = simParams->berendsenPressureFreq;
   if ( ! (berendsenPressure_count % freq ) ) {
    berendsenPressure_count = 0;
    FullAtom *a = patch->atom.begin();
    int numAtoms = patch->numAtoms;
    // Blocking receive for the updated lattice scaling factor.
    Tensor factor = broadcast->positionRescaleFactor.get(step);
    patch->lattice.rescale(factor);
    if ( simParams->useGroupPressure )
    {
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
       int j;
       hgs = a[i].hydrogenGroupSize;
       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
         for ( j = i; j < (i+hgs); ++j ) {
           a[j].position = patch->lattice.apply_transform(
                                 a[j].fixedPosition,a[j].transform);
         }
         continue;
       }
       BigReal m_cm = 0;
       Position x_cm(0,0,0);
       for ( j = i; j < (i+hgs); ++j ) {
         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
         m_cm += a[j].mass;
         x_cm += a[j].mass * a[j].position;
       }
       x_cm /= m_cm;
       Position new_x_cm = x_cm;
       patch->lattice.rescale(new_x_cm,factor);
       Position delta_x_cm = new_x_cm - x_cm;
       for ( j = i; j < (i+hgs); ++j ) {
         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
           a[j].position = patch->lattice.apply_transform(
                                 a[j].fixedPosition,a[j].transform);
           continue;
         }
         a[j].position += delta_x_cm;
       }
     }
    }
    else
    {
     for ( int i = 0; i < numAtoms; ++i )
     {
       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
         a[i].position = patch->lattice.apply_transform(
                                 a[i].fixedPosition,a[i].transform);
         continue;
       }
       patch->lattice.rescale(a[i].position,factor);
     }
    }
   }
   } else {
     berendsenPressure_count = 0;
   }
 }

◆ berendsenPressure_SOA()

void Sequencer::berendsenPressure_SOA ( int step )

protected

Definition at line 3442 of file Sequencer.C.

References berendsenPressure_count, SimParameters::berendsenPressureFreq, broadcast, SimpleBroadcastObject< T >::get(), PatchDataSOA::hydrogenGroupSize, Patch::lattice, PatchDataSOA::mass, namd_reciprocal, PatchDataSOA::numAtoms, Lattice::origin(), patch, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, SimParameters::useGroupPressure, Vector::x, 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_SOA().

 {
 #ifdef SOA_SIMPLIFY_PARAMS
   const int    * __restrict hydrogenGroupSize = patch->patchDataSOA.hydrogenGroupSize;
   const float  * __restrict mass = patch->patchDataSOA.mass;
   double       * __restrict pos_x = patch->patchDataSOA.pos_x;
   double       * __restrict pos_y = patch->patchDataSOA.pos_y;
   double       * __restrict pos_z = patch->patchDataSOA.pos_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
 
   //
   // Loops below simplify if we lift out special cases of fixed atoms
   // and pressure excluded atoms and make them their own branch.
   //
 
   ++berendsenPressure_count;
   if ( berendsenPressure_count == simParams->berendsenPressureFreq ) {
     berendsenPressure_count = 0;
     // Blocking receive for the updated lattice scaling factor.
     Tensor factor = broadcast->positionRescaleFactor.get(step);
     patch->lattice.rescale(factor);
     Vector origin = patch->lattice.origin();
 
     if ( simParams->useGroupPressure ) {
       int hgs;
       for (int i = 0;  i < numAtoms;  i += hgs) {
         int j;
         hgs = hydrogenGroupSize[i];
         // missing fixed atoms implementation
         BigReal m_cm = 0;
         BigReal r_cm_x = 0;
         BigReal r_cm_y = 0;
         BigReal r_cm_z = 0;
         // calculate the center of mass
         for ( j = i; j < (i+hgs); ++j ) {
           m_cm += mass[j];
           r_cm_x += mass[j] * pos_x[j];
           r_cm_y += mass[j] * pos_y[j];
           r_cm_z += mass[j] * pos_z[j];
         }
         BigReal inv_m_cm = namd_reciprocal(m_cm);
         r_cm_x *= inv_m_cm;
         r_cm_y *= inv_m_cm;
         r_cm_z *= inv_m_cm;
         // scale the center of mass with factor
         // shift to origin
         double tx = r_cm_x - origin.x;
         double ty = r_cm_y - origin.y;
         double tz = r_cm_z - origin.z;
         // apply transformation 
         double new_r_cm_x = factor.xx*tx + factor.xy*ty + factor.xz*tz;
         double new_r_cm_y = factor.yx*tx + factor.yy*ty + factor.yz*tz;
         double new_r_cm_z = factor.zx*tx + factor.zy*ty + factor.zz*tz;
         // shift back
         new_r_cm_x += origin.x;
         new_r_cm_y += origin.y;
         new_r_cm_z += origin.z;
         // translation vector from old COM and new COM
         double delta_r_cm_x = new_r_cm_x - r_cm_x;
         double delta_r_cm_y = new_r_cm_y - r_cm_y;
         double delta_r_cm_z = new_r_cm_z - r_cm_z;
         // shift the hydrogen group with translation vector
         for (j = i;  j < (i+hgs);  ++j) {
           pos_x[j] += delta_r_cm_x;
           pos_y[j] += delta_r_cm_y;
           pos_z[j] += delta_r_cm_z;
         }
       }
     } else {
       for (int i = 0;  i < numAtoms;  ++i) {
         // missing fixed atoms implementation
         // scale the coordinates with factor
         // shift to origin
         double tx = pos_x[i] - origin.x;
         double ty = pos_y[i] - origin.y;
         double tz = pos_z[i] - origin.z;
         // apply transformation
         double ftx = factor.xx*tx + factor.xy*ty + factor.xz*tz;
         double fty = factor.yx*tx + factor.yy*ty + factor.yz*tz;
         double ftz = factor.zx*tx + factor.zy*ty + factor.zz*tz;
         // shift back
         pos_x[i] = ftx + origin.x;
         pos_y[i] = fty + origin.y;
         pos_z[i] = ftz + origin.z;
       }
     }
   }
 }

◆ calcFixVirial()

void Sequencer::calcFixVirial	(	Tensor &	fixVirialNormal,
		Tensor &	fixVirialNbond,
		Tensor &	fixVirialSlow,
		Vector &	fixForceNormal,
		Vector &	fixForceNbond,
		Vector &	fixForceSlow
	)

protected

Definition at line 5978 of file Sequencer.C.

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

                                                                        {
 
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
 
   for ( int j = 0; j < numAtoms; j++ ) {
     if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
       Vector dx = a[j].fixedPosition;
       // all negative because fixed atoms cancels these forces
       fixVirialNormal.outerAdd(-1.0, patch->f[Results::normal][j], dx);
       fixVirialNbond.outerAdd(-1.0, patch->f[Results::nbond][j], dx);
       fixVirialSlow.outerAdd(-1.0, patch->f[Results::slow][j], dx);
       fixForceNormal -= patch->f[Results::normal][j];
       fixForceNbond -= patch->f[Results::nbond][j];
       fixForceSlow -= patch->f[Results::slow][j];
     }
   }
 }

◆ calcKineticEnergy()

BigReal Sequencer::calcKineticEnergy ( )

protected

Definition at line 4969 of file Sequencer.C.

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

Referenced by multigratorTemperature().

                                      {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   BigReal kineticEnergy = 0.0;
   if ( simParams->pairInteractionOn ) {
     if ( simParams->pairInteractionSelf ) {
       for (int i = 0; i < numAtoms; ++i ) {
         if ( a[i].partition != 1 ) continue;
         kineticEnergy += a[i].mass * a[i].velocity.length2();
       }
     }
   } else {
     for (int i = 0; i < numAtoms; ++i ) {
       kineticEnergy += a[i].mass * a[i].velocity.length2();
     }
   }
   kineticEnergy *= 0.5;
   return kineticEnergy;
 }

◆ clearDevicePatchMap()

void Sequencer::clearDevicePatchMap ( )

protected

Referenced by ~Sequencer().

◆ constructDevicePatchMap()

void Sequencer::constructDevicePatchMap ( )

protected

Referenced by Sequencer().

◆ correctMomentum()

void Sequencer::correctMomentum	(	int	step,
		BigReal	drifttime
	)

protected

Definition at line 4758 of file Sequencer.C.

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

Referenced by integrate().

                                                            {
 
   //
   // DJH: This test should be done in SimParameters.
   //
   if ( simParams->fixedAtomsOn )
     NAMD_die("Cannot zero momentum when fixed atoms are present.");
 
   // Blocking receive for the momentum correction vector.
   const Vector dv = broadcast->momentumCorrection.get(step);
 
   const Vector dx = dv * ( drifttime / TIMEFACTOR );
 
   FullAtom *a = patch->atom.begin();
   const int numAtoms = patch->numAtoms;
 
 if ( simParams->zeroMomentumAlt ) {
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].velocity += dv * a[i].recipMass;
     a[i].position += dx * a[i].recipMass;
   }
 } else {
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].velocity += dv;
     a[i].position += dx;
   }
 }
 
 }

◆ cycleBarrier()

void Sequencer::cycleBarrier	(	int	doBarrier,
		int	step
	)

protected

Definition at line 6661 of file Sequencer.C.

References broadcast.

Referenced by integrate().

                                                     {
 #if USE_BARRIER
         if (doBarrier)
           // Blocking receive for the cycle barrier.
           broadcast->cycleBarrier.get(step);
 #endif
 }

◆ doMigrationGPU()

void Sequencer::doMigrationGPU	(	const int	startup,
		const int	doGlobal,
		const int	updatePatchMap
	)

protected

◆ hardWallDrude()

void Sequencer::hardWallDrude	(	BigReal	dt,
		int	pressure
	)

protected

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

 {
   if ( simParams->drudeHardWallOn ) {
     Tensor virial;
     Tensor *vp = ( pressure ? &virial : 0 );
     if ( patch->hardWallDrude(dt, vp, pressureProfileReduction) ) {
       iout << iERROR << "Constraint failure in HardWallDrude(); "
         << "simulation may become unstable.\n" << endi;
       Node::Object()->enableEarlyExit();
       terminate();
     }
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
   }
 }

◆ initialize_integrate_CUDA_SOA()

void Sequencer::initialize_integrate_CUDA_SOA	(	int	scriptTask,
		int	step,
		BigReal	timestep,
		int	numberOfSteps,
		int	nbondstep,
		int	slowstep,
		int	maxForceUsed
	)

protected

◆ integrate()

void Sequencer::integrate ( int scriptTask )

protected

Definition at line 3912 of file Sequencer.C.

References SimParameters::accelMDdihe, SimParameters::accelMDdual, SimParameters::accelMDOn, SimParameters::adaptTempOn, adaptTempT, adaptTempUpdate(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), SimParameters::alchOn, SimParameters::alchOutFreq, Results::amdf, Patch::atomMapper, ResizeArray< Elem >::begin(), berendsenPressure(), SimParameters::colvarsOn, SimParameters::commOnly, SimParameters::computeEnergies, ComputeMgr::computeGlobalObject, Node::computeMgr, correctMomentum(), cycleBarrier(), D_MSG, DebugM, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, doKineticEnergy, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, doMomenta, Flags::doNonbonded, Flags::doVirial, SimParameters::dt, ResizeArray< Elem >::end(), endi(), eventEndOfTimeStep, 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_gcd(), NAMD_PROFILE_START, NAMD_PROFILE_STOP, NamdProfileEventStr, Results::nbond, newtonianVelocities(), SimParameters::nonbondedFrequency, Results::normal, SimParameters::numTraceSteps, Node::Object(), 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, SimParameters::stepsPerCycle, SimParameters::stochRescaleOn, 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().

                                         {
     char traceNote[24];
     char tracePrefix[20];
     sprintf(tracePrefix, "p:%d,s:",patch->patchID);
 //    patch->write_tip4_props();
 
     //
     // DJH: Copy all data into SOA (structure of arrays)
     // from AOS (array of structures) data structure.
     //
     //patch->copy_all_to_SOA();
 
 #ifdef TIMER_COLLECTION
     TimerSet& t = patch->timerSet;
 #endif
     TIMER_INIT_WIDTH(t, KICK, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, MAXMOVE, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, DRIFT, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, PISTON, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, SUBMITHALF, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, VELBBK1, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, VELBBK2, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, RATTLE1, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, SUBMITFULL, simParams->timerBinWidth);
     TIMER_INIT_WIDTH(t, SUBMITCOLLECT, simParams->timerBinWidth);
 
     int &step = patch->flags.step;
     step = simParams->firstTimestep;
 
     // drag switches
     const Bool rotDragOn = simParams->rotDragOn;
     const Bool movDragOn = simParams->movDragOn;
 
     const int commOnly = simParams->commOnly;
 
     int &maxForceUsed = patch->flags.maxForceUsed;
     int &maxForceMerged = patch->flags.maxForceMerged;
     maxForceUsed = Results::normal;
     maxForceMerged = Results::normal;
 
     const int numberOfSteps = simParams->N;
     const int stepsPerCycle = simParams->stepsPerCycle;
     const BigReal timestep = simParams->dt;
 
     // what MTS method?
     const int staleForces = ( simParams->MTSAlgorithm == NAIVE );
 
     const int nonbondedFrequency = simParams->nonbondedFrequency;
     slowFreq = nonbondedFrequency;
     const BigReal nbondstep = timestep * (staleForces?1:nonbondedFrequency);
     int &doNonbonded = patch->flags.doNonbonded;
     doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
     if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
     if ( doNonbonded ) maxForceUsed = Results::nbond;
 
     // Do we do full electrostatics?
     const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
     const int fullElectFrequency = simParams->fullElectFrequency;
     if ( dofull ) slowFreq = fullElectFrequency;
     const BigReal slowstep = timestep * (staleForces?1:fullElectFrequency);
     int &doFullElectrostatics = patch->flags.doFullElectrostatics;
     doFullElectrostatics = (dofull && ((step >= numberOfSteps) || !(step%fullElectFrequency)));
     if ( dofull && (fullElectFrequency == 1) && !(simParams->mollyOn) )
                                         maxForceMerged = Results::slow;
     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
 
 //#ifndef UPPER_BOUND
     const Bool accelMDOn = simParams->accelMDOn;
     const Bool accelMDdihe = simParams->accelMDdihe;
     const Bool accelMDdual = simParams->accelMDdual;
     if ( accelMDOn && (accelMDdihe || accelMDdual)) maxForceUsed = Results::amdf;
 
     // Is adaptive tempering on?
     const Bool adaptTempOn = simParams->adaptTempOn;
     adaptTempT = simParams->initialTemp;
     if (simParams->langevinOn)
         adaptTempT = simParams->langevinTemp;
     else if (simParams->rescaleFreq > 0)
         adaptTempT = simParams->rescaleTemp;
 
 
     int &doMolly = patch->flags.doMolly;
     doMolly = simParams->mollyOn && doFullElectrostatics;
     // BEGIN LA
     int &doLoweAndersen = patch->flags.doLoweAndersen;
     doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
     // END LA
 
     int &doGBIS = patch->flags.doGBIS;
     doGBIS = simParams->GBISOn;
 
     int &doLCPO = patch->flags.doLCPO;
     doLCPO = simParams->LCPOOn;
 
     int zeroMomentum = simParams->zeroMomentum;
 
     // Do we need to return forces to TCL script or Colvar module?
     int doTcl = simParams->tclForcesOn;
     int doColvars = simParams->colvarsOn;
 //#endif
     int doGlobal = doTcl || doColvars;
     ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
 
     // Bother to calculate energies?
     int &doEnergy = patch->flags.doEnergy;
     int energyFrequency = simParams->computeEnergies;
 #if defined(NAMD_CUDA) || defined(NAMD_HIP)
     if(simParams->alchOn) energyFrequency = NAMD_gcd(energyFrequency, simParams->alchOutFreq);
 #endif
 #ifndef UPPER_BOUND
     const int reassignFreq = simParams->reassignFreq;
 #endif
 
     int &doVirial = patch->flags.doVirial;
     doVirial = 1;
 
   if ( scriptTask == SCRIPT_RUN ) {
 
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 
 #ifndef UPPER_BOUND
 //    printf("Doing initial rattle\n");
 #ifndef UPPER_BOUND
 D_MSG("rattle1()");
     TIMER_START(t, RATTLE1);
     rattle1(0.,0);  // enforce rigid bond constraints on initial positions
     TIMER_STOP(t, RATTLE1);
 #endif
 
     if (simParams->lonepairs || simParams->singleTopology) {
       patch->atomMapper->registerIDsFullAtom(
                 patch->atom.begin(),patch->atom.end());
     }
 
     if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
        reassignVelocities(timestep,step);
     }
 #endif
 
     doEnergy = ! ( step % energyFrequency );
 #ifndef UPPER_BOUND
     if ( accelMDOn && !accelMDdihe ) doEnergy=1;
     //Update energy every timestep for adaptive tempering
     if ( adaptTempOn ) doEnergy=1;
 #endif
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 D_MSG("runComputeObjects()");
     runComputeObjects(1,step<numberOfSteps); // must migrate here!
 #ifndef UPPER_BOUND
     rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
     adaptTempUpdate(step); // update adaptive tempering temperature
 #endif
 
 #ifndef UPPER_BOUND
     if ( staleForces || doGlobal ) {
       if ( doNonbonded ) saveForce(Results::nbond);
       if ( doFullElectrostatics ) saveForce(Results::slow);
     }
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
     if ( ! commOnly ) {
 D_MSG("newtonianVelocities()");
       TIMER_START(t, KICK);
       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
       TIMER_STOP(t, KICK);
     }
     minimizationQuenchVelocity();
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 #ifndef UPPER_BOUND
 D_MSG("rattle1()");
     TIMER_START(t, RATTLE1);
     rattle1(-timestep,0);
     TIMER_STOP(t, RATTLE1);
 #endif
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 D_MSG("submitHalfstep()");
     TIMER_START(t, SUBMITHALF);
     submitHalfstep(step);
     TIMER_STOP(t, SUBMITHALF);
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
     if ( ! commOnly ) {
 D_MSG("newtonianVelocities()");
       TIMER_START(t, KICK);
       newtonianVelocities(1.0,timestep,nbondstep,slowstep,0,1,1);
       TIMER_STOP(t, KICK);
     }
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 D_MSG("rattle1()");
     TIMER_START(t, RATTLE1);
     rattle1(timestep,1);
     TIMER_STOP(t, RATTLE1);
     if (doGlobal)  // include constraint forces
       computeGlobal->saveTotalForces(patch);
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 D_MSG("submitHalfstep()");
     TIMER_START(t, SUBMITHALF);
     submitHalfstep(step);
     TIMER_STOP(t, SUBMITHALF);
     if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
     if ( ! commOnly ) {
 D_MSG("newtonianVelocities()");
       TIMER_START(t, KICK);
       newtonianVelocities(-0.5,timestep,nbondstep,slowstep,0,1,1);
       TIMER_STOP(t, KICK);
     }
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 #endif
 D_MSG("submitReductions()");
     TIMER_START(t, SUBMITFULL);
     submitReductions(step);
     TIMER_STOP(t, SUBMITFULL);
 //    print_vel_AOS(patch->atom.begin(), 0, patch->numAtoms);
 #ifndef UPPER_BOUND
     if(0){ // if(traceIsOn()){
         traceUserEvent(eventEndOfTimeStep);
         sprintf(traceNote, "%s%d",tracePrefix,step);
         traceUserSuppliedNote(traceNote);
     }
 #endif
     rebalanceLoad(step);
 
   } // scriptTask == SCRIPT_RUN
 
 #ifndef UPPER_BOUND
   bool doMultigratorRattle = false;
 #endif
 
   //
   // DJH: There are a lot of mod operations below and elsewhere to
   // test step number against the frequency of something happening.
   // Mod and integer division are expensive!
   // Might be better to replace with counters and test equality.
   //
 #if 0
     for(int i = 0; i < NamdProfileEvent::EventsCount; i++)
         CkPrintf("-------------- [%d] %s -------------\n", i, NamdProfileEventStr[i]);
 #endif
 
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
   int& eon = patch->flags.event_on;
   int epid = (simParams->beginEventPatchID <= patch->getPatchID()
       && patch->getPatchID() <= simParams->endEventPatchID);
   int beginStep = simParams->beginEventStep;
   int endStep = simParams->endEventStep;
   bool controlProfiling = patch->getPatchID() == 0;
 #endif
 
     for ( ++step; step <= numberOfSteps; ++step )
     {
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED)  || defined(NAMD_ROCTX_ENABLED)
       eon = epid && (beginStep < step && step <= endStep);
 
       if (controlProfiling && step == beginStep) {
         NAMD_PROFILE_START();
       }
       if (controlProfiling && step == endStep) {
         NAMD_PROFILE_STOP();
       }
       char buf[32];
       sprintf(buf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::INTEGRATE_1], patch->getPatchID());
       NAMD_EVENT_START_EX(eon, NamdProfileEvent::INTEGRATE_1, buf);
 #endif
       DebugM(3,"for step "<<step<< " dGlobal " << doGlobal<<"\n"<<endi);
 #ifndef UPPER_BOUND
       rescaleVelocities(step);
       tcoupleVelocities(timestep,step);
       if ( simParams->stochRescaleOn ) {
         stochRescaleVelocities(step);
       }
       berendsenPressure(step);
 
       if ( ! commOnly ) {
         TIMER_START(t, KICK);
         newtonianVelocities(0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
         TIMER_STOP(t, KICK);
       }
 
       // We do RATTLE here if multigrator thermostat was applied in the previous step
       if (doMultigratorRattle) rattle1(timestep, 1);
 
       /* reassignment based on half-step velocities
          if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
          addVelocityToPosition(0.5*timestep);
          reassignVelocities(timestep,step);
          addVelocityToPosition(0.5*timestep);
          rattle1(0.,0);
          rattle1(-timestep,0);
          addVelocityToPosition(-1.0*timestep);
          rattle1(timestep,0);
          } */
 
       TIMER_START(t, MAXMOVE);
       maximumMove(timestep);
       TIMER_STOP(t, MAXMOVE);
 
       NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_1);  // integrate 1
 
       if ( simParams->langevinPistonOn || (simParams->langevinOn && simParams->langevin_useBAOAB) ) {
         if ( ! commOnly ) {
           TIMER_START(t, DRIFT);
           addVelocityToPosition(0.5*timestep);
           TIMER_STOP(t, DRIFT);
         }
         // We add an Ornstein-Uhlenbeck integration step for the case of BAOAB (Langevin)
         langevinVelocities(timestep);
 
         // There is a blocking receive inside of langevinPiston()
         // that might suspend the current thread of execution,
         // so split profiling around this conditional block.
         langevinPiston(step);
 
         if ( ! commOnly ) {
           TIMER_START(t, DRIFT);
           addVelocityToPosition(0.5*timestep);
           TIMER_STOP(t, DRIFT);
         }
       } else {
         // If Langevin is not used, take full time step directly instread of two half steps
         if ( ! commOnly ) {
           TIMER_START(t, DRIFT);
           addVelocityToPosition(timestep);
           TIMER_STOP(t, DRIFT);
         }
       }
 
       NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_2);
 
       // impose hard wall potential for Drude bond length
       hardWallDrude(timestep, 1);
 
       minimizationQuenchVelocity();
 #endif // UPPER_BOUND
 
       doNonbonded = !(step%nonbondedFrequency);
       doFullElectrostatics = (dofull && !(step%fullElectFrequency));
 
 #ifndef UPPER_BOUND
       if ( zeroMomentum && doFullElectrostatics ) {
         // There is a blocking receive inside of correctMomentum().
         correctMomentum(step,slowstep);
       }
 
       // There are NO sends in submitHalfstep() just local summation
       // into the Reduction struct.
       TIMER_START(t, SUBMITHALF);
       submitHalfstep(step);
       TIMER_STOP(t, SUBMITHALF);
 
       doMolly = simParams->mollyOn && doFullElectrostatics;
       // BEGIN LA
       doLoweAndersen = simParams->loweAndersenOn && doNonbonded;
       // END LA
 
       maxForceUsed = Results::normal;
       if ( doNonbonded ) maxForceUsed = Results::nbond;
       if ( doFullElectrostatics ) maxForceUsed = Results::slow;
       if ( accelMDOn && (accelMDdihe || accelMDdual))  maxForceUsed = Results::amdf;
 
       // Migrate Atoms on stepsPerCycle
       doEnergy = ! ( step % energyFrequency );
       if ( accelMDOn && !accelMDdihe ) doEnergy=1;
       if ( adaptTempOn ) doEnergy=1;
 
       // Multigrator
       if (simParams->multigratorOn) {
         doVirial = (!(step % energyFrequency) || ((simParams->outputPressure > 0) && !(step % simParams->outputPressure))
           || !(step % simParams->multigratorPressureFreq));
         doKineticEnergy = (!(step % energyFrequency) || !(step % simParams->multigratorTemperatureFreq));
         doMomenta = (simParams->outputMomenta > 0) && !(step % simParams->outputMomenta);
       } else {
         doVirial = 1;
         doKineticEnergy = 1;
         doMomenta = 1;
       }
 #endif
       NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_2);  // integrate 2
 
       // The current thread of execution will suspend in runComputeObjects().
       runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
 
       NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_3);
 
 #ifndef UPPER_BOUND
       rescaleaccelMD(step, doNonbonded, doFullElectrostatics); // for accelMD
 
       if ( staleForces || doGlobal ) {
         if ( doNonbonded ) saveForce(Results::nbond);
         if ( doFullElectrostatics ) saveForce(Results::slow);
       }
 
       // reassignment based on full-step velocities
       if ( !commOnly && ( reassignFreq>0 ) && ! (step%reassignFreq) ) {
         reassignVelocities(timestep,step);
         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
         rattle1(-timestep,0);
       }
 
       if ( ! commOnly ) {
         TIMER_START(t, VELBBK1);
         langevinVelocitiesBBK1(timestep);
         TIMER_STOP(t, VELBBK1);
         TIMER_START(t, KICK);
         newtonianVelocities(1.0,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
         TIMER_STOP(t, KICK);
         TIMER_START(t, VELBBK2);
         langevinVelocitiesBBK2(timestep);
         TIMER_STOP(t, VELBBK2);
       }
 
       // add drag to each atom's positions
       if ( ! commOnly && movDragOn ) addMovDragToPosition(timestep);
       if ( ! commOnly && rotDragOn ) addRotDragToPosition(timestep);
 
       TIMER_START(t, RATTLE1);
       rattle1(timestep,1);
       TIMER_STOP(t, RATTLE1);
       if (doGlobal)  // include constraint forces
         computeGlobal->saveTotalForces(patch);
 
       TIMER_START(t, SUBMITHALF);
       submitHalfstep(step);
       TIMER_STOP(t, SUBMITHALF);
       if ( zeroMomentum && doFullElectrostatics ) submitMomentum(step);
 
       if ( ! commOnly ) {
         TIMER_START(t, KICK);
         newtonianVelocities(-0.5,timestep,nbondstep,slowstep,staleForces,doNonbonded,doFullElectrostatics);
         TIMER_STOP(t, KICK);
       }
 
         // rattle2(timestep,step);
 #endif
 
         TIMER_START(t, SUBMITFULL);
         submitReductions(step);
         TIMER_STOP(t, SUBMITFULL);
         TIMER_START(t, SUBMITCOLLECT);
         submitCollections(step);
         TIMER_STOP(t, SUBMITCOLLECT);
 #ifndef UPPER_BOUND
        //Update adaptive tempering temperature
         adaptTempUpdate(step);
 
       // Multigrator temperature and pressure steps
       multigratorTemperature(step, 1);
       multigratorPressure(step, 1);
       multigratorPressure(step, 2);
       multigratorTemperature(step, 2);
       doMultigratorRattle = (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq));
 
       NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_3); // integrate 3
 #endif
 
 #if CYCLE_BARRIER
         cycleBarrier(!((step+1) % stepsPerCycle), step);
 #elif PME_BARRIER
         cycleBarrier(doFullElectrostatics, step);
 #elif  STEP_BARRIER
         cycleBarrier(1, step);
 #endif
 
 #ifndef UPPER_BOUND
          if(Node::Object()->specialTracing || simParams->statsOn){
                  int bstep = simParams->traceStartStep;
                  int estep = bstep + simParams->numTraceSteps;
                  if(step == bstep || step == estep){
                          traceBarrier(step);
                  }
          }
 
 #ifdef MEASURE_NAMD_WITH_PAPI
          if(simParams->papiMeasure) {
                  int bstep = simParams->papiMeasureStartStep;
                  int estep = bstep + simParams->numPapiMeasureSteps;
                  if(step == bstep || step==estep) {
                          papiMeasureBarrier(step);
                  }
          }
 #endif
 
         if(0){ // if(traceIsOn()){
             traceUserEvent(eventEndOfTimeStep);
             sprintf(traceNote, "%s%d",tracePrefix,step);
             traceUserSuppliedNote(traceNote);
         }
 #endif // UPPER_BOUND
         rebalanceLoad(step);
 
 #if PME_BARRIER
         // a step before PME
         cycleBarrier(dofull && !((step+1)%fullElectFrequency),step);
 #endif
 
 #if USE_HPM
         if(step == START_HPM_STEP)
           (CProxy_Node(CkpvAccess(BOCclass_group).node)).startHPM();
 
         if(step == STOP_HPM_STEP)
           (CProxy_Node(CkpvAccess(BOCclass_group).node)).stopHPM();
 #endif
 
     }
 
   TIMER_DONE(t);
 #ifdef TIMER_COLLECTION
   if (patch->patchID == SPECIAL_PATCH_ID) {
     printf("Timer collection reporting in microseconds for "
         "Patch %d\n", patch->patchID);
     TIMER_REPORT(t);
   }
 #endif // TIMER_COLLECTION
     //
     // DJH: Copy updates of SOA back into AOS.
     //
     //patch->copy_updates_to_AOS();
 }

◆ integrate_CUDA_SOA()

void Sequencer::integrate_CUDA_SOA ( int scriptTask )

protected

Referenced by algorithm().

◆ integrate_SOA()

void Sequencer::integrate_SOA ( int scriptTask )

protected

Definition at line 2032 of file Sequencer.C.

References addForceToMomentum_SOA(), addVelocityToPosition_SOA(), SimParameters::alchOn, SimParameters::alchOutFreq, berendsenPressure_SOA(), SimParameters::berendsenPressureOn, CheckStep::check(), SimParameters::colvarsOn, SimParameters::computeEnergies, ComputeMgr::computeGlobalObject, Node::computeMgr, SimParameters::cutoff, SimParameters::dcdFrequency, Molecule::dcdSelectionParams, DebugM, Flags::doEnergy, Flags::doFullElectrostatics, Flags::doGBIS, doKineticEnergy, Flags::doLCPO, Flags::doLoweAndersen, Flags::doMolly, doMomenta, Flags::doNonbonded, Flags::doVirial, SimParameters::dt, endi(), PatchDataSOA::f_nbond_x, PatchDataSOA::f_nbond_y, PatchDataSOA::f_nbond_z, PatchDataSOA::f_normal_x, PatchDataSOA::f_normal_y, PatchDataSOA::f_normal_z, PatchDataSOA::f_slow_x, PatchDataSOA::f_slow_y, PatchDataSOA::f_slow_z, SimParameters::firstTimestep, Patch::flags, SimParameters::forceDcdFrequency, dcd_params::frequency, SimParameters::fullElectFrequency, PatchDataSOA::gaussrand_x, PatchDataSOA::gaussrand_y, PatchDataSOA::gaussrand_z, Patch::getPatchID(), PatchDataSOA::hydrogenGroupSize, PatchDataSOA::id, SimParameters::IMDfreq, CheckStep::init(), PatchDataSOA::langevinParam, langevinPiston_SOA(), SimParameters::langevinPistonOn, langevinVelocitiesBBK1_SOA(), langevinVelocitiesBBK2_SOA(), PatchDataSOA::langScalRandBBK2, PatchDataSOA::langScalVelBBK2, PatchDataSOA::mass, Flags::maxForceMerged, Flags::maxForceUsed, maximumMove_SOA(), Node::molecule, SimParameters::N, NAMD_EVENT_START, NAMD_EVENT_START_EX, NAMD_EVENT_STOP, NAMD_gcd(), NAMD_PROFILE_START, NAMD_PROFILE_STOP, NamdProfileEventStr, Results::nbond, SimParameters::nonbondedFrequency, Results::normal, PatchDataSOA::numAtoms, Node::Object(), patch, Patch::patchID, CheckStep::period, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, rattle1_SOA(), rebalanceLoad(), RECIP_TIMEFACTOR, PatchDataSOA::recipMass, SimParameters::restartFrequency, runComputeObjects_SOA(), ComputeGlobal::saveTotalForces(), SCRIPT_RUN, simParams, Results::slow, slowFreq, SOA_SIMPLIFY_PARAMS, SPECIAL_PATCH_ID, Flags::step, SimParameters::stepsPerCycle, SimParameters::stochRescaleOn, stochRescaleVelocities_SOA(), submitCollections_SOA(), submitHalfstep_SOA(), submitReductions_SOA(), SimParameters::tclForcesOn, TIMER_DONE, TIMER_INIT_WIDTH, TIMER_REPORT, TIMER_START, TIMER_STOP, PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, SimParameters::velDcdFrequency, Vector::x, Vector::y, and Vector::z.

Referenced by algorithm().

                                             {
   //
   // Below when accessing the array buffers for position, velocity, force,
   // note that we don't want to set up pointers directly to the buffers
   // because the allocations might get resized after atom migration.
   //
 
 #ifdef TIMER_COLLECTION
   TimerSet& t = patch->timerSet;
 #endif
   TIMER_INIT_WIDTH(t, KICK, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, MAXMOVE, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, DRIFT, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, PISTON, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, SUBMITHALF, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, VELBBK1, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, VELBBK2, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, RATTLE1, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, SUBMITFULL, simParams->timerBinWidth);
   TIMER_INIT_WIDTH(t, SUBMITCOLLECT, simParams->timerBinWidth);
 
   // Keep track of the step number.
   int &step = patch->flags.step;
   step = simParams->firstTimestep;
 
   // For multiple time stepping, which force boxes are used?
   int &maxForceUsed = patch->flags.maxForceUsed;
   int &maxForceMerged = patch->flags.maxForceMerged;
   maxForceUsed = Results::normal;
   maxForceMerged = Results::normal;
 
   // Keep track of total steps and steps per cycle.
   const int numberOfSteps = simParams->N;
   //const int stepsPerCycle = simParams->stepsPerCycle;
   CheckStep stepsPerCycle;
   stepsPerCycle.init(step, simParams->stepsPerCycle);
   // The fundamental time step, get the scaling right for velocity units.
   const BigReal timestep = simParams->dt * RECIP_TIMEFACTOR;
 
   //const int nonbondedFrequency = simParams->nonbondedFrequency;
   //slowFreq = nonbondedFrequency;
   CheckStep nonbondedFrequency;
   slowFreq = simParams->nonbondedFrequency;
   // The step size for short-range nonbonded forces.
   const BigReal nbondstep = timestep * simParams->nonbondedFrequency;
   int &doNonbonded = patch->flags.doNonbonded;
   //doNonbonded = (step >= numberOfSteps) || !(step%nonbondedFrequency);
   doNonbonded = (step >= numberOfSteps) ||
     nonbondedFrequency.init(step, simParams->nonbondedFrequency);
   //if ( nonbondedFrequency == 1 ) maxForceMerged = Results::nbond;
   if ( nonbondedFrequency.period == 1 ) maxForceMerged = Results::nbond;
   if ( doNonbonded ) maxForceUsed = Results::nbond;
 
   // Do we do full electrostatics?
   const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
   //const int fullElectFrequency = simParams->fullElectFrequency;
   //if ( dofull ) slowFreq = fullElectFrequency;
   CheckStep fullElectFrequency;
   if ( dofull ) slowFreq = simParams->fullElectFrequency;
   // The step size for long-range electrostatics.
   const BigReal slowstep = timestep * simParams->fullElectFrequency;
   int &doFullElectrostatics = patch->flags.doFullElectrostatics;
   //doFullElectrostatics = (dofull &&
   //    ((step >= numberOfSteps) || !(step%fullElectFrequency)));
   doFullElectrostatics = (dofull &&
       ((step >= numberOfSteps) ||
        fullElectFrequency.init(step, simParams->fullElectFrequency)));
   //if ( dofull && fullElectFrequency == 1 ) maxForceMerged = Results::slow;
   if ( dofull && fullElectFrequency.period == 1 ) maxForceMerged = Results::slow;
   if ( doFullElectrostatics ) maxForceUsed = Results::slow;
 
   // Bother to calculate energies?
   int &doEnergy = patch->flags.doEnergy;
   //int energyFrequency = simParams->outputEnergies;
   CheckStep energyFrequency;
   int newComputeEnergies = simParams->computeEnergies;
   if(simParams->alchOn) newComputeEnergies = NAMD_gcd(newComputeEnergies, simParams->alchOutFreq);
   doEnergy = energyFrequency.init(step, newComputeEnergies);
 
   // Do we need to return forces to TCL script or Colvar module?
   int doTcl = simParams->tclForcesOn;
   int doColvars = simParams->colvarsOn;
   int doGlobal = doTcl || doColvars;
   ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
   int &doVirial = patch->flags.doVirial;
   doVirial = 1;
 
   // The following flags have to be explicitly disabled in Patch object.
   patch->flags.doMolly = 0;
   patch->flags.doLoweAndersen = 0;
   patch->flags.doGBIS = 0;
   patch->flags.doLCPO = 0;
 
   // Square of maximum velocity for simulation safety check
   const BigReal maxvel2 =
     (simParams->cutoff * simParams->cutoff) / (timestep * timestep);
 
   // check for Langevin piston
   // set period beyond numberOfSteps to disable
   CheckStep langevinPistonFrequency;
   langevinPistonFrequency.init(step,
       (simParams->langevinPistonOn ? slowFreq : numberOfSteps+1 ),
       (simParams->langevinPistonOn ? -1-slowFreq/2 : 0) /* = delta */);
 
   // check for output
   // set period beyond numberOfSteps to disable
   CheckStep restartFrequency;
   restartFrequency.init(step, (simParams->restartFrequency ?
         simParams->restartFrequency : numberOfSteps+1) );
   CheckStep dcdFrequency;
   dcdFrequency.init(step, (simParams->dcdFrequency ?
         simParams->dcdFrequency : numberOfSteps+1) );
   CheckStep velDcdFrequency;
   velDcdFrequency.init(step, (simParams->velDcdFrequency ?
         simParams->velDcdFrequency : numberOfSteps+1) );
   CheckStep forceDcdFrequency;
   forceDcdFrequency.init(step, (simParams->forceDcdFrequency ?
         simParams->forceDcdFrequency : numberOfSteps+1) );
   CheckStep imdFrequency;
   imdFrequency.init(step, (simParams->IMDfreq ?
         simParams->IMDfreq : numberOfSteps+1) );
 
   if ( scriptTask == SCRIPT_RUN ) {
     // enforce rigid bond constraints on initial positions
     TIMER_START(t, RATTLE1);
     rattle1_SOA(0., 0);
     TIMER_STOP(t, RATTLE1);
 
     // must migrate here!
     int natoms = patch->patchDataSOA.numAtoms;
     runComputeObjects_SOA(1, step<numberOfSteps, step);
     // kick -0.5
     TIMER_START(t, KICK);
     addForceToMomentum_SOA(-0.5, timestep, nbondstep, slowstep,
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.recipMass,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms,
 #endif
         maxForceUsed
         );
     TIMER_STOP(t, KICK);
 
     TIMER_START(t, RATTLE1);
     rattle1_SOA(-timestep, 0);
     TIMER_STOP(t, RATTLE1);
 
     TIMER_START(t, SUBMITHALF);
     submitHalfstep_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, SUBMITHALF);
 
     // kick 1.0
     TIMER_START(t, KICK);
     addForceToMomentum_SOA(1.0, timestep, nbondstep, slowstep,
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.recipMass,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms,
 #endif
         maxForceUsed
         );
     TIMER_STOP(t, KICK);
 
     TIMER_START(t, RATTLE1);
     rattle1_SOA(timestep, 1);
     TIMER_STOP(t, RATTLE1);
 
     // save total force in computeGlobal
     if (doGlobal) {
       computeGlobal->saveTotalForces(patch);
     }
 
     TIMER_START(t, SUBMITHALF);
     submitHalfstep_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, SUBMITHALF);
 
     // kick -0.5
     TIMER_START(t, KICK);
     addForceToMomentum_SOA(-0.5, timestep, nbondstep, slowstep,
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.recipMass,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms,
 #endif
         maxForceUsed
         );
     TIMER_STOP(t, KICK);
 
     TIMER_START(t, SUBMITFULL);
     submitReductions_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.pos_x,
         patch->patchDataSOA.pos_y,
         patch->patchDataSOA.pos_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, SUBMITFULL);
 
     rebalanceLoad(step);
   } // scriptTask == SCRIPT_RUN
 
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
   int& eon = patch->flags.event_on;
   int epid = (simParams->beginEventPatchID <= patch->getPatchID()
       && patch->getPatchID() <= simParams->endEventPatchID);
   int beginStep = simParams->beginEventStep;
   int endStep = simParams->endEventStep;
   bool controlProfiling = patch->getPatchID() == 0;
 #endif
 
   for ( ++step; step <= numberOfSteps; ++step ) {
     int dcdSelectionChecks=0;
     Molecule *molecule = Node::Object()->molecule;
     for(int dcdindex=0; dcdindex<16;++dcdindex)
       {
         int dcdSelectionFrequency = molecule->dcdSelectionParams[dcdindex].frequency;
         if(dcdSelectionFrequency && step % dcdSelectionFrequency)
           dcdSelectionChecks++;
       }
     const int isCollection = restartFrequency.check(step) +
       dcdFrequency.check(step) + velDcdFrequency.check(step) +
       forceDcdFrequency.check(step) + imdFrequency.check(step) +
       dcdSelectionChecks;
     const int isMigration = stepsPerCycle.check(step);
     doEnergy = energyFrequency.check(step);
     DebugM(3,"doGlobal now "<< doGlobal<<"\n"<<endi);
 
 #if defined(NAMD_NVTX_ENABLED) || defined(NAMD_CMK_TRACE_ENABLED) || defined(NAMD_ROCTX_ENABLED)
     eon = epid && (beginStep < step && step <= endStep);
 
     if (controlProfiling && step == beginStep) {
       NAMD_PROFILE_START();
     }
     if (controlProfiling && step == endStep) {
       NAMD_PROFILE_STOP();
     }
 //    NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_SOA_1);
     char buf[32];
     sprintf(buf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::INTEGRATE_SOA_1], patch->getPatchID());
     NAMD_EVENT_START_EX(eon, NamdProfileEvent::INTEGRATE_SOA_1, buf);
 #endif
 
     if ( simParams->stochRescaleOn ) {
       stochRescaleVelocities_SOA(step);
     }
 
     if ( simParams->berendsenPressureOn ) {
       berendsenPressure_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.pos_x,
         patch->patchDataSOA.pos_y,
         patch->patchDataSOA.pos_z,
         patch->patchDataSOA.numAtoms,
 #endif
         step);
     }
 
     // kick 0.5
     TIMER_START(t, KICK);
     addForceToMomentum_SOA(0.5, timestep, nbondstep, slowstep,
 #ifndef SOA_SIMPLIFY_PARAMS
                            patch->patchDataSOA.recipMass,
                            patch->patchDataSOA.f_normal_x,
                            patch->patchDataSOA.f_normal_y,
                            patch->patchDataSOA.f_normal_z,
                            patch->patchDataSOA.f_nbond_x,
                            patch->patchDataSOA.f_nbond_y,
                            patch->patchDataSOA.f_nbond_z,
                            patch->patchDataSOA.f_slow_x,
                            patch->patchDataSOA.f_slow_y,
                            patch->patchDataSOA.f_slow_z,
                            patch->patchDataSOA.vel_x,
                            patch->patchDataSOA.vel_y,
                            patch->patchDataSOA.vel_z,
                            patch->patchDataSOA.numAtoms,
 #endif
                            maxForceUsed
       );
     TIMER_STOP(t, KICK);
 
     // maximumMove checks velocity bound on atoms
     TIMER_START(t, MAXMOVE);
     maximumMove_SOA(timestep, maxvel2
 #ifndef SOA_SIMPLIFY_PARAMS
                     ,
                     patch->patchDataSOA.vel_x,
                     patch->patchDataSOA.vel_y,
                     patch->patchDataSOA.vel_z,
                     patch->patchDataSOA.numAtoms
 #endif
       );
     TIMER_STOP(t, MAXMOVE);
 
 
     NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_SOA_1);
 
     // Check to see if Langevin piston is enabled this step:
     //   ! ((step-1-slowFreq/2) % slowFreq)
     if ( langevinPistonFrequency.check(step) ) {
       // if (langevinPistonStep) {
       // drift 0.5
       TIMER_START(t, DRIFT);
       addVelocityToPosition_SOA(0.5*timestep
 #ifndef SOA_SIMPLIFY_PARAMS
                                 ,
                                 patch->patchDataSOA.vel_x,
                                 patch->patchDataSOA.vel_y,
                                 patch->patchDataSOA.vel_z,
                                 patch->patchDataSOA.pos_x,
                                 patch->patchDataSOA.pos_y,
                                 patch->patchDataSOA.pos_z,
                                 patch->patchDataSOA.numAtoms
 #endif
         );
       TIMER_STOP(t, DRIFT);
       // There is a blocking receive inside of langevinPiston()
       // that might suspend the current thread of execution,
       // so split profiling around this conditional block.
       langevinPiston_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.pos_x,
         patch->patchDataSOA.pos_y,
         patch->patchDataSOA.pos_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms,
 #endif
         step
         );
 
       // drift 0.5
       TIMER_START(t, DRIFT);
       addVelocityToPosition_SOA(0.5*timestep
 #ifndef SOA_SIMPLIFY_PARAMS
                                 ,
                                 patch->patchDataSOA.vel_x,
                                 patch->patchDataSOA.vel_y,
                                 patch->patchDataSOA.vel_z,
                                 patch->patchDataSOA.pos_x,
                                 patch->patchDataSOA.pos_y,
                                 patch->patchDataSOA.pos_z,
                                 patch->patchDataSOA.numAtoms
 #endif
         );
       TIMER_STOP(t, DRIFT);
     }
     else {
       // drift 1.0
       TIMER_START(t, DRIFT);
       addVelocityToPosition_SOA(timestep
 #ifndef SOA_SIMPLIFY_PARAMS
                                 ,
                                 patch->patchDataSOA.vel_x,
                                 patch->patchDataSOA.vel_y,
                                 patch->patchDataSOA.vel_z,
                                 patch->patchDataSOA.pos_x,
                                 patch->patchDataSOA.pos_y,
                                 patch->patchDataSOA.pos_z,
                                 patch->patchDataSOA.numAtoms
 #endif
         );
       TIMER_STOP(t, DRIFT);
     }
 
     //NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_SOA_2);
 
     // There are NO sends in submitHalfstep() just local summation
     // into the Reduction struct.
     TIMER_START(t, SUBMITHALF);
     submitHalfstep_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
       patch->patchDataSOA.hydrogenGroupSize,
       patch->patchDataSOA.mass,
       patch->patchDataSOA.vel_x,
       patch->patchDataSOA.vel_y,
       patch->patchDataSOA.vel_z,
       patch->patchDataSOA.numAtoms
 #endif
       );
     TIMER_STOP(t, SUBMITHALF);
 
     //doNonbonded = !(step%nonbondedFrequency);
     doNonbonded = nonbondedFrequency.check(step);
     //doFullElectrostatics = (dofull && !(step%fullElectFrequency));
     doFullElectrostatics = (dofull && fullElectFrequency.check(step));
 
     maxForceUsed = Results::normal;
     if ( doNonbonded ) maxForceUsed = Results::nbond;
     if ( doFullElectrostatics ) maxForceUsed = Results::slow;
 
     // Migrate Atoms on stepsPerCycle
     // Check to see if this is energy evaluation step:
     //   doEnergy = ! ( step % energyFrequency );
     doVirial = 1;
     doKineticEnergy = 1;
     doMomenta = 1;
 
     //NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_SOA_2); // integrate_SOA 2
 
     // The current thread of execution will suspend in runComputeObjects().
     // Check to see if we are at a migration step:
     //   runComputeObjects_SOA(!(step%stepsPerCycle), step<numberOfSteps);
     runComputeObjects_SOA(isMigration, step<numberOfSteps, step);
 
     NAMD_EVENT_START(eon, NamdProfileEvent::INTEGRATE_SOA_3);
 
     TIMER_START(t, VELBBK1);
     langevinVelocitiesBBK1_SOA(
         timestep
 #ifndef SOA_SIMPLIFY_PARAMS
         ,
         patch->patchDataSOA.langevinParam,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, VELBBK1);
 
     // kick 1.0
     TIMER_START(t, KICK);
     addForceToMomentum_SOA(1.0, timestep, nbondstep, slowstep,
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.recipMass,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms,
 #endif
         maxForceUsed
         );
     TIMER_STOP(t, KICK);
 
     TIMER_START(t, VELBBK2);
     langevinVelocitiesBBK2_SOA(
         timestep
 #ifndef SOA_SIMPLIFY_PARAMS
         ,
         patch->patchDataSOA.langevinParam,
         patch->patchDataSOA.langScalVelBBK2,
         patch->patchDataSOA.langScalRandBBK2,
         patch->patchDataSOA.gaussrand_x,
         patch->patchDataSOA.gaussrand_y,
         patch->patchDataSOA.gaussrand_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, VELBBK2);
     
     TIMER_START(t, RATTLE1);
     rattle1_SOA(timestep, 1);
     TIMER_STOP(t, RATTLE1);
 
     // save total force in computeGlobal
     if (doGlobal) {
       computeGlobal->saveTotalForces(patch);
     }
 
     TIMER_START(t, SUBMITHALF);
     submitHalfstep_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, SUBMITHALF);
 
     // kick -0.5
     TIMER_START(t, KICK);
     addForceToMomentum_SOA(-0.5, timestep, nbondstep, slowstep,
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.recipMass,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.numAtoms,
 #endif
         maxForceUsed
         );
     TIMER_STOP(t, KICK);
 
     // XXX rattle2_SOA(timestep,step);
 
     TIMER_START(t, SUBMITFULL);
     submitReductions_SOA(
 #ifndef SOA_SIMPLIFY_PARAMS
         patch->patchDataSOA.hydrogenGroupSize,
         patch->patchDataSOA.mass,
         patch->patchDataSOA.pos_x,
         patch->patchDataSOA.pos_y,
         patch->patchDataSOA.pos_z,
         patch->patchDataSOA.vel_x,
         patch->patchDataSOA.vel_y,
         patch->patchDataSOA.vel_z,
         patch->patchDataSOA.f_normal_x,
         patch->patchDataSOA.f_normal_y,
         patch->patchDataSOA.f_normal_z,
         patch->patchDataSOA.f_nbond_x,
         patch->patchDataSOA.f_nbond_y,
         patch->patchDataSOA.f_nbond_z,
         patch->patchDataSOA.f_slow_x,
         patch->patchDataSOA.f_slow_y,
         patch->patchDataSOA.f_slow_z,
         patch->patchDataSOA.numAtoms
 #endif
         );
     TIMER_STOP(t, SUBMITFULL);
 #ifdef TESTPID
     if (1) {
       int pid = TESTPID;
       if (patch->patchID == pid) {
         const PatchDataSOA& p = patch->patchDataSOA;
         int n = p.numAtoms;
 #if 0
         fprintf(stderr, "Patch %d has %d atoms\n", pid, n);
         fprintf(stderr, "%3s  %8s  %12s  %12s  %12s\n",
             "", "id", "fnormal_x", "fnbond_x", "fslow_x");
         for (int i=0;  i < n;  i++) {
           int index = p.id[i];
           fprintf(stderr, "%3d  %8d  %12.8f  %12.8f  %12.8f\n",
               i, index, p.f_normal_x[i], p.f_nbond_x[i], p.f_slow_x[i]);
         }
 #else
         Vector *f_normal = new Vector[n];
         Vector *f_nbond = new Vector[n];
         Vector *f_slow = new Vector[n];
         for (int i=0;  i < n;  i++) {
           f_normal[i].x = p.f_normal_x[i];
           f_normal[i].y = p.f_normal_y[i];
           f_normal[i].z = p.f_normal_z[i];
           f_nbond[i].x = p.f_nbond_x[i];
           f_nbond[i].y = p.f_nbond_y[i];
           f_nbond[i].z = p.f_nbond_z[i];
           f_slow[i].x = p.f_slow_x[i];
           f_slow[i].y = p.f_slow_y[i];
           f_slow[i].z = p.f_slow_z[i];
         }
         TestArray_write<double>(
             "f_normal_good.bin", "f_normal good", (double*)f_normal, 3*n);
         TestArray_write<double>(
             "f_nbond_good.bin", "f_nbond good", (double*)f_nbond, 3*n);
         TestArray_write<double>(
             "f_slow_good.bin", "f_slow good", (double*)f_slow, 3*n);
         delete [] f_normal;
         delete [] f_nbond;
         delete [] f_slow;
 #endif
       }
     }
 #endif
 
     // Do collections if any checks below are "on."
     // We add because we can't short-circuit.
     TIMER_START(t, SUBMITCOLLECT);
     if (isCollection) {
       submitCollections_SOA(step);
     }
     TIMER_STOP(t, SUBMITCOLLECT);
 
     NAMD_EVENT_STOP(eon, NamdProfileEvent::INTEGRATE_SOA_3); // integrate_SOA 3
 
     rebalanceLoad(step);
   }
 
   patch->copy_updates_to_AOS();
 
   TIMER_DONE(t);
   if (patch->patchID == SPECIAL_PATCH_ID) {
     printf("Timer collection reporting in microseconds for "
         "Patch %d\n", patch->patchID);
     TIMER_REPORT(t);
   }
 }

◆ langevinPiston()

void Sequencer::langevinPiston ( int step )

protected

Definition at line 5340 of file Sequencer.C.

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

 {
   if ( simParams->langevinPistonOn && ! ( (step-1-slowFreq/2) % slowFreq ) )
   {
     //
     // DJH: Loops below simplify if we lift out special cases of fixed atoms
     // and pressure excluded atoms and make them their own branch.
     //
    FullAtom *a = patch->atom.begin();
    int numAtoms = patch->numAtoms;
    // Blocking receive for the updated lattice scaling factor.
    Tensor factor = broadcast->positionRescaleFactor.get(step);
    TIMER_START(patch->timerSet, PISTON);
    // JCP FIX THIS!!!
    Vector velFactor(1/factor.xx,1/factor.yy,1/factor.zz);
    patch->lattice.rescale(factor);
    Molecule *mol = Node::Object()->molecule;
    if ( simParams->useGroupPressure )
    {
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
       int j;
       hgs = a[i].hydrogenGroupSize;
       if ( simParams->fixedAtomsOn && a[i].groupFixed ) {
         for ( j = i; j < (i+hgs); ++j ) {
           a[j].position = patch->lattice.apply_transform(
                                 a[j].fixedPosition,a[j].transform);
         }
         continue;
       }
       BigReal m_cm = 0;
       Position x_cm(0,0,0);
       Velocity v_cm(0,0,0);
       for ( j = i; j < (i+hgs); ++j ) {
         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
         m_cm += a[j].mass;
         x_cm += a[j].mass * a[j].position;
         v_cm += a[j].mass * a[j].velocity;
       }
       x_cm /= m_cm;
       Position new_x_cm = x_cm;
       patch->lattice.rescale(new_x_cm,factor);
       Position delta_x_cm = new_x_cm - x_cm;
       v_cm /= m_cm;
       Velocity delta_v_cm;
       delta_v_cm.x = ( velFactor.x - 1 ) * v_cm.x;
       delta_v_cm.y = ( velFactor.y - 1 ) * v_cm.y;
       delta_v_cm.z = ( velFactor.z - 1 ) * v_cm.z;
       for ( j = i; j < (i+hgs); ++j ) {
         if ( simParams->fixedAtomsOn && a[j].atomFixed ) {
           a[j].position = patch->lattice.apply_transform(
                                 a[j].fixedPosition,a[j].transform);
           continue;
         }
         if ( mol->is_atom_exPressure(a[j].id) ) continue;
         a[j].position += delta_x_cm;
         a[j].velocity += delta_v_cm;
       }
     }
    }
    else
    {
     for ( int i = 0; i < numAtoms; ++i )
     {
       if ( simParams->fixedAtomsOn && a[i].atomFixed ) {
         a[i].position = patch->lattice.apply_transform(
                                 a[i].fixedPosition,a[i].transform);
         continue;
       }
       if ( mol->is_atom_exPressure(a[i].id) ) continue;
       patch->lattice.rescale(a[i].position,factor);
       a[i].velocity.x *= velFactor.x;
       a[i].velocity.y *= velFactor.y;
       a[i].velocity.z *= velFactor.z;
     }
    }
    TIMER_STOP(patch->timerSet, PISTON);
   }
 }

◆ langevinPiston_SOA()

void Sequencer::langevinPiston_SOA ( int step )

protected

Definition at line 3541 of file Sequencer.C.

References broadcast, SimpleBroadcastObject< T >::get(), PatchDataSOA::hydrogenGroupSize, Patch::lattice, PatchDataSOA::mass, namd_reciprocal, PatchDataSOA::numAtoms, Lattice::origin(), patch, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, ControllerBroadcasts::positionRescaleFactor, Lattice::rescale(), simParams, TIMER_START, TIMER_STOP, SimParameters::useGroupPressure, PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, Vector::x, 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_SOA().

 {
 #ifdef SOA_SIMPLIFY_PARAMS
   const int    * __restrict hydrogenGroupSize = patch->patchDataSOA.hydrogenGroupSize;
   const float  * __restrict mass = patch->patchDataSOA.mass;
   double       * __restrict pos_x = patch->patchDataSOA.pos_x;
   double       * __restrict pos_y = patch->patchDataSOA.pos_y;
   double       * __restrict pos_z = patch->patchDataSOA.pos_z;
   double       * __restrict vel_x = patch->patchDataSOA.vel_x;
   double       * __restrict vel_y = patch->patchDataSOA.vel_y;
   double       * __restrict vel_z = patch->patchDataSOA.vel_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
 
   //
   // Loops below simplify if we lift out special cases of fixed atoms
   // and pressure excluded atoms and make them their own branch.
   //
 
   // Blocking receive for the updated lattice scaling factor.
 
   Tensor factor = broadcast->positionRescaleFactor.get(step);
 
   TIMER_START(patch->timerSet, PISTON);
   // JCP FIX THIS!!!
   double velFactor_x = namd_reciprocal(factor.xx);
   double velFactor_y = namd_reciprocal(factor.yy);
   double velFactor_z = namd_reciprocal(factor.zz);
   patch->lattice.rescale(factor);
   Vector origin = patch->lattice.origin();
   if ( simParams->useGroupPressure ) {
     int hgs;
     for (int i=0;  i < numAtoms;  i += hgs) {
       int j;
       hgs = hydrogenGroupSize[i];
       // missing fixed atoms
       BigReal m_cm = 0;
       BigReal r_cm_x = 0;
       BigReal r_cm_y = 0;
       BigReal r_cm_z = 0;
       BigReal v_cm_x = 0;
       BigReal v_cm_y = 0;
       BigReal v_cm_z = 0;
       for ( j = i; j < (i+hgs); ++j ) {
         m_cm += mass[j];
         r_cm_x += mass[j] * pos_x[j];
         r_cm_y += mass[j] * pos_y[j];
         r_cm_z += mass[j] * pos_z[j];
         v_cm_x += mass[j] * vel_x[j];
         v_cm_y += mass[j] * vel_y[j];
         v_cm_z += mass[j] * vel_z[j];
       }
       BigReal inv_m_cm = namd_reciprocal(m_cm);
       r_cm_x *= inv_m_cm;
       r_cm_y *= inv_m_cm;
       r_cm_z *= inv_m_cm;
 
       double tx = r_cm_x - origin.x;
       double ty = r_cm_y - origin.y;
       double tz = r_cm_z - origin.z;
       double new_r_cm_x = factor.xx*tx + factor.xy*ty + factor.xz*tz;
       double new_r_cm_y = factor.yx*tx + factor.yy*ty + factor.yz*tz;
       double new_r_cm_z = factor.zx*tx + factor.zy*ty + factor.zz*tz;
       new_r_cm_x += origin.x;
       new_r_cm_y += origin.y;
       new_r_cm_z += origin.z;
 
       double delta_r_cm_x = new_r_cm_x - r_cm_x;
       double delta_r_cm_y = new_r_cm_y - r_cm_y;
       double delta_r_cm_z = new_r_cm_z - r_cm_z;
       v_cm_x *= inv_m_cm;
       v_cm_y *= inv_m_cm;
       v_cm_z *= inv_m_cm;
       double delta_v_cm_x = ( velFactor_x - 1 ) * v_cm_x;
       double delta_v_cm_y = ( velFactor_y - 1 ) * v_cm_y;
       double delta_v_cm_z = ( velFactor_z - 1 ) * v_cm_z;
       for (j = i;  j < (i+hgs);  j++) {
         pos_x[j] += delta_r_cm_x;
         pos_y[j] += delta_r_cm_y;
         pos_z[j] += delta_r_cm_z;
         vel_x[j] += delta_v_cm_x;
         vel_y[j] += delta_v_cm_y;
         vel_z[j] += delta_v_cm_z;
       }
       // if (i < 10)
       //   printf("cpu: %d, %f, %f, %f, %f, %f, %f\n", i,
       //          pos_x[i], pos_y[i], pos_z[i],
       //          vel_x[i], vel_y[i], vel_z[i]);
     }
   }
   else {
     for (int i=0;  i < numAtoms;  i++) {
       double tx = pos_x[i] - origin.x;
       double ty = pos_y[i] - origin.y;
       double tz = pos_z[i] - origin.z;
       double ftx = factor.xx*tx + factor.xy*ty + factor.xz*tz;
       double fty = factor.yx*tx + factor.yy*ty + factor.yz*tz;
       double ftz = factor.zx*tx + factor.zy*ty + factor.zz*tz;
       pos_x[i] = ftx + origin.x;
       pos_y[i] = fty + origin.y;
       pos_z[i] = ftz + origin.z;
       vel_x[i] *= velFactor_x;
       vel_y[i] *= velFactor_y;
       vel_z[i] *= velFactor_z;
       // if (i < 10)
       //   printf("cpu: %d, %f, %f, %f, %f, %f, %f\n", i,
       //          pos_x[i], pos_y[i], pos_z[i],
       //          vel_x[i], vel_y[i], vel_z[i]);
     }
   }
   TIMER_STOP(patch->timerSet, PISTON);
   // exit(0);
 }

◆ langevinVelocities()

void Sequencer::langevinVelocities ( BigReal dt_fs )

protected

Definition at line 5059 of file Sequencer.C.

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

 {
 // This routine is used for the BAOAB integrator,
 // Ornstein-Uhlenbeck exact solve for the O-part.
 // See B. Leimkuhler and C. Matthews, AMRX (2012)
 // Routine originally written by JPhillips, with fresh errors by CMatthews June2012
 
   if ( simParams->langevinOn && simParams->langevin_useBAOAB )
   {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     Molecule *molecule = Node::Object()->molecule;
     BigReal dt = dt_fs * 0.001;  // convert to ps
     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
     {
         kbT = BOLTZMANN*adaptTempT;
     }
 
     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
 
     for ( int i = 0; i < numAtoms; ++i )
     {
       BigReal dt_gamma = dt * a[i].langevinParam;
       if ( ! dt_gamma ) continue;
 
       BigReal f1 = exp( -dt_gamma );
       BigReal f2 = sqrt( ( 1. - f1*f1 ) * kbT *
                          ( a[i].partition ? tempFactor : 1.0 ) *
                          a[i].recipMass );
       a[i].velocity *= f1;
       a[i].velocity += f2 * random->gaussian_vector();
     }
   }
 }

◆ langevinVelocitiesBBK1()

void Sequencer::langevinVelocitiesBBK1 ( BigReal dt_fs )

protected

Definition at line 5096 of file Sequencer.C.

References ResizeArray< Elem >::begin(), SimParameters::drudeOn, 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().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1);
   if ( simParams->langevinOn && !simParams->langevin_useBAOAB )
   {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     Molecule *molecule = Node::Object()->molecule;
     BigReal dt = dt_fs * 0.001;  // convert to ps
     int i;
 
     if (simParams->drudeOn) {
       for (i = 0;  i < numAtoms;  i++) {
 
         if (i < numAtoms-1 &&
             a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
           //printf("*** Found Drude particle %d\n", a[i+1].id);
           // i+1 is a Drude particle with parent i
 
           // convert from Cartesian coordinates to (COM,bond) coordinates
           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
           BigReal dt_gamma;
 
           // use Langevin damping factor i for v_com
           dt_gamma = dt * a[i].langevinParam;
           if (dt_gamma != 0.0) {
             v_com *= ( 1. - 0.5 * dt_gamma );
           }
 
           // use Langevin damping factor i+1 for v_bnd
           dt_gamma = dt * a[i+1].langevinParam;
           if (dt_gamma != 0.0) {
             v_bnd *= ( 1. - 0.5 * dt_gamma );
           }
 
           // convert back
           a[i].velocity = v_com - m * v_bnd;
           a[i+1].velocity = v_bnd + a[i].velocity;
 
           i++;  // +1 from loop, we've updated both particles
         }
         else {
           BigReal dt_gamma = dt * a[i].langevinParam;
           if ( ! dt_gamma ) continue;
 
           a[i].velocity *= ( 1. - 0.5 * dt_gamma );
         }
 
       } // end for
     } // end if drudeOn
     else {
 
       //
       // DJH: The conditional inside loop prevents vectorization and doesn't
       // avoid much work since addition and multiplication are cheap.
       //
       for ( i = 0; i < numAtoms; ++i )
       {
         BigReal dt_gamma = dt * a[i].langevinParam;
         if ( ! dt_gamma ) continue;
 
         a[i].velocity *= ( 1. - 0.5 * dt_gamma );
       }
 
     } // end else
 
   } // end if langevinOn
 }

◆ langevinVelocitiesBBK1_SOA()

void Sequencer::langevinVelocitiesBBK1_SOA ( BigReal timestep )

protected

Definition at line 3313 of file Sequencer.C.

References Patch::flags, SimParameters::langevinOn, PatchDataSOA::langevinParam, NAMD_EVENT_RANGE_2, PatchDataSOA::numAtoms, patch, simParams, TIMEFACTOR, PatchDataSOA::vel_x, PatchDataSOA::vel_y, and PatchDataSOA::vel_z.

Referenced by integrate_SOA().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::LANGEVIN_VELOCITIES_BBK1_SOA);
 #ifdef SOA_SIMPLIFY_PARAMS
   const float * __restrict langevinParam = patch->patchDataSOA.langevinParam;
   double      * __restrict vel_x = patch->patchDataSOA.vel_x;
   double      * __restrict vel_y = patch->patchDataSOA.vel_y;
   double      * __restrict vel_z = patch->patchDataSOA.vel_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
   if ( simParams->langevinOn /* && !simParams->langevin_useBAOAB */ )
   {
     // scale by TIMEFACTOR to convert to fs and then by 0.001 to ps
     // multiply by the Langevin damping coefficient, units 1/ps
     // XXX we could instead store time-scaled Langevin parameters
     BigReal dt = timestep * (0.001 * TIMEFACTOR);
 
     // XXX missing Drude
 
     //
     // The conditional inside loop prevents vectorization and doesn't
     // avoid much work since addition and multiplication are cheap.
     //
     for (int i=0;  i < numAtoms;  i++) {
       BigReal dt_gamma = dt * langevinParam[i];
       //if ( ! dt_gamma ) continue;
 
       BigReal scaling = 1. - 0.5 * dt_gamma;
       vel_x[i] *= scaling;
       vel_y[i] *= scaling;
       vel_z[i] *= scaling;
     }
   } // end if langevinOn
 }

◆ langevinVelocitiesBBK2()

void Sequencer::langevinVelocitiesBBK2 ( BigReal dt_fs )

protected

Definition at line 5169 of file Sequencer.C.

References SimParameters::adaptTempLangevin, SimParameters::adaptTempOn, adaptTempT, ResizeArray< Elem >::begin(), BOLTZMANN, SimParameters::drudeOn, SimParameters::drudeTemp, 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().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2);
   if ( simParams->langevinOn && !simParams->langevin_useBAOAB )
   {
     //
     // DJH: This call is expensive. Avoid calling when gammas don't differ.
     // Set flag in SimParameters and make this call conditional.
     //
     TIMER_START(patch->timerSet, RATTLE1);
     rattle1(dt_fs,1);  // conserve momentum if gammas differ
     TIMER_STOP(patch->timerSet, RATTLE1);
 
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     Molecule *molecule = Node::Object()->molecule;
     BigReal dt = dt_fs * 0.001;  // convert to ps
     BigReal kbT = BOLTZMANN*(simParams->langevinTemp);
     if (simParams->adaptTempOn && simParams->adaptTempLangevin)
     {
         kbT = BOLTZMANN*adaptTempT;
     }
     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
     int i;
 
     if (simParams->drudeOn) {
       BigReal kbT_bnd = BOLTZMANN*(simParams->drudeTemp);  // drude bond Temp
 
       for (i = 0;  i < numAtoms;  i++) {
 
         if (i < numAtoms-1 &&
             a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
           //printf("*** Found Drude particle %d\n", a[i+1].id);
           // i+1 is a Drude particle with parent i
 
           // convert from Cartesian coordinates to (COM,bond) coordinates
           BigReal m = a[i+1].mass / (a[i].mass + a[i+1].mass);  // mass ratio
           Vector v_bnd = a[i+1].velocity - a[i].velocity;  // vel of bond
           Vector v_com = a[i].velocity + m * v_bnd;  // vel of COM
           BigReal dt_gamma;
 
           // use Langevin damping factor i for v_com
           dt_gamma = dt * a[i].langevinParam;
           if (dt_gamma != 0.0) {
             BigReal mass = a[i].mass + a[i+1].mass;
             v_com += random->gaussian_vector() *
               sqrt( 2 * dt_gamma * kbT *
                   ( a[i].partition ? tempFactor : 1.0 ) / mass );
             v_com /= ( 1. + 0.5 * dt_gamma );
           }
 
           // use Langevin damping factor i+1 for v_bnd
           dt_gamma = dt * a[i+1].langevinParam;
           if (dt_gamma != 0.0) {
             BigReal mass = a[i+1].mass * (1. - m);
             v_bnd += random->gaussian_vector() *
               sqrt( 2 * dt_gamma * kbT_bnd *
                   ( a[i+1].partition ? tempFactor : 1.0 ) / mass );
             v_bnd /= ( 1. + 0.5 * dt_gamma );
           }
 
           // convert back
           a[i].velocity = v_com - m * v_bnd;
           a[i+1].velocity = v_bnd + a[i].velocity;
 
           i++;  // +1 from loop, we've updated both particles
         }
         else {
           BigReal dt_gamma = dt * a[i].langevinParam;
           if ( ! dt_gamma ) continue;
 
           a[i].velocity += random->gaussian_vector() *
             sqrt( 2 * dt_gamma * kbT *
                 ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
           a[i].velocity /= ( 1. + 0.5 * dt_gamma );
         }
 
       } // end for
     } // end if drudeOn
     else {
 
       //
       // DJH: For case using same gamma (the Langevin parameter),
       // no partitions (e.g. FEP), and no adaptive tempering (adaptTempMD),
       // we can precompute constants. Then by lifting the RNG from the
       // loop (filling up an array of random numbers), we can vectorize
       // loop and simplify arithmetic to just addition and multiplication.
       //
       for ( i = 0; i < numAtoms; ++i )
       {
         BigReal dt_gamma = dt * a[i].langevinParam;
         if ( ! dt_gamma ) continue;
 
         a[i].velocity += random->gaussian_vector() *
           sqrt( 2 * dt_gamma * kbT *
               ( a[i].partition ? tempFactor : 1.0 ) * a[i].recipMass );
         a[i].velocity /= ( 1. + 0.5 * dt_gamma );
       }
 
     } // end else
 
   } // end if langevinOn
 }

◆ langevinVelocitiesBBK2_SOA()

void Sequencer::langevinVelocitiesBBK2_SOA ( BigReal timestep )

protected

Definition at line 3359 of file Sequencer.C.

References Patch::flags, Random::gaussian_array_f(), Random::gaussian_vector(), PatchDataSOA::gaussrand_x, PatchDataSOA::gaussrand_y, PatchDataSOA::gaussrand_z, SimParameters::langevinGammasDiffer, SimParameters::langevinOn, PatchDataSOA::langevinParam, PatchDataSOA::langScalRandBBK2, PatchDataSOA::langScalVelBBK2, NAMD_EVENT_RANGE_2, PatchDataSOA::numAtoms, patch, random, rattle1_SOA(), simParams, TIMEFACTOR, TIMER_START, TIMER_STOP, PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, Vector::x, Vector::y, and Vector::z.

Referenced by integrate_SOA().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::LANGEVIN_VELOCITIES_BBK2_SOA);
 #ifdef SOA_SIMPLIFY_PARAMS
   const float * __restrict langevinParam = patch->patchDataSOA.langevinParam;
   const float * __restrict langScalVelBBK2 = patch->patchDataSOA.langScalVelBBK2;
   const float * __restrict langScalRandBBK2 = patch->patchDataSOA.langScalRandBBK2;
   float       * __restrict gaussrand_x = patch->patchDataSOA.gaussrand_x;
   float       * __restrict gaussrand_y = patch->patchDataSOA.gaussrand_y;
   float       * __restrict gaussrand_z = patch->patchDataSOA.gaussrand_z;
   double      * __restrict vel_x = patch->patchDataSOA.vel_x;
   double      * __restrict vel_y = patch->patchDataSOA.vel_y;
   double      * __restrict vel_z = patch->patchDataSOA.vel_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
   if ( simParams->langevinOn /* && !simParams->langevin_useBAOAB */ )
   {
     // XXX missing Drude
 
     // Scale by TIMEFACTOR to convert to fs and then by 0.001 to ps
     // multiply by the Langevin damping coefficient, units 1/ps.
     // XXX we could instead store time-scaled Langevin parameters
     BigReal dt = timestep * (0.001 * TIMEFACTOR);
     // Buffer the Gaussian random numbers
     if ( simParams->langevinGammasDiffer ) {
       // Must re-satisfy constraints if Langevin gammas differ.
       // (conserve momentum?)
       TIMER_START(patch->timerSet, RATTLE1);
       rattle1_SOA(timestep, 1);
       TIMER_STOP(patch->timerSet, RATTLE1);
       //
       // We don't need random numbers for atoms such that gamma=0.
       // If gammas differ, the likely case is that we aren't applying
       // Langevin damping to hydrogen, making those langevinParam=0,
       // in which case we need only numAtoms/3 random vectors.
       //
       // XXX can refine code below, count in advance how many
       // random numbers we need to use Random array filling routine
       //
       // XXX Loop does not vectorize!
       for (int i=0;  i < numAtoms;  i++) {
         Vector rg;  // = 0
         if (langevinParam[i] != 0)  rg = random->gaussian_vector();
         gaussrand_x[i] = float(rg.x);
         gaussrand_y[i] = float(rg.y);
         gaussrand_z[i] = float(rg.z);
       }
     }
     else {
       // Need to completely fill random number arrays.
       random->gaussian_array_f(gaussrand_x, numAtoms);
       random->gaussian_array_f(gaussrand_y, numAtoms);
       random->gaussian_array_f(gaussrand_z, numAtoms);
     }
 
     // do the velocity updates
     for (int i=0;  i < numAtoms;  i++) {
       vel_x[i] += gaussrand_x[i] * langScalRandBBK2[i];
       vel_y[i] += gaussrand_y[i] * langScalRandBBK2[i];
       vel_z[i] += gaussrand_z[i] * langScalRandBBK2[i];
       vel_x[i] *= langScalVelBBK2[i];
       vel_y[i] *= langScalVelBBK2[i];
       vel_z[i] *= langScalVelBBK2[i];
     }
   } // end if langevinOn
 }

◆ maximumMove()

void Sequencer::maximumMove ( BigReal timestep )

protected

Definition at line 5764 of file Sequencer.C.

References ResizeArray< Elem >::begin(), SimParameters::cutoff, 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().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::MAXIMUM_MOVE);
 
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   if ( simParams->maximumMove ) {
     const BigReal dt = timestep / TIMEFACTOR;
     const BigReal maxvel = simParams->maximumMove / dt;
     const BigReal maxvel2 = maxvel * maxvel;
     for ( int i=0; i<numAtoms; ++i ) {
       if ( a[i].velocity.length2() > maxvel2 ) {
         a[i].velocity *= ( maxvel / a[i].velocity.length() );
       }
     }
   } else {
     const BigReal dt = timestep / TIMEFACTOR;
     const BigReal maxvel = simParams->cutoff / dt;
     const BigReal maxvel2 = maxvel * maxvel;
     int killme = 0;
     for ( int i=0; i<numAtoms; ++i ) {
       killme = killme || ( a[i].velocity.length2() > maxvel2 );
     }
     if ( killme ) {
       killme = 0;
       for ( int i=0; i<numAtoms; ++i ) {
         if ( a[i].velocity.length2() > maxvel2 ) {
           ++killme;
           iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
             << ( PDBVELFACTOR * a[i].velocity ) << " (limit is "
             << ( PDBVELFACTOR * maxvel ) << ", atom "
             << i << " of " << numAtoms << " on patch "
             << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
         }
       }
       iout << iERROR <<
         "Atoms moving too fast; simulation has become unstable ("
         << killme << " atoms on patch " << patch->patchID
         << " pe " << CkMyPe() << ").\n" << endi;
       Node::Object()->enableEarlyExit();
       terminate();
     }
   }
 }

◆ maximumMove_SOA()

void Sequencer::maximumMove_SOA	(	const double	dt,
		const double	maxvel2
	)

protected

Parameters

dt	scaled timestep
maxvel2	square of bound on velocity

Definition at line 3255 of file Sequencer.C.

References ResizeArray< Elem >::begin(), Node::enableEarlyExit(), endi(), Patch::flags, CompAtomExt::id, iERROR(), iout, NAMD_EVENT_RANGE_2, PatchDataSOA::numAtoms, Node::Object(), patch, Patch::patchID, PDBVELFACTOR, terminate(), PatchDataSOA::vel_x, PatchDataSOA::vel_y, and PatchDataSOA::vel_z.

Referenced by integrate_SOA().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::MAXIMUM_MOVE_SOA);
 #ifdef SOA_SIMPLIFY_PARAMS
   const double * __restrict vel_x = patch->patchDataSOA.vel_x;
   const double * __restrict vel_y = patch->patchDataSOA.vel_y;
   const double * __restrict vel_z = patch->patchDataSOA.vel_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
 
   // XXX missing maximum move
 
   // Loop vectorizes when replacing logical OR with summing.
   int killme = 0;
   for (int i=0;  i < numAtoms;  i++) {
     BigReal vel2 =
       vel_x[i] * vel_x[i] + vel_y[i] * vel_y[i] + vel_z[i] * vel_z[i];
     killme = killme + ( vel2 > maxvel2 );
   }
   if (killme) {
     // Found at least one atom that is moving too fast.
     // Terminating, so loop performance below doesn't matter.
     // Loop does not vectorize.
     killme = 0;
     for (int i=0;  i < numAtoms;  i++) {
       BigReal vel2 =
         vel_x[i] * vel_x[i] + vel_y[i] * vel_y[i] + vel_z[i] * vel_z[i];
       if (vel2 > maxvel2) {
         const FullAtom *a = patch->atom.begin();
         const Vector vel(vel_x[i], vel_y[i], vel_z[i]);
         const BigReal maxvel = sqrt(maxvel2);
         ++killme;
         iout << iERROR << "Atom " << (a[i].id + 1) << " velocity is "
           << ( PDBVELFACTOR * vel ) << " (limit is "
           << ( PDBVELFACTOR * maxvel ) << ", atom "
           << i << " of " << numAtoms << " on patch "
           << patch->patchID << " pe " << CkMyPe() << ")\n" << endi;
       }
     }
     iout << iERROR <<
       "Atoms moving too fast; simulation has become unstable ("
       << killme << " atoms on patch " << patch->patchID
       << " pe " << CkMyPe() << ").\n" << endi;
     Node::Object()->enableEarlyExit();
     terminate();
   }
 }

◆ minimizationQuenchVelocity()

void Sequencer::minimizationQuenchVelocity ( void )

protected

Definition at line 5809 of file Sequencer.C.

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

Referenced by integrate().

 {
   if ( simParams->minimizeOn ) {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     for ( int i=0; i<numAtoms; ++i ) {
       a[i].velocity = 0.;
     }
   }
 }

◆ minimize()

void Sequencer::minimize ( )

protected

Definition at line 4471 of file Sequencer.C.

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

Referenced by algorithm().

                          {
     //
     // DJH: Copy all data into SOA (structure of arrays)
     // from AOS (array of structures) data structure.
     //
     //patch->copy_all_to_SOA();
 
   const int numberOfSteps = simParams->N;
   const int stepsPerCycle = simParams->stepsPerCycle;
 #if 0 && defined(NODEGROUP_FORCE_REGISTER)
   // XXX DJH: This is a hack that is found to get GPU nonbonded 
   //   force calculation right for --with-single-node-cuda builds
   const int stepsPerCycle_save = stepsPerCycle;
   simParams->stepsPerCycle = 1;
 #endif
   int &step = patch->flags.step;
   step = simParams->firstTimestep;
 
   int &maxForceUsed = patch->flags.maxForceUsed;
   int &maxForceMerged = patch->flags.maxForceMerged;
   maxForceUsed = Results::normal;
   maxForceMerged = Results::normal;
   int &doNonbonded = patch->flags.doNonbonded;
   doNonbonded = 1;
   maxForceUsed = Results::nbond;
   maxForceMerged = Results::nbond;
   const int dofull = ( simParams->fullElectFrequency ? 1 : 0 );
   int &doFullElectrostatics = patch->flags.doFullElectrostatics;
   doFullElectrostatics = dofull;
   if ( dofull ) {
     maxForceMerged = Results::slow;
     maxForceUsed = Results::slow;
   }
   int &doMolly = patch->flags.doMolly;
   doMolly = simParams->mollyOn && doFullElectrostatics;
   int &doMinimize = patch->flags.doMinimize;
   doMinimize = 1;
   // BEGIN LA
   int &doLoweAndersen = patch->flags.doLoweAndersen;
   doLoweAndersen = 0;
   // END LA
 
   int &doGBIS = patch->flags.doGBIS;
   doGBIS = simParams->GBISOn;
 
   int &doLCPO = patch->flags.doLCPO;
   doLCPO = simParams->LCPOOn;
 
   int doTcl = simParams->tclForcesOn;
   int doColvars = simParams->colvarsOn;
   int doGlobal = doTcl || doColvars;
   ComputeGlobal *computeGlobal = Node::Object()->computeMgr->computeGlobalObject;
 
   int &doEnergy = patch->flags.doEnergy;
   doEnergy = 1;
 
   // Do this to stabilize the minimizer, whether or not the user
   // wants rigid bond constraints enabled for dynamics.
   // In order to enforce, we have to call HomePatch::rattle1() directly.
   patch->rattle1(0.,0,0);  // enforce rigid bond constraints on initial positions
 
   if (simParams->lonepairs || simParams->singleTopology) {
     patch->atomMapper->registerIDsFullAtom(
                 patch->atom.begin(),patch->atom.end());
   }
 
   runComputeObjects(1,step<numberOfSteps); // must migrate here!
 
   if ( doGlobal ) {
 #ifdef DEBUG_MINIMIZE
     printf("doTcl = %d   doColvars = %d\n", doTcl, doColvars);
 #endif
     if ( doNonbonded ) saveForce(Results::nbond);
     if ( doFullElectrostatics ) saveForce(Results::slow);
     computeGlobal->saveTotalForces(patch);
   }
 #ifdef DEBUG_MINIMIZE
   else { printf("No computeGlobal\n"); }
 #endif
 
   BigReal fmax2 = TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
 
   submitMinimizeReductions(step,fmax2);
   rebalanceLoad(step);
 
   int downhill = 1;  // start out just fixing bad contacts
   int minSeq = 0;
   for ( ++step; step <= numberOfSteps; ++step ) {
    // Blocking receive for the minimization coefficient.
    BigReal c = broadcast->minimizeCoefficient.get(minSeq++);
 
    if ( downhill ) {
     if ( c ) minimizeMoveDownhill(fmax2);
     else {
       downhill = 0;
       fmax2 *= 10000.;
     }
    }
    if ( ! downhill ) {
     if ( ! c ) {  // new direction
 
       // Blocking receive for the minimization coefficient.
       c = broadcast->minimizeCoefficient.get(minSeq++);
 
       newMinimizeDirection(c);  // v = c * v + f
 
       // Blocking receive for the minimization coefficient.
       c = broadcast->minimizeCoefficient.get(minSeq++);
 
     }  // same direction
     newMinimizePosition(c);  // x = x + c * v
    }
 
     runComputeObjects(!(step%stepsPerCycle),step<numberOfSteps);
     if ( doGlobal ) {
       if ( doNonbonded ) saveForce(Results::nbond);
       if ( doFullElectrostatics ) saveForce(Results::slow);
       computeGlobal->saveTotalForces(patch);
     }
     submitMinimizeReductions(step,fmax2);
     submitCollections(step, 1);  // write out zeros for velocities
     rebalanceLoad(step);
   }
   quenchVelocities();  // zero out bogus velocity
 
   doMinimize = 0;
 
 #if 0
   // when using CUDASOAintegrate, need to update SOA data structures
   if (simParams->CUDASOAintegrateMode && !simParams->useDeviceMigration) {
     patch->copy_atoms_to_SOA();
   }
 #endif
 
 #if 0 && defined(NODEGROUP_FORCE_REGISTER)
   // XXX DJH: all patches in a PE are writing into simParams
   //   so this hack needs a guard
   simParams->stepsPerCycle = stepsPerCycle_save;
 #endif
     //
     // DJH: Copy updates of SOA back into AOS.
     //
     //patch->copy_updates_to_AOS();
 }

◆ minimizeMoveDownhill()

void Sequencer::minimizeMoveDownhill ( BigReal fmax2 )

protected

Definition at line 4617 of file Sequencer.C.

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

                                                   {
 
   FullAtom *a = patch->atom.begin();
   Force *f1 = patch->f[Results::normal].begin();  // includes nbond and slow
   int numAtoms = patch->numAtoms;
 
   for ( int i = 0; i < numAtoms; ++i ) {
     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
     Force f = f1[i];
     if ( f.length2() > fmax2 ) {
       a[i].position += ( 0.1 * f.unit() );
       int hgs = a[i].hydrogenGroupSize;  // 0 if not parent
       for ( int j=1; j<hgs; ++j ) {
         a[++i].position += ( 0.1 * f.unit() );
       }
     }
   }
 
   patch->rattle1(0.,0,0);
 }

◆ monteCarloPressureControl()

void Sequencer::monteCarloPressureControl	(	const int	step,
		const int	doMigration,
		const int	doEnergy,
		const int	doVirial,
		const int	maxForceNumber,
		const int	doGlobal
	)

protected

◆ multigratorPressure()

void Sequencer::multigratorPressure	(	int	step,
		int	callNumber
	)

protected

(stepstepsPerCycle)

Definition at line 4826 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, ResizeArray< Elem >::begin(), broadcast, calcFixVirial(), Flags::doFullElectrostatics, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, SimpleBroadcastObject< T >::get(), CompAtom::hydrogenGroupSize, Tensor::identity(), SubmitReduction::item(), 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().

                                                             {
 // Calculate new positions, momenta, and volume using positionRescaleFactor and
 // velocityRescaleTensor values returned from Controller::multigratorPressureCalcScale()
   if (simParams->multigratorOn && !(step % simParams->multigratorPressureFreq)) {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
 
     // Blocking receive (get) scaling factors from Controller
     Tensor scaleTensor    = (callNumber == 1) ? broadcast->positionRescaleFactor.get(step) : broadcast->positionRescaleFactor2.get(step);
     Tensor velScaleTensor = (callNumber == 1) ? broadcast->velocityRescaleTensor.get(step) : broadcast->velocityRescaleTensor2.get(step);
     Tensor posScaleTensor = scaleTensor;
     posScaleTensor -= Tensor::identity();
     if (simParams->useGroupPressure) {
       velScaleTensor -= Tensor::identity();
     }
 
     // Scale volume
     patch->lattice.rescale(scaleTensor);
     // Scale positions and velocities
     scalePositionsVelocities(posScaleTensor, velScaleTensor);
 
     if (!patch->flags.doFullElectrostatics) NAMD_bug("Sequencer::multigratorPressure, doFullElectrostatics must be true");
 
     // Calculate new forces
     // NOTE: We should not need to migrate here since any migration should have happened in the
     // previous call to runComputeObjects inside the MD loop in Sequencer::integrate()
     const int numberOfSteps = simParams->N;
     const int stepsPerCycle = simParams->stepsPerCycle;
     runComputeObjects(0 , step<numberOfSteps, 1);
 
     reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
     reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
 
     // Virials etc.
     Tensor virialNormal;
     Tensor momentumSqrSum;
     BigReal kineticEnergy = 0;
     if ( simParams->pairInteractionOn ) {
       if ( simParams->pairInteractionSelf ) {
         for ( int i = 0; i < numAtoms; ++i ) {
           if ( a[i].partition != 1 ) continue;
           kineticEnergy += a[i].mass * a[i].velocity.length2();
           virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
         }
       }
     } else {
       for ( int i = 0; i < numAtoms; ++i ) {
         if (a[i].mass < 0.01) continue;
         kineticEnergy += a[i].mass * a[i].velocity.length2();
         virialNormal.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
       }
     }
     if (!simParams->useGroupPressure) momentumSqrSum = virialNormal;
     kineticEnergy *= 0.5;
     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
     ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virialNormal);
 
     if ( simParams->fixedAtomsOn ) {
       Tensor fixVirialNormal;
       Tensor fixVirialNbond;
       Tensor fixVirialSlow;
       Vector fixForceNormal = 0;
       Vector fixForceNbond = 0;
       Vector fixForceSlow = 0;
 
       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
 
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, fixVirialNormal);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, fixVirialNbond);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, fixVirialSlow);
       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NORMAL, fixForceNormal);
       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_NBOND, fixForceNbond);
       ADD_VECTOR_OBJECT(reduction, REDUCTION_EXT_FORCE_SLOW, fixForceSlow);
     }
 
     // Internal virial and group momentum
     Tensor intVirialNormal;
     Tensor intVirialNormal2;
     Tensor intVirialNbond;
     Tensor intVirialSlow;
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
       hgs = a[i].hydrogenGroupSize;
       int j;
       BigReal m_cm = 0;
       Position x_cm(0,0,0);
       Velocity v_cm(0,0,0);
       for ( j = i; j < (i+hgs); ++j ) {
         m_cm += a[j].mass;
         x_cm += a[j].mass * a[j].position;
         v_cm += a[j].mass * a[j].velocity;
       }
       if (simParams->useGroupPressure) momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
       x_cm /= m_cm;
       v_cm /= m_cm;
       if (simParams->fixedAtomsOn) NAMD_bug("Sequencer::multigratorPressure, simParams->fixedAtomsOn not implemented yet");
       if ( simParams->pairInteractionOn ) {
         if ( simParams->pairInteractionSelf ) {
           NAMD_bug("Sequencer::multigratorPressure, this part needs to be implemented correctly");
           for ( j = i; j < (i+hgs); ++j ) {
             if ( a[j].partition != 1 ) continue;
             BigReal mass = a[j].mass;
             Vector v = a[j].velocity;
             Vector dv = v - v_cm;
             intVirialNormal2.outerAdd (mass, v, dv);
             Vector dx = a[j].position - x_cm;
             intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
             intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
             intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
           }
         }
       } else {
         for ( j = i; j < (i+hgs); ++j ) {
           BigReal mass = a[j].mass;
           Vector v = a[j].velocity;
           Vector dv = v - v_cm;
           intVirialNormal2.outerAdd(mass, v, dv);
           Vector dx = a[j].position - x_cm;
           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
         }
       }
     }
 
     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal2);
     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NBOND, intVirialNbond);
     ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_SLOW, intVirialSlow);
     ADD_TENSOR_OBJECT(reduction, REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
 
     reduction->submit();
   }
 }

◆ multigratorTemperature()

void Sequencer::multigratorTemperature	(	int	step,
		int	callNumber
	)

protected

Definition at line 4989 of file Sequencer.C.

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

                                                                {
   if (simParams->multigratorOn && !(step % simParams->multigratorTemperatureFreq)) {
     // Blocking receive (get) velocity scaling factor.
     BigReal velScale = (callNumber == 1) ? broadcast->velocityRescaleFactor.get(step) : broadcast->velocityRescaleFactor2.get(step);
     scaleVelocities(velScale);
     // Calculate new kineticEnergy
     BigReal kineticEnergy = calcKineticEnergy();
     multigratorReduction->item(MULTIGRATOR_REDUCTION_KINETIC_ENERGY) += kineticEnergy;
     if (callNumber == 1 && !(step % simParams->multigratorPressureFreq)) {
       // If this is a pressure cycle, calculate new momentum squared sum
       FullAtom *a = patch->atom.begin();
       int numAtoms = patch->numAtoms;
       Tensor momentumSqrSum;
       if (simParams->useGroupPressure) {
         int hgs;
         for ( int i = 0; i < numAtoms; i += hgs ) {
           hgs = a[i].hydrogenGroupSize;
           int j;
           BigReal m_cm = 0;
           Position x_cm(0,0,0);
           Velocity v_cm(0,0,0);
           for ( j = i; j < (i+hgs); ++j ) {
             m_cm += a[j].mass;
             x_cm += a[j].mass * a[j].position;
             v_cm += a[j].mass * a[j].velocity;
           }
           momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
         }
       } else {
         for ( int i = 0; i < numAtoms; i++) {
           momentumSqrSum.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
         }
       }
       ADD_TENSOR_OBJECT(multigratorReduction, MULTIGRATOR_REDUCTION_MOMENTUM_SQUARED, momentumSqrSum);
     }
     // Submit reductions (kineticEnergy and, if applicable, momentumSqrSum)
     multigratorReduction->submit();
 
   }
 }

◆ newMinimizeDirection()

void Sequencer::newMinimizeDirection ( BigReal c )

protected

Definition at line 4639 of file Sequencer.C.

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

                                               {
   FullAtom *a = patch->atom.begin();
   Force *f1 = patch->f[Results::normal].begin(); // includes nbond and slow
   const bool fixedAtomsOn = simParams->fixedAtomsOn;
   const bool drudeHardWallOn = simParams->drudeHardWallOn;
   int numAtoms = patch->numAtoms;
   BigReal maxv2 = 0.;
 
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].velocity *= c;
     a[i].velocity += f1[i];
     if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
       a[i].velocity = a[i-1].velocity;
     }
     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
     BigReal v2 = a[i].velocity.length2();
     if ( v2 > maxv2 ) maxv2 = v2;
   }
 
   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial); }
 
   maxv2 = 0.;
   for ( int i = 0; i < numAtoms; ++i ) {
     if ( drudeHardWallOn && i && (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
       a[i].velocity = a[i-1].velocity;
     }
     if ( fixedAtomsOn && a[i].atomFixed ) a[i].velocity = 0;
     BigReal v2 = a[i].velocity.length2();
     if ( v2 > maxv2 ) maxv2 = v2;
   }
 
   min_reduction->max(0,maxv2);
   min_reduction->submit();
 
   // prevent hydrogens from being left behind
   BigReal fmax2 = 0.01 * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR;
   // int adjustCount = 0;
   int hgs;
   for ( int i = 0; i < numAtoms; i += hgs ) {
     hgs = a[i].hydrogenGroupSize;
     BigReal minChildVel = a[i].velocity.length2();
     if ( minChildVel < fmax2 ) continue;
     int adjustChildren = 1;
     for ( int j = i+1; j < (i+hgs); ++j ) {
       if ( a[j].velocity.length2() > minChildVel ) adjustChildren = 0;
     }
     if ( adjustChildren ) {
       // if ( hgs > 1 ) ++adjustCount;
       for ( int j = i+1; j < (i+hgs); ++j ) {
         if (a[i].mass < 0.01) continue;  // lone pair
         a[j].velocity = a[i].velocity;
       }
     }
   }
   // if (adjustCount) CkPrintf("Adjusting %d hydrogen groups\n", adjustCount);
 
 }

◆ newMinimizePosition()

void Sequencer::newMinimizePosition ( BigReal c )

protected

Definition at line 4698 of file Sequencer.C.

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

Referenced by minimize().

                                              {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
 
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].position += c * a[i].velocity;
   }
 
   if ( simParams->drudeHardWallOn ) {
     for ( int i = 1; i < numAtoms; ++i ) {
       if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
         a[i].position -= a[i-1].position;
       }
     }
   }
 
   patch->rattle1(0.,0,0);
 
   if ( simParams->drudeHardWallOn ) {
     for ( int i = 1; i < numAtoms; ++i ) {
       if ( (0.05 < a[i].mass) && ((a[i].mass < 1.0)) ) { // drude particle
         a[i].position += a[i-1].position;
       }
     }
   }
 }

◆ newtonianVelocities()

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

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

Referenced by integrate().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::NEWTONIAN_VELOCITIES);
 
   // Deterministic velocity update, account for multigrator
   if (staleForces || (doNonbonded && doFullElectrostatics)) {
     addForceToMomentum3(stepscale*timestep, Results::normal, 0,
                         stepscale*nbondstep, Results::nbond, staleForces,
                         stepscale*slowstep, Results::slow, staleForces);
   } else {
     addForceToMomentum(stepscale*timestep);
     if (staleForces || doNonbonded)
       addForceToMomentum(stepscale*nbondstep, Results::nbond, staleForces);
     if (staleForces || doFullElectrostatics)
       addForceToMomentum(stepscale*slowstep, Results::slow, staleForces);
   }
 }

◆ printDevicePatchMap()

void Sequencer::printDevicePatchMap ( )

protected

◆ quenchVelocities()

void Sequencer::quenchVelocities ( )

protected

Definition at line 4726 of file Sequencer.C.

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

Referenced by minimize().

                                  {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
 
   for ( int i = 0; i < numAtoms; ++i ) {
     a[i].velocity = 0;
   }
 }

◆ rattle1()

void Sequencer::rattle1	(	BigReal	dt,
		int	pressure
	)

protected

Definition at line 5707 of file Sequencer.C.

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

Referenced by integrate(), and langevinVelocitiesBBK2().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::RATTLE1);
   if ( simParams->rigidBonds != RIGID_NONE ) {
     Tensor virial;
     Tensor *vp = ( pressure ? &virial : 0 );
     if ( patch->rattle1(dt, vp, pressureProfileReduction) ) {
       iout << iERROR <<
         "Constraint failure; simulation has become unstable.\n" << endi;
       Node::Object()->enableEarlyExit();
       terminate();
     }
 #if 0
     printf("virial = %g %g %g  %g %g %g  %g %g %g\n",
         virial.xx, virial.xy, virial.xz,
         virial.yx, virial.yy, virial.yz,
         virial.zx, virial.zy, virial.zz);
 #endif
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
 #if 0
     {
       const FullAtom *a = patch->atom.const_begin();
       for (int n=0;  n < patch->numAtoms;  n++) {
         printf("pos[%d] =  %g %g %g\n", n,
             a[n].position.x, a[n].position.y, a[n].position.z);
       }
       for (int n=0;  n < patch->numAtoms;  n++) {
         printf("vel[%d] =  %g %g %g\n", n,
             a[n].velocity.x, a[n].velocity.y, a[n].velocity.z);
       }
       if (pressure) {
         for (int n=0;  n < patch->numAtoms;  n++) {
           printf("force[%d] =  %g %g %g\n", n,
               patch->f[Results::normal][n].x,
               patch->f[Results::normal][n].y,
               patch->f[Results::normal][n].z);
         }
       }
     }
 #endif
   }
 }

◆ rattle1_SOA()

void Sequencer::rattle1_SOA	(	BigReal	timestep,
		int	pressure
	)

protected

Definition at line 3670 of file Sequencer.C.

References ADD_TENSOR_OBJECT, SimParameters::CUDASOAintegrate, Node::enableEarlyExit(), endi(), Patch::flags, iERROR(), iout, NAMD_EVENT_RANGE_2, Node::Object(), patch, pressureProfileReduction, HomePatch::rattle1_SOA(), reduction, RIGID_NONE, SimParameters::rigidBonds, simParams, and terminate().

Referenced by integrate_SOA(), and langevinVelocitiesBBK2_SOA().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::RATTLE1_SOA);
   if ( simParams->rigidBonds != RIGID_NONE ) {
     Tensor virial;
     Tensor *vp = ( pressure ? &virial : 0 );
     // XXX pressureProfileReduction == NULL?
     if ( patch->rattle1_SOA(timestep, vp, pressureProfileReduction) ) {
       iout << iERROR <<
         "Constraint failure; simulation has become unstable.\n" << endi;
       Node::Object()->enableEarlyExit();
       terminate();
     }
 #ifdef NODEGROUP_FORCE_REGISTER
     if(simParams->CUDASOAintegrate){
       ADD_TENSOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_VIRIAL_NORMAL,virial);
     }
     else{
 #endif
       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
 #ifdef NODEGROUP_FORCE_REGISTER
     }
 #endif
   }
 }

◆ reassignVelocities()

void Sequencer::reassignVelocities	(	BigReal	timestep,
		int	step
	)

protected

Definition at line 5498 of file Sequencer.C.

References CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), BOLTZMANN, SimParameters::fixedAtomsOn, Random::gaussian_vector(), 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().

 {
   const int reassignFreq = simParams->reassignFreq;
   if ( ( reassignFreq > 0 ) && ! ( step % reassignFreq ) ) {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     BigReal newTemp = simParams->reassignTemp;
     newTemp += ( step / reassignFreq ) * simParams->reassignIncr;
     if ( simParams->reassignIncr > 0.0 ) {
       if ( newTemp > simParams->reassignHold && simParams->reassignHold > 0.0 )
         newTemp = simParams->reassignHold;
     } else {
       if ( newTemp < simParams->reassignHold )
         newTemp = simParams->reassignHold;
     }
     BigReal kbT = BOLTZMANN * newTemp;
 
     int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
     BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
 
     for ( int i = 0; i < numAtoms; ++i )
     {
       a[i].velocity = ( ( simParams->fixedAtomsOn &&
             a[i].atomFixed && a[i].mass > 0.) ? Vector(0,0,0) :
           sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) *
           random->gaussian_vector() );
     }
   } else {
     NAMD_bug("Sequencer::reassignVelocities called improperly!");
   }
 }

◆ rebalanceLoad()

void Sequencer::rebalanceLoad ( int timestep )

protected

Definition at line 6650 of file Sequencer.C.

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

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

                                           {
   if ( ! ldbSteps ) {
     ldbSteps = LdbCoordinator::Object()->getNumStepsToRun();
   }
   if ( ! --ldbSteps ) {
     patch->submitLoadStats(timestep);
     ldbCoordinator->rebalance(this,patch->getPatchID());
     pairlistsAreValid = 0;
   }
 }

◆ reinitVelocities()

void Sequencer::reinitVelocities ( void )

protected

Definition at line 5530 of file Sequencer.C.

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

Referenced by algorithm().

 {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   BigReal newTemp = simParams->initialTemp;
   BigReal kbT = BOLTZMANN * newTemp;
 
   int lesReduceTemp = simParams->lesOn && simParams->lesReduceTemp;
   BigReal tempFactor = lesReduceTemp ? 1.0 / simParams->lesFactor : 1.0;
 
   for ( int i = 0; i < numAtoms; ++i )
   {
     a[i].velocity = ( ( (simParams->fixedAtomsOn && a[i].atomFixed) ||
           a[i].mass <= 0.) ? Vector(0,0,0) :
         sqrt(kbT * (a[i].partition ? tempFactor : 1.0) * a[i].recipMass) *
         random->gaussian_vector() );
     if ( simParams->drudeOn && i+1 < numAtoms && a[i+1].mass < 1.0 && a[i+1].mass > 0.05 ) {
       a[i+1].velocity = a[i].velocity;  // zero is good enough
       ++i;
     }
   }
 }

◆ reloadCharges()

void Sequencer::reloadCharges ( )

protected

Definition at line 5563 of file Sequencer.C.

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

Referenced by algorithm().

 {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   Molecule *molecule = Node::Object()->molecule;
   for ( int i = 0; i < numAtoms; ++i )
   {
     a[i].charge = molecule->atomcharge(a[i].id);
   }
 }

◆ rescaleaccelMD()

void Sequencer::rescaleaccelMD	(	int	step,
		int	doNonbonded,
		int	doFullElectrostatics
	)

protected

Definition at line 5439 of file Sequencer.C.

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

Referenced by integrate().

 {
    if (!simParams->accelMDOn) return;
    if ((step < simParams->accelMDFirstStep) || ( simParams->accelMDLastStep >0 && step > simParams->accelMDLastStep)) return;
 
    // Blocking receive for the Accelerated MD scaling factors.
    Vector accelMDfactor = broadcast->accelMDRescaleFactor.get(step);
    const BigReal factor_dihe = accelMDfactor[0];
    const BigReal factor_tot  = accelMDfactor[1];
    const int numAtoms = patch->numAtoms;
 
    if (simParams->accelMDdihe && factor_tot <1 )
        NAMD_die("accelMD broadcasting error!\n");
    if (!simParams->accelMDdihe && !simParams->accelMDdual && factor_dihe <1 )
        NAMD_die("accelMD broadcasting error!\n");
 
    if (simParams->accelMDdihe && factor_dihe < 1) {
         for (int i = 0; i < numAtoms; ++i)
           if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
               patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - 1);
    }
 
    if ( !simParams->accelMDdihe && factor_tot < 1) {
         for (int i = 0; i < numAtoms; ++i)
             patch->f[Results::normal][i] *= factor_tot;
         if (doNonbonded) {
             for (int i = 0; i < numAtoms; ++i)
                  patch->f[Results::nbond][i] *= factor_tot;
         }
         if (doFullElectrostatics) {
             for (int i = 0; i < numAtoms; ++i)
                  patch->f[Results::slow][i] *= factor_tot;
         }
    }
 
    if (simParams->accelMDdual && factor_dihe < 1) {
         for (int i = 0; i < numAtoms; ++i)
            if (patch->f[Results::amdf][i][0] || patch->f[Results::amdf][i][1] || patch->f[Results::amdf][i][2])
                patch->f[Results::normal][i] += patch->f[Results::amdf][i]*(factor_dihe - factor_tot);
    }
 
 }

◆ rescaleSoluteCharges()

void Sequencer::rescaleSoluteCharges ( BigReal factor )

protected

Definition at line 5575 of file Sequencer.C.

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

Referenced by algorithm(), and Sequencer().

 {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   Molecule *molecule = Node::Object()->molecule;
   BigReal sqrt_factor = sqrt(factor);
   // apply scaling to the original charge (stored in molecule)
   // of just the marked solute atoms
   for ( int i = 0; i < numAtoms; ++i ) {
     if (molecule->get_ss_type(a[i].id)) {
       a[i].charge = sqrt_factor * molecule->atomcharge(a[i].id);
       if (simParams->SOAintegrateOn) patch->patchDataSOA.charge[i] = a[i].charge;
     }
   }
 }

◆ rescaleVelocities()

void Sequencer::rescaleVelocities ( int step )

protected

Definition at line 5420 of file Sequencer.C.

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

Referenced by integrate().

 {
   const int rescaleFreq = simParams->rescaleFreq;
   if ( rescaleFreq > 0 ) {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     ++rescaleVelocities_numTemps;
     if ( rescaleVelocities_numTemps == rescaleFreq ) {
       // Blocking receive for the velcity scaling factor.
       BigReal factor = broadcast->velocityRescaleFactor.get(step);
       for ( int i = 0; i < numAtoms; ++i )
       {
         a[i].velocity *= factor;
       }
       rescaleVelocities_numTemps = 0;
     }
   }
 }

◆ rescaleVelocitiesByFactor()

void Sequencer::rescaleVelocitiesByFactor ( BigReal factor )

protected

Definition at line 5553 of file Sequencer.C.

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

Referenced by algorithm().

 {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   for ( int i = 0; i < numAtoms; ++i )
   {
     a[i].velocity *= factor;
   }
 }

◆ run()

void Sequencer::run ( void )

Definition at line 257 of file Sequencer.C.

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

Referenced by HomePatch::runSequencer().

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

◆ runComputeObjects()

void Sequencer::runComputeObjects	(	int	migration = `1`,
		int	pairlists = `0`,
		int	pressureStep = `0`
	)

protected

Definition at line 6436 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ResizeArray< Elem >::begin(), COMPUTE_HOME_PRIORITY, Flags::doFullElectrostatics, Flags::doGBIS, Flags::doLoweAndersen, Flags::doMolly, Flags::doNonbonded, Patch::f, Patch::flags, SimParameters::fullElectFrequency, GB1_COMPUTE_HOME_PRIORITY, GB2_COMPUTE_HOME_PRIORITY, HomePatch::gbisComputeAfterP1(), HomePatch::gbisComputeAfterP2(), Patch::getPatchID(), SimParameters::lonepairs, HomePatch::loweAndersenFinish(), HomePatch::mollyMollify(), Results::nbond, Results::normal, Patch::numAtoms, Patch::p, pairlistsAge, pairlistsAgeLimit, pairlistsAreValid, patch, PATCH_PRIORITY, Patch::pExt, HomePatch::positionsReady(), reduction, Flags::savePairlists, Flags::sequence, simParams, SimParameters::singleTopology, Results::slow, suspend(), SWM4, TIP4, Flags::usePairlists, SimParameters::usePairlists, and SimParameters::watmodel.

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

 {
   if ( migration ) pairlistsAreValid = 0;
 #if defined(NAMD_CUDA) || defined(NAMD_HIP) || defined(NAMD_MIC)
   if ( pairlistsAreValid &&
        ( patch->flags.doFullElectrostatics || ! simParams->fullElectFrequency )
                          && ( pairlistsAge > pairlistsAgeLimit ) ) {
     pairlistsAreValid = 0;
   }
 #else
   if ( pairlistsAreValid && ( pairlistsAge > pairlistsAgeLimit ) ) {
     pairlistsAreValid = 0;
   }
 #endif
   if ( ! simParams->usePairlists ) pairlists = 0;
   patch->flags.usePairlists = pairlists || pairlistsAreValid;
   patch->flags.savePairlists =
         pairlists && ! pairlistsAreValid;
 
   if ( simParams->singleTopology ) patch->reposition_all_alchpairs();
   if ( simParams->lonepairs ) patch->reposition_all_lonepairs();
 
   //
   // DJH: Copy updates of SOA back into AOS.
   // The positionsReady() routine starts force computation and atom migration.
   //
   // We could reduce amount of copying here by checking migration status
   // and copying velocities only when migrating. Some types of simulation
   // always require velocities, such as Lowe-Anderson.
   //
   //patch->copy_updates_to_AOS();
 
   patch->positionsReady(migration);  // updates flags.sequence
 
   int seq = patch->flags.sequence;
   int basePriority = ( (seq & 0xffff) << 15 )
                      + PATCH_PRIORITY(patch->getPatchID());
   if ( patch->flags.doGBIS && patch->flags.doNonbonded) {
     priority = basePriority + GB1_COMPUTE_HOME_PRIORITY;
     suspend(); // until all deposit boxes close
     patch->gbisComputeAfterP1();
     priority = basePriority + GB2_COMPUTE_HOME_PRIORITY;
     suspend();
     patch->gbisComputeAfterP2();
     priority = basePriority + COMPUTE_HOME_PRIORITY;
     suspend();
   } else {
     priority = basePriority + COMPUTE_HOME_PRIORITY;
     suspend(); // until all deposit boxes close
   }
 
   //
   // DJH: Copy all data into SOA from AOS.
   //
   // We need everything copied after atom migration.
   // When doing force computation without atom migration,
   // all data except forces will already be up-to-date in SOA
   // (except maybe for some special types of simulation).
   //
   //patch->copy_all_to_SOA();
 
   //
   // DJH: Copy forces to SOA.
   // Force available after suspend() has returned.
   //
   //patch->copy_forces_to_SOA();
 
   if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
     pairlistsAreValid = 1;
     pairlistsAge = 0;
   }
   // For multigrator, do not age pairlist during pressure step
   // NOTE: for non-multigrator pressureStep = 0 always
   if ( pairlistsAreValid && !pressureStep ) ++pairlistsAge;
 
   if (simParams->lonepairs) {
     {
       Tensor virial;
       patch->redistrib_lonepair_forces(Results::normal, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
     }
     if (patch->flags.doNonbonded) {
       Tensor virial;
       patch->redistrib_lonepair_forces(Results::nbond, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
     }
     if (patch->flags.doFullElectrostatics) {
       Tensor virial;
       patch->redistrib_lonepair_forces(Results::slow, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
     }
   } else if (simParams->watmodel == WaterModel::TIP4) {
     {
       Tensor virial;
       patch->redistrib_tip4p_forces(Results::normal, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
     }
     if (patch->flags.doNonbonded) {
       Tensor virial;
       patch->redistrib_tip4p_forces(Results::nbond, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
     }
     if (patch->flags.doFullElectrostatics) {
       Tensor virial;
       patch->redistrib_tip4p_forces(Results::slow, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
     }
   } else if (simParams->watmodel == WaterModel::SWM4) {
     {
       Tensor virial;
       patch->redistrib_swm4_forces(Results::normal, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NORMAL, virial);
     }
     if (patch->flags.doNonbonded) {
       Tensor virial;
       patch->redistrib_swm4_forces(Results::nbond, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_NBOND, virial);
     }
     if (patch->flags.doFullElectrostatics) {
       Tensor virial;
       patch->redistrib_swm4_forces(Results::slow, &virial);
       ADD_TENSOR_OBJECT(reduction, REDUCTION_VIRIAL_SLOW, virial);
     }
   }
 
   if (simParams->singleTopology) {
     patch->redistrib_alchpair_forces(Results::normal);
     if (patch->flags.doNonbonded) {
       patch->redistrib_alchpair_forces(Results::nbond);
     }
     if (patch->flags.doFullElectrostatics) {
       patch->redistrib_alchpair_forces(Results::slow);
     }
   }
 
   if ( patch->flags.doMolly ) {
     Tensor virial;
     patch->mollyMollify(&virial);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,virial);
   }
 
 
   // BEGIN LA
   if (patch->flags.doLoweAndersen) {
       patch->loweAndersenFinish();
   }
   // END LA
 //TODO:HIP check if this applies to HIP
 #ifdef NAMD_CUDA_XXX
   int numAtoms = patch->numAtoms;
   FullAtom *a = patch->atom.begin();
   for ( int i=0; i<numAtoms; ++i ) {
     CkPrintf("%d %g %g %g\n", a[i].id,
         patch->f[Results::normal][i].x +
         patch->f[Results::nbond][i].x +
         patch->f[Results::slow][i].x,
         patch->f[Results::normal][i].y +
         patch->f[Results::nbond][i].y +
         patch->f[Results::slow][i].y,
         patch->f[Results::normal][i].z +
         patch->f[Results::nbond][i].z +
         patch->f[Results::slow][i].z);
     CkPrintf("%d %g %g %g\n", a[i].id,
         patch->f[Results::normal][i].x,
         patch->f[Results::nbond][i].x,
         patch->f[Results::slow][i].x);
     CkPrintf("%d %g %g %g\n", a[i].id,
         patch->f[Results::normal][i].y,
         patch->f[Results::nbond][i].y,
         patch->f[Results::slow][i].y);
     CkPrintf("%d %g %g %g\n", a[i].id,
         patch->f[Results::normal][i].z,
         patch->f[Results::nbond][i].z,
         patch->f[Results::slow][i].z);
   }
 #endif
 
 //#undef PRINT_FORCES
 //#define PRINT_FORCES 1
 #if PRINT_FORCES
   int numAtoms = patch->numAtoms;
   FullAtom *a = patch->atom.begin();
   for ( int i=0; i<numAtoms; ++i ) {
     float fxNo = patch->f[Results::normal][i].x;
     float fxNb = patch->f[Results::nbond][i].x;
     float fxSl = patch->f[Results::slow][i].x;
     float fyNo = patch->f[Results::normal][i].y;
     float fyNb = patch->f[Results::nbond][i].y;
     float fySl = patch->f[Results::slow][i].y;
     float fzNo = patch->f[Results::normal][i].z;
     float fzNb = patch->f[Results::nbond][i].z;
     float fzSl = patch->f[Results::slow][i].z;
     float fx = fxNo+fxNb+fxSl;
     float fy = fyNo+fyNb+fySl;
     float fz = fzNo+fzNb+fzSl;
 
                 float f = sqrt(fx*fx+fy*fy+fz*fz);
     int id = patch->pExt[i].id;
     int seq = patch->flags.sequence;
     float x = patch->p[i].position.x;
     float y = patch->p[i].position.y;
     float z = patch->p[i].position.z;
     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <<% .4e, % .4e, % .4e>>\n", seq,id,
     CkPrintf("FORCE(%04i)[%04i] = % .9e % .9e % .9e\n", seq,id,
     //CkPrintf("FORCE(%04i)[%04i] = <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e> <% .4e, % .4e, % .4e>\n", seq,id,
 //fxNo,fyNo,fzNo,
 fxNb,fyNb,fzNb
 //fxSl,fySl,fzSl,
 //fx,fy,fz
 );
         }
 #endif
 }

◆ runComputeObjects_SOA()

void Sequencer::runComputeObjects_SOA	(	int	migration,
		int	pairlists,
		int	step
	)

protected

Definition at line 3696 of file Sequencer.C.

References PatchDataSOA::charge, COMPUTE_HOME_PRIORITY, SimParameters::CUDASOAintegrate, Flags::doFullElectrostatics, Flags::doNonbonded, PatchDataSOA::f_nbond_x, PatchDataSOA::f_nbond_y, PatchDataSOA::f_nbond_z, PatchDataSOA::f_normal_x, PatchDataSOA::f_normal_y, PatchDataSOA::f_normal_z, PatchDataSOA::f_slow_x, PatchDataSOA::f_slow_y, PatchDataSOA::f_slow_z, Patch::flags, SimParameters::fullElectFrequency, Patch::getPatchID(), Patch::numAtoms, PatchDataSOA::numAtoms, pairlistsAge, pairlistsAgeLimit, pairlistsAreValid, patch, PATCH_PRIORITY, Patch::patchID, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, HomePatch::positionsReady_SOA(), Flags::savePairlists, Flags::sequence, simParams, Flags::step, suspend(), Flags::usePairlists, and SimParameters::usePairlists.

Referenced by integrate_SOA().

 {
   if ( migration ) pairlistsAreValid = 0;
 #if (defined(NAMD_CUDA) || defined(NAMD_HIP)) || defined(NAMD_MIC)
   if ( pairlistsAreValid &&
        ( patch->flags.doFullElectrostatics || ! simParams->fullElectFrequency )
                          && ( pairlistsAge > pairlistsAgeLimit ) ) {
     pairlistsAreValid = 0;
   }
 #else
   if ( pairlistsAreValid && ( pairlistsAge > pairlistsAgeLimit ) ) {
     pairlistsAreValid = 0;
   }
 #endif
   if ( ! simParams->usePairlists ) pairlists = 0;
   patch->flags.usePairlists = pairlists || pairlistsAreValid;
   patch->flags.savePairlists = pairlists && ! pairlistsAreValid;
 
 #if defined(NTESTPID)
   if (1 && patch->patchID == NTESTPID) {
     int step = patch->flags.step;
     int numAtoms = patch->numAtoms;
     double *xyzq = new double[4*numAtoms];
     double *x = patch->patchDataSOA.pos_x;
     double *y = patch->patchDataSOA.pos_y;
     double *z = patch->patchDataSOA.pos_z;
     float *q = patch->patchDataSOA.charge;
     for (int i=0;  i < numAtoms;  i++) {
       xyzq[4*i  ] = x[i];
       xyzq[4*i+1] = y[i];
       xyzq[4*i+2] = z[i];
       xyzq[4*i+3] = q[i];
     }
     char fname[128], remark[128];
     sprintf(fname, "xyzq_soa_pid%d_step%d.bin", NTESTPID, step);
     sprintf(remark, "SOA xyzq, patch %d, step %d", NTESTPID, step);
     TestArray_write<double>(fname, remark, xyzq, 4*numAtoms);
     delete[] xyzq;
   }
 #endif
   // Zero all SOA global forces before computing force 
   patch->zero_global_forces_SOA();
   patch->positionsReady_SOA(migration);  // updates flags.sequence
 
   int seq = patch->flags.sequence;
   int basePriority = ( (seq & 0xffff) << 15 )
                      + PATCH_PRIORITY(patch->getPatchID());
 
   // XXX missing GBIS
   priority = basePriority + COMPUTE_HOME_PRIORITY;
   //char prbuf[32];
   //sprintf(prbuf, "%s: %d", NamdProfileEventStr[NamdProfileEvent::SEQ_SUSPEND], patch->getPatchID());
   //NAMD_EVENT_START_EX(1, NamdProfileEvent::SEQ_SUSPEND, prbuf);
   suspend(); // until all deposit boxes close
   //NAMD_EVENT_STOP(1, NamdProfileEvent::SEQ_SUSPEND);
 
 #ifdef NODEGROUP_FORCE_REGISTER
   if(!simParams->CUDASOAintegrate || migration){
      patch->copy_forces_to_SOA();
   }
 #else
   patch->copy_forces_to_SOA();
 #endif
 
 #if defined(NTESTPID)
   if (1 && patch->patchID == NTESTPID) {
     int step = patch->flags.step;
     int numAtoms = patch->numAtoms;
     char fname[128];
     char remark[128];
     double *fxyz = new double[3*numAtoms];
     double *fx = patch->patchDataSOA.f_normal_x;
     double *fy = patch->patchDataSOA.f_normal_y;
     double *fz = patch->patchDataSOA.f_normal_z;
     for (int i=0;  i < numAtoms;  i++) {
       fxyz[3*i  ] = fx[i];
       fxyz[3*i+1] = fy[i];
       fxyz[3*i+2] = fz[i];
     }
     sprintf(fname, "fxyz_normal_soa_pid%d_step%d.bin", NTESTPID, step);
     sprintf(remark, "SOA fxyz normal, patch %d, step %d", NTESTPID, step);
     TestArray_write<double>(fname, remark, fxyz, 3*numAtoms);
     fx = patch->patchDataSOA.f_nbond_x;
     fy = patch->patchDataSOA.f_nbond_y;
     fz = patch->patchDataSOA.f_nbond_z;
     for (int i=0;  i < numAtoms;  i++) {
       fxyz[3*i  ] = fx[i];
       fxyz[3*i+1] = fy[i];
       fxyz[3*i+2] = fz[i];
     }
     sprintf(fname, "fxyz_nbond_soa_pid%d_step%d.bin", NTESTPID, step);
     sprintf(remark, "SOA fxyz nonbonded, patch %d, step %d", NTESTPID, step);
     TestArray_write<double>(fname, remark, fxyz, 3*numAtoms);
     fx = patch->patchDataSOA.f_slow_x;
     fy = patch->patchDataSOA.f_slow_y;
     fz = patch->patchDataSOA.f_slow_z;
     for (int i=0;  i < numAtoms;  i++) {
       fxyz[3*i  ] = fx[i];
       fxyz[3*i+1] = fy[i];
       fxyz[3*i+2] = fz[i];
     }
     sprintf(fname, "fxyz_slow_soa_pid%d_step%d.bin", NTESTPID, step);
     sprintf(remark, "SOA fxyz slow, patch %d, step %d", NTESTPID, step);
     TestArray_write<double>(fname, remark, fxyz, 3*numAtoms);
     delete[] fxyz;
   }
 #endif
 
 #if 0
   if (1 && patch->patchID == 0) {
     int numAtoms = patch->numAtoms;
     double *fxyz = new double[3*numAtoms];
     double *fx, *fy, *fz;
     char fname[64], remark[128];
     int step = patch->flags.step;
 
     fx = patch->patchDataSOA.f_slow_x;
     fy = patch->patchDataSOA.f_slow_y;
     fz = patch->patchDataSOA.f_slow_z;
     for (int i=0;  i < numAtoms;  i++) {
       fxyz[3*i  ] = fx[i];
       fxyz[3*i+1] = fy[i];
       fxyz[3*i+2] = fz[i];
     }
     sprintf(fname, "fslow_soa_%d.bin", step);
     sprintf(remark, "SOA slow forces, step %d\n", step);
     TestArray_write<double>(fname, remark, fxyz, 3*numAtoms);
 
     fx = patch->patchDataSOA.f_nbond_x;
     fy = patch->patchDataSOA.f_nbond_y;
     fz = patch->patchDataSOA.f_nbond_z;
     for (int i=0;  i < numAtoms;  i++) {
       fxyz[3*i  ] = fx[i];
       fxyz[3*i+1] = fy[i];
       fxyz[3*i+2] = fz[i];
     }
     sprintf(fname, "fnbond_soa_%d.bin", step);
     sprintf(remark, "SOA nonbonded forces, step %d\n", step);
     TestArray_write<double>(fname, remark, fxyz, 3*numAtoms);
 
     fx = patch->patchDataSOA.f_normal_x;
     fy = patch->patchDataSOA.f_normal_y;
     fz = patch->patchDataSOA.f_normal_z;
     for (int i=0;  i < numAtoms;  i++) {
       fxyz[3*i  ] = fx[i];
       fxyz[3*i+1] = fy[i];
       fxyz[3*i+2] = fz[i];
     }
     sprintf(fname, "fnormal_soa_%d.bin", step);
     sprintf(remark, "SOA normal forces, step %d\n", step);
     TestArray_write<double>(fname, remark, fxyz, 3*numAtoms);
 
     delete[] fxyz;
   }
 #endif
 
 #if 0
      //Will print forces here after runComputeObjects
      if(nstep == 1){
         fprintf(stderr, "CPU force arrays for alanin\n" );
         for(int i = 0; i < patch->patchDataSOA.numAtoms; i++){
           fprintf(stderr, "f[%i] = %lf %lf %lf | %lf %lf %lf | %lf %lf %lf\n", i,
              patch->patchDataSOA.f_normal_x[i], patch->patchDataSOA.f_normal_y[i], patch->patchDataSOA.f_normal_z[i],
              patch->patchDataSOA.f_nbond_x[i],  patch->patchDataSOA.f_nbond_y[i],  patch->patchDataSOA.f_nbond_z[i],
              patch->patchDataSOA.f_slow_x[i],   patch->patchDataSOA.f_slow_y[i],   patch->patchDataSOA.f_slow_z[i]);
         }
      }
 #endif
 
   if ( patch->flags.savePairlists && patch->flags.doNonbonded ) {
     pairlistsAreValid = 1;
     pairlistsAge = 0;
   }
   // For multigrator, do not age pairlist during pressure step
   // NOTE: for non-multigrator pressureStep = 0 always
   if ( pairlistsAreValid /* && !pressureStep */ ) ++pairlistsAge;
 
   // XXX missing lonepairs
   // XXX missing Molly
   // XXX missing Lowe-Andersen
 }

◆ runComputeObjectsCUDA()

void Sequencer::runComputeObjectsCUDA	(	int	doMigration,
		int	doGlobal,
		int	pairlists,
		int	nstep,
		int	startup
	)

protected

◆ saveForce()

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

protected

Definition at line 5630 of file Sequencer.C.

References patch, and HomePatch::saveForce().

Referenced by integrate(), and minimize().

 {
   patch->saveForce(ftag);
 }

◆ scalePositionsVelocities()

void Sequencer::scalePositionsVelocities	(	const Tensor &	posScale,
		const Tensor &	velScale
	)

protected

Definition at line 4789 of file Sequencer.C.

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

                                                                                        {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   Position origin = patch->lattice.origin();
   if ( simParams->fixedAtomsOn ) {
     NAMD_bug("Sequencer::scalePositionsVelocities, fixed atoms not implemented");
   }
   if ( simParams->useGroupPressure ) {
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
       hgs = a[i].hydrogenGroupSize;
       Position pos_cm(0.0, 0.0, 0.0);
       Velocity vel_cm(0.0, 0.0, 0.0);
       BigReal m_cm = 0.0;
       for (int j=0;j < hgs;++j) {
         m_cm += a[i+j].mass;
         pos_cm += a[i+j].mass*a[i+j].position;
         vel_cm += a[i+j].mass*a[i+j].velocity;
       }
       pos_cm /= m_cm;
       vel_cm /= m_cm;
       pos_cm -= origin;
       Position dpos = posScale*pos_cm;
       Velocity dvel = velScale*vel_cm;
       for (int j=0;j < hgs;++j) {
         a[i+j].position += dpos;
         a[i+j].velocity += dvel;
       }
     }
   } else {
     for ( int i = 0; i < numAtoms; i++) {
       a[i].position += posScale*(a[i].position-origin);
       a[i].velocity = velScale*a[i].velocity;
     }
   }
 }

◆ scaleVelocities()

void Sequencer::scaleVelocities ( const BigReal velScale )

protected

Definition at line 4961 of file Sequencer.C.

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

Referenced by multigratorTemperature().

                                                       {
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
   for ( int i = 0; i < numAtoms; i++) {
     a[i].velocity *= velScale;
   }
 }

◆ stochRescaleVelocities()

void Sequencer::stochRescaleVelocities ( 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

step	The current timestep

Definition at line 5614 of file Sequencer.C.

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

Referenced by integrate().

 {
   ++stochRescale_count;
   if ( stochRescale_count == simParams->stochRescaleFreq ) {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     // Blocking receive for the temperature coupling coefficient.
     BigReal velrescaling = broadcast->stochRescaleCoefficient.get(step);
     DebugM(4, "stochastically rescaling velocities at step " << step << " by " << velrescaling << "\n");
     for ( int i = 0; i < numAtoms; ++i ) {
       a[i].velocity *= velrescaling;
     }
     stochRescale_count = 0;
   }
 }

◆ stochRescaleVelocities_SOA()

void Sequencer::stochRescaleVelocities_SOA ( int step )

protected

Rescale velocities with the scale factor sent from the Controller.

Parameters

step	The current timestep

Definition at line 3882 of file Sequencer.C.

References broadcast, DebugM, SimpleBroadcastObject< T >::get(), PatchDataSOA::numAtoms, patch, simParams, stochRescale_count, ControllerBroadcasts::stochRescaleCoefficient, SimParameters::stochRescaleFreq, PatchDataSOA::vel_x, PatchDataSOA::vel_y, and PatchDataSOA::vel_z.

Referenced by integrate_SOA().

 {
   ++stochRescale_count;
   if ( stochRescale_count == simParams->stochRescaleFreq ) {
     double * __restrict vel_x = patch->patchDataSOA.vel_x;
     double * __restrict vel_y = patch->patchDataSOA.vel_y;
     double * __restrict vel_z = patch->patchDataSOA.vel_z;
     int numAtoms = patch->patchDataSOA.numAtoms;
     // Blocking receive for the temperature coupling coefficient.
     BigReal velrescaling = broadcast->stochRescaleCoefficient.get(step);
     DebugM(4, "stochastically rescaling velocities at step " << step << " by " << velrescaling << "\n");
     for ( int i = 0; i < numAtoms; ++i ) {
       vel_x[i] *= velrescaling;
       vel_y[i] *= velrescaling;
       vel_z[i] *= velrescaling;
     }
     stochRescale_count = 0;
   }
 }

◆ submitCollections()

void Sequencer::submitCollections	(	int	step,
		int	zeroVel = `0`
	)

protected

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

 {
   //
   // DJH: Copy updates of SOA back into AOS.
   // Do we need to update everything or is it safe to just update
   // positions and velocities separately, as needed?
   //
   //patch->copy_updates_to_AOS();
 
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::SUBMIT_COLLECTIONS);
   int prec;
   int dcdSelectionIndex;
   std::tie (prec, dcdSelectionIndex) = Output::coordinateNeeded(step);
   if ( prec ) {
 #ifndef MEM_OPT_VERSION
     collection->submitPositions(step,patch->atom,patch->lattice,prec,dcdSelectionIndex);
 #else
     collection->submitPositions(step,patch->atom,patch->lattice,prec);
 #endif
   }
   if ( Output::velocityNeeded(step) ) {
     collection->submitVelocities(step,zeroVel,patch->atom);
   }
   if ( Output::forceNeeded(step) ) {
     int maxForceUsed = patch->flags.maxForceUsed;
     if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
     collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
   }
 }

◆ submitCollections_SOA()

void Sequencer::submitCollections_SOA	(	int	step,
		int	zeroVel = `0`
	)

protected

Definition at line 3207 of file Sequencer.C.

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

Referenced by integrate_SOA().

 {
   //
   // Copy updates of SOA back into AOS for collections.
   //
   // XXX Could update positions and velocities separately.
   //
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::SUBMIT_COLLECTIONS_SOA);
   //
   // XXX Poor implementation here!
   // The selector functions called below in Output.C are
   // doing several tests and in an average use case calculating
   // at least two mod functions.
   //
   // However, most steps are NOT output steps!
   //
   int is_pos_needed;
   int dcdIndex;
   std::tie(is_pos_needed, dcdIndex)= Output::coordinateNeeded(step);
   int is_vel_needed = Output::velocityNeeded(step);
   int is_f_needed = Output::forceNeeded(step);
   if (!simParams->useDeviceMigration) {  // This is already done for GPU migration
     if ( is_pos_needed || is_vel_needed ) {
       patch->copy_updates_to_AOS();
     }
     if (is_f_needed) {
       patch->copy_forces_to_AOS();
     }
   }
   if ( is_pos_needed ) {
 #ifndef MEM_OPT_VERSION
     collection->submitPositions(step,patch->atom,patch->lattice,is_pos_needed,dcdIndex);
 #else
     collection->submitPositions(step,patch->atom,patch->lattice,is_pos_needed);
 #endif
   }
   if ( is_vel_needed ) {
     collection->submitVelocities(step,zeroVel,patch->atom);
   }
   if ( is_f_needed ) {
     int maxForceUsed = patch->flags.maxForceUsed;
     if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
     collection->submitForces(step,patch->atom,maxForceUsed,patch->f);
   }
 }

◆ submitHalfstep()

void Sequencer::submitHalfstep ( int step )

protected

Definition at line 5820 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ResizeArray< Elem >::begin(), Lattice::c(), doKineticEnergy, Flags::doVirial, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), Patch::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, and Vector::z.

Referenced by integrate().

 {
   NAMD_EVENT_RANGE_2(patch->flags.event_on, NamdProfileEvent::SUBMIT_HALFSTEP);
 
   // velocity-dependent quantities *** ONLY ***
   // positions are not at half-step when called
   FullAtom *a = patch->atom.begin();
   int numAtoms = patch->numAtoms;
 
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter (0);
 #endif
 
   // For non-Multigrator doKineticEnergy = 1 always
   Tensor momentumSqrSum;
   if (doKineticEnergy || patch->flags.doVirial)
   {
     BigReal kineticEnergy = 0;
     Tensor virial;
     if ( simParams->pairInteractionOn ) {
       if ( simParams->pairInteractionSelf ) {
         for ( int i = 0; i < numAtoms; ++i ) {
           if ( a[i].partition != 1 ) continue;
           kineticEnergy += a[i].mass * a[i].velocity.length2();
           virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
         }
       }
     } else {
       for ( int i = 0; i < numAtoms; ++i ) {
         if (a[i].mass < 0.01) continue;
         kineticEnergy += a[i].mass * a[i].velocity.length2();
         virial.outerAdd(a[i].mass, a[i].velocity, a[i].velocity);
       }
     }
 
     if (simParams->multigratorOn && !simParams->useGroupPressure) {
       momentumSqrSum = virial;
     }
     kineticEnergy *= 0.5 * 0.5;
     reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
     virial *= 0.5;
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
 #ifdef ALTVIRIAL
     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,virial);
 #endif
   }
 
   if (pressureProfileReduction) {
     int nslabs = simParams->pressureProfileSlabs;
     const Lattice &lattice = patch->lattice;
     BigReal idz = nslabs/lattice.c().z;
     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
     int useGroupPressure = simParams->useGroupPressure;
 
     // Compute kinetic energy partition, possibly subtracting off
     // internal kinetic energy if group pressure is enabled.
     // Since the regular pressure is 1/2 mvv and the internal kinetic
     // term that is subtracted off for the group pressure is
     // 1/2 mv (v-v_cm), the group pressure kinetic contribution is
     // 1/2 m * v * v_cm.  The factor of 1/2 is because submitHalfstep
     // gets called twice per timestep.
     int hgs;
     for (int i=0; i<numAtoms; i += hgs) {
       int j, ppoffset;
       hgs = a[i].hydrogenGroupSize;
       int partition = a[i].partition;
 
       BigReal m_cm = 0;
       Velocity v_cm(0,0,0);
       for (j=i; j< i+hgs; ++j) {
         m_cm += a[j].mass;
         v_cm += a[j].mass * a[j].velocity;
       }
       v_cm /= m_cm;
       for (j=i; j < i+hgs; ++j) {
         BigReal mass = a[j].mass;
         if (! (useGroupPressure && j != i)) {
           BigReal z = a[j].position.z;
           int slab = (int)floor((z-zmin)*idz);
           if (slab < 0) slab += nslabs;
           else if (slab >= nslabs) slab -= nslabs;
           ppoffset = 3*(slab + partition*nslabs);
         }
         BigReal wxx, wyy, wzz;
         if (useGroupPressure) {
           wxx = 0.5*mass * a[j].velocity.x * v_cm.x;
           wyy = 0.5*mass * a[j].velocity.y * v_cm.y;
           wzz = 0.5*mass * a[j].velocity.z * v_cm.z;
         } else {
           wxx = 0.5*mass * a[j].velocity.x * a[j].velocity.x;
           wyy = 0.5*mass * a[j].velocity.y * a[j].velocity.y;
           wzz = 0.5*mass * a[j].velocity.z * a[j].velocity.z;
         }
         pressureProfileReduction->item(ppoffset  ) += wxx;
         pressureProfileReduction->item(ppoffset+1) += wyy;
         pressureProfileReduction->item(ppoffset+2) += wzz;
       }
     }
   }
 
   // For non-Multigrator doKineticEnergy = 1 always
   if (doKineticEnergy || patch->flags.doVirial)
   {
     BigReal intKineticEnergy = 0;
     Tensor intVirialNormal;
 
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
 
 #if CMK_BLUEGENEL
       CmiNetworkProgress ();
 #endif
 
       hgs = a[i].hydrogenGroupSize;
       int j;
       BigReal m_cm = 0;
       Velocity v_cm(0,0,0);
       for ( j = i; j < (i+hgs); ++j ) {
         m_cm += a[j].mass;
         v_cm += a[j].mass * a[j].velocity;
       }
       if (simParams->multigratorOn && simParams->useGroupPressure) {
         momentumSqrSum.outerAdd(1.0/m_cm, v_cm, v_cm);
       }
       v_cm /= m_cm;
       if ( simParams->pairInteractionOn ) {
         if ( simParams->pairInteractionSelf ) {
           for ( j = i; j < (i+hgs); ++j ) {
             if ( a[j].partition != 1 ) continue;
             BigReal mass = a[j].mass;
             Vector v = a[j].velocity;
             Vector dv = v - v_cm;
             intKineticEnergy += mass * (v * dv);
             intVirialNormal.outerAdd (mass, v, dv);
           }
         }
       } else {
         for ( j = i; j < (i+hgs); ++j ) {
           BigReal mass = a[j].mass;
           Vector v = a[j].velocity;
           Vector dv = v - v_cm;
           intKineticEnergy += mass * (v * dv);
           intVirialNormal.outerAdd(mass, v, dv);
         }
       }
     }
     intKineticEnergy *= 0.5 * 0.5;
     reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
     intVirialNormal *= 0.5;
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
     if ( simParams->multigratorOn) {
       momentumSqrSum *= 0.5;
       ADD_TENSOR_OBJECT(reduction,REDUCTION_MOMENTUM_SQUARED,momentumSqrSum);
     }
   }
 
 }

◆ submitHalfstep_SOA()

void Sequencer::submitHalfstep_SOA ( )

protected

Definition at line 2849 of file Sequencer.C.

References ADD_TENSOR_OBJECT, SimParameters::CUDASOAintegrate, Patch::flags, PatchDataSOA::hydrogenGroupSize, SubmitReduction::item(), PatchDataSOA::mass, NAMD_EVENT_RANGE_2, PatchDataSOA::numAtoms, patch, reduction, REDUCTION_HALFSTEP_KINETIC_ENERGY, REDUCTION_INT_HALFSTEP_KINETIC_ENERGY, simParams, PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, Tensor::xx, Tensor::xy, Tensor::xz, Tensor::yx, Tensor::yy, Tensor::yz, Tensor::zx, Tensor::zy, and Tensor::zz.

Referenced by integrate_SOA().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::SUBMIT_HALFSTEP_SOA);
 #ifdef SOA_SIMPLIFY_PARAMS
   const int    * __restrict hydrogenGroupSize = patch->patchDataSOA.hydrogenGroupSize;
   const float  * __restrict mass = patch->patchDataSOA.mass;
   const double * __restrict vel_x = patch->patchDataSOA.vel_x;
   const double * __restrict vel_y = patch->patchDataSOA.vel_y;
   const double * __restrict vel_z = patch->patchDataSOA.vel_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
   if ( 1 /* doKineticEnergy || patch->flags.doVirial */ ) {
     BigReal kineticEnergy = 0;
     Tensor virial;
     for (int i=0;  i < numAtoms;  i++) {
       // scalar kineticEnergy += mass[i] * vel[i]^2
       kineticEnergy += mass[i] *
         (vel_x[i]*vel_x[i] + vel_y[i]*vel_y[i] + vel_z[i]*vel_z[i]);
       // tensor virial += mass[i] * outer_product(vel[i], vel[i])
       virial.xx += mass[i] * vel_x[i] * vel_x[i];
       virial.xy += mass[i] * vel_x[i] * vel_y[i];
       virial.xz += mass[i] * vel_x[i] * vel_z[i];
       virial.yx += mass[i] * vel_y[i] * vel_x[i];
       virial.yy += mass[i] * vel_y[i] * vel_y[i];
       virial.yz += mass[i] * vel_y[i] * vel_z[i];
       virial.zx += mass[i] * vel_z[i] * vel_x[i];
       virial.zy += mass[i] * vel_z[i] * vel_y[i];
       virial.zz += mass[i] * vel_z[i] * vel_z[i];
     }
     kineticEnergy *= 0.5 * 0.5;
     virial *= 0.5;
 
 #ifdef NODEGROUP_FORCE_REGISTER
     if(simParams->CUDASOAintegrate){
       CUDASequencer->patchData->reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
       ADD_TENSOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_VIRIAL_NORMAL,virial);
     }
     else{
 #endif
       reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
 #ifdef NODEGROUP_FORCE_REGISTER
     }
 #endif
   }
 
   if ( 1 /* doKineticEnergy || patch->flags.doVirial */ ) {
     BigReal intKineticEnergy = 0;
     Tensor intVirialNormal;
     int hgs;
     for (int i=0;  i < numAtoms;  i += hgs) {
       // find velocity of center-of-mass of hydrogen group
       // calculate mass-weighted velocity
       hgs = hydrogenGroupSize[i];
       BigReal m_cm = 0;
       BigReal v_cm_x = 0;
       BigReal v_cm_y = 0;
       BigReal v_cm_z = 0;
       for (int j = i;  j < (i+hgs);  j++) {
         m_cm += mass[j];
         v_cm_x += mass[j] * vel_x[j];
         v_cm_y += mass[j] * vel_y[j];
         v_cm_z += mass[j] * vel_z[j];
       }
       BigReal recip_m_cm = 1.0 / m_cm;
       v_cm_x *= recip_m_cm;
       v_cm_y *= recip_m_cm;
       v_cm_z *= recip_m_cm;
       // sum virial contributions wrt vel center-of-mass
       for (int j = i;  j < (i+hgs);  j++) {
         BigReal dv_x = vel_x[j] - v_cm_x;
         BigReal dv_y = vel_y[j] - v_cm_y;
         BigReal dv_z = vel_z[j] - v_cm_z;
         // scalar intKineticEnergy += mass[j] * dot_product(vel[j], dv)
         intKineticEnergy += mass[j] *
           (vel_x[j] * dv_x + vel_y[j] * dv_y + vel_z[j] * dv_z);
         // tensor intVirialNormal += mass[j] * outer_product(vel[j], dv)
         intVirialNormal.xx += mass[j] * vel_x[j] * dv_x;
         intVirialNormal.xy += mass[j] * vel_x[j] * dv_y;
         intVirialNormal.xz += mass[j] * vel_x[j] * dv_z;
         intVirialNormal.yx += mass[j] * vel_y[j] * dv_x;
         intVirialNormal.yy += mass[j] * vel_y[j] * dv_y;
         intVirialNormal.yz += mass[j] * vel_y[j] * dv_z;
         intVirialNormal.zx += mass[j] * vel_z[j] * dv_x;
         intVirialNormal.zy += mass[j] * vel_z[j] * dv_y;
         intVirialNormal.zz += mass[j] * vel_z[j] * dv_z;
       }
     }
     intKineticEnergy *= 0.5 * 0.5;
     intVirialNormal *= 0.5;
 #ifdef NODEGROUP_FORCE_REGISTER
     if(simParams->CUDASOAintegrate){
       CUDASequencer->patchData->reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY) += intKineticEnergy;
       ADD_TENSOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_INT_VIRIAL_NORMAL, intVirialNormal);
     }
     else{
 #endif
       reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY)
         += intKineticEnergy;
       ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL,
         intVirialNormal);
 #ifdef NODEGROUP_FORCE_REGISTER
     }
 #endif
   }
 }

◆ submitMinimizeReductions()

void Sequencer::submitMinimizeReductions	(	int	step,
		BigReal	fmax2
	)

protected

Definition at line 6256 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, CompAtomExt::atomFixed, ResizeArray< Elem >::begin(), calcFixVirial(), SimParameters::drudeBondLen, SimParameters::drudeHardWallOn, Patch::f, SimParameters::fixedAtomsOn, Patch::flags, CompAtom::hydrogenGroupSize, SubmitReduction::item(), 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(), FullAtom::velocity, Vector::y, and Vector::z.

Referenced by minimize().

 {
   FullAtom *a = patch->atom.begin();
   Force *f1 = patch->f[Results::normal].begin();
   Force *f2 = patch->f[Results::nbond].begin();
   Force *f3 = patch->f[Results::slow].begin();
   const bool fixedAtomsOn = simParams->fixedAtomsOn;
   const bool drudeHardWallOn = simParams->drudeHardWallOn;
   const double drudeBondLen = simParams->drudeBondLen;
   const double drudeBondLen2 = drudeBondLen * drudeBondLen;
   const double drudeStep = 0.1/(TIMEFACTOR*TIMEFACTOR);
   const double drudeMove = 0.01;
   const double drudeStep2 = drudeStep * drudeStep;
   const double drudeMove2 = drudeMove * drudeMove;
   int numAtoms = patch->numAtoms;
 
   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
 
   for ( int i = 0; i < numAtoms; ++i ) {
 #if 0
     printf("ap[%2d]=  %f  %f  %f\n", i, a[i].position.x, a[i].position.y, a[i].position.z);
     printf("f1[%2d]=  %f  %f  %f\n", i, f1[i].x, f1[i].y, f1[i].z);
     printf("f2[%2d]=  %f  %f  %f\n", i, f2[i].x, f2[i].y, f2[i].z);
     //printf("f3[%2d]=  %f  %f  %f\n", i, f3[i].x, f3[i].y, f3[i].z);
 #endif
     f1[i] += f2[i] + f3[i];  // add all forces
     if ( drudeHardWallOn && i && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) { // drude particle
       if ( ! fixedAtomsOn || ! a[i].atomFixed ) {
         if ( drudeStep2 * f1[i].length2() > drudeMove2 ) {
           a[i].position += drudeMove * f1[i].unit();
         } else {
           a[i].position += drudeStep * f1[i];
         }
         if ( (a[i].position - a[i-1].position).length2() > drudeBondLen2 ) {
           a[i].position = a[i-1].position + drudeBondLen * (a[i].position - a[i-1].position).unit();
         }
       }
       Vector netf = f1[i-1] + f1[i];
       if ( fixedAtomsOn && a[i-1].atomFixed ) netf = 0;
       f1[i-1] = netf;
       f1[i] = 0.;
     }
     if ( fixedAtomsOn && a[i].atomFixed ) f1[i] = 0;
   }
 
   f2 = f3 = 0;  // included in f1
 
   BigReal maxv2 = 0.;
 
   for ( int i = 0; i < numAtoms; ++i ) {
     BigReal v2 = a[i].velocity.length2();
     if ( v2 > 0. ) {
       if ( v2 > maxv2 ) maxv2 = v2;
     } else {
       v2 = f1[i].length2();
       if ( v2 > maxv2 ) maxv2 = v2;
     }
   }
 
   if ( fmax2 > 10. * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR * TIMEFACTOR )
   { Tensor virial; patch->minimize_rattle2( 0.1 * TIMEFACTOR / sqrt(maxv2), &virial, true /* forces */); }
 
   BigReal fdotf = 0;
   BigReal fdotv = 0;
   BigReal vdotv = 0;
   int numHuge = 0;
   for ( int i = 0; i < numAtoms; ++i ) {
     if ( simParams->fixedAtomsOn && a[i].atomFixed ) continue;
     if ( drudeHardWallOn && (a[i].mass > 0.05) && ((a[i].mass < 1.0)) ) continue; // drude particle
     Force f = f1[i];
     BigReal ff = f * f;
     if ( ff > fmax2 ) {
       if (simParams->printBadContacts) {
         CkPrintf("STEP(%i) MIN_HUGE[%i] f=%e kcal/mol/A\n",patch->flags.sequence,patch->pExt[i].id,ff);
       }
       ++numHuge;
       // pad scaling so minimizeMoveDownhill() doesn't miss them
       BigReal fmult = 1.01 * sqrt(fmax2/ff);
       f *= fmult;  ff = f * f;
       f1[i] *= fmult;
     }
     fdotf += ff;
     fdotv += f * a[i].velocity;
     vdotv += a[i].velocity * a[i].velocity;
   }
 
 #if 0
   printf("fdotf = %f\n", fdotf);
   printf("fdotv = %f\n", fdotv);
   printf("vdotv = %f\n", vdotv);
 #endif
   reduction->item(REDUCTION_MIN_F_DOT_F) += fdotf;
   reduction->item(REDUCTION_MIN_F_DOT_V) += fdotv;
   reduction->item(REDUCTION_MIN_V_DOT_V) += vdotv;
   reduction->item(REDUCTION_MIN_HUGE_COUNT) += numHuge;
 
   {
     Tensor intVirialNormal;
     Tensor intVirialNbond;
     Tensor intVirialSlow;
 
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
       hgs = a[i].hydrogenGroupSize;
       int j;
       BigReal m_cm = 0;
       Position x_cm(0,0,0);
       for ( j = i; j < (i+hgs); ++j ) {
         m_cm += a[j].mass;
         x_cm += a[j].mass * a[j].position;
       }
       x_cm /= m_cm;
       for ( j = i; j < (i+hgs); ++j ) {
         BigReal mass = a[j].mass;
         // net force treated as zero for fixed atoms
         if ( simParams->fixedAtomsOn && a[j].atomFixed ) continue;
         Vector dx = a[j].position - x_cm;
         intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
         intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
         intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
       }
     }
 
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
   }
 
   if ( simParams->fixedAtomsOn ) {
     Tensor fixVirialNormal;
     Tensor fixVirialNbond;
     Tensor fixVirialSlow;
     Vector fixForceNormal = 0;
     Vector fixForceNbond = 0;
     Vector fixForceSlow = 0;
 
     calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
 
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
     ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
   }
 
   reduction->submit();
 }

◆ submitMomentum()

void Sequencer::submitMomentum ( int step )

protected

Definition at line 4735 of file Sequencer.C.

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

Referenced by integrate().

                                        {
 
   FullAtom *a = patch->atom.begin();
   const int numAtoms = patch->numAtoms;
 
   Vector momentum = 0;
   BigReal mass = 0;
 if ( simParams->zeroMomentumAlt ) {
   for ( int i = 0; i < numAtoms; ++i ) {
     momentum += a[i].mass * a[i].velocity;
     mass += 1.;
   }
 } else {
   for ( int i = 0; i < numAtoms; ++i ) {
     momentum += a[i].mass * a[i].velocity;
     mass += a[i].mass;
   }
 }
 
   ADD_VECTOR_OBJECT(reduction,REDUCTION_HALFSTEP_MOMENTUM,momentum);
   reduction->item(REDUCTION_MOMENTUM_MASS) += mass;
 }

◆ submitReductions()

void Sequencer::submitReductions ( int step )

protected

Definition at line 5998 of file Sequencer.C.

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

 {
 #ifndef UPPER_BOUND
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::SUBMIT_REDUCTIONS);
   FullAtom *a = patch->atom.begin();
 #endif
   int numAtoms = patch->numAtoms;
 
 #if CMK_BLUEGENEL
   CmiNetworkProgressAfter(0);
 #endif
 
   reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
   reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
 
 #ifndef UPPER_BOUND
   // For non-Multigrator doKineticEnergy = 1 always
   if (doKineticEnergy || doMomenta || patch->flags.doVirial)
   {
     BigReal kineticEnergy = 0;
     Vector momentum = 0;
     Vector angularMomentum = 0;
     Vector o = patch->lattice.origin();
     int i;
     if ( simParams->pairInteractionOn ) {
       if ( simParams->pairInteractionSelf ) {
         for (i = 0; i < numAtoms; ++i ) {
           if ( a[i].partition != 1 ) continue;
           kineticEnergy += a[i].mass * a[i].velocity.length2();
           momentum += a[i].mass * a[i].velocity;
           angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
         }
       }
     } else {
       for (i = 0; i < numAtoms; ++i ) {
         kineticEnergy += a[i].mass * a[i].velocity.length2();
         momentum += a[i].mass * a[i].velocity;
         angularMomentum += cross(a[i].mass,a[i].position-o,a[i].velocity);
       }
       if (simParams->drudeOn) {
         BigReal drudeComKE = 0.;
         BigReal drudeBondKE = 0.;
 
         for (i = 0;  i < numAtoms;  i++) {
           if (i < numAtoms-1 &&
               a[i+1].mass < 1.0 && a[i+1].mass > 0.05) {
             // i+1 is a Drude particle with parent i
 
             // convert from Cartesian coordinates to (COM,bond) coordinates
             BigReal m_com = (a[i].mass + a[i+1].mass);  // mass of COM
             BigReal m = a[i+1].mass / m_com;  // mass ratio
             BigReal m_bond = a[i+1].mass * (1. - m);  // mass of bond
             Vector v_bond = a[i+1].velocity - a[i].velocity;  // vel of bond
             Vector v_com = a[i].velocity + m * v_bond;  // vel of COM
 
             drudeComKE += m_com * v_com.length2();
             drudeBondKE += m_bond * v_bond.length2();
 
             i++;  // +1 from loop, we've updated both particles
           }
           else {
             drudeComKE += a[i].mass * a[i].velocity.length2();
           }
         } // end for
 
         drudeComKE *= 0.5;
         drudeBondKE *= 0.5;
         reduction->item(REDUCTION_DRUDECOM_CENTERED_KINETIC_ENERGY)
           += drudeComKE;
         reduction->item(REDUCTION_DRUDEBOND_CENTERED_KINETIC_ENERGY)
           += drudeBondKE;
       } // end drudeOn
 
     } // end else
 
     kineticEnergy *= 0.5;
     reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
     ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
     ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
   }
 
 #ifdef ALTVIRIAL
   // THIS IS NOT CORRECTED FOR PAIR INTERACTIONS
   {
     Tensor altVirial;
     for ( int i = 0; i < numAtoms; ++i ) {
       altVirial.outerAdd(1.0, patch->f[Results::normal][i], a[i].position);
     }
     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NORMAL,altVirial);
   }
   {
     Tensor altVirial;
     for ( int i = 0; i < numAtoms; ++i ) {
       altVirial.outerAdd(1.0, patch->f[Results::nbond][i], a[i].position);
     }
     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_NBOND,altVirial);
   }
   {
     Tensor altVirial;
     for ( int i = 0; i < numAtoms; ++i ) {
       altVirial.outerAdd(1.0, patch->f[Results::slow][i], a[i].position);
     }
     ADD_TENSOR_OBJECT(reduction,REDUCTION_ALT_VIRIAL_SLOW,altVirial);
   }
 #endif
 
   // For non-Multigrator doKineticEnergy = 1 always
   if (doKineticEnergy || patch->flags.doVirial)
   {
     BigReal intKineticEnergy = 0;
     Tensor intVirialNormal;
     Tensor intVirialNbond;
     Tensor intVirialSlow;
 
     int hgs;
     for ( int i = 0; i < numAtoms; i += hgs ) {
 #if CMK_BLUEGENEL
       CmiNetworkProgress();
 #endif
       hgs = a[i].hydrogenGroupSize;
       int j;
       BigReal m_cm = 0;
       Position x_cm(0,0,0);
       Velocity v_cm(0,0,0);
       for ( j = i; j < (i+hgs); ++j ) {
         m_cm += a[j].mass;
         x_cm += a[j].mass * a[j].position;
         v_cm += a[j].mass * a[j].velocity;
       }
       x_cm /= m_cm;
       v_cm /= m_cm;
       int fixedAtomsOn = simParams->fixedAtomsOn;
       if ( simParams->pairInteractionOn ) {
         int pairInteractionSelf = simParams->pairInteractionSelf;
         for ( j = i; j < (i+hgs); ++j ) {
           if ( a[j].partition != 1 &&
                ( pairInteractionSelf || a[j].partition != 2 ) ) continue;
           // net force treated as zero for fixed atoms
           if ( fixedAtomsOn && a[j].atomFixed ) continue;
           BigReal mass = a[j].mass;
           Vector v = a[j].velocity;
           Vector dv = v - v_cm;
           intKineticEnergy += mass * (v * dv);
           Vector dx = a[j].position - x_cm;
           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
         }
       } else {
         for ( j = i; j < (i+hgs); ++j ) {
           // net force treated as zero for fixed atoms
           if ( fixedAtomsOn && a[j].atomFixed ) continue;
           BigReal mass = a[j].mass;
           Vector v = a[j].velocity;
           Vector dv = v - v_cm;
           intKineticEnergy += mass * (v * dv);
           Vector dx = a[j].position - x_cm;
           intVirialNormal.outerAdd(1.0, patch->f[Results::normal][j], dx);
           intVirialNbond.outerAdd(1.0, patch->f[Results::nbond][j], dx);
           intVirialSlow.outerAdd(1.0, patch->f[Results::slow][j], dx);
         }
       }
     }
 
     intKineticEnergy *= 0.5;
     reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
     ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
   }
 
   if (pressureProfileReduction && simParams->useGroupPressure) {
     // subtract off internal virial term, calculated as for intVirial.
     int nslabs = simParams->pressureProfileSlabs;
     const Lattice &lattice = patch->lattice;
     BigReal idz = nslabs/lattice.c().z;
     BigReal zmin = lattice.origin().z - 0.5*lattice.c().z;
     int useGroupPressure = simParams->useGroupPressure;
 
     int hgs;
     for (int i=0; i<numAtoms; i += hgs) {
       int j;
       hgs = a[i].hydrogenGroupSize;
       BigReal m_cm = 0;
       Position x_cm(0,0,0);
       for (j=i; j< i+hgs; ++j) {
         m_cm += a[j].mass;
         x_cm += a[j].mass * a[j].position;
       }
       x_cm /= m_cm;
 
       BigReal z = a[i].position.z;
       int slab = (int)floor((z-zmin)*idz);
       if (slab < 0) slab += nslabs;
       else if (slab >= nslabs) slab -= nslabs;
       int partition = a[i].partition;
       int ppoffset = 3*(slab + nslabs*partition);
       for (j=i; j < i+hgs; ++j) {
         BigReal mass = a[j].mass;
         Vector dx = a[j].position - x_cm;
         const Vector &fnormal = patch->f[Results::normal][j];
         const Vector &fnbond  = patch->f[Results::nbond][j];
         const Vector &fslow   = patch->f[Results::slow][j];
         BigReal wxx = (fnormal.x + fnbond.x + fslow.x) * dx.x;
         BigReal wyy = (fnormal.y + fnbond.y + fslow.y) * dx.y;
         BigReal wzz = (fnormal.z + fnbond.z + fslow.z) * dx.z;
         pressureProfileReduction->item(ppoffset  ) -= wxx;
         pressureProfileReduction->item(ppoffset+1) -= wyy;
         pressureProfileReduction->item(ppoffset+2) -= wzz;
       }
     }
   }
 
   // For non-Multigrator doVirial = 1 always
   if (patch->flags.doVirial)
   {
     if ( simParams->fixedAtomsOn ) {
       Tensor fixVirialNormal;
       Tensor fixVirialNbond;
       Tensor fixVirialSlow;
       Vector fixForceNormal = 0;
       Vector fixForceNbond = 0;
       Vector fixForceSlow = 0;
 
       calcFixVirial(fixVirialNormal, fixVirialNbond, fixVirialSlow, fixForceNormal, fixForceNbond, fixForceSlow);
 
 #if 0
       auto printTensor = [](const Tensor& t, const std::string& name){
         CkPrintf("%s", name.c_str());
         CkPrintf("\n%12.5lf %12.5lf %12.5lf\n"
                  "%12.5lf %12.5lf %12.5lf\n"
                  "%12.5lf %12.5lf %12.5lf\n",
                  t.xx, t.xy, t.xz,
                  t.yx, t.yy, t.yz,
                  t.zx, t.zy, t.zz);
       };
       printTensor(fixVirialNormal, "fixVirialNormal = ");
       printTensor(fixVirialNbond, "fixVirialNbond = ");
       printTensor(fixVirialSlow, "fixVirialSlow = ");
 #endif
 
       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,fixVirialNormal);
       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NBOND,fixVirialNbond);
       ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_SLOW,fixVirialSlow);
       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,fixForceNormal);
       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NBOND,fixForceNbond);
       ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_SLOW,fixForceSlow);
     }
   }
 #endif // UPPER_BOUND
 
   reduction->submit();
 #ifndef UPPER_BOUND
   if (pressureProfileReduction) pressureProfileReduction->submit();
 #endif
 }

◆ submitReductions_SOA()

void Sequencer::submitReductions_SOA ( )

protected

Definition at line 2969 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, SimParameters::CUDASOAintegrate, PatchDataSOA::f_nbond_x, PatchDataSOA::f_nbond_y, PatchDataSOA::f_nbond_z, PatchDataSOA::f_normal_x, PatchDataSOA::f_normal_y, PatchDataSOA::f_normal_z, PatchDataSOA::f_slow_x, PatchDataSOA::f_slow_y, PatchDataSOA::f_slow_z, Patch::flags, PatchDataSOA::hydrogenGroupSize, SubmitReduction::item(), Patch::lattice, HomePatch::marginViolations, PatchDataSOA::mass, NAMD_EVENT_RANGE_2, namd_reciprocal, PatchDataSOA::numAtoms, Lattice::origin(), patch, PatchDataSOA::pos_x, PatchDataSOA::pos_y, PatchDataSOA::pos_z, reduction, REDUCTION_ATOM_CHECKSUM, REDUCTION_CENTERED_KINETIC_ENERGY, REDUCTION_INT_CENTERED_KINETIC_ENERGY, REDUCTION_MARGIN_VIOLATIONS, simParams, SubmitReduction::submit(), PatchDataSOA::vel_x, PatchDataSOA::vel_y, PatchDataSOA::vel_z, Vector::x, 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_SOA().

       {
   NAMD_EVENT_RANGE_2(patch->flags.event_on,
       NamdProfileEvent::SUBMIT_REDUCTIONS_SOA);
 #ifdef SOA_SIMPLIFY_PARAMS
   const int    * __restrict hydrogenGroupSize = patch->patchDataSOA.hydrogenGroupSize;
   const float  * __restrict mass = patch->patchDataSOA.mass;
   const double * __restrict pos_x = patch->patchDataSOA.pos_x;
   const double * __restrict pos_y = patch->patchDataSOA.pos_y;
   const double * __restrict pos_z = patch->patchDataSOA.pos_z;
   const double * __restrict vel_x = patch->patchDataSOA.vel_x;
   const double * __restrict vel_y = patch->patchDataSOA.vel_y;
   const double * __restrict vel_z = patch->patchDataSOA.vel_z;
   const double * __restrict f_normal_x = patch->patchDataSOA.f_normal_x;
   const double * __restrict f_normal_y = patch->patchDataSOA.f_normal_y;
   const double * __restrict f_normal_z = patch->patchDataSOA.f_normal_z;
   const double * __restrict f_nbond_x = patch->patchDataSOA.f_nbond_x;
   const double * __restrict f_nbond_y = patch->patchDataSOA.f_nbond_y;
   const double * __restrict f_nbond_z = patch->patchDataSOA.f_nbond_z;
   const double * __restrict f_slow_x = patch->patchDataSOA.f_slow_x;
   const double * __restrict f_slow_y = patch->patchDataSOA.f_slow_y;
   const double * __restrict f_slow_z = patch->patchDataSOA.f_slow_z;
   int numAtoms = patch->patchDataSOA.numAtoms;
 #endif
 
 #ifdef NODEGROUP_FORCE_REGISTER
   if(simParams->CUDASOAintegrate){
     CUDASequencer->patchData->reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
     CUDASequencer->patchData->reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
   }else{
 #endif
     reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
     reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;
 #ifdef NODEGROUP_FORCE_REGISTER
   }
 #endif
 
   if ( 1 /* doKineticEnergy || doMomenta || patch->flags.doVirial */ ) {
     BigReal kineticEnergy = 0;
     BigReal momentum_x = 0;
     BigReal momentum_y = 0;
     BigReal momentum_z = 0;
     BigReal angularMomentum_x = 0;
     BigReal angularMomentum_y = 0;
     BigReal angularMomentum_z = 0;
     BigReal origin_x = patch->lattice.origin().x;
     BigReal origin_y = patch->lattice.origin().y;
     BigReal origin_z = patch->lattice.origin().z;
 
     // XXX pairInteraction
 
     for (int i=0;  i < numAtoms;  i++) {
 
       // scalar kineticEnergy += mass[i] * dot_product(vel[i], vel[i])
       kineticEnergy += mass[i] *
         (vel_x[i]*vel_x[i] + vel_y[i]*vel_y[i] + vel_z[i]*vel_z[i]);
 
       // vector momentum += mass[i] * vel[i]
       momentum_x += mass[i] * vel_x[i];
       momentum_y += mass[i] * vel_y[i];
       momentum_z += mass[i] * vel_z[i];
 
       // vector dpos = pos[i] - origin
       BigReal dpos_x = pos_x[i] - origin_x;
       BigReal dpos_y = pos_y[i] - origin_y;
       BigReal dpos_z = pos_z[i] - origin_z;
 
       // vector angularMomentum += mass[i] * cross_product(dpos, vel[i])
       angularMomentum_x += mass[i] * (dpos_y*vel_z[i] - dpos_z*vel_y[i]);
       angularMomentum_y += mass[i] * (dpos_z*vel_x[i] - dpos_x*vel_z[i]);
       angularMomentum_z += mass[i] * (dpos_x*vel_y[i] - dpos_y*vel_x[i]);
     }
 
     // XXX missing Drude
 
     kineticEnergy *= 0.5;
     Vector momentum(momentum_x, momentum_y, momentum_z);
     Vector angularMomentum(angularMomentum_x, angularMomentum_y,
         angularMomentum_z);
 
 #ifdef NODEGROUP_FORCE_REGISTER
     if(simParams->CUDASOAintegrate){
       CUDASequencer->patchData->reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
       ADD_VECTOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_MOMENTUM,momentum);
       ADD_VECTOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
     }else{
 #endif
       reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
       ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
       ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
 #ifdef NODEGROUP_FORCE_REGISTER
     }
 #endif
   }
   // For non-Multigrator doKineticEnergy = 1 always
   if ( 1 /* doKineticEnergy || patch->flags.doVirial */ ) {
     BigReal intKineticEnergy = 0;
     Tensor intVirialNormal;
     Tensor intVirialNbond;
     Tensor intVirialSlow;
 
     int hgs = 1;  // hydrogen group size
     for (int i=0;  i < numAtoms;  i += hgs) {
       hgs = hydrogenGroupSize[i];
       int j;
       BigReal m_cm = 0;
       BigReal r_cm_x = 0;
       BigReal r_cm_y = 0;
       BigReal r_cm_z = 0;
       BigReal v_cm_x = 0;
       BigReal v_cm_y = 0;
       BigReal v_cm_z = 0;
       for ( j = i; j < (i+hgs); ++j ) {
         m_cm += mass[j];
         r_cm_x += mass[j] * pos_x[j];
         r_cm_y += mass[j] * pos_y[j];
         r_cm_z += mass[j] * pos_z[j];
         v_cm_x += mass[j] * vel_x[j];
         v_cm_y += mass[j] * vel_y[j];
         v_cm_z += mass[j] * vel_z[j];
       }
       BigReal inv_m_cm = namd_reciprocal(m_cm);
       r_cm_x *= inv_m_cm;
       r_cm_y *= inv_m_cm;
       r_cm_z *= inv_m_cm;
       v_cm_x *= inv_m_cm;
       v_cm_y *= inv_m_cm;
       v_cm_z *= inv_m_cm;
 
       // XXX removed pairInteraction
       for ( j = i; j < (i+hgs); ++j ) {
         // XXX removed fixed atoms
 
         // vector vel[j] used twice below
         BigReal v_x = vel_x[j];
         BigReal v_y = vel_y[j];
         BigReal v_z = vel_z[j];
 
         // vector dv = vel[j] - v_cm
         BigReal dv_x = v_x - v_cm_x;
         BigReal dv_y = v_y - v_cm_y;
         BigReal dv_z = v_z - v_cm_z;
 
         // scalar intKineticEnergy += mass[j] * dot_product(v, dv)
         intKineticEnergy += mass[j] *
           (v_x * dv_x + v_y * dv_y + v_z * dv_z);
 
         // vector dr = pos[j] - r_cm
         BigReal dr_x = pos_x[j] - r_cm_x;
         BigReal dr_y = pos_y[j] - r_cm_y;
         BigReal dr_z = pos_z[j] - r_cm_z;
 
         // tensor intVirialNormal += outer_product(f_normal[j], dr)
         intVirialNormal.xx += f_normal_x[j] * dr_x;
         intVirialNormal.xy += f_normal_x[j] * dr_y;
         intVirialNormal.xz += f_normal_x[j] * dr_z;
         intVirialNormal.yx += f_normal_y[j] * dr_x;
         intVirialNormal.yy += f_normal_y[j] * dr_y;
         intVirialNormal.yz += f_normal_y[j] * dr_z;
         intVirialNormal.zx += f_normal_z[j] * dr_x;
         intVirialNormal.zy += f_normal_z[j] * dr_y;
         intVirialNormal.zz += f_normal_z[j] * dr_z;
 
         // tensor intVirialNbond += outer_product(f_nbond[j], dr)
         intVirialNbond.xx += f_nbond_x[j] * dr_x;
         intVirialNbond.xy += f_nbond_x[j] * dr_y;
         intVirialNbond.xz += f_nbond_x[j] * dr_z;
         intVirialNbond.yx += f_nbond_y[j] * dr_x;
         intVirialNbond.yy += f_nbond_y[j] * dr_y;
         intVirialNbond.yz += f_nbond_y[j] * dr_z;
         intVirialNbond.zx += f_nbond_z[j] * dr_x;
         intVirialNbond.zy += f_nbond_z[j] * dr_y;
         intVirialNbond.zz += f_nbond_z[j] * dr_z;
 
         // tensor intVirialSlow += outer_product(f_slow[j], dr)
         intVirialSlow.xx += f_slow_x[j] * dr_x;
         intVirialSlow.xy += f_slow_x[j] * dr_y;
         intVirialSlow.xz += f_slow_x[j] * dr_z;
         intVirialSlow.yx += f_slow_y[j] * dr_x;
         intVirialSlow.yy += f_slow_y[j] * dr_y;
         intVirialSlow.yz += f_slow_y[j] * dr_z;
         intVirialSlow.zx += f_slow_z[j] * dr_x;
         intVirialSlow.zy += f_slow_z[j] * dr_y;
         intVirialSlow.zz += f_slow_z[j] * dr_z;
       }
     }
 
     intKineticEnergy *= 0.5;
 
 #ifdef NODEGROUP_FORCE_REGISTER
     if(simParams->CUDASOAintegrate){
       // JM: Every PE will have its own copy of CUDASequencer, since it's a 
       //     group. However, they all share the same nodegroup pointer, 
       //     which means the reduction object is the same across all PEs
       CUDASequencer->patchData->reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
       ADD_TENSOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
       ADD_TENSOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
       ADD_TENSOR_OBJECT(CUDASequencer->patchData->reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
     }else{
 #endif
       reduction->item(REDUCTION_INT_CENTERED_KINETIC_ENERGY) += intKineticEnergy;
       ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NORMAL,intVirialNormal);
       ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_NBOND,intVirialNbond);
       ADD_TENSOR_OBJECT(reduction,REDUCTION_INT_VIRIAL_SLOW,intVirialSlow);
 #ifdef NODEGROUP_FORCE_REGISTER
      }
 #endif
   }
   // XXX removed pressure profile
 
   // XXX removed fixed atoms
 
   if(!simParams->CUDASOAintegrate) reduction->submit();
 
   // XXX removed pressure profile reduction
 }

◆ suspend()

void Sequencer::suspend ( void )

Definition at line 267 of file Sequencer.C.

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

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

 {
     LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
     CthSuspend();
     LdbCoordinator::Object()->startWork(patch->ldObjHandle);
 }

◆ suspendULTs()

void Sequencer::suspendULTs ( )

protected

◆ tcoupleVelocities()

void Sequencer::tcoupleVelocities	(	BigReal	dt_fs,
		int	step
	)

protected

Definition at line 5591 of file Sequencer.C.

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

Referenced by integrate().

 {
   if ( simParams->tCoupleOn )
   {
     FullAtom *a = patch->atom.begin();
     int numAtoms = patch->numAtoms;
     // Blocking receive for the temperature coupling coefficient.
     BigReal coefficient = broadcast->tcoupleCoefficient.get(step);
     Molecule *molecule = Node::Object()->molecule;
     BigReal dt = dt_fs * 0.001;  // convert to ps
     coefficient *= dt;
     for ( int i = 0; i < numAtoms; ++i )
     {
       BigReal f1 = exp( coefficient * a[i].langevinParam );
       a[i].velocity *= f1;
     }
   }
 }

◆ terminate()

void Sequencer::terminate ( void )

protected

Definition at line 6681 of file Sequencer.C.

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

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

                           {
   LdbCoordinator::Object()->pauseWork(patch->ldObjHandle);
   CthFree(thread);
   CthSuspend();
 }

◆ traceBarrier()

void Sequencer::traceBarrier ( int step )

protected

Definition at line 6669 of file Sequencer.C.

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

Referenced by integrate().

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

◆ updateDeviceData()

void Sequencer::updateDeviceData	(	const int	startup,
		const int	maxForceUsed,
		const int	doGlobal
	)

protected

◆ updateDevicePatchMap()

void Sequencer::updateDevicePatchMap ( int startup )

protected

◆ wakeULTs()

void Sequencer::wakeULTs ( )

protected

Friends And Related Function Documentation

◆ HomePatch

friend class HomePatch

friend

Definition at line 45 of file Sequencer.h.

◆ SequencerCUDA

friend class SequencerCUDA

friend

Definition at line 47 of file Sequencer.h.

Referenced by Sequencer().

Member Data Documentation

◆ adaptTempT

BigReal Sequencer::adaptTempT

protected

Definition at line 270 of file Sequencer.h.

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

◆ berendsenPressure_count

int Sequencer::berendsenPressure_count

protected

Definition at line 292 of file Sequencer.h.

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

◆ broadcast

ControllerBroadcasts* Sequencer::broadcast

protected

◆ checkpoint_berendsenPressure_count

int Sequencer::checkpoint_berendsenPressure_count

protected

Definition at line 293 of file Sequencer.h.

Referenced by algorithm().

◆ collection

CollectionMgr* const Sequencer::collection

protected

Definition at line 325 of file Sequencer.h.

Referenced by submitCollections(), and submitCollections_SOA().

◆ doKineticEnergy

int Sequencer::doKineticEnergy

protected

Definition at line 308 of file Sequencer.h.

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

◆ doMomenta

int Sequencer::doMomenta

protected

Definition at line 309 of file Sequencer.h.

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

◆ ldbSteps

int Sequencer::ldbSteps

protected

Definition at line 328 of file Sequencer.h.

Referenced by rebalanceLoad().

◆ masterThread

bool Sequencer::masterThread

protected

Definition at line 329 of file Sequencer.h.

Referenced by Sequencer().

◆ min_reduction

SubmitReduction* Sequencer::min_reduction

protected

Definition at line 226 of file Sequencer.h.

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

◆ multigratorReduction

SubmitReduction* Sequencer::multigratorReduction

protected

Definition at line 307 of file Sequencer.h.

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

◆ pairlistsAge

int Sequencer::pairlistsAge

protected

Definition at line 230 of file Sequencer.h.

Referenced by runComputeObjects(), and runComputeObjects_SOA().

◆ pairlistsAgeLimit

int Sequencer::pairlistsAgeLimit

protected

Definition at line 231 of file Sequencer.h.

Referenced by runComputeObjects(), runComputeObjects_SOA(), and Sequencer().

◆ pairlistsAreValid

int Sequencer::pairlistsAreValid

protected

Definition at line 229 of file Sequencer.h.

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

◆ patch

HomePatch* const Sequencer::patch

protected

Definition at line 321 of file Sequencer.h.

Referenced by addForceToMomentum(), addForceToMomentum3(), addForceToMomentum_SOA(), addMovDragToPosition(), addRotDragToPosition(), addVelocityToPosition(), addVelocityToPosition_SOA(), algorithm(), berendsenPressure(), berendsenPressure_SOA(), calcFixVirial(), calcKineticEnergy(), correctMomentum(), hardWallDrude(), integrate(), integrate_SOA(), langevinPiston(), langevinPiston_SOA(), langevinVelocities(), langevinVelocitiesBBK1(), langevinVelocitiesBBK1_SOA(), langevinVelocitiesBBK2(), langevinVelocitiesBBK2_SOA(), maximumMove(), maximumMove_SOA(), minimizationQuenchVelocity(), minimize(), minimizeMoveDownhill(), multigratorPressure(), multigratorTemperature(), newMinimizeDirection(), newMinimizePosition(), newtonianVelocities(), quenchVelocities(), rattle1(), rattle1_SOA(), reassignVelocities(), rebalanceLoad(), reinitVelocities(), reloadCharges(), rescaleaccelMD(), rescaleSoluteCharges(), rescaleVelocities(), rescaleVelocitiesByFactor(), run(), runComputeObjects(), runComputeObjects_SOA(), saveForce(), scalePositionsVelocities(), scaleVelocities(), Sequencer(), stochRescaleVelocities(), stochRescaleVelocities_SOA(), submitCollections(), submitCollections_SOA(), submitHalfstep(), submitHalfstep_SOA(), submitMinimizeReductions(), submitMomentum(), submitReductions(), submitReductions_SOA(), suspend(), tcoupleVelocities(), and terminate().

◆ pressureProfileReduction

SubmitReduction* Sequencer::pressureProfileReduction

protected

Definition at line 323 of file Sequencer.h.

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

◆ random

Random* Sequencer::random

protected

Definition at line 319 of file Sequencer.h.

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

◆ reduction

SubmitReduction* Sequencer::reduction

protected

Definition at line 322 of file Sequencer.h.

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

◆ rescaleVelocities_numTemps

int Sequencer::rescaleVelocities_numTemps

protected

Definition at line 275 of file Sequencer.h.

Referenced by rescaleVelocities(), and Sequencer().

◆ simParams

SimParameters* const Sequencer::simParams

protected

Definition at line 320 of file Sequencer.h.

Referenced by adaptTempUpdate(), addMovDragToPosition(), addRotDragToPosition(), algorithm(), berendsenPressure(), berendsenPressure_SOA(), calcFixVirial(), calcKineticEnergy(), correctMomentum(), hardWallDrude(), integrate(), integrate_SOA(), langevinPiston(), langevinPiston_SOA(), langevinVelocities(), langevinVelocitiesBBK1(), langevinVelocitiesBBK1_SOA(), langevinVelocitiesBBK2(), langevinVelocitiesBBK2_SOA(), maximumMove(), minimizationQuenchVelocity(), minimize(), minimizeMoveDownhill(), multigratorPressure(), multigratorTemperature(), newMinimizeDirection(), newMinimizePosition(), rattle1(), rattle1_SOA(), reassignVelocities(), reinitVelocities(), rescaleaccelMD(), rescaleSoluteCharges(), rescaleVelocities(), runComputeObjects(), runComputeObjects_SOA(), scalePositionsVelocities(), Sequencer(), stochRescaleVelocities(), stochRescaleVelocities_SOA(), submitCollections_SOA(), submitHalfstep(), submitHalfstep_SOA(), submitMinimizeReductions(), submitMomentum(), submitReductions(), submitReductions_SOA(), tcoupleVelocities(), and ~Sequencer().

◆ slowFreq

int Sequencer::slowFreq

protected

Definition at line 295 of file Sequencer.h.

Referenced by integrate(), integrate_SOA(), and langevinPiston().

◆ stochRescale_count

int Sequencer::stochRescale_count

protected

Count time steps until next stochastic velocity rescaling.

Definition at line 288 of file Sequencer.h.

Referenced by Sequencer(), stochRescaleVelocities(), and stochRescaleVelocities_SOA().

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

Public Member Functions

Protected Member Functions

Protected Attributes

Friends

Detailed Description

Constructor & Destructor Documentation

◆ Sequencer()

◆ ~Sequencer()

Member Function Documentation

◆ adaptTempUpdate()

◆ addForceToMomentum()

◆ addForceToMomentum3()

◆ addForceToMomentum_SOA()

◆ addMovDragToPosition()

◆ addRotDragToPosition()

◆ addVelocityToPosition()

◆ addVelocityToPosition_SOA()

◆ algorithm()

◆ awaken()

◆ berendsenPressure()

◆ berendsenPressure_SOA()

◆ calcFixVirial()

◆ calcKineticEnergy()

◆ clearDevicePatchMap()

◆ constructDevicePatchMap()

◆ correctMomentum()

◆ cycleBarrier()

◆ doMigrationGPU()

◆ hardWallDrude()

◆ initialize_integrate_CUDA_SOA()

◆ integrate()

◆ integrate_CUDA_SOA()

◆ integrate_SOA()

◆ langevinPiston()

◆ langevinPiston_SOA()

◆ langevinVelocities()

◆ langevinVelocitiesBBK1()

◆ langevinVelocitiesBBK1_SOA()

◆ langevinVelocitiesBBK2()

◆ langevinVelocitiesBBK2_SOA()

◆ maximumMove()

◆ maximumMove_SOA()

◆ minimizationQuenchVelocity()

◆ minimize()

◆ minimizeMoveDownhill()

◆ monteCarloPressureControl()

◆ multigratorPressure()

◆ multigratorTemperature()

◆ newMinimizeDirection()

◆ newMinimizePosition()

◆ newtonianVelocities()

◆ printDevicePatchMap()

◆ quenchVelocities()

◆ rattle1()

◆ rattle1_SOA()

◆ reassignVelocities()

◆ rebalanceLoad()

◆ reinitVelocities()

◆ reloadCharges()

◆ rescaleaccelMD()

◆ rescaleSoluteCharges()

◆ rescaleVelocities()

◆ rescaleVelocitiesByFactor()

◆ run()

◆ runComputeObjects()

◆ runComputeObjects_SOA()

◆ runComputeObjectsCUDA()

◆ saveForce()

◆ scalePositionsVelocities()

◆ scaleVelocities()

◆ stochRescaleVelocities()

◆ stochRescaleVelocities_SOA()

◆ submitCollections()

◆ submitCollections_SOA()

◆ submitHalfstep()

◆ submitHalfstep_SOA()

◆ submitMinimizeReductions()

◆ submitMomentum()

◆ submitReductions()

◆ submitReductions_SOA()