framework/integrators/ShadowHMCIntegrator.cpp

#include "Report.h"
#include "mathutilities.h"
#include "pmconstants.h"
#include "Vector3DBlock.h"
#include "ScalarStructure.h"
#include "topologyutilities.h"
#include "GenericTopology.h"
#include "ModifierUpdateBeta.h"
#include "ModifierUpdateBetaAndPush.h"
#include "ShadowHMCIntegrator.h"

#include <iomanip>
using namespace std;

using namespace ProtoMol::Report;
using std::vector;
using std::string;
using std::deque;

namespace ProtoMol {

  //____________________________________________________________ ShadowHMCIntegrator

  const string ShadowHMCIntegrator::keyword("ShadowHMC");

  ShadowHMCIntegrator::ShadowHMCIntegrator()
    : MCIntegrator(),
      myOrder(0),
      myShadowK(0),
      myShadowKover2(0),
      myC(0.0),
      myBeta(0.0),
      sumTotalEnergy(0.0),
      myPrevVelocities(NULL),
      myPrevPositions(NULL),
      myPrevBeta(NULL),
      myModifier1(NULL),
      myModifier2(NULL),
      myOldForces(NULL){}

  //  ---------------------------------------------------------------------  //

  ShadowHMCIntegrator::ShadowHMCIntegrator(int cycles,
					   bool randomCycLen,
					   Real initialTemperature,
					   unsigned int shadowOrder,
					   Real c,
					   ForceGroup *overloadedForces,
					   StandardIntegrator *nextIntegrator)

    : MCIntegrator(cycles,randomCycLen,initialTemperature,overloadedForces,nextIntegrator),
      myOrder(shadowOrder),

      //   This only works for (order == 2k), where k is even.  If/when we
      //   implement k odd, this will need to be changed.
      myShadowK(myOrder >> 1),
      myShadowKover2(myOrder >> 2),

      myC(c),
      myBeta(0.0),
      sumTotalEnergy(0.0),
      myPrevVelocities(new deque<Vector3DBlock>()),
      myPrevPositions(new deque<Vector3DBlock>()),
      myPrevBeta(new deque<Real>()),
      myModifier1(NULL),
      myModifier2(NULL),
      myOldForces(new Vector3DBlock()) {
  }

  //  ---------------------------------------------------------------------  //

