Optimizing probability of barrier crossing with differentiable simulators

Martin Šı́pka 1 2 3 Johannes C. B. Dietschreit 1 Michal Pavelka 2 Lukáš Grajciar 3 Rafael Gómez-Bombarelli 1

Abstract

Finding reaction mechanisms is commonly divided into two
sub-tasks. First, a dimensionality reduction from all DoFs
down to the so-called collective variables (CVs) and second,enhanced sampling along those CVs (Torrie et al., 1977;
Darve & Pohorille, 2001; Laio & Parrinello, 2002; Abrams& Bussi, 2013; Spiwok et al., 2015; Valsson et al., 2016).
Recently, machine learning has partially automatized thetask of CV identification. (Sultan & Pande, 2018; Mendels
et al., 2018; Wehmeyer & Noé, 2018; Wang et al., 2019;
Bonati et al., 2020; Wang & Tiwary, 2021; Sun et al., 2022;Šı́pka et al., 2022). The ideal CVs represent those DoFs,
which fully describe the rare transition event, and thus areassociated with the slowest motions. However, identification of the ideal CVs requires knowledge of the transition
path that one is trying to discover, i.e., one still ends up withthe proverbial ”chicken-and-egg problem” (Rohrdanz et al.,2013). A problem not solved by previously proposed machine learning based tools. Furthermore, once the CVs arechosen it is very difficult to correct them on-the-fly. Onceequipped with a good, low-dimensional representation of
the chemical reaction, enhanced sampling algorithms usually introduce a biasing potential such that the occurrenceof reactive events is increased significantly. Subsequent
analysis can yield the reaction mechanism, reaction barrier,and rate constant.

Simulating events that involve some energy barrier often requires us to promote the barrier crossing in order to increase the probability of the event.
One example of such a system can be a chemical
reaction which we propose to explore using differentiable simulations. Transition path discovery
and estimation of the reaction barrier are merged
into a single end-to-end problem that is solved
by path-integral optimization. We show how the
probability of transition can be formulated in adifferentiable way and increase it by introducing
a trainable position-dependent bias function. We
also introduce improvements over standard methods making DiffSim training stable and efficient.

1. Introduction
A chemical reaction is a transition from one depression(reactant, R) to another (product, P) on the potential energy surface (PES), this includes but is not limited to therearrangement of chemical bonds, the exchange between
conformers, and phase transitions. The reaction mechanism
is defined by the most likely paths connecting the two states.
Determining the reaction path is challenging due to the highdimensionality of molecular configuration space, which canhave thousands of degrees of freedom (DoF). Standard unbiased sampling algorithms like molecular dynamics (MD)or Monte-Carlo (MC) often get trapped in stable regions,making it inefficient to explore candidate paths that haveto cross low probability regions in phase space (Chipot &Pohorille, 2007; Chipot, 2014). Hence, it is necessary tobias the exploration.

Differentiable Simulations (DiffSims) have been developed
for optimization, control, and learning of motion (Degrave
et al., 2016; de Avila Belbute-Peres et al., 2018; Hu et al.,2019; 2020), but also for quantities of interest in molecular dynamics (Wang et al., 2020; Ingraham et al., 2019;
Greener & Jones, 2021). Differentiating along trajectories
follows naturally from optimizing path-dependent quantities.
DiffSim results are often promising, but it is well known(Metz et al., 2021) that naı̈vely backpropagated gradientsmay vanish or explode. This problem is associated withthe spectrum of the system’s Jacobian (Metz et al., 2021;
Galimberti et al., 2021). Therefore, in order to use DiffSim
successfully, gradients have to become controllable. Here,the loss gradient behaviour is thoroughly studied, and mechanisms to control fluctuations and magnitude are proposed.
Employing the improved DiffSim, we define a differentiable
loss function that, when minimized, results in the robust
training of a biasing potential, which enhances the samplingof reactive transitions without prior determination of CVs.

