Petr Plecháč

Giorgos Arampatzis, Markos A. Katsoulakis, Petr Plechac

In this paper we study from a numerical analysis perspective the Fractional Step Kinetic Monte Carlo (FS-KMC) algorithms proposed in [1] for the parallel simulation of spatially distributed particle systems on a lattice. FS-KMC are fractional step algorithms with a time-stepping window $Δt$, and as such they are inherently partially asynchronous since there is no processor communication during the period $Δt$. In this contribution we primarily focus on the error analysis of FS-KMC algorithms as approximations of conventional, serial kinetic Monte Carlo (KMC). A key aspect of our analysis relies on emphasising a goal-oriented approach for suitably defined macroscopic observables (e.g., density, energy, correlations, surface roughness), rather than focusing on strong topology estimates for individual trajectories. One of the key implications of our error analysis is that it allows us to address systematically the processor communication of different parallelization strategies for KMC by comparing their (partial) asynchrony, which in turn is measured by their respective fractional time step $Δt$ for a prescribed error tolerance.

Evangelia Kalligiannaki, Markos A. Katsoulakis, Petr Plechac

We propose a hierarchy of multi-level kinetic Monte Carlo methods for sampling high-dimensional, stochastic lattice particle dynamics with complex interactions. The method is based on the efficient coupling of different spatial resolution levels, taking advantage of the low sampling cost in a coarse space and by developing local reconstruction strategies from coarse-grained dynamics. Microscopic reconstruction corrects possibly significant errors introduced through coarse-graining, leading to the controlled-error approximation of the sampled stochastic process. In this manner, the proposed multi-level algorithm overcomes known shortcomings of coarse-graining of particle systems with complex interactions such as combined long and short-range particle interactions and/or complex lattice geometries. Specifically, we provide error analysis for the approximation of long-time stationary dynamics in terms of relative entropy and prove that information loss in the multi-level methods is growing linearly in time, which in turn implies that an appropriate observable in the stationary regime is the information loss of the path measures per unit time. We show that this observable can be either estimated a priori, or it can be tracked computationally a posteriori in the course of a simulation. The stationary regime is of critical importance to molecular simulations as it is relevant to long-time sampling, obtaining phase diagrams and in studying metastability properties of high-dimensional complex systems. Finally, the multi-level nature of the method provides flexibility in combining rejection-free and null-event implementations, generating a hierarchy of algorithms with an adjustable number of rejections that includes well-known rejection-free and null-event algorithms.

Christian Bayer, Håkon Hoel, Petr Plecháč et al.

Born-Oppenheimer dynamics is shown to provide an accurate approximation of time-independent Schrödinger observables for a molecular system with an electron spectral gap, in the limit of large ratio of nuclei and electron masses, without assuming that the nuclei are localized to vanishing domains. The derivation, based on a Hamiltonian system interpretation of the Schrödinger equation and stability of the corresponding Hamilton-Jacobi equation, bypasses the usual separation of nuclei and electron wave functions, includes caustic states and gives a different perspective on the Born-Oppenheimer approximation, Schrödinger Hamiltonian systems and numerical simulation in molecular dynamics modeling at constant energy microcanonical ensembles.

Evangelia Kalligiannaki, Markos A. Katsoulakis, Petr Plechac

We propose an efficient Markov Chain Monte Carlo method for sampling equilibrium distributions for stochastic lattice models, capable of handling correctly long and short-range particle interactions. The proposed method is a Metropolis-type algorithm with the proposal probability transition matrix based on the coarse-grained approximating measures introduced in a series of works of M. Katsoulakis, A. Majda, D. Vlachos and P. Plechac, L. Rey-Bellet and D.Tsagkarogiannis,. We prove that the proposed algorithm reduces the computational cost due to energy differences and has comparable mixing properties with the classical microscopic Metropolis algorithm, controlled by the level of coarsening and reconstruction procedure. The properties and effectiveness of the algorithm are demonstrated with an exactly solvable example of a one dimensional Ising-type model, comparing efficiency of the single spin-flip Metropolis dynamics and the proposed coupled Metropolis algorithm.

M. A. Katsoulakis, P. Plechac, L. Rey-Bellet et al.