  ShadowHMCIntegrator::~ShadowHMCIntegrator(){
    if(myPrevVelocities != NULL)
      delete myPrevVelocities;
    if(myPrevPositions != NULL)
      delete myPrevPositions;
    if(myPrevBeta != NULL)
      delete myPrevBeta;
    if(myOldForces != NULL)
      delete myOldForces;
  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::initialize( GenericTopology *topo,
				        Vector3DBlock   *positions,
				        Vector3DBlock   *velocities,
				        ScalarStructure *energies ) {

    MCIntegrator::initialize(topo,positions,velocities,energies);
    Report::report << debug(10) <<"[ShadowHMCIntegrator::initialize]"<<Report::endr;

    // Add callbacks
    myModifier1 = new ModifierUpdateBeta(this);
    next()->adoptPreDriftOrNextModifier(myModifier1);

    //  This one is important because it calls the function to save the shadow
    //  data.
    myModifier2 = new ModifierUpdateBetaAndPush(this);
    next()->adoptPostStepModifier(myModifier2);

    //  FIXME  This needs to be hardcoded eventually.  FIXME
    //  Number of samples to estimate myC.
//     int numSamples = 10;
//
//     myC = -0.1 * ( myEnergies->potentialEnergy() +
//             kineticEnergy(myTopo,myVelocities) );
//
//     report << debug(1) << "myC = " << myC << endr;
//
//     run(numSamples);
//
//     myC = 0.90 * ( sumTotalEnergy / numSamples );
//
//     report << debug(1) << "myC = " << myC << endr;

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::run(int numTimesteps) {

    Real initTotEnergy    = 0.0,
         initShadowEnergy = 0.0,
         finShadowEnergy  = 0.0,
         finTotEnergy     = 0.0,
         MDfinEner        = 0.0,
         MDinitEner       = 0.0,
         momentumEnerDiff = 0.0;

    bool acceptMomenta = false;

    //  See note in header file.
    unsigned int localCycleLength = ( myRandomCycLen ? static_cast<int>(
                myCycleLength * ( 1.3 - 0.6 * randomNumber() ) ) : myCycleLength );

    //  Check that cyclelength is long enough.  It is possible to compute the
    //  shadow with a short cyclelength, but it requires a lot of uncessesary
    //  checks.  In general, cyclelength should be much larger than
    //  myShadowK.
    if(localCycleLength < myShadowK) {
      report << warning << "***  The cyclelength for this step is below the "
	     << "number of steps required for accurate computation of the "
	     << "shadow energy.  Adjusting to minimal cycle legnth.  ***" << endr;
      localCycleLength = myShadowK;
    }

    //  Start algorithm.
    for(int i = 0; i < numTimesteps; i++) {

      preStepModify();

      //  Save current positions, forces, and energies in case we reject.
      saveValues();

      //  Repeat until new momenta are accepted.
      acceptMomenta = false;

      while( !acceptMomenta ) {

	//  Clear the history queues and reset beta.
	resetHistory();

	//  Calculate new random velocities.
	perturbSystem();

	//  Save the initial total energy for the acceptance step.
	initTotEnergy = myEnergies->potentialEnergy() + kineticEnergy(myTopo,myVelocities);

        //  Run steps -k/2 to -1.
        runPreSteps( myShadowKover2 );

        //  Save data at step 0.
        pushShadowHistory();

	//  Run steps 1 to k/2.
        next()->run( myShadowKover2 );

        //  We now have data for steps -k/2 -> k/2, so we can calculate the
        //  shadow at step 0 and save it as the initial shadow.
	calculateShadow();

	initShadowEnergy = (*myEnergies)[ScalarStructure::SHADOW];

        //  Accept these momenta with probability proportional to:
        //  min{ 1, exp^(-\beta(H_[2k] - c - H)) }.
	momentumEnerDiff = initShadowEnergy - initTotEnergy - myC;

        //  Decide whether or not to accept new momenta.
	if( metropolisTest( momentumEnerDiff, 0.0 ) ) {

	  acceptMomenta = true;

	  report.precision(8);
          report << debug(1) << "Take momenta (diff = " << momentumEnerDiff <<
                  ")" << endr;

	}
	else {

	  //  Restore saved values from step 0.
	  restoreValues();

	  report.precision(8);
          report << debug(1) << "Deny momenta (diff = " << momentumEnerDiff <<
                  ")" << endr;

	}

      }

      //  Choose max{ H_[2k] - c, H } to be the initial 'energy'.
      MDinitEner = std::max(initShadowEnergy - myC, initTotEnergy);

      sumTotalEnergy += initShadowEnergy - initTotEnergy;

      //  In between the beginning and end of the cycle, we don't need to
      //  save previous values or calculate the shadow so run using LF.  Run
      //  from step k/2 + 1 -> n - k/2 - 1.
      disableBetaUpdate();
      next()->run(localCycleLength - myShadowK - 1);
      enableBetaUpdate();

      //  Finish by running steps n - k/2 -> n.
      next()->run( myShadowKover2 + 1 );

      //  Run k/2 more steps so we can compute the shadow at step n.
      runPostSteps( myShadowKover2 );

      //  Make note of the final total and shadow energies.
      finTotEnergy = myEnergies->potentialEnergy() + kineticEnergy(myTopo,myVelocities);
      finShadowEnergy = (*myEnergies)[ScalarStructure::SHADOW];


      //  FIXME  Temporary.  Used for figuring out a good value of c.

//       cout.setf( ios::fixed );
//       cout << endl << "DEBUG:  dSE: " << setw(10) << setprecision(4) <<
//               finShadowEnergy;
//       cout << " - " << setw(10) << setprecision(4) << initShadowEnergy;
//       cout << " = " << setw(10) << setprecision(4) << finShadowEnergy -
//               initShadowEnergy << endl;
//
//       cout << "DEBUG:  dTE: " << setw(10) << setprecision(4) << finTotEnergy;
//       cout << " - " << setw(10) << setprecision(4) << initTotEnergy;
//       cout << " = " << setw(10) << setprecision(4) << finTotEnergy -
//               initTotEnergy << endl;
//
//       cout << "DEBUG:   dE: " << setw(10) << setprecision(5) << finShadowEnergy
//               - finTotEnergy;
//       cout << " - " << setw(10) << setprecision(5) << initShadowEnergy -
//               initTotEnergy;
//       cout << " = " << setw(10) << setprecision(5) << finShadowEnergy -
//               finTotEnergy - initShadowEnergy + initTotEnergy << endl << endl;

      //  Choose max{ H_[2k] - c, H } to be the final 'energy'.
      MDfinEner = std::max(finShadowEnergy - myC, finTotEnergy);

      //  Based on the difference in the max energies, decide to accept
      //  the new positions or not.
      if( metropolisTest( MDfinEner, MDinitEner ) ) {

	//  Save the new positions.
	saveValues();

	//  Print debug information if desired.
	report << debug(1) << "Accepted (diff = " << MDfinEner - MDinitEner <<")"<< endr;

      }
      else {

	//  If we rejected, the we have a lot of clean up to do.
	restoreValues();

	report << debug(1) << "Rejected (diff = " << MDfinEner - MDinitEner <<")" << endr;

      }

      postStepModify();

    }

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::runPreSteps(int numTimesteps){

    //  Save all values for current timestep (probably time 0).
    Vector3DBlock   positions(*myPositions);
    Vector3DBlock   velocities(*myVelocities);
    ScalarStructure energies(*myEnergies);
    Vector3DBlock   forces(*(next()->getForces()));

    //  Flip clock to run backwards.
    backward();

    next()->run( numTimesteps );

    //  Flip clock to run forward.
    forward();

    //  Restore original values.
    myPositions->intoAssign(positions);
    myVelocities->intoAssign(velocities);
    (*myEnergies) = energies;
    next()->getForces()->intoAssign(forces);
    myBeta = 0.0;

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::runPostSteps(int numTimesteps){

    //  We need to save the current state of the system because the extra steps
    //  we are running are only to compute the shadow.  All other values must
    //  remain at step n.

    Vector3DBlock   positions(*myPositions);
    Vector3DBlock   velocities(*myVelocities);
    ScalarStructure energies(*myEnergies);
    Vector3DBlock   forces(*(next()->getForces()));

    //  Run the extra steps needed to compute the shadow at step n.
    next()->run( numTimesteps );
    calculateShadow();

    //  Restore the energies, forces, positions, velocities and shadow.
    Real finalShEnergy = (*myEnergies)[ScalarStructure::SHADOW];
    myPositions->intoAssign(positions);
    myVelocities->intoAssign(velocities);
    (*myEnergies) = energies;
    next()->getForces()->intoAssign(forces);
    (*myEnergies)[ScalarStructure::SHADOW] = finalShEnergy;

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::updateBeta() {
    Report::report << debug(10) <<"[ShadowHMCIntegrator::updateBeta]"<<Report::endr;
    //  Compute the timestep for a 'half' kick.
    Real h = 0.5 * next()->getTimestep() * Constant::INV_TIMEFACTOR;

    //  Compute the dot product of the positions and forces.
    Real posDotF = 0.0;
    const Vector3DBlock& forces(*next()->getForces());
    for(unsigned int i = 0; i < myPositions->size(); i++) {
      posDotF += (*myPositions)[i].dot(forces[i]);
    }

    //  Update beta.
    myBeta -= h * ( posDotF + 2.0 * myEnergies->potentialEnergy() );

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::pushShadowHistory(){

    Report::report<< debug(10)  <<"[ShadowHMCIntegrator::pushShadowHistory]"<<Report::endr;

    //  When simulating forward in time, add new values to end of queue.
    if(getTimestep() > 0){

      myPrevBeta->push_back( myBeta );
      myPrevVelocities->push_back( *myVelocities );
      myPrevPositions->push_back( *myPositions );

    }
    //  When simulating backward in time, push new values to front of queue.
    else {

      myPrevBeta->push_front( myBeta );
      myPrevVelocities->push_front( *myVelocities );
      myPrevPositions->push_front( *myPositions );

    }

    //  For an order (2k) shadow Hamiltonian, we only need to store (k)+1    //
    //  values.  Release old values.                                         //
    if( myPrevBeta->size() > ( myShadowK + 1 ) ) {

      //  The most recent values are added to the end of the queue, so we remove
      //  from the front.
      if(getTimestep() > 0){

	myPrevBeta->pop_front();
	myPrevVelocities->pop_front();
	myPrevPositions->pop_front();

      }
      else {

	myPrevBeta->pop_back();
	myPrevVelocities->pop_back();
	myPrevPositions->pop_back();

      }

    }

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::calculateShadow() {

    //  Make sure that we have enough saved values to compute the shadow.
    if( myPrevBeta->size() != ( myShadowK + 1 ) ) {
        report << warning << "Called calculateShadow() with less than " <<
                ( myShadowK + 1 ) << " values." << endr;
        return;
    }

    //  Return the appropriate shadow.
    switch( myOrder ) {
    case(4) :
      (*myEnergies)[ScalarStructure::SHADOW] = calcShadow4();
      break;
    case(8) :
      (*myEnergies)[ScalarStructure::SHADOW] = calcShadow8();
      break;
    default :
      return;
    }
  }

  //  ---------------------------------------------------------------------  //

  Real ShadowHMCIntegrator::calcShadow4() {

    Real h = next()->getTimestep() * Constant::INV_TIMEFACTOR;
    Real pos1DotVel0 = 0.0;
    Real pos1DotVel2 = 0.0;
    Real vel1DotPos0 = 0.0;
    Real vel1DotPos2 = 0.0;
    Real a10 = 0.0;
    Real a12 = 0.0;
    const Real c = 1. / (2 * h);
    Real mass_div_2 = 0.;

    //  The shadow Hamiltonian is a symmetric method about a timestep 'n'.  The
    //  way that the data is stored and the fact that myShadowK is 0.5 * the
    //  order of the shadow allows me to access the data in a form consistent
    //  with the algorithm.  i.e.  data[n-2], data[n-1] ...
    int n = myShadowKover2;

    //  Compute the shadow on my range.
    for( unsigned int i = 0; i < myPositions->size(); i++ ) {

      //  Get the mass for this atom.
      mass_div_2 = 0.5 * myTopo->atoms[i].scaledMass;

      //  Compute the dot products.
      pos1DotVel0 += mass_div_2 * (((*myPrevPositions)[n+1])[i] -
				   ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n])[i]);

      vel1DotPos0 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] -
				   ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n])[i]);

      pos1DotVel2 += mass_div_2 * (((*myPrevPositions)[n+1])[i] -
				   ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n+1])[i] -
								     ((*myPrevVelocities)[n])[i] * 2.0 +
								     ((*myPrevVelocities)[n-1])[i]);

      vel1DotPos2 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] -
				   ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n+1])[i] -
								      ((*myPrevPositions)[n])[i] * 2.0 +
								      ((*myPrevPositions)[n-1])[i]);

    }

    //  Compute the a_i_j.
    a10 = pos1DotVel0 - vel1DotPos0 - 0.5 * ( (*myPrevBeta)[n+1] - (*myPrevBeta)[n-1] );
    a12 = pos1DotVel2 - vel1DotPos2;

    report << debug(2) << "calcShadow4() = " << ( c * ( a10 - a12 / 6.0 ) ) <<
            endr;

    //  Calculate the 4th order shadow.
    return( c * ( a10 - a12 / 6.0 ) );

  }

  //  ---------------------------------------------------------------------  //

  Real ShadowHMCIntegrator::calcShadow8() {
    Real h = next()->getTimestep() * Constant::INV_TIMEFACTOR;
    Real pos1DotVel0 = 0.0;
    Real pos1DotVel2 = 0.0;
    Real pos1DotVel4 = 0.0;
    Real vel1DotPos0 = 0.0;
    Real vel1DotPos2 = 0.0;
    Real vel1DotPos4 = 0.0;
    Real pos3DotVel0 = 0.0;
    Real pos3DotVel2 = 0.0;
    Real pos3DotVel4 = 0.0;
    Real vel3DotPos0 = 0.0;
    Real vel3DotPos2 = 0.0;
    Real vel3DotPos4 = 0.;
    Real a10 = 0.0;
    Real a12 = 0.0;
    Real a14 = 0.0;
    Real a30 = 0.0;
    Real a32 = 0.0;
    Real a34 = 0.;
    const Real c = 1. / (2 * h);
    Real mass_div_2 = 0.;

    //  The shadow Hamiltonian is a symmetric method about a timestep 'n'.  The
    //  way that the data is stored and the fact that myShadowK is 0.5 * the
    //  order of the shadow allows me to access the data in a form consistent
    //  with the algorithm.  i.e.  data[n-2], data[n-1] ...
    int n = myShadowKover2;

    //  Compute the shadow on my range.
    for( unsigned int i = 0; i < myPositions->size(); i++ ) {

      mass_div_2 = 0.5 * myTopo->atoms[i].scaledMass;

      //  Compute the dot products.
      pos1DotVel0 += mass_div_2 * (((*myPrevPositions)[n+1])[i] - ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n])[i]);
      vel1DotPos0 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] - ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n])[i]);
      pos1DotVel2 += mass_div_2 * (((*myPrevPositions)[n+1])[i] -
				   ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n+1])[i] -
								     ((*myPrevVelocities)[n])[i] * 2.0 +
								     ((*myPrevVelocities)[n-1])[i]);
      vel1DotPos2 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] -
				   ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n+1])[i] -
								      ((*myPrevPositions)[n])[i] * 2.0 +
								      ((*myPrevPositions)[n-1])[i]);
      pos1DotVel4 += mass_div_2 * (((*myPrevPositions)[n+1])[i] -
				   ((*myPrevPositions)[n-1])[i]).dot(((*myPrevVelocities)[n+2])[i] -
								     ((*myPrevVelocities)[n+1])[i] * 4.0 +
								     ((*myPrevVelocities)[n])[i] * 6.0 -
								     ((*myPrevVelocities)[n-1])[i] * 4.0 +
								     ((*myPrevVelocities)[n-2])[i]);
      vel1DotPos4 += mass_div_2 * (((*myPrevVelocities)[n+1])[i] -
				   ((*myPrevVelocities)[n-1])[i]).dot(((*myPrevPositions)[n+2])[i] -
								      ((*myPrevPositions)[n+1])[i] * 4.0 +
								      ((*myPrevPositions)[n])[i] * 6.0 -
								      ((*myPrevPositions)[n-1])[i] * 4.0 +
								      ((*myPrevPositions)[n-2])[i]);
      pos3DotVel0 += mass_div_2 * (((*myPrevPositions)[n+2])[i] -
				   ((*myPrevPositions)[n+1])[i] * 2.0 +
				   ((*myPrevPositions)[n-1])[i] * 2.0 -
				   ((*myPrevPositions)[n-2])[i]).dot(((*myPrevVelocities)[n])[i]);
      vel3DotPos0 += mass_div_2 * (((*myPrevVelocities)[n+2])[i] -
				   ((*myPrevVelocities)[n+1])[i] * 2.0 +
				   ((*myPrevVelocities)[n-1])[i] * 2.0 -
				   ((*myPrevVelocities)[n-2])[i]).dot(((*myPrevPositions)[n])[i]);
      pos3DotVel2 += mass_div_2 * (((*myPrevPositions)[n+2])[i] -
				   ((*myPrevPositions)[n+1])[i] * 2.0 +
				   ((*myPrevPositions)[n-1])[i] * 2.0 -
				   ((*myPrevPositions)[n-2])[i]).dot(((*myPrevVelocities)[n+1])[i] -
								     ((*myPrevVelocities)[n])[i] * 2.0 +
								     ((*myPrevVelocities)[n-1])[i]);
      vel3DotPos2 += mass_div_2 * (((*myPrevVelocities)[n+2])[i] -
				   ((*myPrevVelocities)[n+1])[i] * 2.0 +
				   ((*myPrevVelocities)[n-1])[i] * 2.0 -
				   ((*myPrevVelocities)[n-2])[i]).dot(((*myPrevPositions)[n+1])[i] -
								      ((*myPrevPositions)[n])[i] * 2.0 +
								      ((*myPrevPositions)[n-1])[i]);
      pos3DotVel4 += mass_div_2 * (((*myPrevPositions)[n+2])[i] -
				   ((*myPrevPositions)[n+1])[i] * 2.0 +
				   ((*myPrevPositions)[n-1])[i] * 2.0 -
				   ((*myPrevPositions)[n-2])[i]).dot(((*myPrevVelocities)[n+2])[i] -
								     ((*myPrevVelocities)[n+1])[i] * 4.0 +
								     ((*myPrevVelocities)[n])[i] * 6.0 -
								     ((*myPrevVelocities)[n-1])[i] * 4.0 +
								     ((*myPrevVelocities)[n-2])[i]);
      vel3DotPos4 += mass_div_2 * (((*myPrevVelocities)[n+2])[i] -
				   ((*myPrevVelocities)[n+1])[i] * 2.0 +
				   ((*myPrevVelocities)[n-1])[i] * 2.0 -
				   ((*myPrevVelocities)[n-2])[i]).dot(((*myPrevPositions)[n+2])[i] -
								      ((*myPrevPositions)[n+1])[i] * 4.0 +
								      ((*myPrevPositions)[n])[i] * 6.0 -
								      ((*myPrevPositions)[n-1])[i] * 4.0 +
								      ((*myPrevPositions)[n-2])[i]);
    }

    //  Compute the a_i_j.
    a10 = pos1DotVel0 - vel1DotPos0 - 0.5 * ( (*myPrevBeta)[n+1] -  (*myPrevBeta)[n-1] );
    a12 = pos1DotVel2 - vel1DotPos2;
    a14 = pos1DotVel4 - vel1DotPos4;
    a30 = pos3DotVel0 - vel3DotPos0 - 0.5 * (*myPrevBeta)[n+2] + (*myPrevBeta)[n+1] - (*myPrevBeta)[n-1] + 0.5 * (*myPrevBeta)[n-2];
    a32 = pos3DotVel2 - vel3DotPos2;
    a34 = pos3DotVel4 - vel3DotPos4;

    report << debug(2) << "calcShadow8() = " << ( (c / 210.) * (210. * a10 + 25.
                * a30 +  26. * a32 - (60. * a12 + 19. * a14 + 1.5 * a34)) ) <<
            endr;

    //  Calculate the 8th order shadow.  I multiplied the algorithm by 210 in
    //  order to easily remove the 5 fractions.  I now only have to do 1
    //  division instead of 5.  I also grouped the additions and subtractions to
    //  avoid as much round off error as possible.

    return ((c / 210.) * (210. * a10 + 25. * a30 +  26. * a32 - (60. * a12 + 19. * a14 + 1.5 * a34)));
  }

  //  ---------------------------------------------------------------------  //

  MTSIntegrator*  ShadowHMCIntegrator::doMake(string& errMsg, const vector<Value>& values,
					      ForceGroup* fg, StandardIntegrator *nextIntegrator)const{

    unsigned int order = values[3];
    if(order != 4 && order != 8){
      errMsg += " Expecting order 4 or 8, but not "+values[3].getString();
      return NULL;
    }

    return new ShadowHMCIntegrator(values[0],values[1],values[2],values[3],values[4],fg,nextIntegrator);
  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::getParameters(vector< Parameter> &parameters) const {
    MCIntegrator::getParameters(parameters);
    parameters.push_back(Parameter("order", Value(myOrder,ConstraintValueType::Positive())));
    parameters.push_back(Parameter("c", Value(myC)));
  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::perturbSystem() {

    report << debug(10) << "[ShadowHMCIntegrator::perturbSystem]"<<endr;

    randomVelocity(getInitialTemperature(),myTopo,myVelocities);

    buildMolecularMomentum(myVelocities,myTopo);

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::resetHistory() {
    //  Reset the queues and beta term.
    myPrevBeta->clear();
    myPrevVelocities->clear();
    myPrevPositions->clear();
    myBeta = 0.;
  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::disableBetaUpdate(){
    if(myModifier1 != NULL)
      myModifier1->disable();
    if(myModifier2 != NULL)
      myModifier2->disable();
  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::enableBetaUpdate(){
    if(myModifier1 != NULL)
      myModifier1->enable();
    if(myModifier2 != NULL)
      myModifier2->enable();

  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::saveValues() {
    (*myOldForces)     = (*(next()->getForces()));
    MCIntegrator::saveValues();
  }

  //  ---------------------------------------------------------------------  //

  void ShadowHMCIntegrator::restoreValues(){
    next()->getForces()->intoAssign(*myOldForces);
    MCIntegrator::restoreValues();
  }

}