1
Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139, USA 2 Mathematical Institute, Faculty of Mathematics andPhysics, Charles University, Sokolovská 83, 186 75 Prague, Czech
Republic 3 Department of Physical and Macromolecular 853 Chemistry, Faculty of Sciences, Charles University, 128 43 854 Prague
2, Czech Republic. Correspondence to: Rafael Gómez-Bombarelli<rafagb@mit.edu>.

Published at the Differentiable Almost Everything Workshop of the
40 th International Conference on Machine Learning, Honolulu,
Hawaii, USA. July 2023. Copyright 2023 by the author(s).

1DiffSim for Rare Events

2. Problem Definition

to the target. Other more complex metrics can be chosen.
Minimizing this loss function leads to a maximization of
the probability (4) and can be seen as the minimization of apath-dependent integral. Note that the loss function is defined only for one point of the trajectory and is influenced by
the dynamics of every preceding point. In this manuscript,
we only consider transitions between two wells, W−1 andW1 . Additional basins would be handled analogously.

Reaction processes on an atomistic level are explored using
molecular dynamics. Let U0 (x) be the system’s PES, where
the column vector x ∈ RN denotes the mass-weighted
coordinates of the system and p the conjugate momenta.
The PES is modified by a learnable bias term B(x, θ) to
increase the probability of the barrier crossing
U (x, θ) = U0 (x) + B(x, θ).

(1)2.2. Differentiable simulations

Commonly, the system is coupled to a thermostat to generate
trajectories from the canonical (NVT) ensemble (Callen &Scott, 1998). Here, we choose the Langevin thermostat
because of its simplicity and its favorable properties withrespect to differentiating along the computational graph, as
will be shown later (see Section 3.2). The thermostat iscoupled to the system through the friction constant γ (2).
Hence, the biased dynamics evolve according to

The loss in (5) is minimized by optimizing the parameters
θ of the bias potential (1). Instead of saving the computational graph for the entire forward path, as it would be
extremely memory-demanding, the optimization process
can be conveniently reformulated using the adjoint equation.
With the adjoint vectors the system dynamics can be run
backwards, leading to memory-savings and the ability to adjust the extent of backpropagation based on the sought-after
dynamical scale (see Section 3). We employ the frameworkand notation adapted recently for neural networks (Chenet al., 2018) from the original work by (Lev Semenovich
Pontryagin et al., 1962).

ẋ(t) = p(t)ṗ(t) = −p
∂U (x(t))− γp(t) + 2γkB T R(t),∂x(2)

where kb is the Boltzmann constant, T the absolute temperature of the bath, and R(t) a Gaussian process.

Propagating z(t) = (x(t), p(t)) using f (z(t), θ) = ż(t)(the right side of (2)) is called the forward process. Notice that f (z(t)) is only time-dependent through z(t). The
adjoint vectors are

2.1. Differentiable barrier crossing
It is often suitable to use general curvilinear coordinate todescribe reactions. These special coordinates are denotedwith ξ(x) ∈ RM and M ≤ N . The wells Wα of reactant(-1) and product (1), divided by a reaction barrier, are characterized by the set of points Γα (α = −1, 1), corresponding
to equilibrium configurations in each well, i.e., an unbiased
simulation on the PES U0 (x) is likely to stay in these wells.
The indicator function for a well is
1 for x ∈ Wα(3)⊮α (x) =0 for x ∈/ Wα .

a(t) =

∂L
.
∂z(t)

(6)