We develop coarse-graining schemes for stochastic many-particle microscopic models with competing short- and long-range interactions on a d-dimensional lattice. We focus on the coarse-graining of equilibrium Gibbs states and using cluster expansions we analyze the corresponding renormalization group map. We quantify the approximation properties of the coarse-grained terms arising from different types of interactions and present a hierarchy of correction terms. We derive semi-analytical numerical schemes that are accompanied with a posteriori error estimates for coarse-grained lattice systems with short and long-range interactions.

Ting Wang, Petr Plechac, Jaroslaw Knap

We develop a class of data-driven generative models that approximate the solution operator for parameter-dependent partial differential equations (PDE). We propose a novel probabilistic formulation of the operator learning problem based on recently developed generative denoising diffusion probabilistic models (DDPM) in order to learn the input-to-output mapping between problem parameters and solutions of the PDE. To achieve this goal we modify DDPM to supervised learning in which the solution operator for the PDE is represented by a class of conditional distributions. The probabilistic formulation combined with DDPM allows for an automatic quantification of confidence intervals for the learned solutions. Furthermore, the framework is directly applicable for learning from a noisy data set. We compare computational performance of the developed method with the Fourier Network Operators (FNO). Our results show that our method achieves comparable accuracy and recovers the noise magnitude when applied to data sets with outputs corrupted by additive noise.

Petr Plechac, Mathias Rousset

We propose and analyze an implicit mass-matrix penalization (IMMP) technique which enables efficient and exact sampling of the (Boltzmann/Gibbs) canonical distribution associated to Hamiltonian systems with fast degrees of freedom (fDOFs). The penalty parameters enable arbitrary tuning of the timescale for the selected fDOFs, and the method is interpreted as an interpolation between the exact Hamiltonian dynamics and the dynamics with infinitely slow fDOFs (equivalent to geometrically corrected rigid constraints). This property translates in the associated numerical methods into a tunable trade-off between stability and dynamical modification. The penalization is based on an extended Hamiltonian with artificial constraints associated with each fDOF. By construction, the resulting dynamics is statistically exact with respect to the canonical distribution in position variables. The algorithms can be easily implemented with standard geometric integrators with algebraic constraints given by the expected fDOFs, and has no additional complexity in terms of enforcing the constraint and force evaluations. The method is demonstrated on a high dimensional system with non-convex interactions. Prescribing the macroscopic dynamical timescale, it is shown that the IMMP method increases the time-step stability region with a gain that grows linearly with the size of the system. The latter property, as well as consistency of the macroscopic dynamics of the IMMP method is proved rigorously for linear interactions. Finally, when a large stiffness parameter is introduced, the IMMP method can be tuned to be asymptotically stable, converging towards the heuristically expected Markovian effective dynamics on the slow manifold.

Petr Plecháč, Erik von Schwerin

We propose a Multi-level Monte Carlo technique to accelerate Monte Carlo sampling for approximation of properties of materials with random defects. The computational efficiency is investigated on test problems given by tight-binding models of a single layer of graphene or of $MoS_2$ where the integrated electron density of states per unit area is taken as a representative quantity of interest. For the chosen test problems the multi-level Monte Carlo estimators significantly reduce the computational time of standard Monte Carlo estimators to obtain a given accuracy.

Ting Wang, Petr Plecháč, David Aristoff

We propose two algorithms for simulating continuous time Markov chains in the presence of metastability. We show that the algorithms correctly estimate, under the ergodicity assumption, stationary averages of the process. Both algorithms, based on the idea of the parallel replica method, use parallel computing in order to explore metastable sets more efficiently. The algorithms require no assumptions on the Markov chains beyond ergodicity and the presence of identifiable metastability. In particular, there is no assumption on reversibility. For simpler illustration of the algorithms, we assume that a synchronous architecture is used throughout of the paper. We present error analyses, as well as numerical simulations on multi-scale stochastic reaction network models in order to demonstrate consistency of the method and its efficiency.

Christian Bayer, Hakon Hoel, Ashraful Kadir et al.

The difference of the values of observables for the time-independent Schroedinger equation, with matrix valued potentials, and the values of observables for ab initio Born-Oppenheimer molecular dynamics, of the ground state, depends on the probability to be in excited states and the electron/nuclei mass ratio. The paper first proves an error estimate (depending on the electron/nuclei mass ratio and the probability to be in excited states) for this difference of microcanonical observables, assuming that molecular dynamics space-time averages converge, with a rate related to the maximal Lyapunov exponent. The error estimate is uniform in the number of particles and the analysis does not assume a uniform lower bound on the spectral gap of the electron operator and consequently the probability to be in excited states can be large. A numerical method to determine the probability to be in excited states is then presented, based on Ehrenfest molecular dynamics and stability analysis of a perturbed eigenvalue problem.

