mici.integrators module#

Symplectic integrators for simulation of Hamiltonian dynamics.

class mici.integrators.BCSSFourStageIntegrator(system, step_size=None)[source]#

Bases: SymmetricCompositionIntegrator

Four-stage symmetric composition integrator due to Blanes, Casas & Sanz-Serna.

Described in equation (6.8) in Blanes, Casas, Sanz-Serna (2014).

Corresponds to specific instance of SymmetricCompositionIntegrator with \(S = 4\) and free coefficients \(a_0 = 0.071353913450279725904\), \(b_1 = 0.191667800000000000000\) and \(a_1 = 0.268548791161230105820\).

References

  1. Blanes, S., Casas, F., & Sanz-Serna, J. M. (2014). Numerical integrators for the Hybrid Monte Carlo method. SIAM Journal on Scientific Computing, 36(4), A1556-A1580.

Parameters:

  • system (TractableFlowSystem) – Hamiltonian system to integrate the dynamics of with tractable Hamiltonian component flows.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.BCSSThreeStageIntegrator(system, step_size=None)[source]#

Bases: SymmetricCompositionIntegrator

Three-stage symmetric composition integrator due to Blanes, Casas & Sanz-Serna.

Described in equation (6.7) in Blanes, Casas, Sanz-Serna (2014).

Corresponds to specific instance of SymmetricCompositionIntegrator with \(S = 3\) and free coefficients \(a_0 = 0.11888010966548\) and \(b_1 = 0.29619504261126\).

References

  1. Blanes, S., Casas, F., & Sanz-Serna, J. M. (2014). Numerical integrators for the Hybrid Monte Carlo method. SIAM Journal on Scientific Computing, 36(4), A1556-A1580.

Parameters:

  • system (TractableFlowSystem) – Hamiltonian system to integrate the dynamics of with tractable Hamiltonian component flows.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.BCSSTwoStageIntegrator(system, step_size=None)[source]#

Bases: SymmetricCompositionIntegrator

Two-stage symmetric composition integrator due to Blanes, Casas & Sanz-Serna.

Described in equation (6.4) in Blanes, Casas, Sanz-Serna (2014).

Corresponds to specific instance of SymmetricCompositionIntegrator with \(S = 2\) and free coefficient \(a_0 = (3 - \sqrt{3}) / 6\).

References

  1. Blanes, S., Casas, F., & Sanz-Serna, J. M. (2014). Numerical integrators for the Hybrid Monte Carlo method. SIAM Journal on Scientific Computing, 36(4), A1556-A1580.

Parameters:

  • system (TractableFlowSystem) – Hamiltonian system to integrate the dynamics of with tractable Hamiltonian component flows.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.ConstrainedLeapfrogIntegrator(system, step_size=None, n_inner_step=1, reverse_check_tol=2e-08, reverse_check_norm=<function maximum_norm>, projection_solver=<function solve_projection_onto_manifold_newton>, projection_solver_kwargs=None)[source]#

Bases: TractableFlowIntegrator

Leapfrog integrator for constrained Hamiltonian systems.

The Hamiltonian function is assumed to be expressible as the sum of two components for which the corresponding (unconstrained) Hamiltonian flows can be exactly simulated. Specifically it is assumed that the Hamiltonian function h takes the form

\[h(q, p) = h_1(q) + h_2(q, p),\]

where \(q\) and \(p\) are the position and momentum variables respectively, and \(h_1\) and \(h_2\) Hamiltonian component functions for which the exact flows, respectively \(\Phi_1\) and \(\Phi_2\), can be computed.

The system is assumed to be additionally subject to a set of holonomic constraints on the position component of the state i.e. that all valid states must satisfy

\[c(q) = 0,\]

for some differentiable and surjective vector constraint function \(c\) and the set of positions satisfying the constraints implicitly defining a manifold the dynamics remain confined to.

The constraints are enforced by introducing a set of Lagrange multipliers \(\lambda\) of dimension equal to number of constraints, and defining a ‘constrained’ Hamiltonian

\[\bar{h}(q, p) = h_1(q) + h_2(q, p) + c(q)^T\lambda ~\text{s.t.}~ c(q) = 0,\]

with corresponding dynamics described by the system of differential algebraic equations

\[\dot{q} = \nabla_2 h(q, p), \quad \dot{p} = -\nabla_1 h(q, p) - \partial c(q)^T\lambda, \quad c(q) = 0.\]

