SHARE
TWEET

Integrator sketch 1.0

a guest Dec 19th, 2015 15 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /* This code sketches a State object, some example integrator objects
  2.  * that own a State (and other pieces of machinery), and use those to
  3.  * implement a wide variety of simulation algorithms where each
  4.  * conforms to the IIntegrator interface. There's copious details
  5.  * missing, multiple TODOs noted, some notes about optimization
  6.  * opportunities, and ways we will later implement stuff with task
  7.  * parallelism. I haven't tried to implement Michael's collectives and
  8.  * elements model, but we can reconsider that if we identify something
  9.  * here that we don't like.
  10.  *
  11.  * The plan here is to have a very flat IIntegrator class hierarchy.
  12.  * See https://en.wikipedia.org/wiki/Composition_over_inheritance (and
  13.  * other places) why this is a good idea. If two integrator classes
  14.  * end up looking like they have enough code duplication to be a
  15.  * problem, then we extract that responsibility to an object of a new
  16.  * class, and each of the integrator classes owns an instance of that
  17.  * object. That is, we don't make sub-classes in order to re-use code
  18.  * from the super-class.
  19.  *
  20.  * Some details have been filled in for a subclass of IIntegrator that
  21.  * does Bussi-style velocity-Verlet unconstrained dynamics. For
  22.  * example, it contains Thermostats and IVelocityUpdater objects that
  23.  * are set up in the factory function to do the right thing according
  24.  * to the inputrec. Some of those objects will get re-used for
  25.  * e.g. Bussi-style leapfrog dynamics, which will also be an immediate
  26.  * subclass of IIntegrator. For example, adding support for
  27.  * constraints could be as straightforward as configuring the object
  28.  * with a NullConstrainer, ShakeConstrainer, or LincsConstrainer,
  29.  * according to the .tpr contents. If we then measured too much
  30.  * overhead for .tprs we care about, we can specialize the
  31.  * integrators accordingly.
  32.  *
  33.  * Once in the MD loop, the integrator just keeps on doing the right
  34.  * thing that it was built to do, rather than continually looking up
  35.  * in the inputrec what it should be doing. This requires a handful of
  36.  * virtual function calls, but passes fewer function parameters and
  37.  * checks fewer flags. (Once we are constructing task DAGs, such
  38.  * concerns move to an entirely different scope, anyway.)
  39.  *
  40.  * Names should be clear and descriptive. If a Chinese grad student
  41.  * can't tell the role of your identifier from the name, then you
  42.  * named it wrongly. :-) We need to write code to be conveniently
  43.  * read, not to be conveniently typed.
  44.  *
  45.  * I haven't sketched how one will implement a multi-stage or MC-MD
  46.  * integrator, but I think the required principles follow from what
  47.  * Michael has planned and I have sketched here.
  48.  */
  49. #include "config.h"
  50.  
  51. #include "gromacs/math/vectypes.h"
  52.  
  53. namespace gmx
  54. {
  55.  
  56. //! Permit zero-cost run-time checks active only for debug builds
  57. const bool isDebugMode = IS_DEBUG_BUILD_TYPE;
  58.  
  59. class StateVector
  60. {
  61.     public:
  62.         // Vector-like interface, including const and non-const
  63.         // data(), and size().
  64.         void setValid(bool validity) { isValid_ = validity; }
  65.         bool valid() const { return isValid_; }
  66.     private:
  67.         bool isValid_;
  68.         // Later, we decide how to minimize any overhead here
  69.         RVec *velocities_;
  70.         rvec *rawVelocities_; // legacy
  71.         int rawVelocitiesAllocSize_; // legacy
  72. };
  73.  
  74. // === State ===
  75.  
  76. class State
  77. {
  78.     public:
  79.         const StateVector *velocities() const
  80.         {
  81.             if (isDebugMode)
  82.             {
  83.                 GMX_RELEASE_ASSERT(velocities.valid(), "Can't read invalid velocities (because ...)");
  84.             }
  85.             return &velocities_;
  86.         }
  87.         // Maybe this should actually return a view object that
  88.         // e.g. can't be re-sized, or even has a destructor that calls
  89.         // validateVelocities()? That means other objects don't have
  90.         // to depend on the main interface of State, but only the view
  91.         // that they need to use.
  92.         StateVector *velocities()
  93.         {
  94.             if (isDebugMode)
  95.             {
  96.                 GMX_RELEASE_ASSERT(velocities.valid(), "Can't modify invalid velocities (because ...)");
  97.                 velocitiesAreValid_ = false;
  98.                 // This forces people to be const correct if they just
  99.                 // want to write code that inspects
  100.                 // velocities. They'll get a run-time error from the
  101.                 // next velocity update if they've used the non-const
  102.                 // getter where they intended read-only access.
  103.             }
  104.             return &velocities_;
  105.         }
  106.         void validateVelocities()
  107.         {
  108.             if (isDebugMode)
  109.             {
  110.                 velocities.setValid(true);
  111.             }
  112.         }
  113.         /* Similar stuff for forces, positions and box (at
  114.            least). Probably virial. Probably not any thermostat or
  115.            barostat stuff - the need for those is particular to
  116.            certain integrators. */
  117.     private:
  118.         // Later, we might extract the concept of validity to FancyRVec type
  119.         StateVector velocities_;
  120. };
  121.  
  122. // === Thermostats ===
  123.  
  124. class IThermostat;
  125.  
  126. IThermostat *buildThermostatForCouplingGroup(t_inputrec ir);
  127.  
  128. //! This is the interface supported by all thermostats
  129. class IThermostat
  130. {
  131.         friend buildThermostatForCouplingGroup;
  132.     public:
  133.         IThermostat(int indexOfCouplingGroup) :
  134.             indexOfCouplingGroup_(indexOfCouplingGroup) {}
  135.         virtual void updateVelocities(real timeIncrement, gmx_int64_t *step, gmx_ekindata_t *ekind, StateVector *velocities) = 0;
  136.     private:
  137.         int indexOfCouplingGroup_; // legacy - needed for using gmx_ekindata_t in current form
  138. };
  139.  
  140. //! Leaves velocity unmodified, e.g. NVE or testing code
  141. class NullThermostat : public IThermostat
  142. {
  143.         // TODO assume single coupling group for now
  144.         friend buildThermostatForCouplingGroup;
  145.     public:
  146.         NullThermostat(int indexOfCouplingGroup) :
  147.             IThermostat(indexOfCouplingGroup) {}
  148.         virtual void updateVelocities(real timeIncrement, gmx_int64_t *step, gmx_ekindata_t *ekind, StateVector *velocities)
  149.         {
  150.             // do nothing, not even re-validate velocities
  151.  
  152.             // (Much) later, calling such functions will spawn tasks
  153.             // into a task DAG, and the job of this function call will
  154.             // be to connect the parent directly to the child, i.e
  155.             // zero overhead when the DAG is later executed.
  156.         }
  157. };
  158.    
  159. /*! \brief Specializes IThermostat for thermostats that do a Bussi-style resample of KE
  160.  *
  161.  * TODO Decide if this approach is reasonable. This extends the
  162.  * IThermostat interface to form a new interface for Bussi
  163.  * thermostats. I'm not sure we need more than one implementation of
  164.  * it, however, if we fix the way ekind stores KE. If we do need
  165.  * multiple implementations, then it would also be reasonable to
  166.  * instead have a KineticEnergyResampler class, and have the different
  167.  * Bussi thermostats subclass IThermostat directory, and call a method
  168.  * on the KineticEnergyResampler instance it contains. */
  169. class IBussiThermostat : public IThermostat
  170. {
  171.         friend buildThermostatForCouplingGroup;
  172.     public:
  173.         IBussiThermostat(int indexOfCouplingGroup /* and other settings from inputrec */) :
  174.             IThermostat(indexOfCouplingGroup) /* and other initializers */ {}
  175.         virtual void updateVelocities(real timeIncrement, gmx_int64_t *step, gmx_ekindata_t *ekind, StateVector *velocities) = 0;
  176.     private:
  177.         real resampleKineticEnergy(real currentKineticEnergy)
  178.         {
  179.             // do actual work, calling RNG, etc.
  180.         }
  181.        
  182.         real tau_;
  183.         int numDegreesOfFreedom_;
  184.         real targetAverageKineticEnergy_;
  185.         // perhaps other pre-computed quantities here
  186. };
  187.  
  188. //! Does Bussi v-rescale thermostat for integrators with on-step KE, e.g VV
  189. class OnStepKineticEnergyBussiThermostat : public IBussiThermostat
  190. {
  191.         friend buildThermostatForCouplingGroup;
  192.     public:
  193.         OnStepKineticEnergyBussiThermostat(int indexOfCouplingGroup /* and other settings from inputrec */) :
  194.             IBussiThermostat(indexOfCouplingGroup /* and other initializers */) {}
  195.         virtual void updateVelocities(real timeIncrement, gmx_int64_t *step, gmx_ekindata_t *ekind, StateVector *velocities)
  196.         {
  197.             // NB no conditional on the integrator or tau or anything
  198.             // - if we get here we just do the right work for this
  199.             // case.
  200.             real newKineticEnergy = IBussiThermostat::resampleKineticEnergy(whateverFunction(ekind->tcstat[indexOfCouplingGroup_].ekinf));
  201.             // Here, do actual rescale of velocity components
  202.         }
  203.     private:
  204. };
  205.    
  206. //! Factory function
  207. IThermostat *buildThermostatForCouplingGroup(t_inputrec *ir, int i)
  208. {
  209.     GMX_ASSERT(ir->etc != etcNO, "Only real thermostats get here");
  210.    
  211.     if (ir->opts->tau_t[i] == 0 || ir->opts->ndrf[i] == 0)
  212.     {
  213.         // There will never be anything to do for this group, so don't do anything!
  214.         return new NullThermostat();
  215.     }
  216.     if (ir->etc == etcVRESCALE)
  217.     {
  218.         if (EI_VV(ir->eI))
  219.         {
  220.             return new OnStepKineticEnergyBussiThermostat(i /* and other stuff from inputrec */);
  221.         }
  222.         else
  223.         {
  224.             // ...
  225.         }
  226.         // Also, bump a counter that will schedule a please_cite()
  227.     }
  228.     else
  229.     {
  230.         // ...
  231.     }
  232. }
  233.  
  234. // TODO use a container that knows it has to call delete on its
  235. // contents (needs to contain pointers to IThermostat for dynamic
  236. // dispatch to work)
  237. typedef std::vector<IThermostat *> Thermostats;
  238.  
  239. //! Factory function
  240. Thermostats buildThermostats(t_inputrec *ir)
  241. {
  242.     std::vector<Thermostat *> thermostats;
  243.     if (ir->etc != etcNO)
  244.     {
  245.         for (int i = 0; i < ir->opts->ngtc; i++)
  246.         {
  247.             thermostats.push_back(buildThermostatForCouplingGroup(ir, i));
  248.         }
  249.     }
  250. }
  251.  
  252. // === Velocity Updaters ===
  253.  
  254. class IVelocityUpdater
  255. {
  256.     public:
  257.         IVelocityUpdater() {};
  258.         virtual ~IVelocityUpdater() {};
  259.         virtual void updater(real timeStep, const StateVector *forces, StateVector *velocities) = 0;
  260. };
  261.  
  262. class NormalVelocityUpdater : public IVelocityUpdater
  263. {
  264.     public:
  265.         NormalVelocityUpdater() {}
  266.         virtual ~NormalVelocityUpdater() {};
  267.         virtual void updater(real timeStep, const StateVector *forces, StateVector *velocities)
  268.         {
  269.             // loop over atoms i
  270.             {
  271.                 // loop over dimensions d
  272.                 {
  273.                     velocities[i][d] += timeStep * forces[i][d] * massFactor_[i];
  274.                 }
  275.             }
  276.         }
  277.     private:
  278.         std::vector<real> massFactor_; // 1/(2 * mass)
  279. };
  280.  
  281. // Other specializations of IVelocityUpdater do freeze groups,
  282. // acceleration, shells, vsites, etc. We also include a fallback "do
  283. // everything supported" implementation, with all the conditionals
  284. // that that requires, and we will unit test that the specializations
  285. // do the same as the fallback.
  286.  
  287. // Position updaters will follow similar style, but are probably simpler.
  288.  
  289. // === Integrators ===
  290.  
  291. class IIntegrator
  292. {
  293.     public:
  294.         IIntegrator(State *state, ForceCalculator *forceCalculator) :
  295.             state_(state), forceCalculator_(forceCalculator) {};
  296.         virtual ~IIntegrator() {};
  297.         virtual void doStep() = 0;
  298.         void doSimulation();
  299.  
  300.     protected:
  301.         State state_;
  302.         //! Configurable object that will compute forces (e.g. normal or shell)
  303.         ForceCalculator *forceCalculator_;
  304.         gmx_int64_t step_;
  305. };
  306.  
  307. /* TODO Later, we will probably specialize a few forms of doStep(),
  308.  * and have doSimulation construct a DAG of e.g. every task to do
  309.  * between successive pair-search steps. */
  310. void IIntegrator::doSimulation()
  311. {
  312.     // Any remaining setup?
  313.     while (true)
  314.     {
  315.         if (/* time to renew pair list */)
  316.         {
  317.             forceCalculator_.updatePairList(state_.positions());
  318.         }
  319.         doStep();
  320.         /* Michael's planning document doesn't address where in the
  321.          * loops we do Hamiltonian updates, but I think this is
  322.          * right. */
  323.         if (/* should update Hamiltonian */)
  324.         {
  325.             forceCalculator_.updateHamiltonian(); // e.g. do replica exchange
  326.         }
  327.         // Handle house-keeping, such as global signals and checkpoint
  328.         // checks
  329.         // ...
  330.         if (checkForLastStep())
  331.         {
  332.             break;
  333.         }
  334.         ++step_;
  335.     }
  336.     writeFinalOutput(step_);
  337. }
  338.  
  339. class GenericIntegrator : public IIntegrator
  340. {
  341.     public:
  342.         GenericIntegrator(State *state, ForceCalculator *forceCalculator, gmx_ekindata_t *ekind, real timeStep) :
  343.             IIntegrator(state, forceCalculator), ekind_(*ekind),
  344.             timeStep_(ir->delta_t), halfTimeStep_(0.5*timeStep_)
  345.             /* probably other stuff here */ {}
  346.         virtual ~GenericIntegrator() {};
  347.         virtual void doStep()
  348.         {
  349.             // This will look much like the current do_md loop, except
  350.             // for the bits handled by doSimulation(). We implement
  351.             // this first (so that we will have code that we can ship
  352.             // that does all the old things). In some cases, this code
  353.             // will be a reference for versions that are specialized
  354.             // by either the kind of integrator or the kind of MD step
  355.             // (e.g. the majority of steps do nothing fancy).
  356.         }
  357.     private:
  358.         // Probably stuff here like generic velocity and position
  359.         // updaters, and initially some of the crud from do_md.
  360. };
  361.  
  362. /*! \brief Implements e.g. Bussi-Parrinello dynamics with unconstrained velocity Verlet
  363.  *
  364.  * ie. B2A|B!X (where X might be the rescaling of KE as in
  365.  * Bussi-Parrinello). */
  366. class ThermostattedVelocityVerletIntegrator : public IIntegrator
  367. {
  368.     public:
  369.         friend buildIntegrator;
  370.  
  371.         ThermostattedVelocityVerletIntegrator(State *state, ForceCalculator *forceCalculator, gmx_ekindata_t *ekind, real timeStep) :
  372.             IIntegrator(state, forceCalculator), ekind_(*ekind),
  373.             timeStep_(timeStep), halfTimeStep_(0.5*timeStep),
  374.             computeGlobalsFlags(CGLO_ENERGY),
  375.             positionUpdater_(nullptr), velocityUpdater_(nullptr), thermostats_() {}
  376.         virtual ~ThermostattedVelocityVerletIntegrator() {};
  377.         virtual void doStep()
  378.         {
  379.             /* Initial conditions: forces are current for x, v is at
  380.                the same time. */
  381.             // TODO As an optimization, we might fuse these first two update stages
  382.             velocityUpdater_->update(halfTimeStep, state_.velocities());
  383.             state_->validateVelocities();
  384.             positionUpdater_->update(timeStep, state_.velocities(), state_.positions());
  385.             state_->validatePositions();
  386.             forceCalculator_->calculate(state_.positions(), state_.forces());
  387.             state_->validateForces();
  388.             velocityUpdater_->update(halfTimeStep, state_.velocities());
  389.             state_->validateVelocities();
  390.             if (/* should communicate */)
  391.             {
  392.                 compute_globals(/* all the arguments */, computeGlobalsFlags | flagsForThisStep);
  393.             }
  394.             if (/* should write trajectories */)
  395.             {
  396.                 writeTrajectory(/* pass required state variables as const */); // Like now, write x on step, v trailing by half a step
  397.             }
  398.             if (/* should write energies */)
  399.             {
  400.                 writeEnergies(); // including dhdl, pull and such
  401.             }
  402.             if (/* should update thermostats */)
  403.             {
  404.                 StateVector *velocities = state_.velocities();
  405.                 for(auto thermostat : Thermostats)
  406.                 {
  407.                     thermostat->updateVelocities(step_, &ekind_, velocities);
  408.                 }
  409.                 state_->validateVelocities();
  410.             }
  411.         }
  412.        
  413.     private:
  414.         gmx_ekindata_t ekind_;
  415.         real timeStep_;
  416.         real halfTimeStep_;
  417.         int computeGlobalsFlags_; // plus any other junk to make compute_globals work. Or do we want a GlobalsComputer object?
  418.         //! Configurable object that will handle dynamical position updates
  419.         IPositionUpdater *positionUpdater_;
  420.         //! Configurable object that will handle dynamical velocity updates
  421.         IVelocityUpdater *velocityUpdater_;
  422.         //! Configurable object that will handle thermostatting
  423.         Thermostats thermostats_;
  424. };
  425.  
  426. //! Factory function for integrators. We pass in stuff read from .tpr and/or .cpt.
  427. IIntegrator *buildIntegrator(t_inputrec *ir, State *state, gmx_ekindata_t *ekind, ForceCalculator *forceCalculator)
  428. {
  429.     IIntegrator *integrator;
  430.  
  431.     /* Make sure that if we don't finish this for GROMACS 2016
  432.      * release, users can't get broken things. This avoids needing to
  433.      * fork the code base and/or have it sit broken for ages. */
  434.     if (!getenv("GMX_USE_NEW_INTEGRATOR_FORMALISM"))
  435.     {
  436.         integrator = new GenericIntegrator(state, forceCalculator, ekind, ir);
  437.     }
  438.     /* This logic is going to get extremely long and complex, but the
  439.      * benefits are that
  440.      * - we currently don't have any good place for mdrun to do run-time
  441.      *   sanity checks for supported algorithm combinations,
  442.      * - this code only runs at setup time,
  443.      * - we will probably have fewer things that can go wrong at run
  444.      *   time,
  445.      * - it's harder to inadvertently support a broken combination of
  446.      *   stuff, and
  447.      * - in our optimized code paths we end up calling a handful of
  448.      *   virtual function calls per step, rather than evaluate scores
  449.      *   of conditionals, pass lots of function arguments on the stack
  450.      *   including many possibly unused ones, and have global
  451.      *   dependencies on stuff like inputrec. */
  452.     else if (EI_VV(ir->eI) && ir->etc != etcNO)
  453.     {
  454.         // Pass in the standard components shared by all integrators
  455.         ThermostattedVelocityVerletIntegrator thisIntegrator =
  456.             new ThermostattedVelocityVerletIntegrator(state, forceCalculator, ekind, ir->delta_t);
  457.         // Build the components that are composed to form the integrator
  458.         thisIntegrator->positionUpdater_ = buildPositionUpdater(ir);
  459.         thisIntegrator->velocityUpdater_ = buildVelocityUpdater(ir);
  460.         thisIntegrator->thermostats_ = buildThermostats(ir);
  461.         integrator = thisIntegrator;
  462.     }
  463.     else if(/* whatever */)
  464.     {
  465.         // ...
  466.     }
  467.     else /* fallback */
  468.     {
  469.         integrator = new GenericIntegrator(state, forceCalculator, ekind, ir);
  470.     }
  471.     return integrator;
  472. }
  473.  
  474. // === High-level calling code ===
  475.  
  476. void mdrunner()
  477. {
  478.     IIntegrator *integrator;
  479.     // Do hardware detection and thread launch and stuff wherever makes sense
  480.     // ...
  481.     {
  482.         State *state, *globalState;
  483.         gmx_ekindata_t *ekind;
  484.         t_inputrec *ir;
  485.         // Read all that stuff from .tpr and/or .cpt files.
  486.         // Do any pre-processing such that we have a consistent state from all starting paths.
  487.         ForceCalculator forceCalculator(/* e.g. whatever init_forcerec gets now */); // Stuff like replica exchange coordinated here
  488.         // ...
  489.         state = DomainDecomposition.doPartition(globalState);
  490.         forceCalculator.searchForPairs(state.positions()); // need initial pair list
  491.         forceCalculator.calculate(state.positions(), state.forces()); // need initial forces
  492.         integrator = buildIntegrator(ir, state, ekind, forceCalculator);
  493.         // Various initialization-time stuff now goes out of scope?
  494.     }
  495.     if (/* rerun */)
  496.     {
  497.         callRerunCode();
  498.     }
  499.     else
  500.     {
  501.         /* Call whatever "integrator," just like now. So we'll also
  502.          * wrap/adapt the EM and TPI code so they fit this
  503.          * formalism. No more crusty switch statement. */
  504.         integrator->doSimulation();
  505.     }
  506.     // Clean up
  507. }
  508.  
  509. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top