Evangelia Kalligiannaki, Vagelis Harmandaris, Markos A. Katsoulakis et al.

Using the probabilistic language of conditional expectations we reformulate the force matching method for coarse-graining of molecular systems as a projection on spaces of coarse observables. A practical outcome of this probabilistic description is the link of the force matching method with thermodynamic integration. This connection provides a way to systematically construct a local mean force in order to optimally approximate the potential of mean force through force matching. We introduce a generalized force matching condition for the local mean force in the sense that allows the approximation of the potential of mean force under both linear and non-linear coarse graining mappings (e.g., reaction coordinates, end-to-end length of chains). Furthermore, we study the equivalence of force matching with relative entropy minimization which we derive for general non-linear coarse graining maps. We present in detail the generalized force matching condition through applications to specific examples in molecular systems.

Petr Plechac, Mathias Rousset

An implicit mass-matrix penalization (IMMP) of Hamiltonian dynamics is proposed, and associated dynamical integrators, as well as sampling Monte-Carlo schemes, are analyzed for systems with multiple time scales. The penalization is based on an extended Hamiltonian with artificial constraints associated with some selected DOFs. The penalty parameters enable arbitrary tuning of timescales for the selected DOFs. The IMMP dynamics is shown to be an interpolation between the exact Hamiltonian dynamics and the dynamics with rigid constraints. This property translates in the associated numerical integrator into a tunable trade-off between stability and dynamical modification. Moreover, a penalty that vanishes with the time-step yields order two convergent schemes for the exact dynamics. Moreover, by construction, the resulting dynamics preserves the canonical equilibrium distribution in position variables, up to a computable geometric correcting potential, leading to Metropolis-like unbiased sampling algorithms. The algorithms can be implemented with a simple modification of standard geometric integrators with algebraic constraints imposed on the selected DOFs, and has no additional complexity in terms of enforcing the constraints and force evaluations. The properties of the IMMP method are demonstrated numerically on the $N$-alkane model, showing that the time-step stability region of integrators and the sampling efficiency can be increased with a gain that grows with the size of the system. This feature is mathematically analyzed for a harmonic atomic chain model. When a large stiffness parameter is introduced, the IMMP method is shown to be asymptotically stable and to converge towards the heuristically expected Markovian effective dynamics on the slow manifold.

Markos A. Katsoulakis, Petr Plechac, Luc Rey-Bellet et al.

The primary objective of this work is to develop coarse-graining schemes for stochastic many-body microscopic models and quantify their effectiveness in terms of a priori and a posteriori error analysis. In this paper we focus on stochastic lattice systems of interacting particles at equilibrium. %such as Ising-type models. The proposed algorithms are derived from an initial coarse-grained approximation that is directly computable by Monte Carlo simulations, and the corresponding numerical error is calculated using the specific relative entropy between the exact and approximate coarse-grained equilibrium measures. Subsequently we carry out a cluster expansion around this first-and often inadequate-approximation and obtain more accurate coarse-graining schemes. The cluster expansions yield also sharp a posteriori error estimates for the coarse-grained approximations that can be used for the construction of adaptive coarse-graining methods. We present a number of numerical examples that demonstrate that the coarse-graining schemes developed here allow for accurate predictions of critical behavior and hysteresis in systems with intermediate and long-range interactions. We also present examples where they substantially improve predictions of earlier coarse-graining schemes for short-range interactions.

Markos A Katsoulakis, Petr Plechac, Alexandros Sopasakis

In this paper we investigate the approximation properties of the coarse-graining procedure applied to kinetic Monte Carlo simulations of lattice stochastic dynamics. We provide both analytical and numerical evidence that the hierarchy of the coarse models is built in a systematic way that allows for error control in both transient and long-time simulations. We demonstrate that the numerical accuracy of the CGMC algorithm as an approximation of stochastic lattice spin flip dynamics is of order two in terms of the coarse-graining ratio and that the natural small parameter is the coarse-graining ratio over the range of particle/particle interactions. The error estimate is shown to hold in the weak convergence sense. We employ the derived analytical results to guide CGMC algorithms and we demonstrate a CPU speed-up in demanding computational regimes that involve nucleation, phase transitions and metastability.