The dynamics implicitly define a set of constraints on the momentum variables, differentiating the constraint equation with respect to time giving that

\[\partial c(q) \nabla_2 h(q, p) = \partial c(q) \nabla_2 h_2(q, p) = 0.\]

The set of momentum variables satisfying the above for given position variables is termed the cotangent space of the manifold (at a position), and the set of position-momentum pairs for which the position is on the constraint manifold and the momentum in the corresponding cotangent space is termed the cotangent bundle.

To define a second-order symmetric integrator which exactly (up to floating point error) preserves these constraints, forming a symplectic map on the cotangent bundle, we follow the approach of Reich (1996).

We first define a map \(\Pi\) parametrised by a vector of Lagrange multipliers \(\lambda\)

\[\Pi(\lambda)(q, p) = (q, p + \partial c(q)^T \lambda),\]

with \(\lambda\) allowed to be an implicitly defined function of \(q\) and \(p\).

We then define a map \(A\) in terms of the \(h_1\) flow map \(\Phi_1\) as

\[A(t) = \Pi(\lambda) \circ \Phi_1(t),\]

with \(\lambda\) implicitly defined such that for \((q', p') = A(t)(q, p)\) we have that \(\partial c(q') \nabla_2 h_2(q', p') = 0\) for any initial state \((q, p)\) in the co-tangent bundle, with \(c(q') = 0\) trivially satisfied as \(\Phi_1\) is an identity map in the position:

\[\Phi_1(t)(q, p) = (q, p - t \nabla h_1(q)).\]

The map \(A(t)\) therefore corresponds to taking an unconstrained step according to the \(h_1\) component flow map \(\Phi_1(t)\) and then projecting the resulting updated momentum back in to the co-tangent space. For the usual case in which \(h\) includes only quadratic terms in the momentum \(p\) such that \(\nabla_2 h(q, p)\) is a linear function of \(p\), then \(\lambda\) can be analytically solved for to give a closed-form expression for the projection into the co-tangent space.

We also define a map \(B\) in terms of the \(h_2\) flow map \(\Phi_2\) as

