Sequencer Class Reference

#include <Sequencer.h>

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

Protected Member Functions
virtual void	algorithm (void)
void	integrate_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 45 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, PatchData::nodeBroadcast, PatchMap::numPatches(), PatchMap::Object(), GlobalGPUMgr::Object(), SynchronousCollectives::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;

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    if (simParams->CUDASOAintegrateMode) {
      CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
      PatchData* patchData = cpdata.ckLocalBranch();
      broadcast = new ControllerBroadcasts(&patch->ldObjHandle, patchData->nodeBroadcast);
    } else
#endif // defined(NAMD_CUDA) || defined(NAMD_HIP)
    {
      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);
      
      syncColl = SynchronousCollectives::Object();
      globalGPUMgr = GlobalGPUMgr::Object();
      constructDevicePatchMap();
#endif
    }
#endif
}

~Sequencer()

Sequencer::~Sequencer ( void  )

virtual

Definition at line 245 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 5467 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 5626 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 5641 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 2743 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 4414 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 4433 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 5665 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 2846 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 289 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 55 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 5261 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 3426 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 5963 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 4954 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 4743 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 6652 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 5677 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 3887 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::doFullDispersion, 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::LJPMEOn, SimParameters::LJPMESerialRealSpaceOn, 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;

    // If doing LJ-PME, track doFullElectrostatics
    int &doFullDispersion = patch->flags.doFullDispersion;
if ( ! simParams->LJPMESerialRealSpaceOn ) {
    doFullDispersion = (simParams->LJPMEOn && doFullElectrostatics);
}

//#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));
if ( ! simParams->LJPMESerialRealSpaceOn ) {
      // XXX in preparation for supporting LJ-PME with MTS
      doFullDispersion = (simParams->LJPMEOn && doFullElectrostatics);
}

#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 2068 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() >= simParams->beginEventPatchID && patch->getPatchID() <= simParams->endEventPatchID);
#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==0)
          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 5325 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 3525 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 5044 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 5081 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 3297 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 5154 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 3343 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 5749 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 3239 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 5794 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 4456 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 4602 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 4811 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 4974 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 4624 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 4683 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 5020 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 4711 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 5692 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 3654 of file Sequencer.C.

References ADD_TENSOR_OBJECT, 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();
    }
    ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
  }
}

reassignVelocities()

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

protected

Definition at line 5483 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 6641 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 5515 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 5548 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 5424 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 5560 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 5405 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 5538 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 269 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 6427 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 3671 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 5615 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 4774 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 4946 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 5599 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 3857 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 6390 of file Sequencer.C.

References broadcast, collection, Output::coordinateNeeded(), IMDSessionInfo::energies_switch, Patch::f, SimParameters::firstTimestep, Patch::flags, Output::forceNeeded(), SimpleBroadcastObject< T >::get(), SimParameters::IMDfreq, SimParameters::IMDon, SimParameters::IMDsendsettings, ControllerBroadcasts::IMDTimeEnergyBarrier, IMDv3, SimParameters::IMDversion, Patch::lattice, Flags::maxForceUsed, NAMD_EVENT_RANGE_2, patch, simParams, Results::slow, CollectionMgr::submitForces(), CollectionMgr::submitPositions(), CollectionMgr::submitVelocities(), IMDSessionInfo::time_switch, 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();

  if (simParams->IMDon &&
      (simParams->IMDversion == IMDversion_t::IMDv3) &&
      !(step % simParams->IMDfreq) &&
      ((simParams->IMDsendsettings.time_switch == 1) || (simParams->IMDsendsettings.energies_switch == 1)) &&
      (step != simParams->firstTimestep)) {
    broadcast->IMDTimeEnergyBarrier.get(step);
  }

  NAMD_EVENT_RANGE_2(patch->flags.event_on,
      NamdProfileEvent::SUBMIT_COLLECTIONS);
  int prec;
  int dcdSelectionIndex;
  std::tie (prec, dcdSelectionIndex) = Output::coordinateNeeded(step);
  if ( prec ) {
    collection->submitPositions(step,patch->atom,patch->lattice,prec,dcdSelectionIndex);
  }
  prec = Output::velocityNeeded(step);
  if ( prec ) {
    collection->submitVelocities(step,zeroVel,patch->atom,prec);
  }
  prec = Output::forceNeeded(step);
  if ( prec ) {
    int maxForceUsed = patch->flags.maxForceUsed;
    if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
    collection->submitForces(step,patch->atom,maxForceUsed,patch->f,prec);
  }
}

submitCollections_SOA()

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

protected

Note

The device migration code path only copies the positions and velocities by copyAoSDataToHost(), so to support force DCD output we have to copy the SOA forces to AOS again.