To solve (6) we introduce the new time τ ∈ (0, τe ) suchthat z(τ = 0) = z(t = te ) and z(τ = τe ) = z(t = 0).
This backward flowing time reflects that the loss is only
influenced by preceding points in time. In backward flowing
time τ the adjoint vectors obey the equations∂L∂z(τ = 0)∂f (z(τ ), θ)ȧ(τ ) = a(τ )T.
∂za(τ = 0) =Under suitable regularity conditions, we can cast the escapeprobability pα of a trajectory X with x(t0 ) ∈ W−α aspα = P sup (⊮α (x(t)) > 0 x(t0 ) ∈ W−α . (4)(7)

The total gradient of the loss function with respect to bias
parameters is then obtained byt<teA soft loss function that is continuous everywhere and differentiable is defined as(0
if ∃ x(t) ∈ Wα
L = Lξ α =,min (ξ(x(t)) − ξ α )2 otherwise

∂L
=
∂θ

Z τe
0

a(τ )T∂f (z(τ ), θ)
dτ .
∂θ

(8)

While solving (7), z(τ ) can be either saved or reconstructedby running dynamics (2) backward, depending on the memory and computational trade-off. The algorithm for running
the adjoint method for the dynamics that includes random
noise is developed and analyzed in (Li et al., 2020).

t0 <t<te(5)
where for each trajectory a random single ξ α ∈ Γα is selected for the loss function. The term (ξ(x(t)) − ξ α )2 isthe distance metric measuring how close a trajectory got
2

DiffSim for Rare Events

3. Challenges and Solutions
3.1. Challenges
Ideally, one would simulate the biased dynamics (2), compute the loss (5), backpropagate by solving (7), and after anumber of training epochs obtain the biasing potential thatenhances transitions. However, currently the use of DiffSim
is hampered by the following issues:
1. Gradient control Past efforts have been devoted tounderstand the behavior of gradients (their explosion orvanishing) that arise while optimizing DiffSims (Suhet al., 2022; Huang et al., 2021; Metz et al., 2021).
A DiffSim can be arbitrarily deep, however, it canbe challenging to controllably backpropagate complex
Hamiltonians. In fact, even simple Hamiltonians canyield exploding gradients (Galimberti et al., 2021).
2. Multiscale Problem The dynamics on U0 (x) can
include very high and very low frequency motions.
However, only the slow dynamics should be controlled
by the trainable bias, as those are associated with therare event. Avoiding fitting fast fluctuations reduces
the noise in the gradients used for training.
3. Chaotic behaviour Small changes in initial conditionsof some Hamiltonian systems can result in exponentially different trajectories (Percival I, 1987).
4. One large parameter update per trajectory DiffSims produce one update per trajectory. Obtaining
a sufficient number of gradient updates is expensive
when long trajectories are required.

Figure 1. Comparison of the adjoint evolution for original andpartially detached graphs for simulations on the 2D Müller-Brown
potential with parameters described in (Šı́pka et al., 2023).

Additionally, the .detach() operator allows for theorems describing the adjoint magnitude and its controllability by thefriction constant γ. Consider a trajectory xt generated by (2)with time t in a possibly infinite time interval I ⊂ (−∞, ∞).
The loss (5) is defined for a point xtL . To optimize B(x, θ),
we need to backpropagate the gradient of this loss throughevery point preceding xtL , using backwards flowing time τ .
Definition 3.1 (Differentiable Trajectory). A Differentiable Trajectory T is defined by the following quadruple(z(t), L(z(tL ), f (z(t)), f˜(z(t))): Let z(t) ∈ Ωx × Ωp
where Ωx ⊂ RN and Ωp ⊂ RN for t ∈ (ti , te ), where−∞ ≤ ti < te ≤ ∞ be the sequence of states generate by
the dynamics f (z(t)) from a certain initial state z(t0 ), t0 ∈[ti , te ]. We define a loss function L(z(tL )) at time tL . The
gradient dynamics of the loss function is guided by the dynamics f˜(z(t)) that includes possible .detach() operators.
The backward dynamics is represented in the reverse flowing
time τ starting from z(tL ) = z(τ = 0) to z(τe ) = z(ti ).

We address all these challenges and show how to efficiently
and controllably learn slow dynamics necessary for the investigation of reactions.

This allows us to formulate a theorem, that is proven in(Šı́pka et al., 2023)

3.2. Partial backpropagation
Stopping the flow of the gradient through x with the .detach() operator introduced, e.g., in Refs. (Foerster et al.,2018; Schulman et al., 2015; Zhang et al., 2019) reduces thecomplexity and level of detail in the equations, the computational graph is pruned such that backpropagation occurs
only in the momenta. The mathematical idea behind theuse of .detach() operator is better explained in (Šı́pka et al.,2023). We demonstrate how the approach suppresses highorder oscillation in Figure 1.

Theorem 3.2 (Converging adjoints). Td be a Differentiable
Trajectory for the molecular dynamics presented in the paper. Let U (x) ∈ C 2 (Ωx ) and denote the spectrum of its
2
Hessian ∂ U∂(x(t))
by λi (x(t)). Define λmin as
2x
λmin = inf min λi (x(τ )).
τ >0

i

(10)

Then for every γ that fulfills: (∆τ λmin + γ) = ϵ > 0, it
holds:
2
2∀τ > 0 : ∥a(τ )∥ ≤ ∥a(0)∥ e−2ϵτ(11)The use of the .detach() operator simplifies the adjoint time
derivative top∂∂U (x(τ ))ȧ(τ ) = aT−
− γp(τ ) + 2γkB T R(τ ) This theorem allows us to backpropagate the dynamics with∂p(τ )∂xout exploding gradients as long as γ is sufficiently large.
2
The
exponential scaling of the adjoints also indicates that
∂ U (x(τ ))
∆τ
−
γa(τ
)
.
= −aT (τ )afteridentifying the point for which the loss is calculated,
∂2x
only
a
few points need to be considered before a(τ ) van(9)3DiffSim for Rare Events

Figure 3. left Bias function (right) log-density obtained from long
biased simulations.

Figure 2. Log-density of simulated points before (left) and after
the training (right) of bias function by differentiable simulations.
The right plot shows how well all important regions are sampled
after training. The background of the Figure is the UMB (x, y), theunderlying Muller-Brown potential.

produces an almost even log-density along the expectedtransition path.
4.2. Cis-trans isomerization of azobenzene

ishes. Further adjoint propagation does not contribute to thegradient update.

Here, we study the cis-trans isomerization of azobenzene in
the electronic ground state (the reaction is shown in the rightside of Figure 3). A trained neural network potentials was
obtained from Ref. (Axelrod et al., 2022) and is used forU0 . For this system we chose two descriptors, the CNNC
dihedral angle and one NNC bond angle. These two angles
describe rotation and inversion mechanisms. See (Axelrod et al., 2022) for more details about different possible
mechanisms.

The theorem also highlights when backpropagation
may lead to exploding gradients. If the expression(∆τ λmin + γ) = ϵ < 0, then the upper bound may not
hold, and gradients can increase exponentially. A strongly
negative λi of the Hessian indicates a concave part in thepotential landscape, generally problematic for control. However, with ∆τ and γ, we have two robust dials to ensure
non-exploding adjoints.

As the isomerization barrier is relatively high and reactant
and product geometries are comparatively stiff, the reaction
investigation is split into two steps. First, we simulate at
300 K with a relatively high learning rate to get an estimate of the reaction barrier. Then, the bias is refined at1000 K with a much lower learning rate (see left part of
Figure 3). The log-density (Figure 3, right side) obtained
from simulations with the fully trained bias clearly shows arotation assisted inversion mechanism. This is analogous to
the findings of (Axelrod et al., 2022).

3.3. Mini-batching the graphTo alleviate challenge 4 we introduce a technique calledgraph mini-batching. The idea is to calculate trajectory
depended gradients first (the adjoints a) in one pass andthen split them to mini-batches. The adjoints are then usedduring training as vectors in Jacobi-vector products (8). The
approach stabilizes learning and allows for much lower
learning rates, better suited for training neural networks. An
example of a use case is more thoroughly discussed in thework (Šı́pka et al., 2023).

5. Conclusion
4. Applications

We present a way for DiffSims to increase the probability ofa barrier crossing by learning a biasing potential. The effectiveness of the framework is demonstrated on a common toy
system used for testing enhanced sampling algorithms andon a challenging cis-trans isomerization. It is shown that aneural network potential to can provide forces for a differentiable simulator and create a setup suitable for accelerators
as dozens of replicas are simulated at the same time, significantly speeding up the algorithm. In the theoretical sectionof the manuscript, the gradient behavior in a differentiable
simulator is discussed and their convergence in a so called
detached and undetached regime is investigated.

4.1. 2D Muller-Brown potential
Details of the Muller-Brown PES (Müller & Brown, 1979)
and the biasing function are given in (Šı́pka et al., 2023).
This PES is difficult as no linear combination of the Cartesian coordinates yields a good CV. After training the bias
potential via DiffSim, the biased dynamics generate increasingly many successful transitions between reactants andproducts along the expected transition path. Comparing
right and left side of Figure 2, the unbiased dynamics wastrapped in the local minima whereas the biased dynamics
4DiffSim for Rare Events

References

Physics for Learning and Control. In Bengio, S.,
Wallach, H., Larochelle, H., Grauman, K., CesaBianchi, N., and Garnett, R. (eds.), Advances in Neural
Information Processing Systems, volume 31. Curran Associates, Inc., 2018. URL https://proceedings.
neurips.cc/paper/2018/file/842424a1d0595b76ec4fa03c46e8d755-Paper.
pdf.

Abrams, C. and Bussi, G. Enhanced Sampling in Molecular
Dynamics Using Metadynamics, Replica-Exchange, andTemperature-Acceleration. Entropy 2014, Vol. 16, Pages
163-199, 16(1):163–199, 12 2013. ISSN 1099-4300.
doi: 10.3390/E16010163.
URL https://www.
mdpi.com/1099-4300/16/1/163/htmhttps://www.mdpi.com/1099-4300/16/1/163.

Degrave, J., Hermans, M., Dambre, J., and wyffels, F. A Differentiable Physics Engine for Deep Learning in Robotics.
arXiv, (1611.01652), 11 2016.

Axelrod, S., Shakhnovich, E., and Gómez-Bombarelli, R.
Excited state non-adiabatic dynamics of large photoswitchable molecules using a chemically transferable
machine learning potential. Nature Communications,13(1):3440, 6 2022. ISSN 2041-1723. doi: 10.1038/
s41467-022-30999-w.

Foerster, J., Farquhar, G., Al-Shedivat, M., Rocktäschel,T., Xing, E. P., and Whiteson, S. DiCE: The Infinitely
Differentiable Monte-Carlo Estimator. 2 2018.

Bonati, L., Rizzi, V., and Parrinello, M. Data-Driven
Collective Variables for Enhanced Sampling. The
journal of physical chemistry letters, 11(8):2998–3004,
4 2020.
ISSN 1948-7185.
doi: 10.1021/ACS.
JPCLETT.0C00535. URL https://pubmed.ncbi.
nlm.nih.gov/32239945/.

Galimberti, C. L., Furieri, L., Xu, L., and Ferrari-Trecate, G.
Hamiltonian Deep Neural Networks Guaranteeing Nonvanishing Gradients by Design. arXiv, (2105.13205), 5
2021.

Callen, H. B. and Scott, H. L. Thermodynamics and an Introduction to Thermostatistics, 2nd ed, volume 66. 1998.
doi: 10.1119/1.19071.

Greener, J. G. and Jones, D. T. Differentiable molecularsimulation can learn all the parameters in a coarse-grained
force field for proteins. PLOS ONE, 16(9):e0256990,
9 2021. ISSN 1932-6203. doi: 10.1371/journal.pone.
0256990.

Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud,
D. Neural ordinary differential equations. In Advancesin Neural Information Processing Systems, volume 2018December, 2018.

Hu, Y., Li, T.-M., Anderson, L., Ragan-Kelley, J., and Durand, F. Taichi. ACM Transactions on Graphics, 38(6):
1–16, 12 2019. ISSN 0730-0301. doi: 10.1145/3355089.
3356506.

Chipot, C. Frontiers in free-energy calculations of biologicalsystems. Wiley Interdisciplinary Reviews: Computational Molecular Science, 4(1):71–89, 1 2014. ISSN
1759-0884. doi: 10.1002/WCMS.1157. URL https://onlinelibrary.wiley.com/doi/full/10.
1002/wcms.1157https://onlinelibrary.
wiley.com/doi/abs/10.1002/wcms.
1157https://wires.onlinelibrary.wiley.
com/doi/10.1002/wcms.1157.

Hu, Y., Anderson, L., Li, T.-M., Sun, Q., Carr, N.,
Ragan-Kelley, J., and Durand, F. DiffTaichi: Differentiable Programming for Physical Simulation. In International Conference on Learning Representations,2020. URL https://openreview.net/forum?
id=B1eB5xSFvr.
Huang, Z., Hu, Y., Du, T., Zhou, S., Su, H., Tenenbaum,
J. B., and Gan, C. PlasticineLab: A Soft-Body Manipulation Benchmark with Differentiable Physics. In
International Conference on Learning Representations,2021. URL https://openreview.net/forum?
id=xCcdBRQEDW.

Chipot, C. and Pohorille, A.
Free Energy Calculations, volume 86 of Springer Series in CHEMICAL PHYSICS.
Springer Berlin Heidelberg,
Berlin, Heidelberg, 2007.
ISBN 978-3-54038447-2. doi: 10.1007/978-3-540-38448-9. URL
http://link.springer.com/10.1007/
978-3-540-38448-9.

Ingraham, J., Riesselman, A., Sander, C., and Marks, D.
Learning Protein Structure with a Differentiable Simulator. In International Conference on Learning Representations, 2019. URL https://openreview.net/
forum?id=Byg3y3C9Km.

Darve, E. and Pohorille, A. Calculating free energies using
average force. 2001. doi: 10.1063/1.1410978. URL
http://jcp.aip.org/jcp/copyright.jsp.

Laio, A. and Parrinello, M.
Escaping freeenergy minima.
Proceedings of the National
Academy of Sciences of the United States ofAmerica, 99(20):12562–12566, 10 2002.
ISSN

de Avila Belbute-Peres, F., Smith, K., Allen, K., Tenenbaum, J., and Kolter, J. Z. End-to-End Differentiable5DiffSim for Rare Events

00278424. doi: 10.1073/PNAS.202427399/ASSET/
2381B9FD-FD4C-4ED9-9DCB-75FD02E8EA7E/ASSETS/GRAPHIC/PQ2024273003.JPEG.
URL
https://www.pnas.org/doi/abs/10.1073/
pnas.202427399.

Spiwok, V., Sucur, Z., and Hosek, P. Enhanced sampling techniques in biomolecular simulations. Biotechnology advances, 33(6 Pt 2):1130–1140, 2015. ISSN
1873-1899. doi: 10.1016/J.BIOTECHADV.2014.11.
011.
URL https://pubmed.ncbi.nlm.nih.
gov/25482668/.

Lev Semenovich Pontryagin, MishGamkrelidze RV,
Boltyanskii VG, and Gamkrelidze RV. The Mathematical
Theory of Optimal Processes. 1962.

Suh, H. J. T., Simchowitz, M., Zhang, K., and Tedrake, R.
Do Differentiable Simulators Give Better Policy Gradients? 2 2022.

Li, X., Wong, T.-K. L., Chen, R. T. Q., and Duvenaud, D. K.
Scalable Gradients and Variational Inference for Stochastic Differential Equations. In Proceedings of The 2nd
Symposium on Advances in Approximate Bayesian Inference, volume 118 of Proceedings of Machine Learning
Research, pp. 1–28. PMLR, 9 2020. URL https://
proceedings.mlr.press/v118/li20a.html.

Sultan, M. M. and Pande, V. S. Automated design ofcollective variables using supervised machine learning.
The Journal of Chemical Physics, 149(9):94106, 2018.
doi: 10.1063/1.5029972. URL https://doi.org/
10.1063/1.5029972.
Sun, L., Vandermause, J., Batzner, S., Xie, Y., Clark, D.,
Chen, W., and Kozinsky, B. Multitask Machine Learning of Collective Variables for Enhanced Sampling ofRare Events. Journal of Chemical Theory and Computation, 18(4):2341–2353, 2022. doi: 10.1021/acs.jctc.
1c00143. URL https://doi.org/10.1021/acs.
jctc.1c00143.

Mendels, D., Piccini, G., and Parrinello, M. Collective Variables from Local Fluctuations. Journal of Physical Chemistry Letters, 9(11):2776–2781, 6 2018. ISSN 19487185.
doi:
10.1021/ACS.JPCLETT.8B00733/ASSET/
IMAGES/LARGE/JZ-2018-00733T{\ }0005.JPEG.
URL https://pubs.acs.org/doi/full/10.
1021/acs.jpclett.8b00733.

Torrie, G. M., Valleau, J. P., Torrie, G. M., andValleau, J. P. Nonphysical Sampling Distributions
in Monte Carlo Free-Energy Estimation: Umbrella
Sampling.
JCoPh, 23(2):187–199, 1977.
ISSN
0021-9991.
doi: 10.1016/0021-9991(77)90121-8.
URL https://ui.adsabs.harvard.edu/abs/
1977JCoPh..23..187T/abstract.

Metz, L., Freeman, C. D., Schoenholz, S. S., and Kachman,
T. Gradients are Not All You Need. arXiv, (2111.05803),
11 2021.
Müller, K. and Brown, L. D. Location of saddle points andminimum energy paths by a constrained simplex optimization procedure. Theoretica chimica acta, 53(1):75–93,
1979. ISSN 1432-2234. doi: 10.1007/BF00547608. URL
https://doi.org/10.1007/BF00547608.

Valsson, O., Tiwary, P., and Parrinello, M.
Enhancing Important Fluctuations:
Rare Eventsand Metadynamics from a Conceptual Viewpoint.
http://dx.doi.org/10.1146/annurev-physchem-040215112229, 67:159–184, 5 2016. ISSN 0066426X. doi:
10.1146/ANNUREV-PHYSCHEM-040215-112229.
URL
https://www.annualreviews.
org/doi/abs/10.1146/
annurev-physchem-040215-112229.

Percival I. Chaos in hamiltonian systems. Proceedingsof the Royal Society of London. A. Mathematical andPhysical Sciences, 413(1844):131–143, 9 1987. ISSN
0080-4630. doi: 10.1098/rspa.1987.0105.
Rohrdanz, M. A., Zheng, W., and Clementi, C. Discovering
Mountain Passes via Torchlight: Methods for the Definition of Reaction Coordinates and Pathways in Complex
Macromolecular Reactions. Annual Review of Physical
Chemistry, 64(1):295–316, 4 2013. ISSN 0066-426X.
doi: 10.1146/annurev-physchem-040412-110006.

Wang, D. and Tiwary, P. State predictive information bottleneck. The Journal of Chemical Physics, 154(13):134111,
4 2021. ISSN 0021-9606. doi: 10.1063/5.0038198.
URL https://aip.scitation.org/doi/abs/
10.1063/5.0038198.

Schulman, J., Heess, N., Weber, T., and Abbeel, P. Gradient
Estimation Using Stochastic Computation Graphs. arXiv,(1506.05254), 6 2015.

Wang, W., Axelrod, S., and Gómez-Bombarelli, R. Differentiable Molecular Simulations for Control and Learning.
arXiv, (2003.00868), 2 2020.

Šı́pka, M., Erlebach, A., and Grajciar, L. Understanding
chemical reactions via variational autoencoder and atomic
representations. arXiv, (2202.00817), 3 2022.

Wang, Y., Ribeiro, J. M. L., and Tiwary, P. Past–future
information bottleneck for sampling molecular reaction coordinate simultaneously with thermodynamics
and kinetics. Nature Communications 2019 10:1, 10Šı́pka, M., Dietschreit, J. C. B., Grajciar, L., and GómezBombarelli, R. Differentiable Simulations for Enhanced
Sampling of Rare Events. 1 2023.
6

DiffSim for Rare Events(1):1–8, 8 2019. ISSN 2041-1723. doi: 10.1038/
s41467-019-11405-4. URL https://www.nature.
com/articles/s41467-019-11405-4.
Wehmeyer, C. and Noé, F. Time-lagged autoencoders:
Deep learning of slow collective variables for molecular kinetics. The Journal of Chemical Physics, 148(24):241703, 2018. doi: 10.1063/1.5011399. URL
https://doi.org/10.1063/1.5011399.
Zhang, S.-X., Wan, Z.-Q., and Yao, H. Automatic Differentiable Monte Carlo: Theory and Application. arXiv,(1911.09117), 11 2019.

7