Predicting the unobserved: A statistical mechanics framework for non-equilibrium material response with quantified uncertainty
Introduction
Predicting the far-from-equilibrium response of materials is of great scientific and practical interest. Any dynamically-loaded system, from an RNA strand being stretched by optical tweezers to steel being cold-rolled, is driven away from equilibrium. Moreover, one may be interested in equilibrium and non-equilibrium quantities like free energies or rheological/flow properties, and wish to infer them from experiments or simulations of potentially different systems or loading conditions. On the simulation side, molecular dynamic simulations have a computational cost and associated time scale limitations that currently precludes the direct exploration of material behavior at low strain rates (Yan and Sharma, 2016). Similarly, it is difficult in many dissipative experimental systems to probe the free energy landscape by evaluating the work at infinitesimally slow loading rates (Collin et al., 2005). Ultimately, extrapolations of material behavior over loading conditions or material systems is crucial for the inverse problem of material design.
Attempting to make predictions about one system based on knowledge of another has a long history: approximating the system of interest as a perturbation on a more easily-soluble one is a standard technique across mechanics and physics (Hashin and Shtrikman, 1963, Zwanzig, 1954), and the idea finds applications in equilibrium free energy calculations with techniques such as thermodynamic integration (Rickman and LeSar, 2002), umbrella sampling (Torrie and Valleau, 1977) or metadynamics (Laio and Parrinello, 2002, Laio and Gervasio, 2008). Other approaches used to infer equilibrium information from non-equilibrium data, coming from either simulations or experiments, include fluctuation theorems, such as the Jarzynski relation (Jarzynski, 1997) or Crooks’ fluctuation theorem (Crooks, 1999). More recently, in the context of simulating non-equilibrium behavior, hyperdynamics methods (Voter, 1997, Kim et al., 2013) have emerged as a means to accelerate the simulation of rare events and alleviate the time-scale bottleneck. These may be abstractly viewed as non-equilibrium analogues of the previous approaches, where the potential is modified to have a faster-evolving system, e.g., by converting a large energy barrier into a smaller one, and the results are then corrected to emulate those of the original system. These approaches are distinct from methods such as transition path sampling (Bolhuis et al., 2002), forward flux sampling (Allen et al., 2009), milestoning (Faradjian and Elber, 2004), parallel replica dynamics (Voter, 1998), and parallel trajectory splicing (Perez et al., 2016), which leave the potential untouched and extend the temporal reach of the simulation via sophisticated algorithms that preserve the statistics of the original system. For a review on accelerated molecular dynamics methods we refer the reader to Perez et al. (2009) and Voter et al. (2002).
In 2007, an exact correspondence between trajectory probabilities in two different systems was derived and exploited (Chen and Horing, 2007, Nummela and Andricioaei, 2007 – see also earlier works with similar ideas, e.g., Zuckerman and Woolf (1999)) to accelerate the sampling of trajectories that overcome both energetic barriers (e.g., surface diffusion in a crystal) and entropic barriers (e.g., nanopore traversal by a polymer; see also Shin et al. (2010)). The fundamental idea in each case is this: add a bias to the system potential , perform simulations of the system governed by , and use the exact relations (discussed in detail below) to reweight the probability of the simulated trajectories so that they are statistically equivalent to those of the original system governed by . The potential can be highly general, and could for example be selected to “fill in” a deep potential well in which the system is trapped, or to point the system towards a particular location of interest in space. This is distinct from other accelerated molecular dynamics formulations as the reweighting is exact for each trajectory, it does not require an a priori knowledge of the reaction coordinates, it could potentially be applied to experimental data, and no time rescaling is performed. There is also no requirement that the system’s evolution can be characterized by long sojourns in deep potential wells, punctuated by occasional transitions between them, meaning that path integral hyperdynamics can be applied even to systems with no well-defined transition state, such as the energetic barriers mentioned above. Furthermore, we are interested not only in accelerating rare event sampling (though this provides an intuitive example), but more generally in the statistical connections between distinct systems. Rather than accelerating the simulation of a system’s evolution, we focus on the computation of ensemble-average quantities, arbitrarily far from equilibrium. The method’s principal shortcoming lies in the reweighting of the trajectories: adding too extreme a bias will destroy the statistics of the quantity one is trying to compute for a fixed computational cost (Ikonen et al., 2011); and this issue limits in practice the approach to systems of relatively small size and short time simulations. In other words, choosing a large bias will greatly accelerate the sampling of rare trajectories, at the cost of introducing great uncertainty in the statistical conclusions. Here, our interest lies not only in accelerated sampling strategies, but also in understanding non-equilibrium material behavior, and its connection to equilibrium properties. We extend the path integral hyperdynamics approach to include time-dependent, out-of-equilibrium boundary conditions, of interest in mechanics, and fully quantify, for the first time, the statistical uncertainty that the biasing procedure introduces.
The paper is organized as follows. In Section 2 below, we give a detailed derivation of the procedure to obtain expectation values of a given statistical observable in system from data generated by stochastic simulations of system . Next, in Section 3, we quantify the uncertainty introduced by the biasing procedure in terms of the difference of interparticle potentials and loading conditions between the target () and simulated () systems/processes. These results are then exemplified in Section 4 over a 1D mass–spring chain subjected to a non-equilibrium, time-dependent boundary condition. While this example is simple, it has proved useful as a prototype for polymeric chains and biological macromolecules, and it is here used to showcase the various possibilities offered by the formalism. In particular, we firstly investigate element transmutation, in which systems and have different interatomic potentials. Secondly, we consider going from equilibrium to non-equilibrium conditions for the same material, and finally, we take the method to its extreme by letting system be a simple Brownian motion, and a nonlinear interacting chain with time-dependent boundary conditions. A second physical system is explored in Section 5, where the caging in two-dimensional glassy systems is predicted from the liquid state. Finally, Section 6 contains our conclusions, and several appendices give further technical details of the calculations for completeness.
Section snippets
Predicting the non-equilibrium behavior of a material/process from that of a reference one
The goal of the present manuscript is to predict the evolution of an observable (e.g., instantaneous stress, work done on the system, diffusion coefficient) for a material and loading conditions , from that of . This general objective thus includes, as particular cases, virtual element transmutation for given loading conditions (e.g., response for different interparticle potentials), as well as different loading conditions for a given material (e.g., predicting the response to an
Uncertainty quantification
From Eq. (4) it is immediate to see that . Although this identity is exactly satisfied for the exact probability density, an empirical average, i.e., , where is the number of realizations, may deviate from 1 if insufficient samples are used. In general, more realizations will be needed the further apart systems/processes and are, as noted in the previous section. Such a distance can be measured by means of the Kullback–Leibler divergence, or relative
Model description and overview of the cases to be examined
The first example that we examine is a prototype for polymers (Doi and Edwards, 1988) and biological macromolecules such as DNA and coiled-coil proteins (Raj and Purohit, 2011). In particular, the model consists of particles following Langevin dynamics, connected through identical springs, as shown in Fig. 4. The first spring is fixed to a wall, while the last one has a prescribed displacement boundary condition of the form , where is a constant pulling velocity. Denoting
Example 2: Caging in two-dimensional glassy systems
One of the most ubiquitous examples of out-of-equilibrium behavior, and one that we are regularly familiar with from everyday experience, is that of glasses (Stillinger and Debenedetti, 2013, Charbonneau et al., 2017). Glassy dynamics is observed in a wide range of length scales (from nanoparticles to grains) spanning a variety of industries from building materials (concrete) to paints (colloidal suspensions) and household goods (foams & gels) (Bonn et al., 2017, Nicolas et al., 2018). Despite
Conclusions
In this paper, we provide a statistical mechanics approach with quantified uncertainty to extrapolate material behavior to distinct loading conditions or material systems. The approach is based on reweighting the probability density for trajectories, building up on the ideas of Chen and Horing (2007), and enables the calculation of ensemble averages of arbitrary observables of system/process from simulations (or, potentially, experimental data) of system/process . The formalism is a priori
Code availability
The codes used to generate the results here presented will be made available upon request to the corresponding author.
CRediT authorship contribution statement
Shenglin Huang: Conceptualization, Methodology, Investigation, Software, Validation, Formal analysis, Writing – original draft, Writing – review & editing. Ian R. Graham: Conceptualization, Investigation, Software, Validation, Writing – original draft, Writing – review & editing. Robert A. Riggleman: Writing – review & editing, Funding acquisition. Paulo E. Arratia: Writing – review & editing, Funding acquisition.. Steve Fitzgerald: Conceptualization, Methodology, Writing – original draft,
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
The authors gratefully acknowledge support from ACS, USA PRF 61793-ND10 Grant (S. H.) and from NSF CAREER Award, CMMI-2047506 (C. R). S. F. acknowledges support from the UK EPSRC , grant no. EP/R005974/1, and helpful discussions with Profs. T. Ala-Nissilä and A. Archer. I. G., R. A. R., and P. E. A. acknowledge support from NSF DMR (Penn MRSEC) 1720530.
References (43)
- et al.
A variational approach to the theory of the elastic behaviour of multiphase materials
J. Mech. Phys. Solids
(1963) - et al.
Harnessing fluctuation theorems to discover free energy and dissipation potentials from non-equilibrium data
J. Mech. Phys. Solids
(2021) - et al.
Exact low-force kinetics from high-force single-molecule unfolding events
Biophys. J.
(2007) - et al.
Accelerated molecular dynamics methods: Introduction and recent developments
Ann. Rep. Comput. Chem.
(2009) - et al.
Phase boundaries as agents of structural change in macromolecules
J. Mech. Phys. Solids
(2011) - et al.
Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling
J. Comput. Phys.
(1977) - et al.
Forward flux sampling for rare event simulations
J. Phys.: Condens. Matter
(2009) - et al.
Theoretical perspective on the glass transition and amorphous materials
Rev. Modern Phys.
(2011) - et al.
Transition path sampling: Throwing ropes over rough mountain passes, in the dark
Ann. Rev. Phys. Chem.
(2002) - et al.
Yield stress materials in soft condensed matter
Rev. Modern Phys.
(2017)