Definition at line 3190 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 ) {
    collection->submitPositions(step,patch->atom,patch->lattice,is_pos_needed,dcdIndex);
  }
  if ( is_vel_needed ) {
    collection->submitVelocities(step,zeroVel,patch->atom,is_vel_needed);
  }
  if ( is_f_needed ) {
    int maxForceUsed = patch->flags.maxForceUsed;
    if ( maxForceUsed > Results::slow ) maxForceUsed = Results::slow;
    collection->submitForces(step,patch->atom,maxForceUsed,patch->f,is_f_needed);
  }
}

submitHalfstep()

void Sequencer::submitHalfstep ( int  step )

protected

Definition at line 5805 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 2885 of file Sequencer.C.

References ADD_TENSOR_OBJECT, 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, 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;

    reduction->item(REDUCTION_HALFSTEP_KINETIC_ENERGY) += kineticEnergy;
    ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
  }

  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;
    reduction->item(REDUCTION_INT_HALFSTEP_KINETIC_ENERGY)
      += intKineticEnergy;
    ADD_TENSOR_OBJECT(reduction, REDUCTION_INT_VIRIAL_NORMAL,
      intVirialNormal);
  }
}

submitMinimizeReductions()

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

protected

Definition at line 6241 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 4720 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 5983 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 2985 of file Sequencer.C.

References ADD_TENSOR_OBJECT, ADD_VECTOR_OBJECT, 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, 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

  reduction->item(REDUCTION_ATOM_CHECKSUM) += numAtoms;
  reduction->item(REDUCTION_MARGIN_VIOLATIONS) += patch->marginViolations;

  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);

    reduction->item(REDUCTION_CENTERED_KINETIC_ENERGY) += kineticEnergy;
    ADD_VECTOR_OBJECT(reduction,REDUCTION_MOMENTUM,momentum);
    ADD_VECTOR_OBJECT(reduction,REDUCTION_ANGULAR_MOMENTUM,angularMomentum);
  }
  // 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;

    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);
  }
  // XXX removed pressure profile

  // XXX removed fixed atoms

  reduction->submit();

  // XXX removed pressure profile reduction
}

suspend()

void Sequencer::suspend

(

void

)

Definition at line 279 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 5576 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 6672 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 6660 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 47 of file Sequencer.h.

SequencerCUDA

friend class SequencerCUDA

friend

Definition at line 49 of file Sequencer.h.

Referenced by Sequencer().

Member Data Documentation

adaptTempT

BigReal Sequencer::adaptTempT

protected

Definition at line 272 of file Sequencer.h.

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

berendsenPressure_count

int Sequencer::berendsenPressure_count

protected

Definition at line 294 of file Sequencer.h.

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

broadcast

ControllerBroadcasts* Sequencer::broadcast

protected

Definition at line 328 of file Sequencer.h.

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

checkpoint_berendsenPressure_count

int Sequencer::checkpoint_berendsenPressure_count

protected

Definition at line 295 of file Sequencer.h.

Referenced by algorithm().

collection

CollectionMgr* const Sequencer::collection

protected

Definition at line 327 of file Sequencer.h.

Referenced by submitCollections(), and submitCollections_SOA().

doKineticEnergy

int Sequencer::doKineticEnergy

protected

Definition at line 310 of file Sequencer.h.

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

doMomenta

int Sequencer::doMomenta

protected

Definition at line 311 of file Sequencer.h.

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

ldbSteps

int Sequencer::ldbSteps

protected

Definition at line 330 of file Sequencer.h.

Referenced by rebalanceLoad().

masterThread

bool Sequencer::masterThread

protected

Definition at line 331 of file Sequencer.h.

Referenced by Sequencer().

min_reduction

SubmitReduction* Sequencer::min_reduction

protected

Definition at line 228 of file Sequencer.h.

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

multigratorReduction

SubmitReduction* Sequencer::multigratorReduction

protected

Definition at line 309 of file Sequencer.h.

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

pairlistsAge

int Sequencer::pairlistsAge

protected

Definition at line 232 of file Sequencer.h.

Referenced by runComputeObjects(), and runComputeObjects_SOA().

pairlistsAgeLimit

int Sequencer::pairlistsAgeLimit

protected

Definition at line 233 of file Sequencer.h.

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

pairlistsAreValid

int Sequencer::pairlistsAreValid

protected

Definition at line 231 of file Sequencer.h.

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

patch

HomePatch* const Sequencer::patch

protected

Definition at line 323 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 325 of file Sequencer.h.

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

random

Random* Sequencer::random

protected

Definition at line 321 of file Sequencer.h.

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

reduction

SubmitReduction* Sequencer::reduction

protected

Definition at line 324 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 277 of file Sequencer.h.

Referenced by rescaleVelocities(), and Sequencer().

simParams

SimParameters* const Sequencer::simParams

protected

Definition at line 322 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(), submitCollections_SOA(), submitHalfstep(), submitMinimizeReductions(), submitMomentum(), submitReductions(), tcoupleVelocities(), and ~Sequencer().

slowFreq

int Sequencer::slowFreq

protected

Definition at line 297 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 290 of file Sequencer.h.

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

---

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

• Sequencer.h
• Sequencer.C