\[B(t) = \Pi(\lambda') \circ \Phi_2(t) \circ \Pi(\lambda),\]

such that for \((q', p') = B(t)(q, p)\), \(\lambda\) is implicitly defined such that \(c(q') = 0\) and \(\lambda'\) is implicitly defined such that \(\partial c(q') \nabla_2 h(q', p') = 0\).

This can be decomposed as first solving for \(\lambda\) such that

\[c((\Phi_2(t) \circ \Pi(\lambda)(q, p))_1) = c((\Phi_2(t)(q, p + \partial c(q)^T \lambda))_1) = 0,\]

i.e. solving for the values of the Lagrange multipliers such that the position component of the output of \(\Phi_2(t) \circ \Pi(\lambda)\) is on the manifold, with this typically a non-linear system of equations that will need to be solved iteratively e.g. using Newton’s method. The momentum output of \(\Phi_2(t) \circ \Pi(\lambda)\) is then projected in to the cotangent space to compute the final state pair, with this projection step as noted above typically having an analytic solution.

For more details see Reich (1996) and section 7.5.1 in Leimkuhler and Reich (2004).

The overall second-order integrator is then defined as the symmetric composition

\[\Psi(t) = A(t / 2) \circ B(t / N)^N \circ A(t / 2),\]

where \(N\) is a positive integer corresponding to the number of ‘inner’ \(h_2\) flow steps, following the ‘geodesic integrator’ formulation proposed by Leimkuhler and Matthews (2016). The additional flexibility introduced by having the possibility of \(N > 1\) is particularly of use when evaluation of \(\Phi_1\) is significantly more expensive than evaluation of \(\Phi_2\); in this case using \(N > 1\) can allow a larger time step \(t\) to be used than may be otherwise possible due to the need to ensure the iterative solver used in \(B\) does not diverge, with a smaller step size \(t / N\) used for the steps involving the iterative solves with the (cheaper) \(\Phi_2\) flow map and a larger step size \(t\) used for the steps involving the (more expensive) \(\Phi_1\) flow map.

This integrator exactly preserves the constraints at all steps, such that if an initial position momentum pair \((q, p)\) are in the cotangent bundle, the corresponding pair after calling the step() method of the integrator will also be in the cotangent bundle, providing the iterative solver converges.

As the iterative solves may fail to converge, or may converge to one of multiple solutions, following the approach proposed by Zappa, Holmes-Cerfon and Goodman (2018), an explicit reversibility check is performed to ensure the overall integrator step is time-reversible; see also Lelievre, Rousset and Stoltz (2019) for an analysis of this approach specifically in the context of Hamiltonian Monte Carlo. If the reversibility check fails or the iterative solver fails to converge an appropriate error is raised (mici.errors.NonReversibleStepError and mici.errors.ConvergenceError respectively).

References

  1. Reich, S. (1996). Symplectic integration of constrained Hamiltonian systems by composition methods. SIAM journal on numerical analysis, 33(2), 475-491.

  2. Leimkuhler, B., & Reich, S. (2004). Simulating Hamiltonian Dynamics (No. 14). Cambridge University Press.

  3. Leimkuhler, B., & Matthews, C. (2016). Efficient molecular dynamics using geodesic integration and solvent-solute splitting. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472(2189), 20160138.

  4. Zappa, E., Holmes-Cerfon, M., & Goodman, J. (2018). Monte Carlo on manifolds: sampling densities and integrating functions. Communications on Pure and Applied Mathematics, 71(12), 2609-2647.

  5. Lelievre, T., Rousset, M., & Stoltz, G. (2019). Hybrid Monte Carlo methods for sampling probability measures on submanifolds. Numerische Mathematik, 143, 379-421.

Parameters:

  • system (ConstrainedTractableFlowSystem) – Hamiltonian system to integrate the constrained dynamics of.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

  • n_inner_step (int) – Positive integer specifying number of ‘inner’ constrained system.h2_flow steps to take within each overall step. As the derivative system.dh1_dpos is not evaluated during the system.h2_flow steps, if this derivative is relatively expensive to compute compared to evaluating system.h2_flow then compared to using n_inner_step = 1 (the default) for a given step_size it can be more computationally efficient to use n_inner_step > 1 in combination within a larger step_size, thus reducing the number of system.dh1_dpos evaluations to simulate forward a given time while still controlling the effective time step used for the constrained system.h2_flow steps which involve solving a non-linear system of equations to retract the position component of the updated state back on to the manifold, with the iterative solver typically diverging if the time step used is too large.

  • reverse_check_tol (float) – Tolerance for check of reversibility of implicit sub-steps which involve iterative solving of a non-linear system of equations. The step is assumed to be reversible if sequentially applying the forward and adjoint updates to a state returns to a state with a position component within a distance (defined by the reverse_check_norm argument) of reverse_check_tol of the original state position component. If this condition is not met a mici.errors.NonReversibleStepError exception is raised.

  • reverse_check_norm (NormFunction) – Norm function accepting a single one-dimensional array input and returning a non-negative floating point value defining the distance to use in the reversibility check. Defaults to mici.solvers.maximum_norm().

  • projection_solver (ProjectionSolver) –

    Function which given two states state and state_prev, floating point time step time_step and a Hamiltonian system object system solves the non-linear system of equations in λ.

    system.constr(
    state.pos
    + dh2_flow_pos_dmom @ system.jacob_constr(state_prev).T @ λ
) == 0

    where dh2_flow_pos_dmom = system.dh2_flow_dmom(time_step)[0] is the derivative of the action of the (linear) system.h2_flow map on the state momentum component with respect to the position component. This is used to project the state position component back on to the manifold after an unconstrained system.h2_flow update. If the solver fails to convege a mici.errors.ConvergenceError exception is raised.

  • projection_solver_kwargs (Optional[dict[str, Any]]) – Dictionary of any keyword arguments to projection_solver.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.ImplicitLeapfrogIntegrator(system, step_size=None, reverse_check_tol=2e-08, reverse_check_norm=<function maximum_norm>, fixed_point_solver=<function solve_fixed_point_direct>, fixed_point_solver_kwargs=None)[source]#

Bases: Integrator

Implicit leapfrog integrator for Hamiltonians with a non-separable component.

Also known as the generalised leapfrog method.

The Hamiltonian function \(h\) is assumed to take the form

\[h(q, p) = h_1(q) + h_2(q, p),\]

where \(q\) and \(p\) are the position and momentum variables respectively, \(h_1\) is a Hamiltonian component function for which the exact flow can be computed and \(h_2\) is a Hamiltonian component function of the position and momentum variables, which may be non-separable and for which exact simulation of the correspond Hamiltonian flow may not be possible.

The overall integrator step \(\Psi\) is defined by the symmetric composition

\[\Psi(t) = A(t/2) \circ B(t/2) \circ C(t/2) \circ C^*(t/2) \circ B^*(t/2) \circ A^*(t/2),\]

where the adjoint of a flow map \(X\) is defined such that \(X^*(t) = X(-t)^{-1}\) and the component maps are defined by

\[\begin{split}A(t)(q, p) = A^*(t)(q, p) = (q, p - t\nabla h_1(q)), \\ B(t)(q, p) = \lbrace (q, p') : p' = p - t \nabla_1 h_2(q, p') \rbrace, \\ B^*(t)(q, p) = (q, p - t \nabla_1 h_2(q, p)), \\ C(t)(q, p) = (q + t \nabla_2 h_2(q, p), p), \\ C^*(t)(q, p) = \lbrace (q', p) : q' = q + t \nabla_2 h_2(q', p) \rbrace.\end{split}\]

The resulting implicit integrator is a symmetric second-order method corresponding to a symplectic partitioned Runge-Kutta method. See Section 6.3.2 in Leimkuhler and Reich (2004) for more details.

Fixed-point iterations are used to solve the non-linear systems of equations in the implicit component updates (\(B\) and \(C^*\)). As the iterative solves may fail to converge, or may converge to one of multiple solutions, following the approach proposed by Zappa, Holmes-Cerfon and Goodman (2018), an explicit reversibility check is performed to ensure the overall integrator step is time-reversible. If the reversibility check fails or the iterative solver fails to converge an appropriate error is raised (mici.errors.NonReversibleStepError and mici.errors.ConvergenceError respectively).

References

  1. Leimkuhler, B., & Reich, S. (2004). Simulating Hamiltonian Dynamics (No. 14). Cambridge University Press.

  2. Zappa, E., Holmes-Cerfon, M., & Goodman, J. (2018). Monte Carlo on manifolds: sampling densities and integrating functions. Communications on Pure and Applied Mathematics, 71(12), 2609-2647.

Parameters:

  • system (System) – Hamiltonian system to integrate the dynamics of.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

  • reverse_check_tol (float) – Tolerance for check of reversibility of implicit sub-steps which involve iterative solving of a non-linear system of equations. The step is assumed to be reversible if sequentially applying the forward and adjoint updates to a state returns to a state with a position component within a distance (defined by the reverse_check_norm argument) of reverse_check_tol of the original state position component. If this condition is not met a mici.errors.NonReversibleStepError exception is raised.

  • reverse_check_norm (NormFunction) – Norm function accepting a single one-dimensional array input and returning a non-negative floating point value defining the distance to use in the reversibility check. Defaults to mici.solvers.maximum_norm().

  • fixed_point_solver (FixedPointSolver) – Function which given a function func and initial guess x0 iteratively solves the fixed point equation func(x) = x initialising the iteration with x0 and returning an array corresponding to the solution if the iteration converges or raising a mici.errors.ConvergenceError otherwise. Defaults to mici.solvers.solve_fixed_point_direct.

  • fixed_point_solver_kwargs (Optional[dict[str, Any]]) – Dictionary of any keyword arguments to fixed_point_solver.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.ImplicitMidpointIntegrator(system, step_size=None, reverse_check_tol=2e-08, reverse_check_norm=<function maximum_norm>, fixed_point_solver=<function solve_fixed_point_direct>, fixed_point_solver_kwargs=None)[source]#

Bases: Integrator

Implicit midpoint integrator for general Hamiltonians.

The Hamiltonian function \(h\) may be a general (non-separable) function of both the position variables \(q\) and momentum variables \(p\).

The Hamiltonian flow \(\Phi\) is approximated with the symmetric composition \(\Psi(t) = A(t/2) \circ A^*(t/2)\) of an implicit Euler half-step \(A(t/2)\) with an explicit Euler half-step \(A^*(t/2)\) (which is adjoint to the implicit Euler step, that is \(A^*(t) = A(-t)^{-1}\)), with the components maps defined by

\[\begin{split}A(t)(q, p) = \lbrace (q', p') : q' = q + t \nabla_2 h(q', p'), p' = p - t \nabla_1 h(q', p')) \rbrace, \\ A^*(t)(q, p) = (q + t \nabla_2 h(q, p), p - t \nabla_1 h(q, p)).\end{split}\]

The resulting implicit integrator is a second-order method corresponding to a symplectic one-stage Runge-Kutta method. See Sections 4.1 and 6.3.1 in Leimkuhler and Reich (2004) for more details.

A fixed-point iteration is used to solve the non-linear system of equations in the implicit Euler step \(A\). As the iterative solve may fail to converge, or may converge to one of multiple solutions, following the approach proposed by Zappa, Holmes-Cerfon and Goodman (2018), an explicit reversibility check is performed to ensure the overall integrator step is time-reversible. If the reversibility check fails or the iterative solver fails to converge an appropriate error is raised (mici.errors.NonReversibleStepError and mici.errors.ConvergenceError respectively).

References

  1. Leimkuhler, B., & Reich, S. (2004). Simulating Hamiltonian Dynamics (No. 14). Cambridge University Press.

  2. Zappa, E., Holmes-Cerfon, M., & Goodman, J. (2018). Monte Carlo on manifolds: sampling densities and integrating functions. Communications on Pure and Applied Mathematics, 71(12), 2609-2647.

Parameters:

  • system (System) – Hamiltonian system to integrate the dynamics of.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step method.

  • reverse_check_tol (float) – Tolerance for check of reversibility of implicit Euler steps which involve iterative solving of a non-linear system of equations. The step is assumed to be reversible if sequentially applying the forward and adjoint updates to a state returns to a state with a position component within a distance (defined by the reverse_check_norm argument) of reverse_check_tol of the original state position component. If this condition is not met a mici.errors.NonReversibleStepError exception is raised.

  • reverse_check_norm (NormFunction) – Norm function accepting a single one-dimensional array input and returning a non-negative floating point value defining the distance to use in the reversibility check. Defaults to mici.solvers.maximum_norm().

  • fixed_point_solver (FixedPointSolver) – Function which given a function func and initial guess x0 iteratively solves the fixed point equation func(x) = x initialising the iteration with x0 and returning an array corresponding to the solution if the iteration converges or raising a mici.errors.ConvergenceError otherwise. Defaults to mici.solvers.solve_fixed_point_direct().

  • fixed_point_solver_kwargs (Optional[dict[str, Any]]) – Dictionary of any keyword arguments to fixed_point_solver.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.Integrator(system, step_size=None)[source]#

Bases: ABC

Base class for integrators for simulating Hamiltonian dynamics.

For a Hamiltonian function \(h\) with position variables \(q\) and momentum variables \(p\), the canonical Hamiltonian dynamic is defined by the ordinary differential equation system

\[\dot{q} = \nabla_2 h(q, p), \qquad \dot{p} = -\nabla_1 h(q, p),\]

with the flow map \(\Phi\) corresponding to the solution of the corresponding initial value problem a time-reversible and symplectic (and by consequence volume-preserving) map.

Derived classes implement a step() method which approximates the flow-map with \(\Psi(t) \approx \Phi(t)\) over some small time interval \(t\), while conserving the properties of being time-reversible and symplectic, with composition of this integrator step method allowing simulation of time-discretised trajectories of the Hamiltonian dynamics.

Parameters:

  • system (System) – Hamiltonian system to integrate the dynamics of.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

step(state)[source]#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.LeapfrogIntegrator(system, step_size=None)[source]#

Bases: TractableFlowIntegrator

Leapfrog integrator for Hamiltonian systems with tractable component flows.

The overall integrator step \(\Psi\) is defined by the symmetric composition

\[\Psi(t) = \Phi_1(t/2) \circ \Phi_2(t) \circ \Phi_1(t/2)\]

where \(\Phi_1\) and \(\Phi_2\) are the exact flow maps associated with the Hamiltonian components \(h_1\) and \(h_2\) respectively.

For separable Hamiltonians of the

\[h(q, p) = h_1(q) + h_2(p),\]

where \(h_1\) is the potential energy and \(h_2\) is the kinetic energy, this integrator corresponds to the classic (position) Störmer-Verlet method.

The integrator can also be applied to the more general Hamiltonian splitting

\[h(q, p) = h_1(q) + h_2(q, p),\]

providing the flows for \(h_1\) and \(h_2\) are both tractable.

For more details see Sections 2.6 and 4.2.2 in Leimkuhler and Reich (2004).

References

  1. Leimkuhler, B., & Reich, S. (2004). Simulating Hamiltonian Dynamics (No. 14). Cambridge University Press.

Parameters:

  • system (TractableFlowSystem) – Hamiltonian system to integrate the dynamics of with tractable Hamiltonian component flows.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.SymmetricCompositionIntegrator(system, free_coefficients, *, step_size=None, initial_h1_flow_step=True)[source]#

Bases: TractableFlowIntegrator

Symmetric composition integrator for Hamiltonians with tractable flows.

The Hamiltonian function is assumed to be expressible as the sum of two analytically tractable components for which the corresponding Hamiltonian flows can be exactly simulated. Specifically it is assumed that the Hamiltonian function \(h\) takes the form

\[h(q, p) = h_1(q) + h_2(q, p),\]

where \(q\) and \(p\) are the position and momentum variables respectively, and \(h_1\) and \(h_2\) are Hamiltonian component functions for which the exact flows, respectively \(\Phi_1\) and \(\Phi_2\), can be computed. An alternating composition can then be formed as

\[\Psi(t) = A(a_S t) \circ B(b_S t) \circ \dots \circ A(a_1 t) \circ B(b_1 t) \circ A(a_0 t),\]

where \(A = \Phi_1\) and \(B = \Phi_2\) or \(A = \Phi_2\) and \(B = \Phi_1\), and \((a_0,\dots, a_S)\) and \((b_1, \dots, b_S)\) are a set of coefficients to be determined with \(S \geq 1\).

To ensure a consistency (i.e. the integrator is at least order one) we require that

\[\sum_{s=0}^S a_s = \sum_{s=1}^S b_s = 1.\]

For symmetric compositions we restrict that

\[a_{S-m} = a_m, \quad b_{S+1-m} = b_m,\]

with symmetric consistent methods of at least order two.

The combination of the symmetry and consistency requirements mean that for each \(S \geq 1\) a symmetric composition method can be described by \(S - 1\) ‘free’ coefficients \((a_0, b_1, \dots, a_{K-1}, b_K)\) with \(K = (S - 1) / 2\) if \(S > 1\) is odd (with no free coefficients for \(S = 1\) case) or \((a_0, b_1, \dots, a_K)\) with \(K = (S - 2) / 2\) if \(S > 2\) is even (with a single free coefficient \(a_0\) for \(S=2\) case).

The Störmer-Verlet ‘leapfrog’ integrator is the special case corresponding to the unique (symmetric and consistent) ‘1-stage’ (\(S = 1\)) integrator.

For more details see Section 6.2 in Leimkuhler and Reich (2004).

References

  1. Leimkuhler, B., & Reich, S. (2004). Simulating Hamiltonian Dynamics (No. 14). Cambridge University Press.

Parameters:

  • system (TractableFlowSystem) – Hamiltonian system to integrate the dynamics of with tractable Hamiltonian component flows.

  • free_coefficients (Sequence[float]) – Sequence of \(S - 1\) scalar values, where \(S\) is the number of stages in the symmetric composition, specifying the free coefficients \((a_0, b_1, a_1, \dots, a_K, b_K)\) with \(K = (S - 1) / 2\) if \(S\) is odd or \((a_0, b_1, a_1, \dots, a_K)\) with \(k = (S - 2) / 2\) if S is even.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

  • initial_h1_flow_step (bool) – Whether the initial \(A\) flow in the composition should correspond to the flow of the h_1 Hamiltonian component (True) or to the flow of the \(h_2\) component (False).

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState

class mici.integrators.TractableFlowIntegrator(system, step_size=None)[source]#

Bases: Integrator

Base class for integrators for Hamiltonian systems with tractable flows.

The Hamiltonian function is assumed to be expressible as the sum of two analytically tractable components for which the corresponding Hamiltonian flows can be exactly simulated. Specifically it is assumed that the Hamiltonian function \(h\) takes the form

\[h(q, p) = h_1(q) + h_2(q, p),\]

where \(q\) and \(p\) are the position and momentum variables respectively, and \(h_1\) and \(h_2\) are Hamiltonian component functions for which the exact flow maps, \(\Phi_1\) and \(\Phi_2\) respectively, can be computed exactly.

Parameters:

  • system (TractableFlowSystem) – Hamiltonian system to integrate the dynamics of with tractable Hamiltonian component flows.

  • step_size (Optional[float]) – Integrator time step. If set to None it is assumed that a step size adapter will be used to set the step size before calling the step() method.

step(state)#

Perform a single integrator step from a supplied state.

Parameters:

state (ChainState) – System state to perform integrator step from.

Returns:

New object corresponding to stepped state.

Return type:

ChainState