• In the classical work by Irving and Zwanzig [Irving J.H. and Zwanzig R.W., J. Chem. Phys. 19 (1951), 1173-1180 ] it has been shown that quantum observables for macroscopic density, momentum and energy satisfy the conservation laws of fluid dynamics. This work derives the corresponding classical molecular dynamics limit by extending Irving and Zwanzig's result to matrix-valued potentials for a general quantum particle system. The matrix formulation provides the semi-classical limit of the quantum observables in the conservation laws, also in the case where the temperature is large compared to the electron eigenvalue gaps. The classical limit of the quantum observables in the conservation laws is useful in order to determine the constitutive relations for the stress tensor and the heat flux by molecular dynamics simulations. The main new steps to obtain the molecular dynamics limit is to: (i) approximate the dynamics of quantum observables accurately by classical dynamics, by diagonalizing the Hamiltonian using a non linear eigenvalue problem, (ii) define the local energy density by partitioning a general potential, applying perturbation analysis of the electron eigenvalue problem, (iii) determine the molecular dynamics stress tensor and heat flux in the case of several excited electron states, and (iv) construct the initial particle phase-space density as a local grand canonical quantum ensemble determined by the initial conservation variables.
  • In this paper, we study Monte Carlo estimators based on the likelihood ratio approach for steady-state sensitivity. We first extend the result of Glynn and Olvera-Cravioto [arXiv:1707.02659] to the setting of continuous time Markov chains (CTMC) with a countable state space which include models such as stochastic reaction kinetics and kinetic Monte Carlo lattice system. Then we show that the variance of the centered LR estimators do not grow in time. This result suggests that the centered estimators should be favored when the mixing time of the CTMC is large. We demonstrate a practical implication of this analysis on a numerical benchmark of two examples for the biochemical reaction networks.
  • It is known that ab initio molecular dynamics based on the electron ground state eigenvalue can be used to approximate quantum observables in the canonical ensemble when the temperature is low compared to the first electron eigenvalue gap. This work proves that a certain weighted average of the different ab initio dynamics, corresponding to each electron eigenvalue, approximates quantum observables for any temperature. The proof uses the semiclassical Weyl law to show that canonical quantum observables of nuclei-electron systems, based on matrix valued Hamiltonian symbols, can be approximated by ab initio molecular dynamics with the error proportional to the electron-nuclei mass ratio. The result covers observables that depend on time-correlations. A combination of the Hilbert-Schmidt inner product for quantum operators and Weyl's law shows that the error estimate holds for observables and Hamiltonian symbols that have three and five bounded derivatives, respectively, provided the electron eigenvalues are distinct for any nuclei position and the observables are in the diagonal form with respect to the electron eigenstates.
  • We introduce a multiscale framework which combines time-dependent nonequilibrium Green function (TD-NEGF) algorithms, scaling linearly in the number of time steps and describing quantum-mechanically conduction electrons in the presence of time-dependent fields of arbitrary strength or frequency, with classical time evolution of local magnetic moments described by the Landau-Lifshitz-Gilbert (LLG) equation. Our TD-NEGF+LLG framework can be applied to a variety of problems where current-driven spin torque induces dynamics of magnetic moments as the key resource for next generation spintronics. Previous approaches to such nonequilibrium many-body system have neglected noncommutativity of a quantum Hamiltonian of conduction electrons at different times and, therefore, the impact of time-dependent magnetic moments on electrons which induce pumping of spin and charge currents that, in turn, can self-consistently affect the dynamics of magnetic moments themselves. Using magnetic domain wall (DW) as an example, we predict that its motion will pump time-dependent spin and charge currents (on the top of unpolarized DC charge current injected through normal metal leads to drive the DW motion), where the latter can be viewed as a realization of nonadiabatic quantum charge pumping due to time-dependence of the Hamiltonian and left-right symmetry breaking of the two-terminal device structure. The conversion of AC components of spin current, whose amplitude increases (decreases) as the DW approaches (distances from) the normal metal lead, into AC voltage via the inverse spin Hall effect offers a tool to precisely track the DW position along magnetic nanowire. We also quantify the DW transient inertial displacement due to its acceleration and deceleration by pulse current and the entailed spin and charge pumping.
  • We review a unified approach for computing: (i) spin-transfer torque in magnetic trilayers like spin-valves and magnetic tunnel junction, where injected charge current flows perpendicularly to interfaces; and (ii) spin-orbit torque in magnetic bilayers of the type ferromagnet/spin-orbit-coupled-material, where injected charge current flows parallel to the interface. Our approach requires to construct the torque operator for a given Hamiltonian of the device and the steady-state nonequilibrium density matrix, where the latter is expressed in terms of the nonequilibrium Green's functions and split into three contributions. Tracing these contributions with the torque operator automatically yields field-like and damping-like components of spin-transfer torque or spin-orbit torque vector, which is particularly advantageous for spin-orbit torque where the direction of these components depends on the unknown-in-advance orientation of the current-driven nonequilibrium spin density in the presence of spin-orbit coupling. We provide illustrative examples by computing spin-transfer torque in a one-dimensional toy model of a magnetic tunnel junction and realistic Co/Cu/Co spin-valve, both of which are described by first-principles Hamiltonians obtained from noncollinear density functional theory calculations; as well as spin-orbit torque in a ferromagnetic layer described by a tight-binding Hamiltonian which includes spin-orbit proximity effect within ferromagnetic monolayers assumed to be generated by the adjacent monolayer transition metal dichalcogenide.
  • Metastable condensed matter typically fluctuates about local energy minima at the femtosecond time scale before transitioning between local minima after nanoseconds or microseconds. This vast scale separation limits the applicability of classical molecular dynamics methods and has spurned the development of a host of approximate algorithms. One recently proposed method is diffusive molecular dynamics which aims to integrate a system of ordinary differential equations describing the likelihood of occupancy by one of two species, in the case of a binary alloy, while quasistatically evolving the locations of the atoms. While diffusive molecular dynamics has shown to be efficient and provide agreement with observations, it is fundamentally a model, with unclear connections to classical molecular dynamics. In this work, we formulate a spin-diffusion stochastic process and show how it can be connected to diffusive molecular dynamics. The spin-diffusion model couples a classical overdamped Langevin equation to a kinetic Monte Carlo model for exchange amongst the species of a binary alloy. Under suitable assumptions and approximations, spin-diffusion can be shown to lead to diffusive molecular molecular dynamics type models. The key assumptions and approximations include a well defined time scale separation, a choice of spin exchange rates, a low temperature approximation, and a mean field type approximation. We derive several models from different assumptions and show their relationship to diffusive molecular dynamics. Differences and similarities amongst the models are explored in a simple test problem.
  • Stochastic reaction networks that exhibit bistability are common in many fields such as systems biology and materials science. Sampling of the stationary distribution is crucial for understanding and characterizing the long term dynamics of bistable stochastic dynamical systems. However, this is normally hindered by the insufficient sampling of the rare transitions between the two metastable regions. In this paper, we apply the parallel replica (ParRep) method for continuous time Markov chain to accelerate the stationary distribution sampling of bistable stochastic reaction networks. The proposed method uses parallel computing to accelerate the sampling of rare transitions and it is very easy to implement. We combine ParRep with the path space information bounds for parametric sensitivity analysis. We demonstrate the efficiency and accuracy of the method by studying the Schl\"{o}gl model and the genetic switches network.
  • 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. 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.
  • 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.
  • In this paper, we discuss information-theoretic tools for obtaining optimized coarse-grained molecular models for both equilibrium and non-equilibrium molecular dynamics. The latter are ubiquitous in physicochemical and biological applications, where they are typically associated with coupling mechanisms, multi-physics and/or boundary conditions. In general the non-equilibrium steady states are not known explicitly as they do not necessarily have a Gibbs structure. The presented approach can compare microscopic behavior of molecular systems to parametric and non-parametric coarse-grained one using the relative entropy between distributions on the path space and setting up a corresponding path space variational inference problem. The methods can become entirely data-driven when the microscopic dynamics are replaced with corresponding correlated data in the form of time series. Furthermore, we present connections and generalizations of force matching methods in coarse-graining with path-space information methods, as well as demonstrate the enhanced transferability of information-based parameterizations to general observables due to information inequalities. We further discuss methodological connections between information-based coarse-graining of molecular systems and variational inference methods primarily developed in the machine learning community. However, we note that the work presented here addresses variational inference for correlated time series due to the focus on dynamics. The applicability of the proposed methods is demonstrated on high-dimensional stochastic processes given by Langevin, overdamped and driven Langevin dynamics of interacting particles.
  • In the presence of multiscale dynamics in a reaction network, direct simulation methods become inefficient as they can only advance the system on the smallest scale. This work presents stochastic averaging techniques to accelerate computations for obtaining estimates of expected values and sensitivities with respect to the steady state distribution. A two-time-scale formulation is used to establish bounds on the bias induced by the averaging method. Further, this formulation provides a framework to create an accelerated `averaged' version of most single-scale sensitivity estimation method. In particular, we propose a new lower-variance ergodic likelihood ratio type estimator for steady-state estimation and show how one can adapt it to accelerated simulations of multiscale systems.Lastly, we develop an adaptive "batch-means" stopping rule for determining when to terminate the micro-equilibration process.
  • Uncertainty quantification is a primary challenge for reliable modeling and simulation of complex stochastic dynamics. Such problems are typically plagued with incomplete information that may enter as uncertainty in the model parameters, or even in the model itself. Furthermore, due to their dynamic nature, we need to assess the impact of these uncertainties on the transient and long-time behavior of the stochastic models and derive corresponding uncertainty bounds for observables of interest. A special class of such challenges is parametric uncertainties in the model and in particular sensitivity analysis along with the corresponding sensitivity bounds for stochastic dynamics. Moreover, sensitivity analysis can be further complicated in models with a high number of parameters that render straightforward approaches, such as gradient methods, impractical. In this paper, we derive uncertainty and sensitivity bounds for path-space observables of stochastic dynamics in terms of new goal-oriented divergences; the latter incorporate both observables and information theory objects such as the relative entropy rate. These bounds are tight, depend on the variance of the particular observable and are computable through Monte Carlo simulation. In the case of sensitivity analysis, the derived sensitivity bounds rely on the path Fisher Information Matrix, hence they depend only on local dynamics and are gradient-free. These features allow for computationally efficient implementation in systems with a high number of parameters, e.g., complex reaction networks and molecular simulations.
  • 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.
  • 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.
  • In this paper we focus on the development of new methods suitable for efficient and reliable coarse-graining of {\it non-equilibrium} molecular systems. In this context, we propose error estimation and controlled-fidelity model reduction methods based on Path-Space Information Theory, and combine it with statistical parametric estimation of rates for non-equilibrium stationary processes. The approach we propose extends the applicability of existing information-based methods for deriving parametrized coarse-grained models to Non-Equilibrium systems with Stationary States (NESS). In the context of coarse-graining it allows for constructing optimal parametrized Markovian coarse-grained dynamics, by minimizing information loss (due to coarse-graining) on the path space. Furthermore, the associated path-space Fisher Information Matrix can provide confidence intervals for the corresponding parameter estimators. We demonstrate the proposed coarse-graining method in a non-equilibrium system with diffusing interacting particles, driven by out-of-equilibrium boundary conditions.
  • 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 $\Delta t$, and as such they are inherently partially asynchronous since there is no processor communication during the period $\Delta 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 $\Delta t$ for a prescribed error tolerance.
  • 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.
  • Born-Oppenheimer dynamics is shown to provide an accurate approximation of time-independent Schr\"odinger 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\"odinger 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\"odinger Hamiltonian systems and numerical simulation in molecular dynamics modeling at constant energy microcanonical ensembles.
  • In this work we propose a hierarchy of Monte Carlo methods for sampling equilibrium properties of stochastic lattice systems with competing short and long range interactions. Each Monte Carlo step is composed by two or more sub - steps efficiently coupling coarse and microscopic state spaces. The method can be designed to sample the exact or controlled-error approximations of the target distribution, providing information on levels of different resolutions, as well as at the microscopic level. In both strategies the method achieves significant reduction of the computational cost compared to conventional Markov Chain Monte Carlo methods. Applications in phase transition and pattern formation problems confirm the efficiency of the proposed methods.
  • We present a mathematical framework for constructing and analyzing parallel algorithms for lattice Kinetic Monte Carlo (KMC) simulations. The resulting algorithms have the capacity to simulate a wide range of spatio-temporal scales in spatially distributed, non-equilibrium physiochemical processes with complex chemistry and transport micro-mechanisms. The algorithms can be tailored to specific hierarchical parallel architectures such as multi-core processors or clusters of Graphical Processing Units (GPUs). The proposed parallel algorithms are controlled-error approximations of kinetic Monte Carlo algorithms, departing from the predominant paradigm of creating parallel KMC algorithms with exactly the same master equation as the serial one. Our methodology relies on a spatial decomposition of the Markov operator underlying the KMC algorithm into a hierarchy of operators corresponding to the processors' structure in the parallel architecture. Based on this operator decomposition, we formulate Fractional Step Approximation schemes by employing the Trotter Theorem and its random variants; these schemes, (a) determine the communication schedule} between processors, and (b) are run independently on each processor through a serial KMC simulation, called a kernel, on each fractional step time-window. Furthermore, the proposed mathematical framework allows us to rigorously justify the numerical and statistical consistency of the proposed algorithms, showing the convergence of our approximating schemes to the original serial KMC. The approach also provides a systematic evaluation of different processor communicating schedules.
  • A new method of deriving reduced models of Hamiltonian dynamical systems is developed using techniques from optimization and statistical estimation. Given a set of resolved variables that define a model reduction, the quasi-equilibrium ensembles associated with the resolved variables are employed as a family of trial probability densities on phase space. The residual that results from submitting these trial densities to the Liouville equation is quantified by an ensemble-averaged cost function related to the information loss rate of the reduction. From an initial nonequilibrium state, the statistical state of the system at any later time is estimated by minimizing the time integral of the cost function over paths of trial densities. Statistical closure of the underresolved dynamics is obtained at the level of the value function, which equals the optimal cost of reduction with respect to the resolved variables, and the evolution of the estimated statistical state is deduced from the Hamilton-Jacobi equation satisfied by the value function. In the near-equilibrium regime, or under a local quadratic approximation in the far-from-equilibrium regime, this best-fit closure is governed by a differential equation for the estimated state vector coupled to a Riccati differential equation for the Hessian matrix of the value function. Since memory effects are not explicitly included in the trial densities, a single adjustable parameter is introduced into the cost function to capture a time-scale ratio between resolved and unresolved motions. Apart from this parameter, the closed equations for the resolved variables are completely determined by the underlying deterministic dynamics.
  • 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.
  • 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.
  • 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.
  • 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.