
In the classical work by Irving and Zwanzig [Irving J.H. and Zwanzig R.W., J.
Chem. Phys. 19 (1951), 11731180 ] 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 matrixvalued
potentials for a general quantum particle system. The matrix formulation
provides the semiclassical 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
phasespace 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 steadystate sensitivity. We first extend the result of Glynn and
OlveraCravioto [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
nucleielectron systems, based on matrix valued Hamiltonian symbols, can be
approximated by ab initio molecular dynamics with the error proportional to the
electronnuclei mass ratio. The result covers observables that depend on
timecorrelations. A combination of the HilbertSchmidt 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 timedependent
nonequilibrium Green function (TDNEGF) algorithms, scaling linearly in the
number of time steps and describing quantummechanically conduction electrons
in the presence of timedependent fields of arbitrary strength or frequency,
with classical time evolution of local magnetic moments described by the
LandauLifshitzGilbert (LLG) equation. Our TDNEGF+LLG framework can be
applied to a variety of problems where currentdriven spin torque induces
dynamics of magnetic moments as the key resource for next generation
spintronics. Previous approaches to such nonequilibrium manybody system have
neglected noncommutativity of a quantum Hamiltonian of conduction electrons at
different times and, therefore, the impact of timedependent magnetic moments
on electrons which induce pumping of spin and charge currents that, in turn,
can selfconsistently affect the dynamics of magnetic moments themselves. Using
magnetic domain wall (DW) as an example, we predict that its motion will pump
timedependent 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 timedependence of the Hamiltonian and leftright symmetry breaking of
the twoterminal 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) spintransfer torque in
magnetic trilayers like spinvalves and magnetic tunnel junction, where
injected charge current flows perpendicularly to interfaces; and (ii)
spinorbit torque in magnetic bilayers of the type
ferromagnet/spinorbitcoupledmaterial, 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 steadystate
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 fieldlike
and dampinglike components of spintransfer torque or spinorbit torque
vector, which is particularly advantageous for spinorbit torque where the
direction of these components depends on the unknowninadvance orientation of
the currentdriven nonequilibrium spin density in the presence of spinorbit
coupling. We provide illustrative examples by computing spintransfer torque in
a onedimensional toy model of a magnetic tunnel junction and realistic
Co/Cu/Co spinvalve, both of which are described by firstprinciples
Hamiltonians obtained from noncollinear density functional theory calculations;
as well as spinorbit torque in a ferromagnetic layer described by a
tightbinding Hamiltonian which includes spinorbit 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 spindiffusion stochastic process and show how
it can be connected to diffusive molecular dynamics. The spindiffusion 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, spindiffusion 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 multiscale stochastic reaction
network models in order to demonstrate consistency of the method and its
efficiency.

We propose a Multilevel 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
tightbinding 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 multilevel Monte Carlo estimators
significantly reduce the computational time of standard Monte Carlo estimators
to obtain a given accuracy.

In this paper, we discuss informationtheoretic tools for obtaining optimized
coarsegrained molecular models for both equilibrium and nonequilibrium
molecular dynamics. The latter are ubiquitous in physicochemical and biological
applications, where they are typically associated with coupling mechanisms,
multiphysics and/or boundary conditions. In general the nonequilibrium 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 nonparametric coarsegrained 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
datadriven 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 coarsegraining with
pathspace information methods, as well as demonstrate the enhanced
transferability of informationbased parameterizations to general observables
due to information inequalities.
We further discuss methodological connections between informationbased
coarsegraining 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 highdimensional 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 twotimescale
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 singlescale sensitivity estimation method. In
particular, we propose a new lowervariance ergodic likelihood ratio type
estimator for steadystate estimation and show how one can adapt it to
accelerated simulations of multiscale systems.Lastly, we develop an adaptive
"batchmeans" stopping rule for determining when to terminate the
microequilibration 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 longtime 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 pathspace observables of stochastic
dynamics in terms of new goaloriented 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 gradientfree. 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 timeindependent
Schroedinger equation, with matrix valued potentials, and the values of
observables for ab initio BornOppenheimer 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 spacetime 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 coarsegraining 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 nonlinear coarse graining
mappings (e.g., reaction coordinates, endtoend length of chains).
Furthermore, we study the equivalence of force matching with relative entropy
minimization which we derive for general nonlinear 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 coarsegraining of {\it nonequilibrium} molecular
systems. In this context, we propose error estimation and controlledfidelity
model reduction methods based on PathSpace Information Theory, and combine it
with statistical parametric estimation of rates for nonequilibrium stationary
processes. The approach we propose extends the applicability of existing
informationbased methods for deriving parametrized coarsegrained models to
NonEquilibrium systems with Stationary States (NESS). In the context of
coarsegraining it allows for constructing optimal parametrized Markovian
coarsegrained dynamics, by minimizing information loss (due to
coarsegraining) on the path space. Furthermore, the associated pathspace
Fisher Information Matrix can provide confidence intervals for the
corresponding parameter estimators. We demonstrate the proposed coarsegraining
method in a nonequilibrium system with diffusing interacting particles, driven
by outofequilibrium boundary conditions.

