Data-driven parameterization of the generalized Langevin equation
We present a data-driven approach to determine the memory kernel and random noise in generalized Langevin equations. To facilitate practical implementations, we parameterize the kernel function in the Laplace domain by a rational function, with coefficients directly linked to the equilibrium statistics of the coarse-grain variables. We show that such an approximation can be constructed to arbitrarily high order and the resulting generalized Langevin dynamics can be embedded in an extended stochastic model without explicit memory. We demonstrate how to introduce the stochastic noise so that the second fluctuation-dissipation theorem is exactly satisfied. Results from several numerical tests are presented to demonstrate the effectiveness of the proposed method.
Generalized Langevin equations (GLEs) have recently re-emerged in the area of molecular modeling as a promising description of reduced-dimension coarse-grained variables. In principle, GLEs can be derived using the Mori-Zwanzig projection formalism Mori (1965); Zwanzig (1973). Examples of such derivations can be found for a variety of applications Berne et al. (1990); Chorin and Stinis (2005); Curtarolo and Ceder (2002); Hijón et al. (2010); Izvekov and Voth (2006); Li (2010); Ma et al. (2016); Givon et al. (2004); e.g., climate modeling Wouters and Lucarini (2013); Majda and Harlim (2012). However, practical implementations of GLEs require specification of the memory function which can be difficult to obtain, even when the full dynamics of the system is known. For example, the memory functions obtained in past studies Li (2010); Ma et al. (2016); Chen et al. (2014); Darve et al. (2009) have involved functions of high-dimensional matrices. Darve et al. proposed a more efficient algorithm Darve et al. (2009) to compute the memory kernel by solving an equation for the orthogonal dynamics derived from the Mori-Zwanzig formalism. However, the orthogonal dynamics equation can be expensive to solve when the original system is large. Furthermore, even when the memory kernel function is available, direct evaluation of the memory term can be costly because it requires the history of the coarse-grained (CG) variables at every time step and the associated numerical integration. Sampling of the random noise is also a challenging component of GLEs: to generate the correct equilibrium statistics for the CG model, the random noise has to obey the second fluctuation-dissipation theorem (FDT) Kubo (1966). The theory of stationary processes Doob (1944) states that the random process is uniquely determined by the correlation function, which is proportional to the memory kernel; however, sampling the random noise is nontrivial in practice. Methods based on matrix factorization are computationally challenging because they require decomposition of a correlation matrix with dimension proportional to the total simulation period. Alternatively, more efficient methods based on Fast Fourier Transforms (FFTs) may create artificial periodicity Li et al. (2015).
In addition to the direct derivation of memory kernels Li (2010); Ma et al. (2016); Chen et al. (2014); Darve et al. (2009), there have been numerous attempts to compute the memory kernel from full molecular dynamics (MD) simulations Oliva et al. (2000); Berkowitz et al. (1981, 1983), especially for systems with zero net mean force. Such analyses lead to integral equations of the first kind, which are numerically unstable without additional regularization. Another approach for estimating the kernel uses Kalman filtering and assumes functions of exponential form Harlim and Li (2015); Fricks et al. (2009) such that the GLE can be embedded in a Markovian dynamics framework. In recent work, Chorin and Lu Chorin and Lu (2015) considered a time-discrete representation, representing the memory effects using the NARMAX (nonlinear autoregression moving average with exogenous input) method. Voth et al. proposed an alternative approach Davtyan et al. (2015) to recover the coarse-grained dynamics by introducing fictitious particles that interact with the coarse-grained variables, effectively introducing an approximation of the kernel function.
In this work, we present a hierarchical approach to obtain GLE kernel functions from simulation data. Such a data-driven approach is particularly useful for complex models (e.g., biomolecules or climate) in which full dynamics models are typically unavailable or inaccessible with finite computing resources. The key idea is parameterization of the kernel function Laplace transform by a rational approximation. The first two approximations in the rational approximation hierarchy correspond to the Markovian approximation Hijón et al. (2010); Lei et al. (2010); Kauzlarić et al. (2011) and approximation of noise by a Ornstein-Uhlenbeck process, which is the ansatz used in some previous works. As we will show, these two approximations are often insufficient to predict dynamics properties; however, our hierarchy can be used to construct arbitrarily high-order models. The parameters in our ansatz can be computed from statistical properties of the coarse-grained variables. Additionally, this ansatz makes it possible to eliminate memory from the GLE by introducing auxiliary variables. In particular, we will show how to introduce inexpensive white noise terms into the extended dynamics to approximate the random noise while satisfying the second FDT is exactly. As a result, no memory term needs to be evaluated, and no colored noise needs to be sampled.
Ii Numerical Method
The generalized Langevin equation (GLE) can be expressed in the following form
where and are the coarse-grained coordinates and momenta, and is the conservative force term determined by the potential of mean force (PMF) . denotes a memory kernel function, which is the main focus of this paper. The noise is a stationary Gaussian process with zero mean, satisfying the second fluctuation-dissipation theorem (FDT) Kubo (1966):
where . Such equations have been derived in previous work Li (2010); Chen et al. (2014) using the Mori-Zwanzig formalism Mori (1965); Zwanzig (1961). In the present work, our main goal is to estimate the memory kernel and noise terms from dynamics simulation data.
We assume that we have a time series data set of and () drawn from an equilibrium simulation such that the time series corresponds to a stationary random process. We right-multiply the second equation in (1) by to obtain
Here we have defined the correlation matrices,
and we have made the assumption that , which was verified in previous work Chen et al. (2014).
Given the correlation functions, Eq. (3) can be regarded as an integral equation from which the memory function can be computed. However, this is an integral equation of the first kind, and it is not well-posed, leading to unreliable solutions. Instead of a determining the kernel function directly in the time-domain, we can instead parameterize its Laplace transform. We define the Laplace transform as
Similarly, we denote the Laplace transforms of and by and , respectively. Taking the Laplace transform of Eq. (3), we arrive at
We utilize the values of at specific time points to construct . By taking , we obtain
It is clear that
which makes it possible to incorporate a Green-Kubo type of formula in the approximation and model accuracy over long time scales.
For short or intermediate time scales, we use the point . Using (6), one can find the limiting values of the kernel and its derivatives as . This calculation amounts to computing and similarly , which is straightforward. For instance, by integrating by parts repeatedly in (5), we find that
For example, we have . In addition, we can find
In the derivations above, we have incorporated the values at and But the ansatz of the rational approximation is quite flexible, and other interpolation points can be used as well. For stationary process of Hamiltonian system, ; therefore, and .
Given limiting values available extracted from the data, we seek a rational function approximation for , in the form of
The coefficients can be determined by matching the limits of . The matching conditions lead to a linear system of equations, which can be solved analytically for small or numerically for large .
The zero-order approximation treats as a constant matrix set to . Accordingly, one gets a Markovian approximation by a Langevin dynamics with damping coefficient given by . We can determine by Eq. (7). In fact, is proportional to the diffusion tensor; i.e., the matching condition recovers the Einstein relation Kubo (1966); Einstein (1905).
For the first order approximation (), we have
In this case, the memory function in the time domain is given by,
Depending on the eigenvalues of , the memory function can exhibit both exponential decay and oscillations.
A computational difficulty in GLE simulation is that the integral has to be evaluated at every step. However, this difficulty can be removed by introducing auxiliary equations based on the rational approximation of the memory function (11). More specifically, we can define Then, Eq. (12) implies that the approximate GLE can be written as,
where we added a white noise term satisfying . We pick the initial state of the auxiliary variable to satisfy .
We can show that this new memory-less dynamics corresponds to an approximation of the GLE (1). The memory function is approximated by the rational function in the frequency space, which is precisely (13). More importantly, with this proper choice of the initial condition for , the random noise is given by which is a stationary Gaussian process that satisfies the second FDT (2) exactly with an invariant distribution given by
The procedure above can be extended to arbitrarily high order, and the extended system can be written as follows,
where is a matrix and is added white noise. For example, for the fourth-order method, is given by:
The matrices and can be determined by matching the Laplace transform of the memory kernel given by with the rational approximation (11). Similar to Eq. (14), we can also show that, by choosing the white noise and the initial conditions of properly; i.e., , , the colored noise generated by the extended dynamics also satisfies the second FDT (2) exactly. Eq. (16) has an invariant distribution given by
Iii Numerical Results
We demonstrate our method through coarse-graining the dynamics of a tagged particle within a one-dimensional harmonic chain. In particular, we consider the first particle (on the free end) as the target particle and treat the remaining particles as the heat bath. As shown in Ref. Adelman and Doll (1974); Li and E (2007), the dynamics of the target particle can be modeled by a GLE with kernel given by , where is the mass of each particle, is the force constant of harmonic interaction and is a Bessel function of the first kind. We calculated the trajectory of the tagged particle in a harmonic chain consisting of particles with , and set to unity. Fig. 1 shows the numerical results for obtained using the rational approximation method described above. As the approximation increases to third order, agrees well with the exact formulation, which is given by .
We also studied a tagged particle immersed in a fluid system governed by pairwise-conservative forces similar to those in dissipative particle dynamics (DPD) simulations Hoogerbrugge and Koelman (1992); Espanol and Warren (1995); i.e., defined by
where , and and is the force magnitude and is the cut-off radius beyond which all interactions vanish. One tagged particle and solvent particles were placed in a cubic box , with a periodic boundary condition imposed in each direction. The mass of both the tagged particle and solvent particle were chosen to be unity. A Nosé-Hoover thermostat was used to equilibrate the system.
Following the Mori-Zwanzig method, the dynamics of the target particle can be modeled by a GLE with zero mean force. We considered two cases: (I) , and (II) , . Based on the velocity correlation function obtained from MD simulation data, we can compute the different orders of kernel approximation by Eq. (11). As shown in Fig. 2, the exact kernel function agrees well with the numerical result directly obtained from as the approximation order increases to second order or above.
Unlike case (I), in case (II) shows a pronounced peak near , indicating significant oscillations in the time domain of . Fig. 3 shows the velocity correlation function and the mean-squared displacement obtained from solving Eq. (1) with kernel constructed by different orders of rational approximations. The zero-order approximation corresponds to the Markovian approximation of the kernel term
In case (I), Eq. (1) can reproduce the dynamics of the system very well for rational approximations of second order (and above). In contrast, for case (II), obtained from the second-order approximation exhibits artificial oscillation and deviations from the dynamics results; the third and fourth order approximations yield much better agreement. The different performance for cases (I) and (II) can be understood by examining the time scale separation of and in the memory term . As shown in Fig. 4, the velocity correlation of case (I) decays much more slowly than for case (II). The plateau region of of case(I) shown in Fig. 2 illustrates the similarities between and . These similarities are why the Markovian approximation by Eq. (20) yields fairly good agreement for case(I). In contrast, there is no apparent time scale separation between and for case(II), explaining the need for higher-order approximations to characterize the coupling between and .
To further demonstrate our method, we simulated the transition rate of a tagged particle in a double-well potential using both by GLE and full molecular dynamics. The double-well potential had the form
where and refer to the depth and width of the potential field. The tagged particle interacted with solvent particles through interactions defined in Eq. (19) with and , interactions which drove the transitions between the two energy minima at and . The instantaneous transition rate was computed by
where is a collective variable defining the dividing iso-surface between the states and is the characteristic function for the reaction trajectory. For the double-well system, and , where is the Heaviside function. is the reaction partition function over the phase space of state . The transition flux correlation function is given by , where is the reaction rate predicted from transition state theory. Fig. 5 shows obtained from the full molecular dynamics the GLE systems with kernels modeled by different orders of rational approximation. The zero-order approximation (corresponding to Langevin dynamics with a constant friction coefficient) yielded a which deviated from the full molecular dynamics result, indicating a pronounced non-Markovian effect. In contrast, obtained from the GLE using the third- and fourth-order rational approximations agree well with the full molecular dynamics results. Moreover, as shown in the inset of Fig. 5, both transition state theory and the zero-order (Langevin dynamics) model overestimate the reaction rate; approaches the full molecular dynamics results only when using the third- and fourth-order approximations where the memory kernel can be more accurately modeled.
We have presented a new data-driven approach to obtain the memory kernel for a GLE through rational-function approximation in the Laplace transform domain. The data-driven nature of this method arises through connection of the rational function coefficients with equilibrium statistics of the coarse-grained variables – statistics that can be calculated through simulation time series data. The zero- and first-order approximations recover the Langevin and Ornstein-Uhlenbeck stochastic processes, respectively. Higher-order approximations have also been systematically derived for systems with significant memory effects. Unlike the time-domain kernel function representation, numerical simulations of the GLE using the rational approximation can be conveniently implemented by introducing auxiliary variables that follow linear stochastic dynamics with no memory, eliminating the need for expensive calculations of history-dependent memory terms. We have also shown that the second FDT (2) can be satisfied automatically using our approach. This method has been tested with simple systems but applicable to much more complicated biological and material systems with pronounced memory effects; e.g., sub-diffusion in single-molecule measurements Kou and Xie (2004); Min et al. (2005) or transition dynamics of chemical and biological reaction systems Schenter et al. (1992); Lange and Grubmüller (2006). The rational function approximations (Eq. (14) and Eq. (16)) of the GLE enable us to analyze the transition dynamics via the extended dynamics in an augmented phase space , where direct analysis via the GLE could be difficult/inaccessible.
Acknowledgements.We thank George Karniadakis, Eric Darve, Panos Stinis, Greg Schenter, Lei Wu, Dave Sept, J. Andrew McCammon, and Zhen Li for helpful discussion. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research as part of the Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4).
- Mori (1965) H. Mori, Prog. Theor. Phys. 33, 423 (1965).
- Zwanzig (1973) R. Zwanzig, J. Stat. Phys. 9, 215 (1973).
- Berne et al. (1990) B. Berne, M. Tuckerman, J. E. Straub, and A. Bug, The Journal of Chemical Physics 93, 5084 (1990).
- Chorin and Stinis (2005) A. J. Chorin and P. Stinis, Comm. Appl. Math. Comp. Sc. 1, 1 (2005).
- Curtarolo and Ceder (2002) S. Curtarolo and G. Ceder, Phys. Rev. Lett. 88 (2002).
- Hijón et al. (2010) C. Hijón, P. Español, E. Vanden-Eijnden, and R. Delgado-Buscalioni, Faraday discuss. 144, 301 (2010).
- Izvekov and Voth (2006) S. Izvekov and G. A. Voth, J. Chem. Phys. 125, 151101 (2006).
- Li (2010) X. Li, Int. J. Numer. Meth. Engng. 83, 986 (2010).
- Ma et al. (2016) L. Ma, X. Li, and C. Liu, Arxiv preprint , 1605.04886 (2016).
- Givon et al. (2004) D. Givon, R. Kupferman, and A. Stuart, Nonlinearity 17, R55 (2004).
- Wouters and Lucarini (2013) J. Wouters and V. Lucarini, Journal of Statistical Physics 151, 850 (2013).
- Majda and Harlim (2012) A. J. Majda and J. Harlim, Nonlinearity 26, 201 (2012).
- Chen et al. (2014) M. Chen, X. Li, and C. Liu, J. Chem. Phys. 141, 064112 (2014).
- Darve et al. (2009) E. Darve, J. Solomon, and A. Kia, Proc. Natl. Acad. Sci. 106, 10884 (2009).
- Kubo (1966) R. Kubo, Rep. Prog. Phys. 29(1), 255 (1966).
- Doob (1944) J. L. Doob, Ann. Math. Stat. 15, 229 (1944).
- Li et al. (2015) Z. Li, X. Bian, X. Li, and G. E. Karniadakis, The Journal of chemical physics 143, 243128 (2015).
- Oliva et al. (2000) B. Oliva, X. Daura, E. Querol, F. X. Avilés, and O. Tapia, Theor. Chem. Acc. 105, 101 (2000).
- Berkowitz et al. (1981) M. Berkowitz, J. D. Morgan, D. J. Kouri, and J. A. McCammon, J. Chem. Phys. 75, 2462 (1981).
- Berkowitz et al. (1983) M. Berkowitz, J. Morgan, and J. A. McCammon, J. Chem. Phys. 78, 3256 (1983).
- Harlim and Li (2015) J. Harlim and X. Li, Physical Review E 91, 053306 (2015).
- Fricks et al. (2009) J. Fricks, L. Yao, T. C. Elston, and M. G. Forest, SIAM journal on applied mathematics 69, 1277 (2009).
- Chorin and Lu (2015) A. J. Chorin and F. Lu, Proceedings of the National Academy of Sciences 112, 9804 (2015).
- Davtyan et al. (2015) A. Davtyan, J. F. Dama, G. A. Voth, and H. C. Andersen, J. Chem. Phys. 142, 154104 (2015).
- Lei et al. (2010) H. Lei, B. Caswell, and G. E. Karniadakis, Phys. Rev. E 81, 026704 (2010).
- Kauzlarić et al. (2011) D. Kauzlarić, J. T. Meier, P. Español, S. Succi, A. Greiner, and J. G. Korvink, J. Chem. Phys. 134, 064106 (2011).
- Zwanzig (1961) R. Zwanzig, Lectures in Theoretical Physics 3, 106 (1961).
- Einstein (1905) A. Einstein, Annalen der Physik 17, 549 (1905).
- Adelman and Doll (1974) S. A. Adelman and J. D. Doll, J. Chem. Phys. 61, 4242 (1974).
- Li and E (2007) X. Li and W. E, Phys. Rev. B 76, 104107 (2007).
- Hoogerbrugge and Koelman (1992) P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992).
- Espanol and Warren (1995) P. Espanol and P. Warren, Europhysics Letters 30, 191 (1995).
- Kou and Xie (2004) S. C. Kou and X. S. Xie, Phys. Rev. Lett. 93, 180603 (2004).
- Min et al. (2005) W. Min, G. Luo, B. J. Cherayil, S. C. Kou, and X. S. Xie, Phys. Rev. Lett. 94, 198302 (2005).
- Schenter et al. (1992) G. K. Schenter, R. P. McRae, and B. C. Garrett, J. Chem. Phys. 97, 9116 (1992).
- Lange and Grubmüller (2006) O. F. Lange and H. Grubmüller, J. Chem. Phys. 124, 214903 (2006).