In this paper we study from a numerical analysis perspective the Fractional
Step Kinetic Monte Carlo (FSKMC) algorithms proposed in [1] for the parallel
simulation of spatially distributed particle systems on a lattice. FSKMC are
fractional step algorithms with a timestepping 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 FSKMC algorithms as approximations of
conventional, serial kinetic Monte Carlo (KMC). A key aspect of our analysis
relies on emphasising a goaloriented 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 multilevel kinetic Monte Carlo methods for
sampling highdimensional, 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
coarsegrained dynamics. Microscopic reconstruction corrects possibly
significant errors introduced through coarsegraining, leading to the
controllederror approximation of the sampled stochastic process. In this
manner, the proposed multilevel algorithm overcomes known shortcomings of
coarsegraining of particle systems with complex interactions such as combined
long and shortrange particle interactions and/or complex lattice geometries.
Specifically, we provide error analysis for the approximation of longtime
stationary dynamics in terms of relative entropy and prove that information
loss in the multilevel 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
longtime sampling, obtaining phase diagrams and in studying metastability
properties of highdimensional complex systems. Finally, the multilevel nature
of the method provides flexibility in combining rejectionfree and nullevent
implementations, generating a hierarchy of algorithms with an adjustable number
of rejections that includes wellknown rejectionfree and nullevent
algorithms.

BornOppenheimer dynamics is shown to provide an accurate approximation of
timeindependent 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 HamiltonJacobi
equation, bypasses the usual separation of nuclei and electron wave functions,
includes caustic states and gives a different perspective on the
BornOppenheimer 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 controllederror 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 spatiotemporal scales
in spatially distributed, nonequilibrium physiochemical processes with complex
chemistry and transport micromechanisms. The algorithms can be tailored to
specific hierarchical parallel architectures such as multicore processors or
clusters of Graphical Processing Units (GPUs). The proposed parallel algorithms
are controllederror 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 timewindow.
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
quasiequilibrium 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 ensembleaveraged 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 HamiltonJacobi equation satisfied by the value function. In
the nearequilibrium regime, or under a local quadratic approximation in the
farfromequilibrium regime, this bestfit 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
timescale 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 shortrange particle interactions. The proposed method is a
Metropolistype algorithm with the proposal probability transition matrix based
on the coarsegrained approximating measures introduced in a series of works of
M. Katsoulakis, A. Majda, D. Vlachos and P. Plechac, L. ReyBellet 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 Isingtype model, comparing efficiency of the
single spinflip Metropolis dynamics and the proposed coupled Metropolis
algorithm.

An implicit massmatrix penalization (IMMP) of Hamiltonian dynamics is
proposed, and associated dynamical integrators, as well as sampling MonteCarlo
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 tradeoff between stability and dynamical
modification. Moreover, a penalty that vanishes with the timestep 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
Metropolislike 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 timestep 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 massmatrix 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 tradeoff 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 nonconvex
interactions. Prescribing the macroscopic dynamical timescale, it is shown that
the IMMP method increases the timestep 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 coarsegraining schemes for
stochastic manybody 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
Isingtype models. The proposed algorithms are derived from an initial
coarsegrained 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 coarsegrained
equilibrium measures. Subsequently we carry out a cluster expansion around this
firstand often inadequateapproximation and obtain more accurate
coarsegraining schemes. The cluster expansions yield also sharp a posteriori
error estimates for the coarsegrained approximations that can be used for the
construction of adaptive coarsegraining methods. We present a number of
numerical examples that demonstrate that the coarsegraining schemes developed
here allow for accurate predictions of critical behavior and hysteresis in
systems with intermediate and longrange interactions. We also present examples
where they substantially improve predictions of earlier coarsegraining schemes
for shortrange interactions.