Path Space Markov Chain Monte Carlo Methods for Molecular Simulation Citation Rosa-Raíces, Jorge Luis (2022) Path Space Markov Chain Monte Carlo Methods for Molecular Simulation. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/10jr-hg67. https://resolver.caltech.edu/CaltechTHESIS:03292022-234015177 Abstract Path space Markov-chain Monte Carlo (McMC) provides a versatile framework for simulating the structure and dynamics of condensed-phase systems aptly described by classical and quantum Boltzmann statistics. This thesis comprises our efforts to design, analyze and improve path space McMC algorithms to achieve numerically advantageous, and physically accurate, simulation of molecular processes across a range of scales. To improve molecular dynamics (MD) simulations of atomically resolved systems exhibiting pronounced nuclear quantum effects, we introduce a family of integrators for non-preconditioned path-integral MD exhibiting dimension-free statistical accuracy and efficiency, and enabling a many-fold increase in time-step stability relative to conventional approaches at no additional computational cost or implementation complexity. The integrators come with robust performance guarantees that are borne out in thermostatted ring-polymer MD simulations of realistic condensed-phase models. Concurrently, toward extending the range of accessible timescales in stochastic MD simulations of mesoscale coarse-grained molecular systems, we introduce a parallel-in-time integrator for the overdamped Langevin equation based on McMC evaluation of a path-integral representation of the many time-step stochastic MD transition kernel. The parallel-in-time integrator achieves simultaneous integration of multiple stochastic MD time-steps at no greater wall-time cost and with no lesser accuracy than a standard Euler--Maruyama integrator does in serial, and thus instantiates new opportunities to accelerate stochastic dynamics simulations on massively parallel computer architectures. Our work along these two methodological avenues extends the utility of path space McMC across applications in molecular simulation and has broader implications in other disciplines that require accurate and efficient simulations of Markov diffusion processes in state spaces or path spaces. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Markov chain Monte Carlo, path integrals, molecular dynamics integrators, parallel-in-time integration, dimension-free accuracy Degree Grantor: California Institute of Technology Division: Chemistry and Chemical Engineering Major Option: Chemistry Thesis Availability: Public (worldwide access) Research Advisor(s): Miller, Thomas F. (advisor) Stoltz, Brian M. (co-advisor) Thesis Committee: Wang, Zhen-Gang (chair) Okumura, Mitchio Miller, Thomas F. Stuart, Andrew M. Stoltz, Brian M. Defense Date: 28 April 2022 Funders: Funding Agency Grant Number Department of Energy (DOE) DE-FOA-0001912 Department of Energy (DOE) DE-AC02-05CH11231 Department of Energy (DOE) DE-SC0019390 Office of Naval Research (ONR) N00014-10-1-0884 NSF DMS-1816378 NIH R01GM125063 Record Number: CaltechTHESIS:03292022-234015177 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:03292022-234015177 DOI: 10.7907/10jr-hg67 Related URLs: URL URL Type Description https://doi.org/10.1063/1.5125455 DOI Path-accelerated molecular dynamics (Published article in J. Chem. Phys.) https://doi.org/10.1063/1.5134810 DOI Dimension-free non-preconditioned path-integral molecular dynamics (Published article in J. Chem. Phys.) https://doi.org/10.1063/5.0036954 DOI Generalized integrators for dimension-free thermostatted ring-polymer molecular dynamics integrators (Published article in J. Chem. Phys.) http://dx.doi.org/10.22002/D1.20111 DOI Supplementary data to Path-accelerated molecular dynamics article (DOI: 10.1063/1.5125455) ORCID: Author ORCID Rosa-Raíces, Jorge Luis 0000-0003-2311-2948 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 14534 Collection: CaltechTHESIS Deposited By: Jorge Rosa Deposited On: 17 May 2022 18:29 Last Modified: 08 Nov 2023 00:39 Thesis Files PDF - Final Version See Usage Policy. 7MB Repository Staff Only: item control page

Path space Markov chain Monte Carlo methods for molecular simulation Thesis by Jorge L. Rosa-Raíces In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Chemistry CALIFORNIA INSTITUTE OF TECHNOLOGY Pasadena, California 2022 Defended April 28, 2022 ii © 2022 Jorge L. Rosa-Raíces ORCID: 0000-0003-2311-2948 All rights reserved except where otherwise noted iii Acknowledgments This thesis is the product of efforts and experiences well beyond academic in scope. Throughout my graduate studies I received copious support and encouragement through multiple avenues, and acknowledgments are due to those responsible. I am thankful to my advisor, Tom Miller, for his inspiring dedication as a scientist, educator, and leader. During my entire time as a graduate student, Tom provided a perfect balance between stalwart direction, kind encouragement, and truly unwavering support for our common and, most remarkably, my own academic and professional goals. I will forever be indebted to him for his patience and understanding, and for offering within his group a safe space where I could tumble while searching for a sense of belonging in the Caltech community. For graciously extending beyond his role as collaborator to provide mentorship and guidance in a similar capacity to Tom, I wholeheartedly thank Professor Nawaf BouRabee at Rutgers–Camden. I hope that our future interactions will continue to yield fruitful scientific exploits and foster a long-standing friendship. During my early years in the Miller group I benefited greatly from the kindness of senior colleagues Matt Welborn, Brett Savoie, Michiel Niesen, Ioan-Bogdan Magdău, and Philip Shushkov, all of whom, intentionally or otherwise, served as role models and offered invaluable career advice I hold dear to this day. Thanks are also due to my cohort of contemporary graduates in the group with whom I was able to grow, comprising Sebastian Lee, Matthew Zimmer, Jeongmin Kim, Xuecheng Tao, Dan Jacobson, Roman Korol, Sherry Cheng, Jiace Sun, Zhuoran Qiao, and Doris Mai; they all provided a pleasant work environment and a forum for interesting discussion of matters scientific and beyond, selflessly stepped in for others in difficult times, and motivated each other to strive further through their impressive accomplishments. It is great to see that our most recent members, namely Jorge Emiliano Deustua, Tomislav Begušić, Rui Li, Vignesh Bhethanabotla, and Bo Li, will carry on the group’s spirit by exhibiting these qualities. To all friends and colleagues, mentioned and otherwise, I wish only the best. Lastly, I am immensely grateful to my family for their love and unconditional support throughout my pursuit of higher education. There is no one more important to me than all of you, and I can only hope to repay any fraction of all you have done for me. Thank you. iv Abstract Path space Markov chain Monte Carlo (McMC) provides a versatile framework for simulating the structure and dynamics of condensed-phase systems aptly described by classical and quantum Boltzmann statistics. This thesis comprises our efforts to design, analyze, and improve path space McMC algorithms to achieve numerically advantageous, and physically accurate, simulation of molecular processes across a range of scales. To improve molecular dynamics (MD) simulations of atomically resolved systems exhibiting pronounced nuclear quantum effects, we introduce a family of integrators for non-preconditioned path-integral MD exhibiting dimensionfree statistical accuracy and efficiency, and enabling a many-fold increase in timestep stability relative to conventional approaches at no additional computational cost or implementation complexity. The integrators come with robust performance guarantees that are borne out in thermostatted ring-polymer MD simulations of realistic condensed-phase models. Concurrently, toward extending the range of accessible timescales in stochastic MD simulations of mesoscale coarse-grained molecular systems, we introduce a parallel-in-time integrator for the overdamped Langevin equation based on McMC evaluation of a path-integral representation of the many time-step stochastic MD transition kernel. The parallel-in-time integrator achieves simultaneous integration of multiple stochastic MD time-steps at no greater walltime cost and with no lesser accuracy than a standard Euler–Maruyama integrator does in serial, and thus instantiates new opportunities to accelerate stochastic dynamics simulations on massively parallel computer architectures. Our work along these two methodological avenues extends the utility of path space McMC across applications in molecular simulation and has broader implications in other disciplines that require accurate and efficient simulations of Markov diffusion processes in state spaces or path spaces. v Published Content and Contributions [3] J. L. Rosa-Raíces, J. Sun, N. Bou-Rabee, and T. F. Miller III, “A generalized class of strongly stable and dimension-free T-RPMD integrators”, The Journal of Chemical Physics 154, 024106 (2021) 10.1063/5.0036954. J.L.R.-R. participated in conceiving the project, participated in the formulation, analysis, and implementation of some of the methods introduced, prepared the simulations to illustrate the methods, and contributed to writing the article. [2] R. Korol, J. L. Rosa-Raíces, N. Bou-Rabee, and T. F. Miller III, “Dimensionfree path-integral molecular dynamics without preconditioning”, The Journal of Chemical Physics 152, 104102 (2020) 10.1063/1.5134810. J.L.R.-R. prepared the simulations to illustrate some of the methods introduced, participated in the formulation, analysis, and implementation of some of the methods introduced, and contributed to writing the article. [1] J. L. Rosa-Raíces, B. Zhang, and T. F. Miller III, “Path-accelerated stochastic molecular dynamics: Parallel-in-time integration using path integrals”, The Journal of Chemical Physics 151, 164120 (2019) 10.1063/1.5125455. J.L.R.-R. participated in designing and implementing the methods introduced, prepared the simulations to illustrate the methods, and contributed to writing the article. vi Table of Contents Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Published Content and Contributions . . . . . . . . . . . . . . . . . . . . . . v Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 I Dimension-free ring-polymer molecular dynamics Chapter 2: Thermostatted ring-polymer molecular dynamics integrators with dimensionality freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Non-preconditioned PIMD . . . . . . . . . . . . . . . . . . . . . . . 2.3 BCOCB avoids pathologies in the infinite-bead limit . . . . . . . . . 2.4 Consequences for the primitive kinetic energy expectation value . . . 2.5 Dimensionality freedom for OBCBO via force mollification . . . . . 2.6 Results for anharmonic oscillator . . . . . . . . . . . . . . . . . . . 2.7 Results for liquid water . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Appendix A: Other splittings . . . . . . . . . . . . . . . . . . . . . Chapter 3: Generalization and optimization of dimension-free integrators for thermostatted ring-polymer molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Appendix A: Necessary and sufficient condition for eigenvalues of a 2 × 2 real matrix to be inside the unit circle . . . . . . . . . . . . . . 3.6 Appendix B: Stability condition for harmonic external potentials . . . 3.7 Appendix C: Dimension-free quantitative contraction rate for harmonic external potentials in the infinite-friction limit . . . . . . . . . 3.8 Appendix D: Total variation bound on the equilibrium accuracy error for harmonic external potentials . . . . . . . . . . . . . . . . . . . . 3.9 Appendix E: Asymptotic variance of kinetic energy observables for harmonic external potentials in the infinite-friction limit . . . . . . . 3.10 Appendix F: Stability interval calibration for liquid water simulations 3.11 Appendix G: Simulation and estimation details . . . . . . . . . . . . 8 9 10 12 14 18 21 26 28 33 35 36 37 38 52 58 60 62 63 65 66 68 69 vii II Parallel-in-time stochastic molecular dynamics 81 Chapter 4: Path-accelerated stochastic molecular dynamics: Parallel-in-time integration using path integrals . . . . . . . . . . . . . . . . . . . . . . . 82 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3 Calculation details . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 viii List of Illustrations Number Page 2.1 Illustration of the consequences of dimensionality freedom, and lack thereof, to the statistical accuracy of non-preconditioned path-integral molecular dynamics (PIMD) simulations. . . . . . . . . . . . . . . . 19 2.2 Convergence with increasing path dimension of kinetic energy expectation values for a weakly anharmonic potential with standard and dimension-free non-preconditioned PIMD. . . . . . . . . . . . . . . 27 2.3 Stability with increasing time-step of kinetic energy expectation values for weakly and strongly anharmonic potentials, using standard and dimension-free non-preconditioned PIMD integrators. . . . . . . 29 2.4 Convergence with increasing path dimension of the expected kinetic energy per proton in a strongly anharmonic model of room-temperature liquid water, as per standard and dimension-free non-preconditioned PIMD integration schemes. . . . . . . . . . . . . . . . . . . . . . . 31 2.5 Stability with increasing time-step of the expected kinetic energy per proton in a strongly anharmonic model of room-temperature liquid water, as per standard and dimension-free non-preconditioned PIMD integration schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.6 Dynamical properties of liquid water within the thermostatted ringpolymer molecular dynamics (T-RPMD) model, computed using the standard, Bussi–Parinello “OBABO” and the dimension-free, Cayleymodified Leimkuhler “BAOAB” integration schemes. . . . . . . . . . 34 3.1 Spectra of generalized symplectic propagators for the jth mode of a free ring polymer. . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2 Spectral properties of the T-RPMD update for the free ring polymer for various choices of θ. . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3 Dimension-free convergence to equilibrium of BAOAB-like T-RPMD schemes with a harmonic external potential. . . . . . . . . . . . . . . 51 3.4 Performance at equilibrium of various BAOAB-like T-RPMD schemes applied to the one-dimensional quantum harmonic oscillator. . . . . . 54 3.5 Performance of various BAOAB-like T-RPMD schemes applied to q-TIP4P/F liquid water at room temperature. . . . . . . . . . . . . . 56 ix 3.6 3.7 3.8 3.9 3.10 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 Molecular dipole autocovariance function, corresponding infrared absorption spectrum, and molecular velocity autocovariance function in room-temperature liquid water for various BAOAB-like T-RPMD schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 General properties of the eigenspectrum of a symplectic ring-polymer mode propagator within a strongly stable T-RPMD scheme. . . . . . 61 Stability interval calibration for q-TIP4P/F room-temperature liquid water simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 Using the automatic windowing method to estimate integrated autocorrelation times of several observables of the one-dimensional harmonic oscillator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Thermalization of various BAOAB-like schemes applied to q-TIP4P/F water at room temperature. . . . . . . . . . . . . . . . . . . . . . . . 74 Illustration of a simple path-based MD integration scheme. . . . . . . 86 Illustration of a path-accelerated molecular dynamics (PAMD). . . . 87 The interplay of path length, shift length and number of sampling operations in determining the accuracy of the PAMD integration scheme. 89 PAMD integrates equilibrium trajectories by relaxing non-equilibrium path segments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Illustration of the multilevel sliding and sampling (S&S) algorithm for MC updates of the PAMD path configuration. . . . . . . . . . . . 95 Benchmarks of PAMD for an overdamped harmonic oscillator. . . . . 99 Hard-sphere schedule employed to precondition MC path sampling in the reported Lennard–Jones PAMD simulations. . . . . . . . . . . 101 Benchmarks of PAMD for a Lennard–Jones colloid. . . . . . . . . . 103 x List of Tables Number Page 4.1 Summary of notation employed to describe the PAMD integration scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.2 Summary of PAMD simulation parameters used for the harmonic oscillator application. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Summary of PAMD simulation parameters used for the application to the Lennard–Jones liquid. . . . . . . . . . . . . . . . . . . . . . . 104 1 Chapter 1 Introduction This thesis addresses the expansive field of molecular simulation, which comprises an ever-growing palette of computational tools for studying the structure and dynamics of molecular systems through the evaluation of principled models. Molecular simulation offers the prospect to observe molecular processes in greater spatiotemporal resolution than is possible with modern experimental techniques, prospectively yielding insights into the structure and dynamics of matter at length-scales ranging from the electronic to the macroscopic [1]; into macromolecular conformational transition rates and their underlying mechanisms [2]; into chemical reactions mediated by complex environments [3]; and into the dynamical response properties of materials in relation to their microstructural fluctuations. Contingent on sufficient modeling accuracy, as increasingly provided by the upscaling of ab initio electronic structure techniques via data-driven and machine-learning avenues [4, 5], molecular simulation can therefore help predict and interpret experimental outcomes, as well as guide the design of materials and technologies toward the betterment of society. Despite impressive progress [6–8], it could be argued that the potential of molecular simulation to deliver on the above prospects remains largely unrealized; this is especially so for molecular processes governed by strong and long-ranged quantummechanical effects [9, 10] and that involve long time-scale, large length-scale, and mechanistically complex conformational and chemical transformations [11, 12]. The state of affairs calls for further development of theoretical models and computational methods, implementable on suitably optimized hardware [13, 14], that advance the accurate molecular simulation of processes exhibiting these traits. The current thesis aims to facilitate the application of path space Markov chain Monte Carlo (McMC) methods toward this overarching goal. Besides its modern applications in Bayesian inference [15–17], mathematical finance [18], and various areas of statistical and machine learning [19], path space McMC has proven to be a versatile framework for molecular simulation as demonstrated by the long-standing success of transition path sampling (TPS) [20–23] and path-integral molecular dynamics (PIMD) [24–28]. Like McMC simulation schemes formulated on the state space (e.g., position space or phase space) of a molec- 2 ular model, such as stochastically thermostatted molecular dynamics (MD) [29], Metropolis–Hastings Monte Carlo (MC) [30], and hybrid MD/MC schemes [31, 32], path space McMC entails the Markovian, step-wise evolution of a path, or an ordered collection of replicas of the molecular model, in accordance with a discretized equation of motion or update rule designed to ensure that admissible path configurations are visited in accordance with prescribed relative probabilities in the long-time limit [21, 33]. The prescribed stationary distribution of path configurations is determined by the features of the underlying molecular model, as illustrated in Ref. [20] in the case of TPS and in Ref. [28] in the case of PIMD. While recasting a state-space McMC sampling problem as a path-space McMC sampling problem might appear computationally disadvantageous due to the intrinsically higher dimensionality of the latter, it can result in enhanced performance as evidenced by the success of (i) non-preconditioned PIMD at estimating the effects of nuclear quantization (as manifest in, e.g., zero-point energy motion and quantum tunneling) on the structural and dynamical properties of condensed-phase systems [3, 10], and of (ii) diffusive TPS at generating mechanistically insightful transition path ensembles between metastable conformational states that would be more computationally demanding to discover via direct diffusive MD simulations [22, 34]. Employing molecular paths as variables to be sampled via McMC nevertheless entails unique methodological challenges, and the current thesis collects some of our efforts to address them. Part I of this thesis concerns non-preconditioned PIMD integration, where pathdiscretization artifacts had previously encumbered the performance and limited the stability of conventional time-discretization schemes [35]. Colleagues recently observed that conventional non-preconditioned PIMD integrators are incompatible with structural features unique to the path space distributions they intend to sample, resulting in a sharp deterioration of their performance and stability under path-discretization refinements [16, 36, 37]. Through the work detailed in Chapters 2 and 3, we introduce a host of dimension-free integration schemes for the non-preconditioned PIMD equations of motion that are free of stability-limiting artifacts in both toymodel and realistic applications. The dimension-free integrators are furthermore designed to yield greater long-time statistical accuracy, higher long-time sampling efficiency as measured by the asymptotic variance of observables of practical interest, and faster convergence to the stationary Boltzmann–Gibbs path distribution than conventional schemes. Having found their way into widely available molecular simulation libraries [38, 39], our dimension-free non-preconditioned PIMD inte- 3 grators are already extending the utility of path space McMC toward elucidating the effects of nuclear quantization in the structure and dynamics of molecular processes modeled at all-atom resolution with ab initio potential energy landscapes [40]. Part II of this thesis explores how path space McMC can enhance long time-scale MD simulation of mesoscale coarse-grained molecular systems. Direct sampling of the trajectory space might appear a convenient avenue to evaluate a long MD trajectory at a lower wall-time cost than with direct integration, but a naïve implementation would be hindered by the slow mixing of McMC schemes when sampling long paths [41, 42]. To effectively apply path space McMC to accelerate the integration of long time-scale stochastic MD, Chapter 4 introduces the path-accelerated molecular dynamics (PAMD) scheme, which enables parallel-in-time integration of a long stochastic MD trajectory through McMC-based evaluation of a relatively short path distribution governing the near-term evolution of a system over the course of several MD time-steps. Compared to competing parallel-in-time integration schemes for stochastic MD [43, 44], which are based on a prediction-correction paradigm [45] and whose performance can be limited by the quality of the prediction employed, the PAMD scheme shows greater potential for speedup due to its McMC-based formulation, which admits implementations through a plethora of path space sampling methods that can efficiently leverage massively parallel computing architectures. The implementation of PAMD outlined in Chapter 4 is numerically shown to substantially outperform a standard serial-in-time integration scheme for overdamped Langevin dynamics in simple systems, and encourages efforts toward refined implementations amenable to real-world applications in molecular simulation and other disciplines. 4 References Cited in Chapter 1 [1] R. Radhakrishnan, “A survey of multiscale modeling: Foundations, historical milestones, current status, and future prospects”, AIChE Journal 67, e17026 (2021). [2] F. Noé and E. Rosta, “Markov models of molecular kinetics”, The Journal of Chemical Physics 151, 190401 (2019). [3] J. S. Kretchmer, N. Boekelheide, J. J. Warren, J. R. Winkler, H. B. Gray, and T. F. Miller III, “Fluctuating hydrogen-bond networks govern anomalous electron transfer kinetics in a blue copper protein”, Proceedings of the National Academy of Sciences 115, 6129 (2018). [4] O. A. von Lilienfeld and K. Burke, “Retrospective on a decade of machine learning for chemical discovery”, Nature Communications 11, 4895 (2020). [5] H. Kulik, T. Hammerschmidt, J. Schmidt, S. Botti, M. A. L. Marques, M. Boley, M. Scheffler, M. Todorović, P. Rinke, C. Oses, A. Smolyanyuk, S. Curtarolo, A. Tkatchenko, A. Bartok, S. Manzhos, M. Ihara, T. Carrington, J. Behler, O. Isayev, M. Veit, A. Grisafi, J. Nigam, M. Ceriotti, K. T. Schütt, J. Westermayr, M. Gastegger, R. Maurer, B. Kalita, K. Burke, R. Nagai, R. Akashi, O. Sugino, J. Hermann, F. Noé, S. Pilati, C. Draxl, M. Kuban, S. Rigamonti, M. Scheidgen, M. Esters, D. Hicks, C. Toher, P. Balachandran, I. Tamblyn, S. Whitelam, C. Bellinger, and L. M. Ghiringhelli, “Roadmap on machine learning in electronic structure”, Electronic Structure, Accepted Manuscript (2022). [6] R. O. Dror, R. M. Dirks, J. P. Grossman, H. Xu, and D. E. Shaw, “Biomolecular simulation: A computational microscope for molecular biology”, Annual Review of Biophysics 41, 429 (2012). [7] N. E. Jackson, M. A. Webb, and J. J. de Pablo, “Recent advances in machine learning towards multiscale soft materials design”, Current Opinion in Chemical Engineering 23, 106 (2019). [8] F. Noé, A. Tkatchenko, K.-R. Müller, and C. Clementi, “Machine learning for molecular simulation”, Annual Review of Physical Chemistry 71, 361 (2020). [9] M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales, and T. E. Markland, “Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges”, Chemical Reviews 116, 7529 (2016). [10] T. E. Markland and M. Ceriotti, “Nuclear quantum effects enter the mainstream”, Nature Reviews Chemistry 2, 0109 (2018). [11] R. Elber, “Perspective: Computer simulations of long time dynamics”, The Journal of Chemical Physics 144, 060901 (2016). 5 [12] D. E. Makarov, “Life in silico: Are we close yet?”, Proceedings of the National Academy of Sciences 118, e2100278118 (2021). [13] D. E. Shaw, R. O. Dror, J. K. Salmon, J. P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. Deneroff, B. Batson, K. J. Bowers, E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson, K. LindorffLarsen, P. Maragakis, M. A. Moraes, S. Piana, Y. Shan, and B. Towles, “Millisecond-scale molecular dynamics simulations on Anton”, in Proceedings of the conference on high performance computing networking, storage and analysis, SC ’09 (2009), pp. 1–11. [14] J. P. Grossman, B. Towles, B. Greskamp, and D. E. Shaw, “Filtering, reductions and synchronization in the Anton 2 network”, Parallel and Distributed Processing Symposium, 860 (2015). [15] A. M. Stuart, “Inverse problems: A Bayesian perspective”, Acta Numerica 19, 451 (2010). [16] S. L. Cotter, G. O. Roberts, A. M. Stuart, and D. White, “MCMC methods for functions: Modifying old algorithms to make them faster”, Statistical Science 28, 424 (2013). [17] M. Dashti and A. M. Stuart, “The Bayesian Approach to Inverse Problems”, in Handbook of uncertainty quantification, edited by R. Ghanem, D. Higdon, and H. Owhadi (Springer International Publishing, 2017), pp. 311–428. [18] V. Linetsky, “The path integral approach to financial modeling and options pricing”, Computational Economics 11, 129 (1997). [19] G. W. Benton, W. J. Maddox, J. P. Salkey, J. Albinati, and A. G. Wilson, “Function-space distributions over kernels”, in Proceedings of the 33rd international conference on neural information processing systems (Curran Associates, Inc., Red Hook, NY, 2019), pp. 1–12. [20] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, “Transition path sampling: Throwing ropes over rough mountain passes, in the dark”, Annual Reviews of Physical Chemistry 53, 291 (2002). [21] M. G. Reznikoff and E. Vanden-Eijnden, “Invariant measures of stochastic partial differential equations and conditioned diffusions”, Comptes Rendus Mathematique 340, 305 (2005). [22] T. F. Miller III and C. Predescu, “Sampling diffusive transition paths”, The Journal of Chemical Physics 126, 144102 (2007). [23] F. J. Pinski and A. M. Stuart, “Transition paths in molecules at finite temperature”, The Journal of Chemical Physics 132, 184104 (2010). [24] D. Chandler and P. G. Wolynes, “Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids”, The Journal of Chemical Physics 74, 4078 (1981). 6 [25] M. Parrinello and A. Rahman, “Study of an F center in molten KCl”, The Journal of Chemical Physics 80, 860 (1984). [26] I. R. Craig and D. E. Manolopoulos, “Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics”, The Journal of Chemical Physics 121, 3368 (2004). [27] T. F. Miller III and D. E. Manolopoulos, “Quantum diffusion in liquid water from ring polymer molecular dynamics”, The Journal of Chemical Physics 123, 154504 (2005). [28] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, “Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space”, Annual Reviews of Physical Chemistry 64, 387 (2013). [29] N. Bou-Rabee, “Time integrators for molecular dynamics”, Entropy 16, 138 (2014). [30] A. Sokal, “Monte Carlo methods in statistical mechanics: Foundations and new algorithms”, in Functional integration: Basics and applications, edited by C. DeWitt-Morette, P. Cartier, and A. Folacci (Springer Science & Business Media, 1997), pp. 131–192. [31] J. P. Nilmeier, G. E. Crooks, D. D. L. Minh, and J. D. Chodera, “Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation”, Proceedings of the National Academy of Sciences 108, E1009 (2011). [32] N. Bou-Rabee and J. M. Sanz-Serna, “Geometric integrators and the Hamiltonian Monte Carlo method”, Acta Numerica 27, 113 (2018). [33] M. Hairer, A. M. Stuart, and J. Voss, “Analysis of SPDEs arising in path sampling part II: The nonlinear case”, The Annals of Applied Probability 17, 1657 (2007). [34] J. Lee, I.-H. Lee, I. Joung, L. Jooyoung, and B. R. Brooks, “Finding multiple reaction pathways via global optimization of action”, Nature Communications 8, 15443 (2017). [35] R. Korol, N. Bou-Rabee, and T. F. Miller III, “Cayley modification for strongly stable path-integral and ring-polymer molecular dynamics”, The Journal of Chemical Physics 151, 124103 (2019). [36] A. Beskos, G. Roberts, A. Stuart, and J. Voss, “MCMC methods for diffusion bridges”, Stochastics and Dynamics 8, 319 (2008). [37] N. Bou-Rabee, “Cayley splitting for second-order Langevin stochastic partial differential equations”, Preprint (2017). 7 [38] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS: a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales”, Computer Physics Communications 271, 108171 (2022). [39] V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, B. A. Helfrecht, P. Juda, S. P. Bienvenue, W. Fang, J. Kessler, I. Poltavsky, S. Vandenbrande, J. Wieme, C. Corminboeuf, T. D. Kühne, D. E. Manolopoulos, T. E. Markland, J. O. Richardson, A. Tkatchenko, G. A. Tribello, V. Van Speybroeck, and M. Ceriotti, “i-PI 2.0: A universal force engine for advanced molecular simulations”, Computer Physics Communications 236, 214 (2019). [40] D. L. Eldridge, R. Korol, M. K. Lloyd, A. C. Turner, M. A. Webb, T. F. Miller III, and D. A. Stolper, “Comparison of experimental vs theoretical abundances of 13 CH3 D and 12 CH2 D2 for isotopically equilibrated systems from 1 to 500 ◦ C”, ACS Earth and Space Chemistry 3, 2747 (2019). [41] M. Hairer, A. M. Stuart, J. Voss, and P. Wiberg, “Analysis of SPDEs arising in path sampling. Part I: The Gaussian case”, Communications in Mathematical Sciences 3, 587 (2005). [42] T. R. Gingrich and P. L. Geissler, “Preserving correlations between trajectories for efficient path sampling”, The Journal of Chemical Physics 142, 234104 (2015). [43] A. L. Blumers, Z. Li, and G. E. Karniadakis, “Supervised parallel-in-time algorithm for long-time Lagrangian simulations of stochastic dynamics: Application to hydrodynamics”, Journal of Computational Physics 393, 214 (2019). [44] F. Legoll, T. Lelièvre, and U. Sharma, “An adaptive Parareal algorithm: Application to the simulation of molecular dynamics trajectories”, SIAM Journal on Scientific Computing 44, B146 (2022). [45] J.-L. Lions, Y. Maday, and G. Turinici, “A parareal discretization in time of PDE’s”, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 332, 661 (2001). Part I DIMENSION-FREE RING-POLYMER MOLECULAR DYNAMICS 9 Chapter 2 Thermostatted ring-polymer molecular dynamics integrators with dimensionality freedom Adapted from R. Korol, J. L. Rosa-Raíces, N. Bou-Rabee, and T. F. Miller III, “Dimension-free path-integral molecular dynamics without preconditioning”, The Journal of Chemical Physics 152, 104102 (2020) DOI: 10.1063/1.5134810. Abstract Convergence with respect to imaginary-time discretization (i.e., the number of ringpolymer beads) is an essential part of any path-integral-based molecular dynamics (MD) calculation. However, an unfortunate property of existing non-preconditioned numerical integration schemes for path-integral molecular dynamics (PIMD), including essentially all existing ring-polymer molecular dynamics (RPMD) and thermostatted RPMD (T-RPMD) methods, is that for a given MD time-step, the overlap between the exact ring-polymer Boltzmann–Gibbs distribution and that sampled using MD becomes zero in the infinite-bead limit. This has clear implications for “Metropolized” MD sampling schemes such as hybrid Monte Carlo, and it also causes the divergence with bead number of the primitive path-integral kinetic energy expectation value when using standard RPMD or T-RPMD. We show that these and other problems can be avoided through the introduction of “dimension-free” numerical integration schemes for which the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite-bead limit for the case of a harmonic potential. Most notably, we introduce the BCOCB integration scheme, which achieves dimension freedom via a particular symmetric splitting of the integration time-step and a novel implementation of the Cayley modification [J. Chem. Phys. 151, 124103 (2019)] for the free ring-polymer half-steps. More generally, we show that dimension freedom can be achieved via mollification of the forces from the external physical potential. The dimension-free path-integral numerical integration schemes introduced here yield finite error bounds for a given MD time-step, 10 even as the number of beads is taken to infinity; these conclusions are proven for the case of a harmonic potential and borne out numerically for anharmonic systems that include liquid water. The numerical results for BCOCB are particularly striking, allowing for nearly three-fold increases in the stable time-step for liquid water with respect to the Bussi-Parrinello (OBABO) and Leimkuhler (BAOAB) integrators while introducing negligible errors in the calculated statistical properties and absorption spectrum. Importantly, the dimension-free, non-preconditioned integration schemes introduced here preserve strong stability, symplecticity, time reversibility, and global second-order accuracy; and they remain simple, black-box methods that avoid additional computational costs, tunable parameters, or system-specific implementations. 2.1 Introduction Considerable effort has been dedicated to the development of numerical integration schemes for imaginary-time path-integral molecular dynamics (PIMD) [1]. In comparison to standard classical molecular dynamics, PIMD numerical integration faces the additional challenge of the highly oscillatory dynamics of the ring-polymer internal modes. Work on numerical integrators for PIMD generally falls into two distinct categories. In the first, the PIMD equations of motion are preconditioned by modifying the ring-polymer mass matrix [2–10]; this approach, which includes the widely used staging algorithms [11], causes the integrated trajectories to differ from those of the ring-polymer molecular dynamics (RPMD) model for real-time dynamics [12, 13], but it can lead to efficient [4–6] sampling of the quantum Boltzmann–Gibbs distribution [14, 15]. In the second category, no modification is made to the ringpolymer mass matrix, i.e., the equations of motion are non-preconditioned [13, 16– 21]. With the aim of providing useful models for real-time quantum dynamics, as well as simple and efficient algorithms for equilibrium thermal sampling, the current chapter focuses on non-preconditioned PIMD numerical integration, notable examples of which include RPMD [12, 13] and its thermostatted variant T-RPMD [20]. Numerical integration schemes for these methods typically employ symmetric factorizations of the time-evolution operator of the form [11, 16–25] ∆t ∆t ∆t ∆t ∆t ∆t e∆tLn ≈ ea 2 On e 2 Bn e 2 An e(1−a)∆tOn e 2 An e 2 Bn ea 2 On with a ∈ {0, 1} , (2.1) 11 where the operator Ln = An + Bn + On includes contributions from the purely harmonic free ring-polymer motion (An ), the external potential gradient (Bn ), and a thermostat (On ), and where n is the (typically large) number of beads or modes comprising the ring polymer. Note that the standard microcanonical RPMD numerical integration scheme is recovered in the limit of zero coupling to the thermostat, and that Eq. (2.1) yields the “OBABO” scheme of Bussi and Parrinello [22] when a = 1 and the “BAOAB” scheme of Leimkuhler [25] when a = 0. In a recent publication [26], it was emphasized that earlier PIMD numerical integration schemes had overlooked a fundamental aspect of the exp((∆t/2)An ) sub-step of the time evolution in Eq. (2.1). Standard practice in these integration schemes has been to exactly evolve the harmonic free ring-polymer dynamics associated with exp((∆t/2)An ) using the uncoupled free ring-polymer normal modes [11, 16–18], which was shown to lack the property of strong stability in the numerical integration, leading to resonance instabilities for microcanonical RPMD and loss of ergodicity for T-RPMD [26]. Use of the Cayley modification to the free ring-polymer motion was shown to impart strong stability to the time-evolution, thereby improving numerical stability for microcanonical RPMD and restoring ergodicity for T-RPMD [26]. In the current study, we focus on the accuracy of both statistical and dynamical properties of the OBABO and BAOAB schemes, as well as the corresponding integrators obtained when the exact free ring-polymer step is replaced by the strongly stable Cayley modification (OBCBO and BCOCB, respectively). Particular attention is paid to the effect of finite-time-step error with these integrators in the limit of large bead numbers. Of these four integrators, it is found that only BCOCB is dimension-free, in the sense that the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite-bead limit for the case of a harmonic potential. Further, it is demonstrated that the OBCBO scheme can be made dimension-free via the technique of force mollification. It is shown that the newly introduced BCOCB integrator yields better accuracy than all other considered non-preconditioned PIMD integrators and allows for substantially larger time-steps in the calculation of both statistical and dynamical properties. Importantly, these gains are made without loss of computational efficiency or algorithmic simplicity. 12 2.2 Non-preconditioned PIMD Consider a one-dimensional quantum particle with the Hamiltonian operator Ĥ = 1 2 p̂ + V (q̂) , 2m (2.2) where m is the particle mass, q̂ and p̂ the position and momentum operators, and V (q̂) a potential energy surface. Ignoring exchange statistics, the properties of this system at thermal equilibrium are encoded in the quantum partition function Q = Tr[e−β Ĥ ] , (2.3) where β = (kB T )−1 , kB is the Boltzmann constant and T the thermodynamic temperature. Using a path-integral discretization (i.e., a Trotter factorization of the Boltzmann–Gibbs operator [27]), Q = limn→∞ Qn can be approximated by the classical partition function Qn of a ring polymer with n beads [14, 15], Qn = mn w n w n −βHn (q,v) d q d ve , (2πℏ)n (2.4) where q = [q0 , . . . , qn−1 ]T is the vector of bead positions and v the corresponding vector of velocities. The ring-polymer Hamiltonian is given by Hn (q, v) = Hn0 (q, v) + Vnext (q) , (2.5) which includes contributions from the physical potential n−1 1X V (qj ) n j=0 (2.6)  mn X  2 vj + ωn2 (qj+1 − qj )2 , 2 j=0 (2.7) Vnext (q) = and the free ring-polymer Hamiltonian n−1 Hn0 (q, v) = where mn = m/n, ωn = n/(ℏβ) and qn = q0 . Non-preconditioned PIMD evolves the phase [q(t)T , v(t)T ]T of the ring polymer as per q̇(t) = v(t) ; v̇(t) = − Ω2 q(t) + m−1 n F (q(t)) (2.8) p 1/2 − Γv(t) + 2β −1 m−1 Γ Ẇ (t) , n 13 which corresponds to the Hamiltonian dynamics of Hn (q, v) with an Ornstein– Uhlenbeck thermostat. In Eq. (3.1) we introduced F (q) = −∇Vnext (q), an ndimensional standard Brownian motion W (t) and the n × n matrices Ω = U diag (0, ω1,n , . . . , ωn−1,n ) U T and Γ = U diag (0, γ1 , . . . , γn−1 ) U T , (2.9) where γj ≥ 0 is the jth friction coefficient, U an n × n real orthogonal discrete Fourier transform matrix such that   −2 1 0 ··· 0 1    1 −2 1  0 · · · 0     . . . . . .   . . . 2 2  , Ω = −ωn  (2.10)  .. .. .. . . .      0 ··· 0  1 −2 1   1 0 ··· 0 1 −2 and the ωj,n , with j = 0, . . . , n−1, are the so-called Matsubara frequencies [28] 
2ωn sin πj if j is even ,  2n (2.11) ωj,n = 2ωn sin π(j+1) else . 2n Observe that the zero-frequency (i.e., centroid) ring-polymer mode is uncoupled from the thermostat, and the coefficients {γj }n−1 j=1 in Eq. (2.9) constitute the friction schedule applied to the non-centroid modes. The standard method for discretizing Eq. (2.8) is to use a symmetric splitting method of the form of Eq. (2.1) that consists of a combination of three types of sub-steps: (i) exact free ring-polymer evolution of time-step τ , " # " # " # q q 0 I ← exp(τA) with A = , (2.12) v v −Ω2 0 where A is the Hamiltonian matrix associated to the free ring polymer, (ii) velocity updates of time-step τ due to forces from the external potential, v ← v + τ m−1 n F (q) , (2.13) and (iii) velocity updates of time-step τ due to the thermostat, v ← exp(−τ Γ) v + p 1/2 β −1 m−1 ξ, n (I − exp(−2τ Γ)) (2.14) 14 where I is the n × n identity matrix and ξ is an n-dimensional vector whose components are independent, standard normal random variables. The acronyms “OBABO” and “BAOAB” indicate the order in which these sub-steps are applied, as specified in Eq. (2.1) with a = 1 or a = 0, respectively. In previous work [26], some of us showed that the matrix exponential for the free ring polymer evolution in Eq. (2.12) is not a strongly stable symplectic matrix, and that as a consequence, the OBABO and BAOAB schemes can display non-ergodicity at time-steps ∆t = kπ/ωj,n for any 1 ≤ j ≤ n and k ≥ 1. We also identified a maximum safe time-step size ∆t⋆ = βℏπ/(2n), below which the matrix exponential is strongly stable. As n → ∞, this maximum safe time-step goes to zero, such that no finite time-step for the scheme in Eq. (2.1) is safe in this limit from non-ergodicity. This non-ergodicity motivates the Cayley modification [26] which consists of approximating the matrix exponential appearing in Eq. (2.12) with the Cayley transform. Specifically, for the Cayley-modified OBABO scheme (called OBCBO), we replace the exact free ring-polymer update of time-step τ = ∆t with cay(∆tA) = (I − (1/2)∆tA)−1 (I + (1/2)∆tA) . (2.15) For the Cayley-modified BAOAB scheme (called BCOCB), we replace the two exact free ring-polymer updates of half-time-step τ = ∆t/2 with cay(∆tA)1/2 . While it might be expected that these half-time-step updates would instead be replaced with cay((∆t/2)A), such a choice leads to a loss of strong stability. Our use of the square root of the Cayley transform preserves strong stability, symplecticity, time reversibility, local third-order accuracy, and guarantees that cay(∆tA)1/2 cay(∆tA)1/2 = cay(∆tA). Furthermore, the square root of the Cayley transform is no more complicated to evaluate than the Cayley transform itself. Both the OBCBO and BCOCB Cayley modifications of Eq. (2.1) are ergodic for a fixed time-step, irrespective of the number of beads; moreover, like Eq. (2.1), the Cayley modified integrators exhibit locally third-order accuracy in the time-step and leave invariant the free ringpolymer Boltzmann–Gibbs distribution in the special case of a constant external potential (V ≡ const.) [26]. 2.3 BCOCB avoids pathologies in the infinite-bead limit In this section, we show that of the OBABO, BAOAB, OBCBO and BCOCB integration schemes, only BCOCB is dimension-free. Although the current section presents 15 analytical results for the specific case of a harmonic external potential, these results are supported by numerical results for anharmonic external potentials in Secs. 2.6 and 2.7. To this end, consider the jth internal ring-polymer mode with frequency ωj,n , in the presence of a harmonic external potential V (q) = (1/2)Λq 2 and a Langevin thermostat with friction γj . Expressed in terms of the normal mode coordinates, obtained from the Cartesian positions and velocities via the orthogonal transformation ϱ = U T q and φ = U T v (2.16) where U is defined in Eq. (3.2), the non-preconditioned PIMD equations of motion for this mode are # " # " # " ϱ̇j (t) ϱj (t) 0 = Kj,n + p 2β −1 m−1 φ̇j (t) φj (t) (2.17) n γj Ẇj (t) Kj,n = Aj,n + B + Oj , where Ẇj is a scalar white-noise and we have introduced the 2 × 2 matrices " # " # " # 0 1 0 0 0 0 Aj,n = , B= and Oj = . 2 −ωj,n 0 −Λ/m 0 0 −γj The solution (ϱj (t), φj (t)) of Eq. (2.17) is a bivariate Gaussian, and in the limit as t → ∞, the probability distribution of (ϱj (t), φj (t)) converges to a centered bivariate normal distribution with covariance matrix " # 2 sj,n 0 1 1 where s2j,n = Σj,n = . (2.18) 2 βmn 0 1 Λ/m + ωj,n For this system, a single time-step of Eq. (2.1) can be compactly written as " # " # " # ϱj (t + ∆t) ϱj (t) ξ0 1/2 = Mj,n,∆t + Rj,n,∆t , (2.19) φj (t + ∆t) φj (t) η0 where ξ0 and η0 are independent standard normal random variables and we have introduced the 2 × 2 matrices ∆t ∆t ∆t ∆t ∆t ∆t Mj,n,∆t = ea 2 Oj e 2 B e 2 Aj,n e(1−a)∆tOj e 2 Aj,n e 2 B ea 2 Oj 1 − e−2(1−a)γj ∆t T Nj,n,∆t P Nj,n,∆t βmn
∆t ∆t 1 − e−aγj ∆t + (Mj,n,∆t e−a 2 Oj )P (Mj,n,∆t e−a 2 Oj )T + P βmn Rj,n,∆t = 16 ∆t ∆t ∆t where P = [ 00 01 ] and Nj,n,∆t = ea 2 Oj e 2 B e 2 Aj,n,∆t . The corresponding step for the Cayley modification is obtained by replacing exp((∆t/2)Aj,n ) in Eq. (2.19) with cay(∆tAj,n )1/2 , which is given by " # s 2 ∆t 1 . (2.20) cay(∆tAj,n )1/2 = 2 2 4 + ωj,n ∆t2 −ωj,n ∆t 2 A sufficient condition1 for ergodicity of Eq. (2.19) is 1 > A2j,n,∆t cosh2 ((∆t/2)γj ) , where Aj,n,∆t = cos(∆tωj,n ) − (2.21) (Λ/m)∆t sin(∆tωj,n ) . 2ωj,n For the Cayley modification of Eq. (2.19), Eq. (2.21) still provides a sufficient condition for ergodicity, except with 2 4 − (ωj,n + 2(Λ/m))∆t2 Aj,n,∆t = . 2 4 + ωj,n ∆t2 Due to the lack of strong stability in the exact free ring-polymer evolution, Eq. (2.19) fails to meet the condition in Eq. (2.21) and becomes non-ergodic whenever ∆t = kπ/ωj,n where k ≥ 1 [26]; no such problem exists for the Cayley modification. Regardless, assuming that the condition in Eq. (2.21) holds, the numerical stationary distribution is a centered Gaussian with 2×2 covariance matrix Σj,n,∆t that satisfies the linear equation T Σj,n,∆t = Mj,n,∆t Σj,n,∆t Mj,n,∆t + Rj,n,∆t , for which the solution is " # s2j,n,∆t 0 1 Σj,n,∆t = . 2 βmn 0 rj,n,∆t (2.22) The variances in the position and velocity marginals of the numerical stationary distribution with the covariance matrix in Eq. (2.22) are, respectively, (βmn )−1 s2j,n,∆t 1 In the special case when Λ = 0, the given condition for OBCBO corrects a sign error in Eq. (37) of Ref. [26]. 17 2 and (βmn )−1 rj,n,∆t with  1    2 Λ∆tωj,n ωj,n + m cot(∆tωj,n ) − ( Λ∆t )2 2m s2j,n,∆t = 1    2 Λ∆tωj,n ω ) ωj,n + 2m cot( ∆t 2 j,n   a=1 1 2 ∆t rj,n,∆t = 2mωj,n − Λ∆t tan( ωj,n ) 2   a=0 2mωj,n a=1 (2.23) a=0 (2.24) For the Cayley modification of Eq. (2.19), s2j,n,∆t = 2 rj,n,∆t = 4m s2 , 4m − a∆t2 Λ j,n 4m − (1 − a)∆t2 Λ . 4m (2.25) (2.26) Note that these numerical stationary distributions are independent of the friction parameter γj , which is a benefit of schemes based on splitting the non-preconditioned PIMD dynamics into Hamiltonian and thermostat parts, and using the exact Ornstein– Uhlenbeck flow in Eq. (2.14) to evolve the thermostat part. Moreover, comparing the exact covariance matrix in Eq. (2.18) with the finite-time-step approximations in Eqs. (2.22)–(2.26), note that in all cases Σj,n = lim∆t→0 Σj,n,∆t . These results have previously been reported for the OBABO (Eqs. (2.23) and (2.24), a = 1) and BAOAB (Eqs. (2.23) and (2.24), a = 0) schemes [8, 29] but not for the OBCBO (Eqs. (2.25) and (2.26), a = 1) or BCOCB (Eqs. (2.25) and (2.26), a = 0) schemes. In normal mode coordinates, the exact and numerical position-marginals can be written as an infinite product of one-dimensional centered normal distributions with variances given by (βmn )−1 s2j,n and (βmn )−1 s2j,n,∆t , respectively. By Kakutani’s theorem [30, 31], these two distributions have a non-zero overlap in the infinitebead limit if and only if the limit series 2 n−1  X sj,n (2.27) lim 1− n→∞ sj,n,∆t j=1 converges. For OBABO and BAOAB, due to the oscillatory cotangent term appearing in sj,n,∆t , the limit limn→∞ (1 − sn−1,n /sn−1,∆t )2 does not exist, and therefore, the series does not converge. For OBCBO, the jth summand of this series is !−2 r ∆t4 Λ2 4m − ∆t2 Λ 1+ , 16m2 4m 18 which more obviously diverges. Therefore, for OBABO, OBCBO, and BAOAB, the numerical stationary distribution has no overlap with the exact stationary distribution in the infinite-bead limit; it is in this sense that these schemes fail to exhibit the property of dimensionality freedom. Remarkably, BCOCB is exact in the positionmarginal and thus exhibits dimensionality freedom. See Sec. 2.9 for a brief summary of the properties of other symmetric splittings that were considered. 2.4 Consequences for the primitive kinetic energy expectation value In the current section, we show that the non-overlap pathology of the OBABO, BAOAB, and OBCBO schemes causes a divergence with increasing bead number of the primitive path-integral kinetic-energy expectation value, an issue that is numerically well known for OBABO and BAOAB [8, 29, 32, 33]. We further show that this divergence is fully eliminated via the BCOCB scheme, as expected. The primitive kinetic energy expectation value is given by [34, 35] n X mn ωn2 n − ⟨(qj − qj−1 )2 ⟩ 2β j=1 2  n−1  2 X mn ωj,n 1 1 2 + − ⟨ϱj ⟩ = 2β j=1 2β 2 ⟨KEpri n ⟩ = (2.28) (2.29) where the first equality involves a sum over the ring-polymer beads in Cartesian coordinates (with qn = q0 ) and the second equality performs the summation in terms of the ring-polymer normal modes. The divergence of this expectation value is numerically illustrated for the simple case of a harmonic oscillator (Figs. 2.1a-d); note that for larger MD time-steps, the OBABO, BAOAB, and OBCBO schemes fail to reach a plateau with increasing bead number and dramatically deviate from the exact result (dashed line). The same divergence for OBABO and BAOAB has been numerically observed in many systems [8, 29, 32, 33], including liquid water which we discuss later. A striking observation from Figs. 2.1a-d is that the BCOCB exhibits no such divergence or error in the primitive kinetic energy expectation value at high bead number, regardless of the employed time-step. Using Eq. (2.18), note that the contribution to the primitive kinetic energy expectation value from the jth ring-polymer mode is ⟨KEpri j,n ⟩ =
1 2 1 − ωj,n s2j,n , 2β 4 3 (a) 2 1 0 t = 0.5 fs OBABO OMCMO BAOAB 4 8 16 OBCBO OmCmO BCOCB 32 Primitive KE (kBT) Primitive KE (kBT) 19 4 3 4 8 3 2 1 0 (c) t = 1.0 fs 4 8 16 32 64 128 256 # of beads 16 32 64 128 256 # of beads Primitive KE (kBT) Primitive KE (kBT) # of beads 4 t = 0.5 fs 1 0 64 128 256 (b) 2 4 3 2 1 0 4 8 16 (d) t = 2.0 fs 32 64 128 256 # of beads Figure 2.1. Primitive kinetic energy expectation values for a harmonic potential V (q) = 1 2 2 Λq with force constant Λ = 256, ℏ = 1, mass m = 1 and reciprocal temperature β = 1; choosing energies to be in units of kB T at room temperature (300 K), βℏ ≈ 25.5 fs and Λ = mω 2 where ω = 3315 cm−1 . (Panels a–d) For various MD time-steps ∆t, the primitive kinetic energy expectation value as a function of the number of ring-polymer beads, with the exact kinetic energy indicated as a dashed gray line. The standard error of all visible data points in each plot is smaller than the symbol size. (Panel e) Per-mode error in the variance of position coordinate of the normal modes for simulations run with 128 ringpolymer beads and a time-step of 1 fs; solid lines are analytic predictions from Eq. (2.33) with (2.23) and (2.25) defining s2j,n,∆t for the different schemes; points indicate the results of numerical PIMD simulations using the various integration schemes. The BCOCB scheme is not shown since it has zero error for all internal modes. The black vertical line indicates the crossover frequency (ωx = 2/∆t) for the error of OBCBO and OMCMO based on the bounds in Eqs. (2.43) and (2.44). 20 such that in the infinite-bead limit, lim n→∞ n−1 X j=1 ⟨KEpri j,n ⟩ = ℏ 4 r Λ coth m βℏ 2 r Λ m ! . (2.30) Similarly using Eq. (2.22), the jth-mode contribution to the kinetic energy from the finite-time-step numerical expectation value is ⟨KEpri j,n,∆t ⟩ =
1 2 s2j,n,∆t . 1 − ωj,n 2β (2.31) Thus, the per-mode error in kinetic energy is pri | ⟨KEpri j,n ⟩ − ⟨KEj,n,∆t ⟩ | = 2 mn ωj,n ρj,n,∆t , 2 (2.32) where the per-mode error in the position-marginal for internal mode j is ρj,n,∆t = 1 s2 − s2j,n,∆t , βmn j,n (2.33) where sj,n,∆t is given by Eq. (2.23) for the cases of OBABO (a = 1) and BAOAB (a = 0) and by Eq. (2.25) for the cases of OBCBO (a = 1) and BCOCB (a = 0). Note that this error vanishes only for the BCOCB scheme, which satisfies ρj,n,∆t = 0 for each mode j, irrespective of the time-step ∆t. Equations (2.32) and (2.33) indicate that the primitive kinetic energy estimator is a sensitive measure of the finite-time-step error in the sampled ring-polymer position distribution associated with the high-frequency modes. Figure 2.1e resolves this per-mode error, ρj,n,∆t , for each internal mode in simulations that employ a total of 128 beads, including results from OBABO (red), BAOAB (magenta), and OBCBO (blue), using a time-step of 1 fs, with the solid lines indicating the analytical predictions in Eq. (2.33) and with the dots indicating the result of numerical simulations. The analytical results are fully reproduced by the simulations. Note that the OBABO per-mode error exhibits dramatic spikes for ωj,n ∆t = kπ where 1 ≤ j ≤ n and for some k ≥ 1, which coincide with the loss of ergodicity of that integration scheme. The BAOAB scheme exhibits these resonance instabilities at even values of k. However, it is the failure of this per-mode error to sufficiently decay as a function of the mode number for all three of OBABO, BAOAB, and OBCBO that gives rise upon summation to the divergence of the primitive kinetic energy expectation 2 value, as seen for this particular time-step value in Fig. 2.1d. Since ωj,n s2j,n → 1 Pn−1 pri as n → ∞, the convergence of limn→∞ j=1 | ⟨KEpri j,n ⟩ − ⟨KEj,n,∆t ⟩ | reduces to P 2 2 the convergence of the series limn→∞ n−1 j=1 sj,n − sj,n,∆t , which diverges for both OBABO and OBCBO due to the same reasons as discussed in Sec. 2.3. 21 2.5 Dimensionality freedom for OBCBO via force mollification Sections 2.3 and 2.4 have demonstrated that whereas the BCOCB integrator exhibits dimensionality freedom, the OBCBO integrator does not. In the current section, we show that this shortcoming of OBCBO can be addressed by the use of force mollification, in which the external potential energy in Eq. (2.6) is replaced by Ṽnext (q) = Vnext (sinc(Ω̃∆t/2)q) , (2.34) where Ω̃ is any positive semi-definite n × n matrix that has the same eigenvectors as Ω (Eq. (3.2)) while possibly having different eigenvalues. Force mollification has not previously been employed for PIMD, although the strategy originates from a variation-of-constants formulation of the solution to Eq. (2.8) [36–39]; specifically, the protocol in Eq. (2.34) is a generalization of the mollified impulse method [36]. The use of force mollification for enhancing non-preconditioned PIMD integration can be motivated on physical grounds: In the absence of a physical potential, four of the considered integration schemes (OBABO, BAOAB, OBCBO, and BCOCB) leave invariant the exact free ring-polymer Boltzmann–Gibbs distribution [26]. Therefore, the loss of any overlap between the exact stationary distribution of the position marginals in the infinite-bead limit for OBABO, BAOAB, and OBCBO must be attributed to the influence of the time-evolution from the external potential in the schemes (i.e., the “B” sub-step) as implemented in Eq. (2.13); the BCOCB scheme does not suffer from this problem. To remove this pathology in the OBCBO scheme, we thus use mollification to “taper down” the external forces on the high-frequency modes, such that the resulting integration correctly reverts to free-particle motion for those modes, which should become decoupled from the external potential as the frequency increases. The specific appearance of the 1/2 factor in the sinc function argument ensures that the sinc function switches from its high-frequency effect to its low-frequency effect when the period of the Matsubara frequency is commensurate with ∆t; the zero-frequency ring-polymer centroid mode is untouched by mollification. Force mollification requires only a small algorithmic modification of the OBCBO integrator. Specifically, the “B” sub-step in Eq. (2.13) is replaced with ∆t −1 (2.35) v←v+ m F̃ (q) , 2 n where the mollified forces are F̃ (q) = sinc(Ω̃∆t/2)F (q̃) = U D∆t U T F (q̃) (2.36) 22 where q̃ = U D∆t U T q are the mollified bead positions and where D∆t is the diagonal matrix of eigenvalues associated with sinc(Ω̃∆t/2), i.e., D∆t = diag (sinc(ω̃0,n ∆t/2), . . . , sinc(ω̃n−1,n ∆t/2)) (2.37) where ω̃j,n is the jth eigenvalue of Ω̃. In practice, the mollified forces are computed in normal mode coordinates as follows: (a) Starting with the ring-polymer bead position in normal mode coordinates, obtain a copy of the mollified bead positions via q̃ = U D∆t ϱ . (2.38) (b) Evaluate the external forces at the mollified ring-polymer bead positions, F (q̃). (c) Apply the remaining mollification to the forces in Eq. (2.36) via U T F̃ (q) = D∆t U T F (q̃) . (2.39) We emphasize that in comparison to the standard force update (Eq. (2.13)) the use of the mollified force update (Eq. (2.35)) introduces neither additional evaluations of the external forces nor n × n matrix multiplies associated with the discrete Fourier transform; it therefore avoids any significant additional computational cost. This mollification scheme preserves reversibility and symplecticity as well as local third-order accuracy of the OBCBO scheme with time-step. We emphasize that the sinc-function-based mollification scheme in Eq. (2.35) is not unique and alternatives can certainly be devised. Even within the functional form of the mollification in Eq. (2.35), flexibility remains with regard to the choice of the matrix Ω̃, which allows for mode-specificity in the way the mollification is applied. A simple choice for this matrix is Ω̃ = Ω, such that mollification is applied to all of the non-zero ringpolymer internal modes. With this choice, we arrive at a fully-specified integration scheme that replaces the original “B” sub-step in Eq. (2.13) with the mollified-force sub-step in Eq. (2.35); we shall refer to this force-mollified version of OBCBO integration scheme as “OMCMO.” In Sec. 2.5.1, we propose a partially mollified choice for Ω̃ that further improves accuracy. For the harmonic external potential, the previously derived relations for OBCBO (most notably Eqs. (2.21), (2.25)–(2.26) and (2.32)–(2.33)) also hold for OMCMO with Λ suitably replaced by Λ̃j = sinc2 (ωj,n ∆t/2)Λ. Note that Λ̃j ≤ Λ, since 23 sinc2 (x) ≤ 1 for all x ≥ 0, making clear that the mollification reduces the effect of the external potential on the higher-frequency internal ring-polymer modes. We now show that mollifying the forces in the B sub-step fixes the pathologies of OBCBO in the infinite-bead limit, by restoring overlap between the sampled and exact stationary distributions. To see this, note that the jth summand in Eq. (2.27) for OMCMO satisfies !2  2 s2j,n sj,n ∆t4 Λ2 1− ≤ 1− 2 ≤ f (ωj ∆t/2) , sj,n,∆t sj,n,∆t 16m2
2 where f (x) = (1 − sinc2 (x))/ x2 + sinc2 (x) and we have used the infinite-bead limit for the ring-polymer internal-mode frequencies  πj  if j is even ,  ℏβ (2.40) ωj = lim ωj,n = π(j + 1) n→∞   else . ℏβ since2 ∞ X j=1 we obtain f (ωj ∆t/2) ≤ 6 n−1  X lim 1− n→∞ j=1 sj sj,n,∆t 2 ℏβ +4, π∆t   4 2 ℏβ ∆t Λ ≤ 6 +4 . π∆t 16m2 (2.41) Again invoking Kakutani’s theorem (Eq. (2.27)), it follows that the numerical stationary distribution has an overlap with the exact stationary distribution. As a byproduct of this analysis, we can also quantify the amount of overlap between the exact 2 This inequality comes from using Eq. (2.40) to write P∞ ⌊ℏβ/(π∆t)⌋ I=2 X ∞ X f (jπ∆t/(ℏβ)) and II = 2 j=1 j=1 f (ωj ∆t/2) = I + II where f (jπ∆t/(ℏβ)) . j=⌈ℏβ/(π∆t)⌉ Then the first term admits the bound I ≤ 2f (1)ℏβ/(π∆t) < 4ℏβ/(π∆t) ; for the second term, we use II ≤ F (1) + ℏβ/(π∆t) w∞ 1 F (x)dx F (x) = 2((1 − sinc2 (x))/x2 + 1/x2 )2 is monotone decreasing on [1, ∞) with F (1) ≤ 4 and rwhere ∞ F (x)dx ≤ 2. 1 24 and numerically sampled stationary distributions, revealing that the total variation distance [40] between these distributions is given by3 s  2 ℏβ ∆t Λ dTV (µ, µ∆t ) ≤ 6 +4 . (2.42) π∆t 2m In summary, the force mollification strategy introduced here provably removes the pathologies due to the “B” sub-step in the case of a harmonic oscillator potential. Moreover, for any finite number of beads, the total variation distance between the exact and numerically sampled stationary distribution can be bounded by Eq. (2.42), and thus OMCMO admits dimension-free error bounds. Before proceeding, we first return to Fig. 2.1 to compare the accuracy of OMCMO with the un-mollified OBCBO scheme for the internal mode position marginal of the harmonic oscillator. As seen in Fig. 2.1e for the results with a time-step of 1 fs, the per-mode error obtained by the mollified scheme (OMCMO, green) decays more rapidly with mode number than does OBCBO. Figure 2.1d further illustrates that upon summation of the per-mode contributions, the OMCMO prediction for the primitive kinetic energy converges to a well-defined asymptote with respect to the number of ring-polymer beads, whereas OBCBO diverges as discussed earlier. Similar behavior is seen for shorter MD time-steps (panels a-c), although the failure of OBCBO becomes less severe with this range of bead numbers as the time-step is reduced. Although it is satisfying that mollification via OMCMO both formally and numerically ameliorates the problems of the OBCBO scheme in the high-bead-number limit, the OMCMO results in Fig. 2.1 are not ideal, since in some cases the OMCMO error is substantially larger than that of OBCBO when a modest number of beads is used (e.g., for 16 beads in panel d). This observation points to a simple and general refinement of the OMCMO scheme, which we discuss in the next subsection. 3 This quantification uses: (i) dTV ≤ 21/2 dH where dTV is the total variation distance and dH is the Hellinger distance; and (ii) subadditivity of the squared Hellinger distance, which implies that d2H (µ, µ∆t ) ≤ lim n→∞ ≤ lim n→∞ n−1 X j=1 n−1 X j=1 d2H (N (0, s2j,n ), N (0, s2j,n,∆t )) (1 − s2j,n /s2j,n,∆t )2 ≤ (3ℏβ/(π∆t) + 2) ∆t4 Λ2 . 8m2 25 2.5.1 Partial mollification Comparison of the per-mode errors from OBCBO and OMCMO in Fig. 2.1e reveals that lower errors for OMCMO are only enjoyed for internal modes that exceed a particular frequency (indicated by the vertical black line). This observation suggests that if a “crossover frequency” could be appropriately defined, then a refinement to OMCMO could be introduced for which mollification is applied only to the ringpolymer internal modes with frequencies that exceed this crossover value. For the case of a harmonic external potential, this crossover frequency ωx can be found by comparing a bound for the per-mode error (Eq. (2.33)) for OBCBO   ∆t2 Λ 1 (2.43) ρj,n,∆t ≤ 2 mn ωj,n β 4m − ∆t2 Λ to that for OMCMO  ρj,n,∆t ≤ g(ωj,n ∆t/2) 1 ∆t2 Λ 2 mn ωj,n β 4m − ∆t2 Λ  , (2.44) where g(x) = (1 − sinc2 (x))/x2 + sinc2 (x). Since g(x) ≥ 1 only when x ≤ 1, we expect better accuracy if mollification is only applied to those ring-polymer internal modes with frequencies ωj,n ≥ ωx , where ωx = 2/∆t. Although this result was derived for the case of a harmonic potential, it does not depend on Λ. We call this resulting partly mollified integration scheme “OmCmO.” This scheme has the nice properties of OMCMO, including strong stability and dimensionality freedom. Implementation of OmCmO is a trivial modification of OMCMO, requiring only that the diagonal elements of D∆t in Eq. (2.37) are evaluated using  1 for ωj,n < ωx sinc(ω̃j,n ∆t/2) = (2.45) sinc(ω ∆t/2) otherwise , j,n where j = 0, . . . , n − 1. In physical terms, the emergence of 2/∆t in the crossover frequency is intuitive, since as was previously mentioned, it corresponds to having the ring-polymer mode undergo a full period per time-step ∆t. Finally, numerical results for the case of a harmonic potential (Figs. 2.1a-d) reveal that the partially modified OmCmO scheme (cyan) achieves both robust convergence of the primitive kinetic energy with increasing bead number, as well as consistently better accuracy than the OBCBO and OMCMO integration schemes, as expected. However, it must be emphasized that for all panels of Fig. 2.1, the BCOCB scheme (which requires no force mollification) is by far the most accurate and stable. 26 2.6 Results for anharmonic oscillator Having numerically characterized the performance of the various non-preconditioned PIMD integrators for the case of the harmonic oscillator external potential in Fig. 2.1, we now turn our attention to anharmonic external potentials. In this section, we consider both a weakly anharmonic (aHO) potential   1 3 1 4 1 2 (2.46) q + q + q V (q) = Λ 2 10 100 and the more strongly anharmonic quartic potential 1 V (q) = q 4 . 4 (2.47) All calculations are performed using ℏ = 1, m = 1 and β = 1. Assuming the system to be at room temperature (300 K), the thermal timescale corresponds to βℏ ≈ 25.5 fs and Λ = mω 2 , where ω = 3315 cm−1 for Λ = 256. The trajectories are performed with the centroid mode uncoupled from the thermostat (i.e., in the manner of T-RPMD); for the remaining n − 1 internal modes, simulations performed with the OBABO and BAOAB schemes use the standard [18, 20] damping schedule of Γ = Ω, and simulations performed using the Cayley modification (i.e., BCOCB, OBCBO, OMCMO, and OmCmO) use friction γj = min{ωj,n , 0.9γjmax (Λ), 0.9γjmax (0)} for the jth mode, where γjmax (Λ) is the friction that saturates the inequality in Eq. (2.21); for the quartic potential, we set Λ = 1 in this calculation of γjmax . Figures 2.2a and 2.2b presents kinetic energy expectation values for the aHO potential corresponding to 3315 cm−1 at room temperature. For the primitive kinetic energy expectation value, the results obtained using the various integration schemes with time-steps of both 0.5 fs (panel a) and 1.0 fs (panel b) are consistent with the observations for the harmonic potential in Fig. 2.1; specifically, the integrators without dimensionality freedom (OBABO, BAOAB, and OBCBO) fail to converge with increasing bead number, while the mollified integrators (OMCMO and OmCmO) smoothly converge with increasing bead number and the partially mollified scheme (OmCmO) is consistently more accurate than OBCBO and OMCMO. However, it is also clear that BCOCB exhibits the best accuracy with increasing bead number, converging to the exact result without perceivable time-step error. Figures 2.2c and 2.2d present the corresponding results for the virial kinetic energy 3 2 1 0 (a) 0.5 fs 4 5.0 Virial KE (kBT) Primitive KE (kBT) 4 8 16 32 64 128 256 OBABO OMCMO BAOAB 4.5 OBCBO OmCmO BCOCB 4.0 3.5 (c) 0.5 fs 16 4 3 2 1 0 # of beads 32 64 128 256 # of beads (b) 1.0 fs 4 8 16 32 64 128 256 # of beads 5.0 Virial KE (kBT) Primitive KE (kBT) 27 4.5 4.0 3.5 (d) 1.0 fs 16 32 64 128 256 # of beads Figure 2.2. Primitive and virial kinetic energy expectation values as a function of bead number for the weakly anharmonic potential corresponding to 3315 cm−1 at room temperature, with results obtained using a time-step of 0.5 fs (panels a and c) and 1.0 fs (panels b and d). The standard error of all visible data points in each plot is smaller than the symbol size. The exact kinetic energy is indicated with a dashed line. expectation value, ⟨KEvir n ⟩ = 1 1 − (q − q̄1n )T F (q) 2β 2 (2.48) where q̄ is the centroid (bead-averaged) position and 1n is the n-dimensional vector with all entries equal to 1. Whereas the virial kinetic energy for all of the strongly stable integration schemes is well behaved, the OBABO and BAOAB schemes perform erratically at large time-steps due to their provable non-ergodicities [26]. Appealingly, the BCOCB scheme is consistently the most accurate for the virial kinetic energy expectation value, as it was for the primitive kinetic energy expectation value. Figure 2.3a-d shows the results of the various numerical integration schemes for the primitive and virial kinetic energy expectation values, as a function of the MD time-step using 64 ring-polymer beads. Results are shown for both the aHO and the strongly anharmonic quartic oscillator. In all cases, the BCOCB scheme is consistently the most accurate across this array of model systems. 28 Finally, Fig. 2.3e illustrates the use of the BCOCB integrator for the calculation of real-time quantum dynamics within the T-RPMD approximation, where it replaces the often-employed OBABO integration scheme. Using 64 beads, the T-RPMD results are plotted for a range of integration time-steps. Strikingly, over the entire range of considered time-steps, BCOCB introduces negligible error in the calculated position time autocorrelation function; it is confirmed that these results are graphically indistinguishable from those obtained using the OBABO integrator in the small-time-step limit. 2.7 Results for liquid water In model systems, Sec. 2.6 demonstrated the strong performance of the BCOCB integrator at evaluating both PIMD statistics and real-time dynamics through the TRPMD model. In the current section, we test the accuracy and stability of the various un-mollified integration schemes (i.e., OBABO, OBCBO, BAOAB, and BCOCB) in liquid water, a high-dimensional and relatively complex system. Specifically, we consider a periodic 32-molecule water box at a temperature of 298 K and a density of 0.998 g/cm3 , as described by the q-TIP4P/F force field [41]. In Fig. 2.4, we compare the accuracy achieved by the different integrators for the average kinetic energy per hydrogen atom as a function of the number of ring-polymer beads. As in previous sections, we consider both the primitive (Eq. (2.29)) and virial (Eq. (2.48)) estimators for the kinetic energy. For each choice of integrator, time-step and bead number, the primitive and virial estimators for the kinetic energy of a single hydrogen atom were averaged over a 1-nanosecond trajectory integrated in the manner of T-RPMD (i.e., with the centroid mode uncoupled from the thermostat); for the remaining n − 1 internal modes, simulations performed with the OBABO and BAOAB schemes use the standard [18, 20] damping schedule of Γ = Ω, and simulations performed using the Cayley modification use friction γj = min{ωj , 0.9γjmax (0)}, where γjmax (Λ/m) saturates the inequality in Eq. (2.21) for the given values of j and Λ/m at the given time-step. Multi-nanosecond staging PIMD [8, 11] simulations at the small time-step of 0.1 fs were performed to obtain a bead-converged reference value for the H-atom kinetic energy, plotted as a dashed line in Figs. 2.4 and 2.5. The primitive kinetic energy expectation values in Figs. 2.4a and 2.4b show similar trends to those seen in Figs. 2.1 and 2.2 for the harmonic and weakly anharmonic Primitive KE (kBT) Virial KE (kBT) 4 2 0 0.125 0.25 0.5 timestep (fs) 0.8 1 0.4 0.2 1 (b) aHO 4 OBABO OBCBO BAOAB 3 0.125 0.25 0.5 timestep (fs) 0.8 (c) Quartic 0.6 0.0 Cqq (t) 5 (a) aHO Virial KE (kBT) Primitive KE (kBT) 29 2 timestep (fs) (d) Quartic 0.6 0.4 OMCMO OmCmO BCOCB 0.2 0.0 4 1 t=0.125 fs t=2.0 fs 0.5 1 2 timestep (fs) 4 t=4.0 fs t=8.0 fs 0.0 0.5 (e) Quartic 0 50 100 150 200 250 time (fs) Figure 2.3. Primitive and virial kinetic energy expectation values as a function of the timestep for the weakly anharmonic potential corresponding to 3315 cm−1 at room temperature (panels a and b), and for the quartic potential (panels c and d). The exact kinetic energy is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size. Also, the position autocorrelation function (panel e) for the quartic oscillator at room temperature computed using T-RPMD with the BCOCB integrator. Results are obtained using 64 ring-polymer beads using time-steps of ∆t = 0.125, 2, 4, and 8 fs. 30 oscillators. For a 0.5-fs time-step (Fig. 2.4a), at which all integrators exhibit strong stability for ring polymers with up to 64 beads at the system temperature [26], the OBABO, BAOAB, and OBCBO primitive kinetic energy estimates diverge from the converged result as the number of beads increases, in agreement with the proven result that the error in the ring-polymer configurational distribution generated with these schemes grows unboundedly with increasing bead number. At the larger, 0.8-fs time-step, (Fig. 2.4b), OBABO and BAOAB formally lose strong stability and their respective primitive kinetic energy estimates dramatically diverge for bead numbers greater than 32; the strongly stable OBCBO scheme also yields a divergent result for the same reason as in Fig. 2.4a. As seen on the HO and aHO model systems, the primitive kinetic energy expectation value from the BCOCB integrator monotonically converges to the reference value with increasing bead number, avoiding any perceptible time-step error. Figures 2.4c and 2.4d show the corresponding virial kinetic energy expectation values. For the smaller time-step of 0.5 fs, which is a common choice for path-integral simulations of water, all of the integrators perform similarly. However, upon increasing the time-step to 0.8 fs, significant differences in the performance of the integrators emerges, with only BCOCB avoiding perceptible time-step error. To further compare the accuracy and stability of the OBABO, BAOAB, OBCBO, and BCOCB integrators, Fig. 2.5 considers the average kinetic energy per hydrogen atom obtained using 64 beads over a wide range of time-steps. These results show that BCOCB remains remarkably accurate for time-steps as large as 1.4 fs for liquid water, which corresponds to the limit of stability for Verlet integration of the centroid mode. In comparison, OBCBO diverges monotonically as the time-step increases, reaching unphysical values for the primitive expectation value and yielding sizable error (20%) for the virial expectation value. The erratic performance of both OBABO and BAOAB is due to the emergence of numerical resonance instabilities at time-steps greater than 0.6 fs at 64 beads; indeed, the largest safe time-step at which OBABO and BAOAB remain strongly stable at this bead number, ∆t⋆ ≈ 0.63 fs, precedes the range of time-steps in Fig. 2.5 where these two integrators vary erratically. Extending beyond statistics, we now consider quantum dynamical properties of liquid water within the T-RPMD approximation. Given the outstanding performance of the BCOCB scheme for the calculated statistical properties in Figs. 2.4 and 2.5, we present results that focus on this scheme in comparison to the most widely used 31 Figure 2.4. Primitive and virial kinetic energy expectation values as a function of the bead number for the hydrogen atom in liquid water at 298 K and 0.998 g/cm3 at time-step ∆t = 0.5 fs (panels a and c) and ∆t = 0.8 fs (panels b and d). The reference kinetic energy, converged along a staging PIMD simulation at time-step ∆t = 0.1 fs and 256 beads, is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size. OBABO scheme. In particular, we consider the liquid water infrared absorption ˜ ˜ where the dipole spectrum I(ω) = spectrum [42], which is proportional to ω 2 I(ω) r −iωt C̃µ·µ (t) is the Fourier transform of the Kubo-transformed dipole autocordt e R relation function C̃µ·µ (t). The latter is approximated in the RPMD model by [16] P C̃µ·µ (t) = N1 N i=1 ⟨µ̄i (t) · µ̄i (0)⟩, where N is the number of molecules in the liquid, µ̄i (t) is the bead-averaged dipole moment of molecule i at time t and the angle brackets denote an average against the ring-polymer thermal distribution. To obtain the time-correlation functions and spectra shown in Fig. 2.6 for the OBABO and BCOCB integration schemes, 12-nanosecond T-RPMD trajectories were simulated for a ring-polymer with 64 beads and time-steps ranging from 0.2 to 1.4 fs, using the same friction schedule as described for Figs. 2.4 and 2.5. Along each trajectory, the velocities of all degrees of freedom in the system were drawn anew from the Maxwell–Boltzmann distribution every 20 ps; the autocorrelation function was evaluated out to 2 ps by averaging over staggered windows of that 32 Figure 2.5. Primitive and virial kinetic energy expectation values as a function of the timestep for the hydrogen atom in liquid water at 298-K and 0.998 g/cm3 as described by a 64-bead ring polymer. The reference kinetic energy, converged along a staging PIMD simulation at time-step ∆t = 0.1 fs and 256 beads, is indicated with a dashed line. The standard error of all visible data points in each plot is smaller than the symbol size. time-length within every 20-picosecond trajectory segment; and exponential-decay extrapolation was used to extend the autocorrelation function before evaluating its numerical Fourier transform to obtain the infrared absorption spectrum. Figures 2.6a and 2.6b present the dipole autocorrelation functions obtained using the OBABO and BCOCB integrators with a range of time-steps. For the OBABO integrator, the calculated correlation function is qualitatively incorrect for time-steps as large as 0.8 fs. For the BCOCB integrator, the resulting correlations functions are far more robust with respect to time-step. Although modest differences are seen in the exponential tail of the correlation function, the dynamics on vibrational timescales (see inset) is largely unchanged as the time-step is varied from 0.2 fs to 1.4 fs. Figure 2.6c further emphasizes this point by showing the absorption spectrum that is obtained from the BCOCB time-correlation functions with the various time-steps. To minimize bias, we avoided any smoothing of the spectra shown in panel c. It is 33 clearly seen that the librational and bending features (below 2500 cm−1 ) are graphically indistinguishable over the entire range of considered time-steps. To clarify the comparison for the stretching region above 3000 cm−1 , we smooth the raw spectra in that region by convolution against a Gaussian kernel with a width of 150 cm−1 (see inset). Again, the robustness of the simulated spectrum over this span of time-steps is excellent, with the only significant effect due to finite-time-step error being a slight blue-shifting of the OH stretching frequency for the results using a 1.4-fs time-step, which is nearly three times larger than the typical value employed for the OBABO scheme for simulations with 64 beads. Taken together, these results indicate that the BCOCB integrator provides an excellent description of both PIMD statistics and T-RPMD dynamics in realistic molecular systems, substantially improving the accuracy and stability of previously employed numerical integrators. 2.8 Summary In previous work [26], some of us showed that several widely used schemes for the non-preconditioned equations of motion of PIMD, including the widely used OBABO scheme, lack strong stability due to the use of exact free ring-polymer time evolution in the “A” sub-step, and we proved that this lack of strong stability gives rise to a lack of ergodicity in the thermostatted trajectories. We further showed that ergodicity can be restored by simply replacing the “A” sub-step with the Cayley transform. In the current chapter, we show that a completely distinct, yet equally important, pathology exists in the “B” sub-step of previously developed non-preconditioned PIMD integrators, due to the outsized effect of the external potential on the dynamics of the high-frequency ring-polymer modes. Specifically, we show that previous integrators (including OBABO, BAOAB, and OBCBO) yield a numerical stationary distribution for which the overlap with the exact stationary distribution vanishes in the infinite-bead limit. We then show that this pathology is completely avoided in the BCOCB scheme, and we further show that the pathology can be eliminated in the OBCBO scheme by suitably mollifying the “B” sub-step, yielding the dimensionfree non-preconditioned PIMD integrators, namely BCOCB, OMCMO, and OmCmO. Implementation of the dimension-free integration schemes involves no significant additional computational cost, no additional parameters, and no increase in algorithmic complexity in comparison to either OBABO or BAOAB. Furthermore, since the integrators considered here are all non-preconditioned, they can immedi- 34 Figure 2.6. Dynamical properties of liquid water computed using T-RPMD with the (panel a) OBABO and (panels b and c) BCOCB integration schemes. Panels a and b present the Kubo-transformed dipole autocorrelation function computed with various time-steps, and panel c presents the absorption spectrum from the BCOCB correlation function at each time-step. The inset to panel c presents the OH stretching region with smoothing. 35 ately be used for computing the equilibrium statistical properties as well as dynamical properties within the T-RPMD approximation. The numerical performance of the BCOCB scheme is particularly striking, yielding results that are markedly better in terms of accuracy and time-step stability than any of the other considered integrators. For liquid water, it is shown that BCOCB allows for time-steps as large as 1.4 fs while exhibiting minimal time-step error in the calculation of both equilibrium expectation values and the dipole absorption spectrum. Acknowledgments. We thank Andreas Eberle and Ondrej Marsalek for helpful discussions. N. B.-R. acknowledges support from the National Science Foundation under Award No. DMS-1816378 and the Alexander von Humboldt foundation. R. K., J. L. R.-R., and T. F. M. acknowledge support from the Department of Energy under Award No. DE-FOA-0001912 and the Office of Naval Research under Award No. N00014-10-1-0884. 2.9 Appendix A: Other splittings There are exactly four locally third-order accurate symmetric splitting schemes that involve one new force evaluation per integration step and that involve splitting the T-RPMD dynamics into Hamiltonian and thermostat sub-steps: OBABO, BAOAB, OABAO, and ABOBA. In Sec. 2.3, we quantified the properties of OBABO, BAOAB, and their Cayley modifications in the case of a harmonic external potential. The corresponding properties of the Cayley modifications of OABAO and ABOBA are given below. • OCBCO is exact in the velocity marginal, but the variance in the position 2 marginal is (βmn )−1 s2j,n,∆t where s2j,n,∆t = (4m − ∆t2 Λ)/(4Λ + 4mωj,n ); • CBOBC is exact in the position marginal, but the variance in the velocity 2 2 marginal is (βmn )−1 rj,n,∆t where rj,n,∆t = 4m/(4m − ∆t2 Λ). Numerical experiments confirmed these properties but did not show significant improvement in accuracy compared with BCOCB. Therefore, we did not include numerical results for these schemes. 36 Chapter 3 Generalization and optimization of dimension-free integrators for thermostatted ring-polymer molecular dynamics Adapted from J. L. Rosa-Raíces, J. Sun, N. Bou-Rabee, and T. F. Miller III, “A generalized class of strongly stable and dimension-free T-RPMD integrators”, The Journal of Chemical Physics 154, 024106 (2021) DOI: 10.1063/5.0036954. Abstract Recent work, some presented in Chapter 2, showed that strong stability and dimensionality freedom are essential for robust numerical integration of thermostatted ring-polymer molecular dynamics (T-RPMD) and path-integral molecular dynamics (PIMD), without which standard integrators exhibit non-ergodicity and other pathologies [J. Chem. Phys. 151, 124103 (2019); J. Chem. Phys. 152, 104102 (2020)]. In particular, the BCOCB scheme, obtained via Cayley modification of the standard BAOAB scheme, features a simple reparametrization of the free ringpolymer sub-step that confers strong stability and dimensionality freedom and has been shown to yield excellent numerical accuracy in condensed-phase systems with large time-steps. Here, we build upon these findings by benchmarking the performance of BCOCB against a broad class of numerical integrators that also exhibit strong stability and dimensionality freedom for commonly used T-RPMD friction schedules. In addition to considering equilibrium accuracy and time-step stability as in Chapter 2, we compare the integrators on the basis of their rates of convergence to equilibrium and their efficiency at evaluating equilibrium expectation values. We find BCOCB to be optimal within this generalized class with respect to asymptotic statistical accuracy and efficiency for various configuration-dependent observables, although other integrators within the generalized class perform better for the evaluation of velocity-dependent quantities and exhibit faster convergence to equilibrium 37 at large bead numbers. Numerical evidence indicates that this result, here proven analytically for the harmonic oscillator, also holds for a strongly anharmonic model of liquid water. Both analytical and numerical results indicate that BCOCB outperforms other strongly stable and dimension-free T-RPMD integrators in terms of configurational accuracy, asymptotic sampling efficiency, and time-step stability across both toy-model and real-world applications. 3.1 Introduction Path-integral molecular dynamics (PIMD) provides a practical and popular tool to simulate condensed-phase systems subject to strong nuclear quantum effects [1, 13, 43]. Based on the ring-polymer correspondence between quantum and classical Boltzmann–Gibbs statistics [14, 15], PIMD exploits the computational methods of molecular dynamics [44–47] to approximate quantum thermodynamics and kinetics through various classical models [12, 48–52]. Applications of PIMD include calculations of chemical reaction rates [53, 54], diffusion coefficients [16, 17], absorption spectra [42, 55], solid and liquid structure [56, 57], and equilibrium isotope effects [58, 59]. Many numerical integration schemes for PIMD are based on a symmetric Trotter splitting [27, 60] of the exact time-evolution operator, and feature a sub-step for free ring-polymer propagation [8, 11, 18]. Due to fast harmonic motions present in the free ring polymer, a strongly stable implementation of this sub-step is essential [61, 62]. Strong stability can be achieved by one of two approaches. The first approach introduces a preconditioned form of the equations of motion by modifying the ringpolymer mass matrix. Preconditioning improves the stability of the exact free ringpolymer update at the expense of consistent dynamics [3, 6, 8–11]. The second approach does not modify the ring-polymer mass matrix, leaving the dynamics nonpreconditioned [7, 18–21], and instead replaces the exact free ring-polymer update with a strongly stable approximation [26]. We apply the latter approach in the current chapter to thermostatted ring-polymer molecular dynamics (T-RPMD) [20], a nonpreconditioned variant of PIMD featuring an Ornstein–Uhlenbeck thermostat that approximately preserves the real-time dynamical accuracy of RPMD for quantum correlation functions of a wide range of observables [63]. In addition to strong stability of the free ring-polymer update, another basic requirement of a numerical integrator for T-RPMD is non-zero overlap between the numer- 38 ically sampled and exact ring-polymer configurational distributions in the limit of an infinite number of ring-polymer beads. As argued in Chapter 2, standard integrators fail to satisfy this requirement at any finite integration time-step, which motivates the introduction of dimension-free T-RPMD schemes that allow for accurate configurational sampling with large time-stepping and arbitrarily many ring-polymer beads. We showed in Chapter 2 that standard integrators can be made dimensionfree through the introduction of a suitable strongly stable ring-polymer update, and the current chapter investigates this finding in much greater generality. To this end, we introduce a function θ that defines the free ring-polymer update and deduce how the choice of θ impacts the properties and performance of the corresponding T-RPMD integrator. The case θ(x) = x, i.e., θ is the identity, corresponds to the exact free ring-polymer update. Therefore, to ensure second-order global accuracy, θ must approximate the identity near the origin, i.e., θ(0) = 0, θ′ (0) = 1 and θ′′ (0) = 0. Moreover, strong stability requires that the range of the function θ is within (0, π) for x > 0, and ergodicity and dimensionality freedom of the corresponding T-RPMD integrator impose additional requirements on θ. There are many choices of θ that fulfill the identified conditions, including θ(x) = 2 arctan(x/2) which leads to the BCOCB scheme introduced in Chapter 2. In fact, we find that this choice of θ may be optimal for the estimation of configurational averages via T-RPMD from the perspectives of accuracy and efficiency in general applications. This chapter is organized as follows. In Sec. 3.2 we recall T-RPMD and its time discretization, present the new function θ that determines the free ring-polymer update, and obtain sufficient conditions on θ to guarantee strong stability and dimensionality freedom of the corresponding T-RPMD integrator. In Sec. 3.3, we compare the performance of various θ in applications to the one-dimensional quantum harmonic oscillator and to a quantum-mechanical model of room-temperature liquid water. Section 3.4 summarizes the chapter, and subsequent appendices provide supporting mathematical proofs and list calculation details. 3.2 Theory 3.2.1 T-RPMD integrator desiderata We review the setting addressed by T-RPMD integrators, recalling details from Chapter 2 as necessary. Consider a one-dimensional molecular system with po- 39 tential energy function V (q), mass m, and Hamiltonian Ĥ as defined in Eq. (2.2). Following the arguments and definitions in Sec. 2.2, the equations of motion for the corresponding n-bead ring polymer held at constant temperature T by an Ornstein– Uhlenbeck thermostat are q̇(t) = v(t) , v̇(t) = − Ω2 q(t) + m−1 n F (q(t)) p 1/2 Ẇ (t) , − Γv(t) + 2β −1 m−1 n Γ (3.1) which reproduces Eq. (2.8). Here, W is an n-dimensional standard Brownian motion, q(t) = (q0 (t), . . . , qn−1 (t))T is the vector of positions for the n ring-polymer beads at time t ≥ 0, v(t) are the corresponding velocities, mn = m/n, β = (kB T )−1 , and F (q) = −∇Vnext (q) where Vnext is the bead-averaged potential defined in (2.6). Moreover, Ω2 is the n × n symmetric positive semi-definite matrix given in Eq. (2.10), diagonalizable by an n × n orthonormal real discrete Fourier transform matrix U as Ω = U diag (0, ω1,n , . . . , ωn−1,n ) U T , (3.2) where ωj,n is the jth Matsubara frequency defined in (2.11). Similarly, Γ in Eq. (2.8) is an n × n symmetric positive semi-definite friction matrix of the form Γ = U diag (0, γ1,n , . . . , γn−1,n ) U T , (3.3) where γj,n is the friction factor in the jth normal mode. In RPMD and T-RPMD calculations, one is often interested in the dynamics of Eq. (3.1) with initial conditions drawn from the stationary distribution with nonnormalized density exp(−βHn (q, v)), where Hn (q, v) is the ring-polymer Hamiltonian defined by Eqs. (2.5)–(2.7). We shall consider this case, as well as the case of initialization from non-equilibrium distributions, later in this chapter. Numerical integrators for Eq. (3.1) typically employ symmetric propagator splittings of the form given in Eq. (2.1), where the operator Ln = An + Bn + On includes contributions from the n-bead free ring-polymer motion (An ), the external potential (Bn ) and the thermostat (On ), and ∆t is a sufficiently small time-step. Recall from Chapter 2 that the standard RPMD integrator is recovered in the limit of zero coupling to the thermostat [18], and that Eq. (2.1) yields the OBABO scheme of Bussi and Parrinello [22] if a = 1 and the BAOAB scheme of Leimkuhler [25] if a = 0. Standard implementations of the T-RPMD splittings in Eq. (2.1) use the exact free ring-polymer propagator exp((∆t/2)An ) to evolve the free ring polymer; however, 40 previous work by some of us [26] showed that such implementations exhibit poor ergodicity if large numbers n of ring-polymer beads are employed in conjunction with large time-steps ∆t, and suggested replacing the exact ring-polymer propagator with its Cayley approximation [64] for improved performance. Subsequent work, presented in Chapter 2, introduced a Cayley-modified BAOAB scheme, denoted BCOCB, and presented numerical evidence that cemented the scheme as an improvement over standard BAOAB due to its superior equilibrium accuracy and time-step stability. Generalizing beyond the Cayley modification, the current chapter studies a family of modified BAOAB schemes that contains BCOCB and introduces others with similar theoretical guarantees. These BAOAB-like schemes are obtained by replacing the exact free ring-polymer update in Eq. (2.1) with approximations that endow the properties listed below. (P1) Strong stability. For a free ring polymer (i.e., for V ≡ const.), the integrator with γj,n = 0 is both strongly stable and second-order accurate in ∆t. (P2) Free ring-polymer ergodicity. For a free ring polymer, the integrator with γj,n > 0 is ergodic with respect to the distribution with density proportional to exp(−βHn0 (q, v)). (P3) Dimension-free stability. For a harmonically confined ring polymer (i.e., for V (q) = (Λ/2) q 2 ), the integrator with γj,n = 0 is stable for any n if ∆t leads to stable integration for n = 1. (P4) Dimension-free ergodicity. For a harmonically confined ring polymer, the integrator with γj,n > 0 and stable ∆t is ergodic with respect to its stationary distribution for any n. (P5) Dimension-free equilibrium accuracy. For a harmonically confined ring polymer, the integrator leaves invariant an accurate approximation of the distri
bution with density proportional to exp −(βmn /2) q T (Λ/m + Ω2 )q , with bounded error for any n. To obtain integrators satisfying properties (P1)–(P5), we introduce a function θ that defines the free ring-polymer update and then place conditions on θ that lead to the 41 desired properties. To this end, let " # iθ(ωj,n ∆t)/2 e 0 1/2 S j,n,∆t = Qj,n (3.4) Q−1 j,n , −iθ(ωj,n ∆t)/2 0 e   where Qj,n = iω1j,n −iω1j,n and essential properties of θ are determined in the sequel. We focus on T-RPMD schemes derived from the BAOAB splitting (i.e., a = 0 1/2 in Eq. (2.1)) with the exact free ring-polymer update replaced by S j,n,∆t . For such schemes, an integration time-step is comprised by the following sequence of substeps: B: Update velocities for half a step: v ← v + ∆t m−1 n F (q). 2 Convert bead Cartesian coordinates to normal modes using ϱ = U T q and φ = U T v . (3.5) A: Evolve the free ring polymer in normal-mode coordinates for half a step: " # " # ϱj ϱj 1/2 ← S j,n,∆t for 0 ≤ j ≤ n − 1 . φj φj O: Perform an Ornstein–Uhlenbeck velocity update for a full time-step: s 1 − e−2γj,n ∆t ξj , φj ← e−γj,n ∆t φj + βmn where ξj are independent standard normal random variables and 0 ≤ j ≤ n − 1. A: Evolve the free ring polymer in normal-mode coordinates for half a step: " # " # ϱj ϱj 1/2 ← S j,n,∆t for 0 ≤ j ≤ n − 1 . φj φj Convert back to bead Cartesian coordinates using the inverse of U , which is just its transpose since U is orthogonal. B: Update velocities for half a step: v ← v + ∆t m−1 n F (q). 2 In the remainder of this section, we identify conditions on the choice of θ that imply properties (P1)–(P5) for the corresponding T-RPMD integrator. Despite our focus on BAOAB-like splittings, we describe how the conditions on θ can be adjusted to construct integrators derived from the OBABO splitting (i.e., a = 1 in Eq. (2.1)) that satisfy properties (P1)–(P5). 42 3.2.2 Strong stability of RPMD with a constant external potential In this section, sufficient conditions on θ are identified to satisfy property (P1) in Sec. 3.2.1. Let V ≡ const. and γj,n = 0 for 1 ≤ j ≤ n − 1, corresponding to the free ring polymer. The jth normal mode [ϱj , φj ]T satisfies " # " # " # ϱ̇j ϱj 0 1 = Aj,n where Aj,n = . (3.6) 2 φ̇j φj −ωj,n 0 In this case, the algorithm from Sec. 3.2.1 reduces to a step of S j,n,∆t ≈ exp(∆tAj,n ), i.e., " # " # ϱj ϱj ← S j,n,∆t for 0 ≤ j ≤ n − 1 , (3.7) φj φj 1/2 1/2 where S j,n,∆t = S j,n,∆t S j,n,∆t follows from Eq. (3.4) and the function θ is such that property (P1) holds. We proceed to identify sufficient conditions on θ such that the corresponding free ring-polymer update satisfies property (P1). First note that for any function θ such 1/2 that θ(−x) = −θ(x) for x > 0, the structure of S j,n,∆t guarantees that the corresponding free ring-polymer update is reversible, symplectic, and preserves the free ring-polymer Hamiltonian Hn0 (q, v). Now, observe that S j,n,∆t is exact if θ(x) = x; therefore, second-order accuracy requires that θ approximates the identity near the origin, i.e., θ(0) = 0, θ′ (0) = 1, and θ′′ (0) = 0. (C1) Moreover, strong stability follows if the eigenvalues exp(±iθ(ωj,n ∆t)) of S j,n,∆t are distinct [26]; to this end we require that 0 < θ(x) < π for x > 0 . (C2) Jointly, conditions (C1) and (C2) guarantee that the update in Eq. (3.7) satisfies property (P1). There are many different choices of θ that obey these conditions, e.g., θ(x) = arctan(x), θ(x) = arccos(sech(x))1 and θ(x) = 2 arctan(x/2). The latter choice leads to the Cayley approximation of the free ring-polymer update, as 1/2 can be verified by substitution in Eq. (3.4) and comparison of the resulting S j,n,∆t with Eq. (17) in Chapter 2. Figure 3.1 compares the eigenvalues of S j,n,∆t with θ(x) = x and several choices of θ that meet conditions (C1) and (C2). 1 The function θ(x) = arccos(sech(x)) is not differentiable at the origin and hence, strictly speaking, does not satisfy condition (C1). Moreover, the function has even symmetry and hence fails to yield a reversible free ring-polymer update. These formal shortcomings can be fixed by multiplying the function by sign(x), which we implicitly do for this and other functions θ with similar features. 43 (a) ✓(x) = x (b) ✓(x) = arctan(x) (c) ✓(x) = 2 arctan(x/2) (d) ✓(x) = arccos(sech(x)) Figure 3.1. Eigenvalues of S j,n,∆t for 50 different time-steps between 0.05 and 5.0 (evenly spaced) and fixed Matsubara frequency ω = 3. The colors go from blue (smallest time-step) through green and yellow to red (largest time-step). In panel a, the eigenvalues rotate around the unit circle several times, which indicates that the corresponding S j,n,∆t is not always strongly stable. In panels b, c, and d, the eigenvalues are distinct and on the unit circle; thus the corresponding S j,n,∆t is strongly stable. 3.2.3 Ergodicity of T-RPMD with a constant external potential In this section, it is shown that condition (C2) implies property (P2) in Sec. 3.2.1. Let V ≡ const. and γj,n > 0 for 1 ≤ j ≤ n − 1, corresponding to the free ring polymer with a Ornstein–Uhlenbeck thermostat. In this case, the jth normal mode satisfies  " # " #  0 ϱ̇j ϱj = (Aj,n + Oj,n ) + q 2γj,n  , (3.8) Ẇj φ̇j φj βmn   where Oj,n = 00 −γ0j,n and Ẇj is a scalar white-noise. The solution [ϱj (t), φj (t)]T of Eq. (3.8) is an ergodic Markov process, and in the limit as t → ∞, its distribution 44 converges to the centered bivariate normal with covariance " # s2j,n 0 1 1 Σj,n = where s2j,n = 2 . βmn 0 1 ωj,n This distribution corresponds to the jth marginal of the free ring-polymer equilibrium distribution with density proportional to exp(−βHn0 (q, v)). The choice of γj,n > 0 in Eq. (3.8) determines the rate at which the associated Markov process converges to its stationary distribution if initialized away from it. When γj,n < 2ωj,n , the process is dominated by the deterministic Hamiltonian dynamics and is characterized as underdamped; on the other hand, when γj,n > 2ωj,n , the process is overdamped; and at the critical value γj,n = 2ωj,n the process is characterized as critically damped and converges to equilibrium fastest [65, 66]. This analytical result motivates the so-called PILE friction schedule [18, 20]. We specialize to this schedule in the remainder of this section and set γj,n = 2ωj,n for 1 ≤ j ≤ n − 1. The BAOAB-like update in Sec. 3.2.1, applied to Eq. (3.8) with PILE friction, can be written compactly as " # " # " # ϱj ϱj ξj 1/2 ← Mj,n,∆t + Rj,n,∆t for 0 ≤ j ≤ n − 1 , (3.9) φj φj ηj where ξj and ηj are independent standard normal random variables and we have introduced the 2 × 2 matrices " # 1 0 1/2 1/2 Mj,n,∆t = S j,n,∆t O j,n,∆t S j,n,∆t , O j,n,∆t = and 0 e−2ωj,n ∆t " # 0 0 1 − e−4ωj,n ∆t 1/2 1/2 Rj,n,∆t = S j,n,∆t (S j,n,∆t )T . βmn 0 1 1/2 Since S j,n,∆t and the Ornstein–Uhlenbeck update are individually preservative irrespective of the chosen θ, Eq. (3.9) exactly preserves the free ring-polymer equilibrium distribution for any choice of θ that satisfies (C1) and (C2). The ergodicity of the integrator specified by Eq. (3.9) depends entirely on the asymptotic stability of Mj,n,∆t , that is, on whether ∥Mkj,n,∆t ∥ → 0 as k → ∞, where ∥ · ∥ is a matrix norm. The matrix Mj,n,∆t is asymptotically stable if its spectral radius (i.e., the modulus of its largest eigenvalue) is smaller than unity [62], which depends on det(Mj,n,∆t ) = e−2ωj,n ∆t and Tr(Mj,n,∆t ) = cos(θ(ωj,n ∆t))(1 + e−2ωj,n ∆t ) . 45 In particular, the eigenvalues of Mj,n,∆t are both inside the unit circle if and only if | Tr(Mj,n,∆t )| < 1 + det(Mj,n,∆t ) < 2 ; a proof of this claim is provided in Appendix 3.5. This inequality reveals that condition (C2) implies property (P2). Moreover, if Tr(Mj,n,∆t )2 − 4 det(Mj,n,∆t ) ≤ 0, p then the spectral radius of Mj,n,∆t is minimal and equal to det(Mj,n,∆t ) = e−ωj,n ∆t ; this occurs when | cos(θ(ωj,n ∆t))| ≤ sech(ωj,n ∆t) for all ωj,n ∆t, which holds if the function θ satisfies arccos(sech(x)) ≤ θ(x) ≤ π − arccos(sech(x)) for x > 0 . (3.10) Any choice of θ that does not satisfy Eq. (3.10) will be overdamped in some modes, in the sense that the corresponding Mj,n,∆t will have a spectral radius strictly larger than e−ωj,n ∆t . The function θ(x) = arccos(sech(x)) saturates the (left) inequality in Eq. (3.10) while satisfying conditions (C1) and (C2), and hence provides a strongly stable and critically damped integrator for the thermostatted free ring polymer. As illustration of this, Fig. 3.2a shows that θ(x) = arctan(x) is overdamped for all modes whereas the Cayley angle θ(x) = 2 arctan(x/2) exhibits mixed damping. In contrast, the function θ(x) = arccos(sech(x)) preserves the critically damped behavior of its continuous counterpart under the PILE friction schedule. Figure 3.2b confirms that the spectral radius of Mj,n,∆t is minimal at θ(x) = arccos(sech(x)) for x > 0; consequently, this choice of θ optimizes the convergence of the integrator to stationarity. Conditions (C1) and (C2) also imply property (P2) for the OBABO-like update 1/2 1/2 associated with a compliant choice of θ, because the matrices S j,n,∆t O j,n,∆t S j,n,∆t 1/2 1/2 and O j,n,∆t S j,n,∆t O j,n,∆t have equal spectral radii. 3.2.4 Dimension-free stability of RPMD with a harmonic external potential In this section, we identify a condition on θ that yields property (P3) in Sec. 3.2.1. Let V (q) = (Λ/2) q 2 and γj,n = 0 for 1 ≤ j ≤ n − 1, corresponding to the nonthermostatted ring polymer with a harmonic external potential. In this case, the jth normal mode satisfies " # " # ϱ̇j ϱj = (Aj,n + B) (3.11) φ̇j φj 46 (x) underdamped 1 2 1 4 overdamped 0 1 spectral radius (a) overdamped 3 4 (b) (x) = arccos(sech(x)) (critical) (x) = 2arctan(x/2) (Cayley) (x) = arctan(x) (overdamped) 0.75 0.5 0.25 0 0 1 2 x 3 4 5 Figure 3.2. Spectral properties of the T-RPMD update for the free ring polymer for various choices of θ. Panel a plots the functions θ(x) = arccos(sech(x)), arctan(x) and arctan(x/2), and regions of overdamping and underdamping with PILE friction, separated at the locus of points where | cos(θ(x)) cosh(x)| = 1. The dynamics is underdamped in the gray region (| cos(θ(x)) cosh(x)| < 1), while in the white region (| cos(θ(x)) cosh(x)| > 1) the dynamics is overdamped. The function θ(x) = arctan(x) lies in the overdamped region for x > 0, whereas θ(x) = 2 arctan(x/2) is in the underdamped region for x ⪅ 2.4 and in the overdamped region otherwise. The function θ(x) = arccos(sech(x)), however, is critically damped for x > 0 and optimizes the convergence rate of the integrator. Panel b plots the spectral radius of Mj,n,∆t corresponding to each choice of θ as a function of x. where B =  0 0 −Λ/m 0  , and conserves the Hamiltonian Hj,n (ρj , φj ) =
mn 2 |φj |2 + (ωj,n + Λ/m)|ϱj |2 . 2 For this system, the BAOAB-like update in Sec. 3.2.1 reduces to " # " # ϱj ϱj ← Mj,n,∆t for 0 ≤ j ≤ n − 1 , φj φj where we have introduced the 2 × 2 matrices " 1 0 Mj,n,∆t = B1/2 S j,n,∆t B1/2 and B1/2 = −∆t(Λ/m)/2 1 (3.12) # . 47 This update may be interpreted as a symplectic perturbation of the free ring-polymer update in Eq. (3.7) due to the harmonic external potential [26], and conserves a modification of Hj,n that depends on the choices of θ and ∆t [67]. The update in Eq. (3.12) is stable if [9] 1 | Tr(Mj,n,∆t )| = max |Aj,n | < 1 , 0≤j≤n−1 2 0≤j≤n−1 max where Aj,n = cos(θ(ωj,n ∆t)) − (3.13) ∆t2 (Λ/m) sin(θ(ωj,n ∆t)) . 2 ωj,n ∆t Moreover, the 0th (i.e., centroid) mode, like the single-bead ring polymer, evolves through the velocity Verlet algorithm, whose stability requires that ∆t2 Λ/m < 4. Combining this requirement with condition (C2) yields a sufficient condition for Eq. (3.13) to hold at any n, 0 < θ(x) ≤ 2 arctan(x/2) for x > 0 . (C3) A proof of this result is provided in Sec. 3.6. The functions θ(x) = 2 arctan(x/2), θ(x) = arctan(x) and θ(x) = arccos(sech(x)) all satisfy condition (C3), which ensures that the corresponding RPMD integrator meets property (P3). 3.2.5 Dimension-free ergodicity and equilibrium accuracy of T-RPMD with a harmonic external potential In this section, it is shown that condition (C3) implies property (P4) in Sec. 3.2.1, and an additional condition is introduced to ensure that property (P5) holds. Let V (q) = (Λ/2) q 2 and γj,n = 2ωj,n for 1 ≤ j ≤ n − 1. In this case, the jth normal mode satisfies  " # " #  0 ϱ̇j ϱj = (Aj,n + B + Oj,n ) (3.14) + q 4ωj,n  . Ẇ φ̇j φj j βmn The solution [ϱj (t), φj (t)]T of Eq. (3.14) is an ergodic Markov process, and its distribution as t → ∞ converges to the centered bivariate normal with covariance matrix " # 2 s 0 1 1 j,n Σj,n = where s2j,n = ; (3.15) 2 βmn 0 1 Λ/m + ωj,n the associated position-marginal is the jth marginal of the ring-polymer configurational distribution with density exp(−(βmn /2) q T (Λ/m + Ω2 )q). 48 For this system, the BAOAB-like update in Sec. 3.2.1 is of the same form as Eq. (3.9) with 1/2 1/2 Mj,n,∆t = B1/2 S j,n,∆t O j,n,∆t S j,n,∆t B1/2 and " # 0 0 1 − e−4ωj,n ∆t 1/2 1/2 1/2 Rj,n,∆t = B S j,n,∆t (B1/2 S j,n,∆t )T . βmn 0 1 (3.16) As in the case of a constant external potential, the ergodicity of this integrator depends on the spectral radius of Mj,n,∆t . By Theorem 1 in Appendix 3.5 and the fact that det(Mj,n,∆t ) = e−2ωj,n ∆t and Tr(Mj,n,∆t ) = Aj,n (1 + e−2ωj,n ∆t ) , it follows that condition (C3) gives a simple and sufficient condition for ergodicity at any bead number n and hence implies property (P4) for the BAOAB-like update specified by Eqs. (3.9) and (3.16). Furthermore, because the matrix Mj,n,∆t of the corresponding OBABO-like update has equal trace and determinant, condition (C3) also guarantees property (P4) in that case.2 If condition (C3) holds, the BAOAB-like update is ergodic with respect to a centered bivariate normal distribution whose covariance matrix Σj,n,∆t satisfies the linear equation Σj,n,∆t = Mj,n,∆t Σj,n,∆t MT (3.17) j,n,∆t + Rj,n,∆t , for which the solution is " # s2j,n,∆t 0 1 Σj,n,∆t = 2 βmn 0 rj,n,∆t (3.18) where the variance in the position- and velocity-marginal is respectively (βmn )−1 s2j,n,∆t 2 and (βmn )−1 rj,n,∆t with  −1 Λ ωj,n ∆t/2 2 2 sj,n,∆t = ωj,n + and m tan (θ(ωj,n ∆t)/2) (3.19) 2 tan (θ(ω ∆t)/2) ∆t Λ j,n 2 . rj,n,∆t =1− 4m ωj,n ∆t/2 Because the tangent function is monotonically increasing on the range of θ specified by condition (C3), we have the correspondence 2 0 < s2j,n,∆t ≤ s2j,n and ∆t2 Λ 2 1− ≤ rj,n,∆t <1 4m (3.20) Condition (C3) may be viewed as a relaxation of the sufficient condition for ergodicity given in Eq. (18) of Chapter 2. Indeed, condition (C3) implies ergodicity irrespective of the Ornstein– Uhlenbeck friction schedule, whereas Eq. (18) in Chapter 2 does not imply ergodicity for friction schedules that lead to overdamped dynamics. 49 between the exact and numerical variances of the jth ring-polymer mode. Equation (3.19) reveals that θ(x) = 2 arctan(x/2) is the unique function that complies with condition (C3) and saturates the inequality s2j,n,∆t ≤ s2j,n in Eq. (3.20); consequently, the corresponding BAOAB-like scheme preserves the exact positionmarginal in all modes and trivially satisfies property (P5). The BCOCB integrator introduced in Chapter 2 corresponds to this choice of θ and thus uniquely provides optimal equilibrium position-marginal accuracy for harmonic external potentials. To identify other BAOAB-like schemes compliant with condition (C3) that satisfy property (P5), we examine the overlap between the numerical stationary positionmarginal distribution µn,∆t and the exact distribution µn where µn = n−1 Y j=1 N  s2j,n 0, βm n  and µn,∆t = n−1 Y j=1  s2  j,n,∆t N 0, βm . n Centroid-mode marginals have been suppressed in the definitions of µn and µn,∆t . A BAOAB-like scheme is dimension-free if it admits an n-independent upper bound on the distance dTV (µn , µn,∆t ) between µn and µn,∆t , where dTV is the total variation metric [40]. In particular, if we require x ≤ θ(x) ≤ 2 arctan(x/2) for x > 0 , 1 + |x| (C4) then we have the dimension-free bound r dTV (µn , µn,∆t ) < 4 ℏβ 3 ∆t ! ∆t2 Λ . m (3.21) A proof of this claim is provided in Sec. 3.8. Condition (C4) ensures that any BAOAB-like integrator with a compliant choice of θ meets property (P5). For OBABO-like schemes, the bound in condition (C4) must be tightened to guarantee non-zero overlap between µn and µn,∆t for arbitrarily large n. In particular, replacing 2 arctan(x/2) with min{2 arctan(x/2), C} for some C ∈ (0, π) in the upper bound of condition (C4) yields a n-independent bound on dTV (µn , µn,∆t ) for all compliant OBABO-like integrators, as can be shown through arguments similar to those in Sec. 3.8. Taken together, conditions (C1)–(C4) specify a family of BAOAB-like schemes with dimension-free stability, ergodicity and equilibrium accuracy for applications with harmonic external potentials. Numerical results in Sec. 3.3 suggest that the integrators exhibit similar properties in a more realistic setting with a strongly anharmonic external potential. 50 3.2.6 Dimension-free convergence to equilibrium of T-RPMD with a harmonic external potential Beyond ensuring ergodicity of the T-RPMD update in Eq. (3.16), condition (C3) leads to explicit dimension-free equilibration rates for compliant schemes. Theorem 4 in Sec. 3.7 proves this result in the infinite-friction limit for ring-polymer modes with arbitrarily high frequency. In detail, the theorem shows that the configurational transition kernel associated with the T-RPMD update of the jth mode in Eq. (3.16) is contractive in the 2-Wasserstein metric [68] and equilibrates any given initial distribution at a rate determined by the function θ, the (external) potential curvature Λ, and the (stable) time-step ∆t if condition (C3) holds. The rate in Theorem 4, though obtained in the infinite-friction limit, holds for finite friction coefficients γj,n leading to spectral radii ρ(Mj,n,∆t ) ≤ |Aj,n |, where Aj,n is defined in the display after Eq. (3.13) and |Aj,n | = limγj,n →∞ ρ(Mj,n,∆t ) is the spectral radius at infinite friction. To illustrate dimension-free convergence, Fig. 3.3 plots the 2-Wasserstein distance between the stationary configurational (i.e., position-marginal) distribution µn,∆t and the distribution µkn,∆t at the kth T-RPMD step evolved from a point mass at the origin using the schemes specified by θ(x) = arccos(sech(x)) (Fig. 3.3a), θ(x) = 2 arctan(x/2) (Fig. 3.3b), and θ(x) = arctan(x) (Fig. 3.3c) for a range of bead numbers n. These choices of θ respectively lead to overdamped, critical, and Cayley evolution of the thermostatted free ring polymer under PILE friction (see Sec. 3.2.3), and are identified accordingly in Fig. 3.3. The ring-polymer system considered in Fig. 3.3 approximates the O–H stretch dynamics in liquid water at room temperature with the parameters listed in Sec. 3.3.1. Velocity-marginals were initialized as in the setting of Theorem 4 (see Sec. 3.7), and the position of the jth ring-polymer mode at time k∆t follows a centered normal distribution with variance (βmn )−1 (skj,n,∆t )2 , where (skj,n,∆t )2 = (Mkj,n,∆t )212 + βmn k−1 X
Mℓj,n,∆t Rj,n,∆t (Mℓj,n,∆t )T 11 for k > 0 . ℓ=0 Here and throughout, the notation ( · )ij chooses the element in the ith row and jth column of the matrix enclosed in the parentheses. The 2-Wasserstein distances in Fig. 3.3 were evaluated using a well-known result for multivariate normal distributions [69]. 51 (a) critical t = 2.0 fs 10 10 6 0 104 10 5 t) / 2( 0, t) 20 (b) Cayley 2( kt, 10 4 103 102 10 3 101 100 30 2-Wasserstein distance time k t (fs) 30 10 2 10 1 (c) overdamped 100 20 10 0 16 64 256 1024 4096 bead number n Figure 3.3. Dimension-free convergence to equilibrium of BAOAB-like T-RPMD schemes with a harmonic external potential. The physical parameters of the ring-polymer system (i.e., Λ, m, and β) are listed in Sec. 3.3.1. Panels a, b, and c plot the normalized 2-Wasserstein distance between the configurational ring-polymer distribution at stationarity and at time k∆t, as evolved via various BAOAB-like schemes from an initial point-mass distribution. Regions with darker color indicate smaller 2-Wasserstein distance to stationarity, and black lines mark iso-distance contours. The contours plateau at some value of n for all tested schemes, which checks that they exhibit dimension-free convergence as predicted by Theorem 4. 52 Figures 3.3a and 3.3c clearly show that the critical and overdamped schemes converge at dimension-free rates, but this is less evident from Fig. 3.3b for the Cayley scheme. The latter scheme nonetheless displays an n-independent, and hence dimension-free, distance to stationarity at all times k∆t > 0, indicated by plateauing of the contour lines towards the right of Fig. 3.3b. The ladder-like pattern that precedes this plateau shows a transition from geometric (i.e., fast) to sub-geometric (i.e., slow) convergence upon introducing higher-frequency modes into the ring polymer. The transition manifests with the Cayley scheme because of its aggressive overdamping of the high-frequency modes, which is absent in the other two schemes (see Fig. 3.2). The example considered in this section illustrates that the equilibration timescale (e.g., the time until the 2-Wasserstein distance decays below 10−6 ) of the Cayley scheme at large n can dramatically exceed that of other BAOAB-like schemes. Although this negative feature may render the scheme impractical for pathological applications, we find in the next section that the Cayley scheme’s superior configurational sampling provides compelling justification for its preferred use in realistic settings. 3.3 Numerical results The current section provides numerical comparisons of the BAOAB-like T-RPMD integrators in Sec. 3.2, on applications featuring harmonic (Sec. 3.3.1) and anharmonic (Sec. 3.3.2) external potentials. Three choices of θ are considered in the numerical comparisons, namely θ(x) = arctan(x), θ(x) = arccos(sech(x)), and θ(x) = 2 arctan(x/2). These choices respectively lead to overdamped, critical, and Cayley evolution of the thermostatted free ring polymer under PILE friction (Sec. 3.2.3), and are identified accordingly throughout this section. It is borne out from the numerical comparisons that the Cayley scheme exhibits superior configurational sampling among the tested schemes in both applications. 3.3.1 One-dimensional quantum harmonic oscillator In the current section, we numerically integrate Eq. (3.1) with the harmonic potential p V (q) = (Λ/2) q 2 using PILE friction (i.e., Γ = 2Ω), m = 0.95 amu, Λ/m = 3886 cm−1 , and T = 298 K. This choice of physical parameters corresponds to a harmonic approximation of the Morse contribution to the O–H bond potential in 53 the q-TIP4P/F force field for water [41], and sets a least upper bound for the Tp RPMD stability interval at ∆tmax = 2/ Λ/m = 2.74 fs. The simulations reported throughout this section employ the time-step ∆t = 0.73 × ∆tmax = 2.00 fs. Figure 3.4 compares the accuracy and efficiency of various BAOAB-like T-RPMD schemes at equilibrium as a function of the bead number n. For a description of the numerical simulation and statistical estimation procedures used to generate the numerical data (filled circles) in Fig. 3.4, the reader is referred to Sec. 3.11. Figures 3.4a and 3.4c report the mean quantum kinetic energy at equilibrium as per the primitive and virial estimators, KEpri n (q) = n−1 X n mn ωn2 − (qj+1 − qj )2 and 2β j=0 2 n−1 1X 1 + KEvir (q) = (qj − q) ∂qj Vnext (q) , n 2β 2 j=0 (3.22) P where q = n1 n−1 j=0 qj is the centroid position of the n-bead ring polymer. For these two observables, Figs. 3.4b and 3.4d quantify the equilibrium sampling efficiency of the schemes in terms of the integrated autocorrelation time (or normalized asymptotic variance) [70–74]   PK−1 1 (k∆t) √ O (ξ ) lim Var n K→∞ k=0 aVar(On ) K
= Var(On ) Var On ∞ X =1+2 Cor(On (ξ (0) ), On (ξ (k∆t) )) , (3.23) k=1 (k∆t) , v (k∆t) )}∞ where On is an n-bead observable, {ξ (k∆t) }∞ k=0 = {(q k=0 a T-RPMD trajectory initialized at stationarity, Var(On ) the variance of On at equilibrium, and Cor(On (ξ (0) ), On (ξ (k∆t) )) the lag-k∆t autocorrelation of On along the T-RPMD trajectory. The integrated autocorrelation time of On is interpreted as the timescale over which adjacent observations along an equilibrium trajectory become statistically uncorrelated [70–74] and is hence a measure of the efficiency of a T-RPMD scheme at estimating the mean of On with respect to the numerically sampled equilibrium distribution. Figures 3.4a-d show that the scheme specified by the Cayley angle (orange) outperforms others in terms of both accuracy and efficiency at estimating the equilibrium average of the quantum kinetic energy observables. From the perspective of configurational accuracy, the optimality of the Cayley angle displayed in Figs. 3.4a and 3.4c is not surprising in light of the findings in Sec. 3.2.5. primitive KE KEpri n 54 equilibrium mean (kBT) 8 (a) 6 integrated autocorrelation time (fs) 3.5 (b) 3.0 4 2.5 2 virial KE KEvir n 0 8 filled: numerical empty: analytical (c) 1.5 3.5 6 3.0 4 2.5 2 2.0 t = 2.0 fs 0 8 1e 1 (e) classical KE KEcla n 2.0 (d) critical Cayley overdamped 1.5 104 (f) 103 6 102 4 2 101 4 16 64 256 1024 bead number n 100 4 16 64 256 1024 bead number n Figure 3.4. Performance at equilibrium of various BAOAB-like T-RPMD schemes applied to the one-dimensional quantum harmonic oscillator with physical parameters listed in Sec. 3.3.1. Panels a, c, and e, respectively, plot the equilibrium mean primitive kinetic energy, virial kinetic energy, and non-centroid classical kinetic energy per mode as a function of bead number n; the corresponding means in the exact infinite-bead limit are plotted as dashed lines. Panels b, d, and f plot the integrated autocorrelation times (Eq. (3.23)) of the respective observables. Exact (resp. numerically estimated) values of the plotted quantities are shown with empty (resp. filled) circles. Numerical estimates were obtained using the protocol described in Sec. 3.11. 55 Less expected are the results in Figs. 3.4b and 3.4d, which suggest that the Cayley angle is also optimal from the standpoint of configurational sampling efficiency for the quantum kinetic energy observables in Eq. (3.22). Sec. 3.9 supports this conjecture with an analytical result for harmonic external potentials. Figure 3.4e plots the mean classical kinetic energy at equilibrium as computed from the non-centroid ring-polymer velocities, KEcla n (v) = n−1 mn X 2 1 (vj − v 2 ) ≈ , 2(n − 1) j=0 2β (3.24) and Fig. 3.4f plots the corresponding integrated autocorrelation time as given by Eq. (3.23). For this observable, the equilibrium accuracy and efficiency of the Cayley scheme are significantly worse than those of the others as n increases. This is a consequence of the strongly overdamped behavior of Cayley T-RPMD at high frequencies (see Fig. 3.2), for which the integrator’s ergodicity degrades as its spectral radius approaches unity. Note that this shortcoming of the Cayley scheme presents no adverse implications to the equilibrium sampling of observables that exclusively depend on the ring-polymer configuration, as confirmed by Figs. 3.4a-d. In summary, Fig. 3.4 establishes that the T-RPMD scheme specified by the Cayley angle provides optimally accurate and efficient configurational sampling at equilibrium. To exploit this remarkable feature in practice, the scheme must manifest rapid converge to equilibrium when initialized away from it, as is necessary in most realistic applications of T-RPMD. Fortunately, Theorem 4 guarantees that any BAOABlike scheme compliant with conditions (C1)–(C4) features a contractive configurational transition kernel for any number of ring-polymer beads, and Fig. 3.3 in Sec. 3.2.6 illustrates this fact for the quantum harmonic oscillator considered in the current section. 3.3.2 Room-temperature liquid water While theoretical analysis and numerical tests of BAOAB-like T-RPMD schemes in Secs. 3.2.1 and 3.3.1 have focused on harmonic external potentials, the current section demonstrates that the resulting insights carry over to a realistic, strongly anharmonic model of room-temperature liquid water. Our test system is a periodic box containing 32 water molecules at a temperature of 298 K and a density of 0.998 g/cm3 , with potential energy described by the q-TIP4P/F force field [41]. As in Sec. 3.3.1, we compare the performance of various BAOAB-like T-RPMD schemes 56 H primitive KE 11 equilibrium mean (kBT) (a) 9 1.5 H virial KE (c) 6.0 O H bond PE (d) 10.0 7.5 5.0 5.5 12 H O H angle PE (b) 2.5 7 5 6.5 H classical KE 3.5 integrated autocorrelation time (fs) (e) 11 75 10 50 9 2.1 (g) 2.0 50 t = 1.4 fs (i) 1.5 102 critical Cayley overdamped 101 1.4 1.3 (h) 100 75 1.9 1.8 1.6 (f) 100 32 64 128 bead number n 100 32 (j) 64 128 bead number n Figure 3.5. Performance of various BAOAB-like T-RPMD schemes applied to q-TIP4P/F liquid water at room temperature. As a function of the bead number n and for a 1.4-fs time-step, panels a and c plot the equilibrium kinetic energy per H atom as per the primitive and virial estimators (Eq. (3.22)), and panels b and d plot the corresponding integrated autocorrelation times. Likewise, panels e and g plot the equilibrium potential energy per H2 O molecule due to the O−H-stretch and H−O−H-bend contributions, as defined in the q-TIP4P/F force field, and the corresponding autocorrelation times are plotted by panels f and h. Finally, panel i plots the classical kinetic energy per H atom computed from the noncentroid velocity estimator (Eq. (3.24)), and panel j plots the corresponding autocorrelation time. The numerical estimates and reference results (dashed lines) were obtained using the protocols described in Sec. 3.11. 57 for integrating the many-dimensional analogue of Eq. (3.1) with PILE friction, using the simulation time-step ∆t = 1.4 fs in all simulations. Numerical tests reported in Sec. 3.10 show that this value of ∆t closely approximates the upper limit of the Verlet (i.e., n = 1) stability interval for q-TIP4P/F liquid water. In agreement with Sec. 3.3.1, the experiments reveal that among the tested T-RPMD schemes, the Cayley scheme offers superior configurational sampling. For details on the numerical simulation and statistical estimation procedures used to generate the data presented in this section, the reader is referred to Sec. 3.11. Figure 3.5 compares the equilibrium accuracy achieved by the tested schemes in terms of the quantum and classical kinetic energy per hydrogen atom (Figs. 3.5a, 3.5c, and 3.5i) and the intramolecular potential energy per water molecule (Figs. 3.5e and 3.5g); also plotted are the respective integrated autocorrelation times as a function of bead number n. The kinetic energy estimates in Figs. 3.5a and 3.5c exhibit similar trends to those seen in Fig. 3.4 for the one-dimensional harmonic oscillator. In particular, the T-RPMD scheme specified by the Cayley angle outperforms others in terms of quantum kinetic energy accuracy as n increases, most outstandingly with a highly accurate primitive kinetic energy estimate despite the large time-step employed. Still in close agreement with the harmonic oscillator results, Figs. 3.5b and 3.5d show that the Cayley scheme displays the shortest integrated autocorrelation time among the tested schemes for the quantum kinetic energy observables. Similar trends manifest in the intramolecular potential energy averages and their autocorrelation times (Figs. 3.5e-h), where the Cayley scheme also achieves superior accuracy and efficiency. Finally, Figs. 3.5i and 3.5j confirm that the relative performance of the compared schemes in terms of velocity-marginal sampling is qualitatively consistent with the harmonic results. Taken together, the results in Fig. 3.5 suggest that the superiority of the Cayley scheme for configurational sampling, proven in the model setting of a harmonic external potential, is also reflected in realistic applications. In a final numerical test, Fig. 3.6 confirms that the sampling advantages of the Cayley T-RPMD scheme are obtained without downside in the estimation of dynamical quantities of typical interest. Specifically, Fig. 3.6b shows (unnormalized) infrared absorption spectra for room-temperature liquid water, computed from the 128-bead T-RPMD trajectories used to generate Fig. 3.5 using linear response theory and the T-RPMD approximation to real-time quantum dynamics [20, 42]. Linear response dictates that the absorption spectrum is proportional to ω 2 Ĩ(ω), where r 58 Ĩ(ω) = R dt e−iωt C̃µµ (t) is the Fourier transform of the quantum mechanical Kubotransformed dipole autocorrelation function C̃µµ (t). The latter is approximated within PNH2 O the T-RPMD framework [12, 16] by C̃µµ (t) ≈ NH1 O i=1 ⟨µi (t) · µi (0)⟩, where 2 NH2 O is the number of molecules in the liquid, µi (t) is the bead-averaged dipole moment vector of molecule i at time t, and the covariance ⟨µi (t) · µi (0)⟩ is estimated from a stationary T-RPMD trajectory as indicated in Sec. 3.11. Figure 3.6a plots the T-RPMD estimates of C̃µµ (t) leading to the absorption spectra in Fig. 3.6b. On the scale in which the absorption spectrum exhibits its key features, the spectra in Fig. 3.6b show very minor qualitative discrepancies. A similar conclusion holds for Fig. 3.6c, where the T-RPMD approximation of the Kubo-transformed velocPNH2 O ⟨v i (t) · v i (0)⟩ is plotted for the ity autocovariance function C̃vv (t) ≈ NH1 O i=1 2 three tested T-RPMD schemes. Collectively, these observations indicate that the accuracy of dynamical properties computed with BAOAB-like schemes is not significantly affected by the particular θ employed if conditions (C1)–(C4) in Sec. 3.2 are met. This result is expected due to the fact that the considered dynamical properties depend on bead-averaged (i.e., centroid-mode) coordinates, whose evolution is largely independent of the choice of θ under weak coupling between the centroid and non-centroid ring-polymer modes. 3.4 Summary Recent works showed that strong stability [26] and dimensionality [75] are essential features of a robust T-RPMD integration scheme that many widespread integrators do not possess. A T-RPMD scheme with these features, denoted BCOCB, was introduced via a simple and inexpensive Cayley modification of the free ring-polymer update (i.e., the “A” sub-step) of the standard BAOAB integrator. The BCOCB scheme was then shown to dramatically outperform BAOAB at estimating static and dynamic properties of various systems with remarkable accuracy at unprecedented time-steps. The current chapter generalizes beyond the Cayley modification by introducing a simple parameterization of the free ring-polymer update and a corresponding family of strongly stable and dimension-free modifications of the BAOAB scheme. Among these schemes lies BCOCB, which is found to exhibit superior configurational sampling despite exhibiting worse accuracy and efficiency for observables that depend on the non-centroid ring-polymer velocities. This conclusion is obtained theoretically via exhaustive analysis of a harmonic model, and numerically via simulation dipole autocovariance C (k t) (Debye2) 1500 velocity autocovariance Cvv(k t) (Å2 ps 2) 2I 6 infrared absorption ( ) (Debye2 cm 1) 59 (a) 5.4 5.0 4.7 4 2 0 0 5 (b) 0.0 0.1 0.2 10 15 lag k t (ps) 20 n = 128, t = 1.4 fs 1000 500 0 0 40 1000 2000 3000 4000 5000 frequency (cm 1) (c) critical Cayley overdamped 20 0 0.0 0.2 0.4 0.6 lag k t (ps) 0.8 1.0 Figure 3.6. Molecular dipole autocovariance function (panel a), corresponding infrared absorption spectrum (panel b), and molecular velocity autocovariance function (panel c) in room-temperature liquid water for various BAOAB-like T-RPMD schemes. The plotted quantities autocovariance exhibit minor qualitative discrepancies across schemes, which suggests that all schemes compliant with conditions (C1)–(C4) exhibit comparable accuracy in the computation of dynamical properties. Numerical estimates of the autocovariance functions were obtained using the protocol described in Sec. 3.11. 60 of a realistic quantum-mechanical model of liquid water at room temperature. In this way, the chapter supports the outstanding utility of BCOCB for accurate and efficient equilibrium simulation of condensed-phase systems with T-RPMD. To conclude, we stress that implementing BCOCB or any of the new dimension-free and strongly-stable schemes adds no cost, no algorithmic parameters and insignificant coding overhead relative to the standard BAOAB integrator. The modified integrators thus provide “turnkey” means to significantly improve the accuracy and stability of existing RPMD and T-RPMD implementations [76, 77]. Acknowledgments. J. L. R.-R. and J. S. contributed equally to the work presented in this chapter. The work was supported in part by the U.S. Department of Energy (DE-SC0019390) and the National Institutes of Health (R01GM125063) N. B.-R. acknowledges support by the Alexander von Humboldt foundation and the National Science Foundation (DMS-1816378). 3.5 Appendix A: Necessary and sufficient condition for eigenvalues of a 2 × 2 real matrix to be inside the unit circle This section provides a proof of the standard result used in Secs. 3.2.3 and 3.2.5 to infer ergodicity of the T-RPMD update for free and harmonically-confined ring polymers. Theorem 1. The spectral radius of a 2 × 2 real matrix M is strictly less than one if and only if | Tr(M )| < 1 + det(M ) < 2 . (3.25) Figure 3.7 plots eigenvalue pairs λ1 , λ2 that satisfy Eq. (3.25) for a fixed value of det(M ) = λ1 λ2 . Note that the spectral radius of M is minimized when λ1 and λ2 p are on the circle with radius r = det(M ). Proof. Let λ1 , λ2 be the (possibly complex) eigenvalues of M . By definition, the spectral radius of M is max(|λ1 |, |λ2 |) =: ρ. Since M is real, both Tr(M ) = λ1 + λ2 and det(M ) = λ1 λ2 are real. Thus, either: 1. λ1 , λ2 are a complex conjugate pair; or, 2. λ1 , λ2 are both real. 61 1.0 1 2 imaginary 0.5 0.0 0.5 1.0 1.0 0.5 0.0 real 0.5 1.0 Figure 3.7. All possible eigenvalue pairs λ1 , λ2 of a matrix M that satisfies Eq. (3.25) with det(M ) = λ1 λ2 = 1/4. The eigenvalue pairs either lie on the circle with radius r = 1/2 or are both real, and in the former case, the spectral radius of M is minimal. In the first case, λ1 = a + ib and λ2 = a − ib for some real numbers a and b with √ b ̸= 0, and hence, det(M ) = λ1 λ2 = a2 +b2 > 0, and ρ = |λ1 | = |λ2 | = a2 + b2 , p √ i.e., the eigenvalues lie on the circle with radius ρ = a2 + b2 = det(M ). In this case, the first inequality in Eq. (3.25) holds since b ̸= 0 implies | Tr(M )| = 2|a| < 2ρ ≤ 1 + ρ2 = 1 + det(M ) . Hence, Eq. (3.25) is equivalent to 1 + det(M ) < 2 or ρ < 1. In the second case, λ1 , λ2 are both real, and the condition | Tr(M )| < 1 + det(M ) is equivalent to 1 + λ1 λ2 + λ1 + λ2 = (1 + λ1 )(1 + λ2 ) > 0 , and 1 + λ1 λ2 − λ1 − λ2 = (1 − λ1 )(1 − λ2 ) > 0 . Together with det(M ) = λ1 λ2 < 1, these conditions imply ρ = max(|λ1 |, |λ2 |) < 1. 62 3.6 Appendix B: Stability condition for harmonic external potentials This section proves that condition (C3) implies property (P3), as claimed in Sec. 3.2.4. For notational brevity, we define A(x) := cos(θ(x)) − ∆t2 (Λ/m) sin(θ(x)) . 2 x Note that A(x) is equal to Aj,n in the display under Eq. (3.13) if x = ωj,n ∆t. Theorem 2. For any α⋆ > 0, (F2) implies (F1). (F1) For all Λ ≥ 0, m > 0 and ∆t > 0 satisfying ∆t2 Λ/m < α⋆ , the function θ satisfies |A(x)| < 1 for x > 0 . (F2) The function θ satisfies: 0 < θ(x) < 2 arctan(2x/α⋆ ) for x > 0 . Proof. Let α = ∆t2 (Λ/m). For notational brevity, define ϕα (x) := arctan(α/(2x)) for x > 0 . By the harmonic addition identity cos(θ) − tan(ϕα ) sin(θ) = cos(θ + ϕα ) , cos(ϕα ) note that (F1) can be rewritten as cos(θ(x) + ϕα (x)) < 1 for x > 0 , 0 < α < α⋆ . cos(ϕα (x)) (3.26) For 0 < θ(x) < π, Eq. (3.26) holds if and only if ϕα (x) < θ(x) + ϕα (x) < π − ϕα (x) , which can be rewritten as 0 < θ(x) < 2 arctan(2x/α) , where we used the identity π − 2 arctan(x) = 2 arctan(1/x) valid for x > 0 . (3.27) 63 Since arctan is monotone increasing, and 0 < α < α⋆ by assumption, we may conclude that 0 < θ(x) < 2 arctan(2x/α⋆ ) < 2 arctan(2x/α) . Thus, if (F1) holds, then Eq. (3.27) holds and therefore (F1) holds. Fix ϵ ∈ (0, 1). Since Theorem 2 is true for arbitrary α⋆ , if we take α⋆ = 4 − ϵ, then the theorem holds with ∆t2 Λ/m < 4 − ϵ in Theorem 2 (F1), and θ(x) < 2 arctan(2x/(4 − ϵ)) in Theorem 2 (F1). Since ϵ > 0 is arbitrary, and arctan is monotone increasing, we can conclude that the theorem holds with ∆t2 Λ/m < 4 and θ(x) ≤ 2 arctan(x/2). Summarizing, Corollary 3. Suppose that the function θ satisfies 0 < θ(x) ≤ 2 arctan(x/2) for x > 0 . Then for all Λ ≥ 0, m > 0 and ∆t > 0 satisfying ∆t2 Λ/m < 4, we have |A(x)| < 1 for x > 0 . 3.7 Appendix C: Dimension-free quantitative contraction rate for harmonic external potentials in the infinite-friction limit In the infinite-friction limit, Eq. (3.16) simplifies to " # 1 0 1/2 1/2 Mj,n,∆t = B1/2 S j,n,∆t S j,n,∆t B1/2 and 0 0 " # 0 0 1 1/2 1/2 B1/2 S j,n,∆t (B1/2 S j,n,∆t )T . Rj,n,∆t = βmn 0 1 The kth step of the corresponding T-RPMD integrator can be written compactly as " (k) # " (k−1) # " (k−1) # ϱj ϱj ξ 1/2 = Mj,n,∆t (k−1) + Rj,n,∆t j(k−1) , (k) φj φj ηj (k−1) (k−1) where ξj and ηj are independent standard normal random variables. Suppose (0) that the initial velocity is drawn from the Maxwell–Boltzmann distribution, φj ∼ N (0, (βmn )−1 ), and the initial position is drawn from an arbitrary distribution µj on (0) R, ϱj ∼ µj . Let pkj,n denote the k-step transition kernel of the position-marginal, (k) (0) i.e., µj pkj,n is the probability distribution of ϱj with ϱj ∼ µj . 64 The next theorem shows that starting from any two initial distributions µj and νj on R, the distance between the distributions µj pkj,n and νj pkj,n is contractive. We quantify the distance between these distributions in terms of the 2-Wasserstein metric. For two probability distributions µ and ν on R, the 2-Wasserstein distance between µ and ν is defined as: W2 (µ, ν) :=  inf E(|X − Y |2 ) X∼µ Y ∼ν 1/2 , where the infimum is taken over all bivariate random variables (X, Y ) such that X ∼ µ and Y ∼ ν [68], and where E(·) denotes a mathematical expectation. Theorem 4. Suppose that the function θ satisfies 0 < θ(x) ≤ 2 arctan(x/2) for x > 0 . Then for all k > 1, Λ ≥ 0, m > 0 and ∆t > 0 satisfying ∆t2 Λ/m < 4, and for all initial distributions µj and νj on R,   A(ωj,n ∆t)k−1 W2 (µj , νj ) if A(ωj,n ∆t) > 0, k k W2 (µj pj,n , νj pj,n ) ≤ (3.28) 1  W2 (µj , νj ) else.  2(k − 1) Proof. In the infinite-friction limit, the eigenvalues of Mj,n,∆t are {0, A(ωj,n ∆t)}, (0) (0) where A(x) is defined in Appendix 3.6. Let ϱj ∼ µj and ϱ̃j ∼ νj be an opti(0) (0) mal coupling of µj and νj , which implies that W2 (µj , νj ) = E(|ϱj − ϱ̃j |2 )1/2 . (0) (0) (k) (k) Conditional on ϱj and ϱ̃j , ϱj and ϱ̃j are Gaussian random variables with equal variances, but different means. By a well-known result for the 2-Wasserstein distance between Gaussian distributions [69], W2 (µj pkj,n , νj pkj,n )2 = |A(ωj,n ∆t)|2(k−1) (Mj,n,∆t )211 W2 (µj , νj )2 2 2(k−1) (1 + A(ωj,n ∆t)) = |A(ωj,n ∆t)| 4 W2 (µj , νj )2 , (3.29) where we used (Mj,n,∆t )11 = (1 + A(ωj,n ∆t))/2. Now we distinguish between two cases. In the case where A(ωj,n ∆t) > 0, we obtain the required result since |A(ωj,n ∆t)| < 1 by Corollary 3, and therefore, (1 + A(ωj,n ∆t))2 ≤1. 4 (3.30) 65 Otherwise, for −1 < A(ωj,n ∆t) ≤ 0 the quantity |A(ωj,n ∆t)|2(k−1) (1+A(ωj,n ∆t))2
2k is maximized at (1 − k)/k (k − 1)−2 , and therefore, |A(ωj,n ∆t)|2(k−1) 1 (1 + A(ωj,n ∆t))2 ≤ . 4 4(k − 1)2 (3.31) Inserting Eq. (3.30) and Eq. (3.31) into Eq. (3.29), and then taking square roots, gives the required result. 3.8 Appendix D: Total variation bound on the equilibrium accuracy error for harmonic external potentials In this section, we show that Eq. (3.21) follows from conditions (C1)–(C4) in the setting of Sec. 3.2.5. It is helpful to recall the quantities  πj  if j is even ,  ℏβ ωj = lim ωj,n = π(j + 1) (3.32) n→∞   else . ℏβ In the following, µj,n,∆t and µj,n respectively denote the jth factor of the product distributions µn,∆t and µn introduced in Sec. 3.2.5. Theorem 5. Suppose that the function θ satisfies conditions (C1)–(C4). Then for all Λ ≥ 0, m > 0 and ∆t > 0 satisfying ∆t2 Λ/m < 4, the total variation distance between µn and µn,∆t is bounded as in Eq. (3.21). Proof. Subadditivity of the total variation distance dTV between product distributions and its equivalence with the Hellinger distance [40] dH lead to the inequalities 2 dTV (µn , µn,∆t ) ≤ ≤ n−1 X n−1 X dTV (µj,n , µj,n,∆t )2 j=1 2 2 dH (µj,n , µj,n,∆t ) ≤ n−1 X 2(sj,n − sj,n,∆t )2 (s2j,n + s2j,n,∆t ) !2 2 X n−1  n−1 X s2j,n sj,n ≤ ≤ 1− 2 , 1− sj,n,∆t sj,n,∆t j=1 j=1 j=1 j=1 (3.33) where the second-to-last step uses Eq. (3.20) and the last step uses the elementary inequality (1 − x2 )2 ≥ (1 − x)2 valid for all x ≥ 0. Since tan increases superlinearly on the interval (0, π), we have θ(x)/2 ≤ tan(θ(x)/2) ≤ x/2 for x > 0 , 66 where the second inequality uses (C3). Consequently, the jth summand in Eq. (3.33) admits the bound !2   2 s2j,n ωj,n ∆t/2 Λ/m 1− 2 −1 = 2 sj,n,∆t ωj,n + Λ/m tan (θ(ωj,n ∆t)/2)  2  2 ωj,n ∆t ∆t Λ/m −1 ≤ (ωj,n ∆t)2 θ(ωj,n ∆t)  2 2 ∆t Λ 1 ≤ , m (ωj,n ∆t)2 where the last line uses the lower bound in (C4). Using that for any even positive integer n  2 X n−1 ∞ ∞ n−1 X X X 1 1 ℏβ 2 1 , < lim = < 2 2 2 2 n→∞ ω ω ω π j j,n j,n j j=1 j=1 j=1 j=1 where we used the definition in Eq. (3.32), the bound in Eq. (3.33) becomes 2 dTV (µn , µn,∆t ) <  ∆t2 Λ m 2  ℏβ π∆t 2 X ∞ 2 . 2 j j=1 Taking square roots and using the Riemann zeta function [78] to evaluate the infinite sum yields Eq. (3.21). 3.9 Appendix E: Asymptotic variance of kinetic energy observables for harmonic external potentials in the infinite-friction limit In Sec. 3.3.1, Figs. 3.4b and 3.4d show that the T-RPMD scheme specified by θ(x) = 2 arctan(x/2), which coincides with the Cayley-modified BAOAB scheme introduced in Chapter 2, provides the smallest integrated autocorrelation time (Eq. (3.23)) for quantum kinetic energy observables (Eq. (3.22)) among several schemes with properties (P1)–(P5). In this section, we show that this scheme minimizes an upper bound (Eq. (3.35)) on the integrated autocorrelation time of the quantum kinetic energy among all dimension-free and strongly-stable BAOAB-like schemes for harmonic external potentials. To this end, note that for a n-bead thermostatted ring polymer with external potential 67 Λ Vnext (q) = 2n |q|2 , Eq. (3.22) can be rewritten KEpri n (ϱ) = n−1 2 X mn ωj,n n − ϱ2j and 2β j=1 2 KEvir n (ϱ) = n−1 X 1 Λ 2 + ϱ 2β j=1 2n j (3.34) where ϱ is defined in Eq. (3.5). In the following, we denote both observables in Eq. (3.34) as KEn and distinguish between the two as needed. To control the integrated autocorrelation time of KEn , we need the stationary auto
correlation Cor KEn (ϱ(0) ), KEn (ϱ(k∆t) ) for k ≥ 0. Note that the distributions of n−1 ϱ(k∆t) and ϱ(0) are equal by stationarity, and that components (ϱj )j=0 are uncorrelated in a harmonic external potential. Thus, n−1 X
(k∆t) 2
(0) | , Cor KEn (ϱ(0) ), KEn (ϱ(k∆t) ) = wj,n Cor |ϱj |2 , |ϱj j=1 where (0) κ2j,n Var |ϱj |2
wj,n = Pn−1 (0) 2
2 i=1 κi,n Var |ϱi | and  2   mn ωj,n for KEpri , n 2 κj,n =  Λ for KEvir n . 2n If the evolution of the ring polymer is governed by the BAOAB-like update in Eq. (3.9), then the jth mode satisfies (k∆t) 2
(0) | Cov |ϱj |2 , |ϱj (0) 2 (k∆t) 2
| = Cor |ϱj | , |ϱj (0) 2
Var |ϱj | = (Mkj,n,∆t )211 , (k∆t) (k∆t) where we used that the phase [ϱj , φj ]T follows a centered Gaussian distribution with covariance given in Eq. (3.18) for all k ≥ 0. Therefore, in the infinitefriction limit where Mj,n,∆t is given in Appendix 3.7, the integrated autocorrelation time of KEn evaluates to n−1 ∞ X X aVar(KEn ) =1+2 wj,n (Mkj,n,∆t )211 Var(KEn ) j=1 k=1 ≤1+ 1 1 + A(ωj,n ∆t) max , 1≤j≤n−1 2 1 − A(ωj,n ∆t) (3.35) 68 where simplification of (Mkj,n,∆t )11 was aided by the Cayley–Hamilton theorem for 2 × 2 matrices [79], A(x) is defined in Appendix 3.6, and in the last line we used P that n−1 j=1 wj,n = 1. Equation (3.35) states that the integrated autocorrelation time of KEn can only be as small as that of the component |ϱj |2 exhibiting the slowest decorrelation at stationarity. Having derived Eq. (3.35), we now prove our claim for this section. Let x := ωj,n ∆t > 0 and α := ∆t2 Λ/m ∈ (0, 4). For fixed x and α, the function A(x) = α cos(θ(x)) − 2x sin(θ(x)) monotonically decreases toward −1 as the angle θ(x) in

creases toward π. Consequently, the function 1 + A(x) / 1 − A(x) decreases (toward 0) as θ(x) increases (toward π), but condition (C3) requires that θ(x) ≤ 2 arctan(x/2) to achieve stable evolution. Therefore, because it yields the largest stable angle, the choice θ(x) = 2 arctan(x/2), corresponding to the Cayley angle, minimizes the upper bound in Eq. (3.35). A similar argument can be made to support the conjecture, suggested by Fig. 3.4f, that the non-centroid velocity estimator for the classical kinetic energy KEcla n in Eq. (3.24), equivalently written KEcla n (φ) = n−1 mn X 2 φ 2(n − 1) j=1 j (3.36) with φ defined in Eq. (3.5), exhibits a maximal integrated autocorrelation time if the Cayley angle θ(x) = 2 arctan(x/2) is used. Indeed, the integrated autocorrelation time of this estimator is bounded by 1 − A(ωj,n ∆t) aVar(KEcla 1 n ) max , (3.37) ≤1+ cla 2 1≤j≤n−1 1 + A(ωj,n ∆t) Var(KEn )

where the function 1 − A(x) / 1 + A(x) monotonically increases as θ(x) approaches the largest stable angle for fixed x and α, which again corresponds to the Cayley angle. We stress that the conclusions of this section hold for arbitrary friction schedules despite our use of the infinite-friction limit in Eqs. (3.35) and (3.37). 3.10 Appendix F: Stability interval calibration for liquid water simulations This section describes the computational procedure used to identify ∆t = 1.4 fs as close to the upper bound of the stability interval of T-RPMD applied to q-TIP4P/F 69 liquid water at 298 K and 0.998 g/cm3 . Starting from an ensemble of 104 thermally distributed initial conditions obtained as per the following paragraph, we integrated an ensemble of 104 T-RPMD trajectories using the algorithm outlined in Sec. 3.2.1 in its single-bead realization (which is identical to velocity Verlet in classical MD [46]). We then counted the fraction of trajectories that did not “blow up,” or exhibit numerical overflow, throughout their duration for each tested time-step. A time-step was deemed “stable” if 99% or more of the corresponding trajectories did not exhibit numerical overflow over a 50-ps time period. A range of time-steps between 1 ps and 2 ps was tested, and the fraction of stable trajectories at each time-step is reported in Fig. 3.8. To avoid initialization bias in the stability interval estimation, thermalized initial phase-points were generated with a Metropolized Markov-chain Monte Carlo sampler targeted at the equilibrium configurational distribution of the liquid. Specifically, a randomized Hamiltonian Monte Carlo [9, 80] (rHMC) simulation of sufficient length was used to thermalize a crystalline configuration of the system at the target density, and 102 configurations were extracted from well-separated points along the rHMC trajectory. Each of these (approximately) independent draws from the equilibrium configurational distribution of the liquid at the target physical conditions was subsequently paired with 102 independent velocities drawn from the corresponding Maxwell–Boltzmann distribution, yielding 104 (approximately) independent draws from the phase space distribution of the liquid at thermal equilibrium. 3.11 Appendix G: Simulation and estimation details This section compiles simulation protocols and statistical estimation methods used to generate Figs. 3.4 for the one-dimensional quantum harmonic oscillator, and Figs. 3.5 and 3.6 for room-temperature liquid water. 3.11.1 One-dimensional quantum harmonic oscillator Numerical equilibrium averages and integrated autocorrelation times for the quantum harmonic oscillator were estimated by averaging over a 10-nanosecond T-RPMD trajectory integrated using the algorithm listed in Sec. 3.2.1, and initialized at an exact sample from the numerical stationary distribution (listed for the jth ring polymer mode in Eq. (3.19)) corresponding to the physical parameters (i.e., Λ, m, and β) and fraction of stable trajectories 70 1.0 0.5 0.0 1.3 1.4 1.5 1.6 1.7 1.8 timestep t (fs) Figure 3.8. Stability interval calibration for q-TIP4P/F room-temperature liquid water simulations. Data points correspond to the fraction of thermally initialized single-bead T-RPMD trajectories that “survived” (i.e., did not exhibit numerical overflow) over a 50-ps time period at the respective integration time-step ∆t. Error bars correspond to the standard error of the fraction of surviving trajectories across initialization points with different configurations. The gray dashed line marks the ≥ 99% trajectory survival threshold beyond which we deem a time-step “stable.” Note that no time-step beyond ∆t = 1.4 fs reaches the survival threshold. simulation parameters (i.e., n, ∆t, and the function θ) listed in Sec. 3.3.1. Specifically, the statistics reported in Fig. 3.4 were obtained by partitioning the T-RPMD trajectory into 10 disjoint blocks, estimating the equilibrium average and autocorrelation time within each block, and computing the sample mean and standard error among the resulting block estimates with 1000 bootstrap resamples. We now describe the formulas and methods used to obtain block estimates for the equilibrium mean and integrated autocorrelation time. The equilibrium average µOn of observable On within each block of the partitioned T-RPMD trajectory was estimated using the standard estimator [81] K−1 1 X (k∆t) µ̂On = O , K k=0 n (3.38) (k∆t) where K is the number of steps in the block (i.e., the block size) and On the value of On at the kth step within the block. Similarly, the lag-k∆t autocovariance COn (k∆t) was estimated using [81] ĈOn (k∆t) = K−k−1 X ℓ=0 (ℓ∆t) On
((ℓ+k)∆t)
− µ̂On On − µ̂On K −k for 0 ≤ k∆t ≤ (K − 1)∆t = 1 ns. The integrated autocorrelation time was 71 subsequently estimated using [71, 81] M X \ aVar ĈOn (k∆t) On (M ) = 1 + 2 , VarOn ĈOn (0) k=1 (3.39) where 0 < M ≤ K is a suitable cutoff. The choice of M is nontrivial, as it carries a trade-off between bias (more pronounced at small M ) and variance (more pronounced at large M ) [71]. To choose M judiciously, we follow the automatic windowing (AW) method described in Appendix C of Ref. [82]. The AW method dictates that M should correspond to the smallest lag that satisfies the inequality M ≥c \ aVar On (M ) , VarOn where the parameter c > 0 dictates the variance-bias trade-off in place of M , and is chosen as large as possible to reduce the bias of the estimator for a given variance threshold. Fig. 3.9 illustrates usage of the AW method for integrated autocorrelation time estimation, using trajectory data generated by the T-RPMD scheme with the choice θ(x) = 2 arctan(x/2) at n = 64 beads and ∆t = 2.0 fs, and focusing on the obvir cla servables KEpri n (black), KEn (red), and KEn (cyan) introduced in Sec. 3.3.1. The estimated integrated autocorrelation times are plotted with solid lines in Fig. 3.9a for various values of c, and the corresponding cutoffs M are plotted in Fig. 3.9b. Exact integrated autocorrelation times are plotted with dashed lines in Fig. 3.9a. Note that as c (and thus M ) increases, the estimates converge to the corresponding exact values at the expense of a larger variance, which can nonetheless be controlled by adjusting the block size K. 3.11.2 Room-temperature liquid water The equilibrium averages and integrated autocorrelation times reported in Fig. 3.5 were obtained by averaging over 10-nanosecond T-RPMD trajectories integrated for each considered bead number n, time-step ∆t, and function θ. All trajectories were initialized at an approximate sample from the corresponding numerical equilibrium distribution, obtained by thermalizing for 20 ps a classical (i.e., n = 1) configuration of the system into the n-bead ring-polymer phase space. Figure 3.10 below checks that this thermalization protocol indeed lead to near-equilibrium initialization of the T-RPMD trajectories. The reference equilibrium averages plotted with dashed lines int. autocorr. time aVarOn/VarOn (fs) 72 16 12 (a) 8 4 AW cutoff M t (fs) 0 100 75 On = KEpri n On = KEpri n On = KEcla n (b) harmonic oscillator n = 64, t = 2.0 fs 50 25 0 1 2 3 4 5 AW parameter c Figure 3.9. Integrated autocorrelation times of several observables of the one-dimensional harmonic oscillator in Sec. 3.3.1, estimated with the AW method. Trajectory data for the estimates was generated using the T-RPMD scheme with θ(x) = 2 arctan(x/2) at n = 64 beads and ∆t = 2.0 fs, and processed as described in the current section. Estimated (resp. vir cla exact) integrated autocorrelation times for observables KEpri n (black), KEn (red), and KEn (cyan) are shown in solid (resp. dashed) lines in panel a as a function of the windowing parameter c. Panel b plots the cutoffs determined by the choice of c for the three observables, where the linear relation between M ∆t and c at large values of the latter corroborates the non-spurious convergence of the autocorrelation time estimates. in Fig. 3.5 were obtained by averaging over a one-nanosecond, 256-bead staging PIMD [11] trajectory integrated at a 0.1-fs time-step with the mass and friction parameters recommended in Ref. [8], and initialized with the same protocol used for the T-RPMD simulations. The observables considered in Fig. 3.5 measure properties per H atom or per H2 O molecule, and thus the reported values are averages over estimates obtained for each simulated moiety. The equilibrium mean and integrated autocorrelation time of observable On for each moiety was estimated by partitioning the trajectory of the moiety into 10 disjoint 1-nanosecond blocks, evaluating Eqs. (3.38) and (3.39) within each block, and determining the sample mean and standard error among the block estimates with 1000 bootstrap resamples. The AW method [82] was applied to choose 73 a cutoff lag M ≤ 1 ns in Eq. (3.39), as illustrated in Fig. 3.9 for the harmonic oscillator application. The T-RPMD trajectories used to generate Fig. 3.5 also yielded Fig. 3.6, where panPNH2 O els a and c plot autocovariance functions of the form NH1 O i=1 Ōi (0)·Ōi (k∆t) , 2 where NH2 O = 32 is the number of simulated H2 O molecules and Ōi (k∆t) is the bead-averaged value of observable O (e.g., the molecular dipole moment or centerof-mass velocity) on the ith molecule at time k∆t along a stationary T-RPMD trajectory. The autocovariance Ōi (0) · Ōi (t) was estimated for the lags k∆t shown in Fig. 3.6 by Ōi (0) · Ōi (k∆t) ≈ K−k−1 X ℓ=0 (ℓ∆t) Ōi ((ℓ+k)∆t) · Ōi K −k , where K∆t = 1 ns is the length of each block in a partitioned 10-nanosecond T-RPMD trajectory. As with the results in Fig. 3.5, statistics for each molecule were obtained from block estimates via bootstrapping, and Figs. 3.6a and 3.6c report molecule-averaged statistics. Figure 3.10 validates the 20-ps thermalization interval used to initialize the trajectories that generated Figs. 3.5 and 3.6. In detail, Figs. 3.10a and 3.10b (resp., Figs. 3.10c and 3.10d) plot the non-equilibrium mean of the primitive and virial quantum kinetic energy per H atom (resp. the mean O − H bond and H − O − H angle potential energy per water molecule) as it approaches the equilibrium value in Figs. 3.5a and 3.5c (resp., Figs. 3.5e and 3.5g) for a 64-bead ring polymer at a 1.4 fs time-step with the considered choices of θ. At each time k∆t within the 20ps time interval, the non-equilibrium mean is estimated by averaging across 1000 independent trajectories initialized at a point-mass distribution on the n-bead ringpolymer phase space centered at the classical (i.e., n = 1) sample used to initialize the reported simulations. Within statistical uncertainty, the non-equilibrium mean for each observable converges to its equilibrium value within the 20-ps interval at visually indistinguishable rates across the tested choices of θ. 74 11.5 (a) 9 O H bond PE mean (kBT) H atom primitive KE mean (kBT) 10 8 7 6 10.5 10.0 5 6.0 5.8 5.6 n = 64, t = 1.4 fs 0.0 0.25 0.5 0.75 1.0 time k t (ps) H O H angle PE mean (kBT) H atom virial KE mean (kBT) 6.2 9.5 9.0 2.1 (c) 6.4 (b) 11.0 (d) 2.0 1.9 1.8 1.7 (x) = arccos(sech(x)) (x) = 2arctan(x/2) (x) = arctan(x) 1.6 0.0 2.5 5.0 7.5 10.0 time k t (ps) Figure 3.10. Convergence to equilibrium of the BAOAB-like schemes considered in Sec. 3.3.2 with n = 64 ring-polymer beads and a ∆t = 1.4 fs time-step. With respect to the non-equilibrium 64-bead configurational distribution evolved from a point mass at a classical (i.e., n = 1) configuration, panels a and c plot the mean kinetic energy per H atom for the n-bead system as per the primitive and virial estimators, respectively, for times up to 1.0 ps. Panels b and d, respectively, plot the non-equilibrium mean O−H-bond and H−O−H-angle potential energy per q-TIP4P/F water molecule, for times up to 10 ps. The lightly shaded interval around each curve corresponds to the standard error of the estimated non-equilibrium mean, computed with 1000 bootstrap resamples from a sample of 1000 independent trajectories. 75 References Cited in Part I [1] M. Parrinello and A. Rahman, “Study of an F center in molten KCl”, The Journal of Chemical Physics 80, 860 (1984). [2] G. J. Martyna, A. Hughes, and M. E. Tuckerman, “Molecular dynamics algorithms for path integrals at constant pressure”, The Journal of Chemical Physics 110, 3275 (1999). [3] P. Minary, G. J. Martyna, and M. E. Tuckerman, “Algorithms and novel applications based on the isokinetic ensemble. I. Biophysical and path integral molecular dynamics”, The Journal of Chemical Physics 118, 2510 (2003). [4] A. Beskos, G. Roberts, A. Stuart, and J. Voss, “MCMC methods for diffusion bridges”, Stochastics and Dynamics 8, 319 (2008). [5] A. Beskos, F. J. Pinski, J. M. Sanz-Serna, and A. M. Stuart, “Hybrid Monte Carlo on Hilbert spaces”, Stochastic Processes and their Applications 121, 2201 (2011). [6] J. Lu, Y. Lu, and Z. Zhou, “Continuum limit and preconditioned Langevin sampling of the path integral molecular dynamics”, Journal of Computational Physics 423, 109788 (2020). [7] Z. Zhang, X. Liu, Z. Chen, H. Zheng, K. Yan, and J. Liu, “A unified thermostat scheme for efficient configurational sampling for classical/quantum canonical ensembles via molecular dynamics”, The Journal of Chemical Physics 147, 034109 (2017). [8] J. Liu, D. Li, and X. Liu, “A simple and accurate algorithm for path integral molecular dynamics with the Langevin thermostat”, The Journal of Chemical Physics 145, 024103 (2016). [9] N. Bou-Rabee and J. M. Sanz-Serna, “Geometric integrators and the Hamiltonian Monte Carlo method”, Acta Numerica 27, 113 (2018). [10] N. Bou-Rabee and A. Eberle, “Two-scale coupling for preconditioned Hamiltonian Monte Carlo in infinite dimensions”, Stochastics and Partial Differential Equations: Analysis and Computations 9, 207 (2021). [11] M. E. Tuckerman, B. J. Berne, G. J. Martyna, and M. L. Klein, “Efficient molecular dynamics and hybrid Monte Carlo algorithms for path integrals”, The Journal of Chemical Physics 99, 2796 (1993). [12] I. R. Craig and D. E. Manolopoulos, “Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics”, The Journal of Chemical Physics 121, 3368 (2004). 76 [13] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, “Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space”, Annual Reviews of Physical Chemistry 64, 387 (2013). [14] R. P. Feynman, Quantum mechanics and path integrals (McGraw-Hill, 1965). [15] D. Chandler and P. G. Wolynes, “Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids”, The Journal of Chemical Physics 74, 4078 (1981). [16] T. F. Miller III and D. E. Manolopoulos, “Quantum diffusion in liquid water from ring polymer molecular dynamics”, The Journal of Chemical Physics 123, 154504 (2005). [17] T. F. Miller III and D. E. Manolopoulos, “Quantum diffusion in liquid parahydrogen from ring-polymer molecular dynamics”, The Journal of Chemical Physics 122, 184503 (2005). [18] M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, “Efficient stochastic thermostatting of path integral molecular dynamics”, The Journal of Chemical Physics 133, 124104 (2010). [19] M. Ceriotti, D. E. Manolopoulos, and M. Parrinello, “Accelerating the convergence of path integral dynamics with a generalized Langevin equation”, The Journal of Chemical Physics 134, 084104 (2011). [20] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, “How to remove the spurious resonances from ring polymer molecular dynamics”, The Journal of Chemical Physics 140, 234116 (2014). [21] M. Rossi, V. Kapil, and M. Ceriotti, “Fine tuning classical and quantum molecular dynamics using a generalized Langevin equation”, Journal of Chemical Physics 148, 102301 (2018). [22] G. Bussi, D. Donadio, and M. Parrinello, “Canonical sampling through velocity rescaling”, The Journal of Chemical Physics 126, 014101 (2007). [23] N. Bou-Rabee and H. Owhadi, “Long-run accuracy of variational integrators in the stochastic context”, SIAM Journal on Numerical Analysis 48, 278 (2010). [24] N. Bou-Rabee and E. Vanden-Eijnden, “Pathwise accuracy and ergodicity of metropolized integrators for SDEs”, Communications on Pure and Applied Mathematics 63, 655 (2010). [25] B. Leimkuhler and C. Matthews, “Rational construction of stochastic numerical methods for molecular sampling”, Applied Mathematics Research eXpress 2013, 34 (2013). 77 [26] R. Korol, N. Bou-Rabee, and T. F. Miller III, “Cayley modification for strongly stable path-integral and ring-polymer molecular dynamics”, The Journal of Chemical Physics 151, 124103 (2019). [27] H. F. Trotter, “On the product of semi-groups of operators”, Proceedings of the American Mathematical Society 10, 545 (1959). [28] T. Matsubara, “A new approach to quantum-statistical mechanics”, Progress of Theoretical Physics 14, 351 (1955). [29] D. Li, X. Han, Y. Chai, C. Wang, Z. Zhang, Z. Chen, J. Liu, and J. Shao, “Stationary state distribution and efficiency analysis of the Langevin equation via real or virtual dynamics”, The Journal of Chemical Physics 147, 184104 (2017). [30] S. Kakutani, “On equivalence of infinite product measures”, Annals of Mathematics 49, 214 (1948). [31] V. I. Bogachev, Gaussian measures, Vol. 62, Mathematical Surveys and Monographs (American Mathematical Society, 1998). [32] A. Pérez and M. E. Tuckerman, “Improving the convergence of closed and open path integral molecular dynamics via higher order Trotter factorization schemes”, The Journal of Chemical Physics 135, 064104 (2011). [33] O. Marsalek, P.-Y. Chen, R. Dupuis, M. Benoit, M. Méheut, Z. Bačić, and M. E. Tuckerman, “Efficient calculation of free energy differences associated with isotopic substitution using path-integral molecular dynamics”, Journal of Chemical Theory and Computation 10, 1440 (2014). [34] J. A. Barker, “A quantum-statistical Monte Carlo method; path integrals with boundary conditions”, The Journal of Chemical Physics 70, 2914 (1979). [35] M. F. Herman, E. J. Bruskin, and B. J. Berne, “On path integral Monte Carlo simulations”, The Journal of Chemical Physics 76, 5150 (1982). [36] B. García-Archilla, J. M. Sanz-Serna, and R. D. Skeel, “Long-time-step methods for oscillatory differential equations”, SIAM Journal on Scientific Computing 20, 930 (1998). [37] E. Hairer and C. Lubich, “Long-time energy conservation of numerical methods for oscillatory differential equations”, SIAM Journal on Numerical Analysis 38, 414 (2000). [38] J. M. Sanz-Serna, “Mollified impulse methods for highly oscillatory differential equations”, SIAM Journal on Scientific Computing 46, 1040 (2008). [39] R. I. McLachlan and A. Stern, “Modified trigonometric integrators”, SIAM Journal on Scientific Computing 52, 1378 (2014). [40] A. L. Gibbs and F. E. Su, “On choosing and bounding probability metrics”, International Statistical Review 70, 419 (2002). 78 [41] S. Habershon, T. E. Markland, and D. E. Manolopoulos, “Competing quantum effects in the dynamics of a flexible water model”, The Journal of Chemical Physics 131, 024501 (2009). [42] S. Habershon, G. S. Fanourgakis, and D. E. Manolopoulos, “Comparison of path integral molecular dynamics methods for the infrared absorption spectrum of liquid water”, The Journal of Chemical Physics 129, 074501 (2008). [43] T. E. Markland and M. Ceriotti, “Nuclear quantum effects enter the mainstream”, Nature Reviews Chemistry 2, 0109 (2018). [44] D. Frenkel and B. Smit, Understanding molecular simulation: From algorithms to applications, 2nd (Academic Press, 2002). [45] D. C. Rapaport, The art of molecular dynamics simulation, 2nd (Cambridge University Press, 2004). [46] B. Leimkuhler and C. Matthews, Molecular dynamics, Vol. 39, Interdisciplinary Applied Mathematics (Springer International Publishing, 2015). [47] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, 2nd edition (Oxford University Press, 2017). [48] J. Cao and G. A. Voth, “The formulation of quantum statistical mechanics based on the Feynman path centroid density. II. Dynamical properties”, The Journal of Chemical Physics 100, 5106 (1994). [49] J. Liu, “Path integral Liouville dynamics for thermal equilibrium systems”, The Journal of Chemical Physics 140, 224107 (2014). [50] T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, “Boltzmannconserving classical dynamics in quantum time-correlation functions: “Matsubara dynamics””, The Journal of Chemical Physics 142, 134103 (2015). [51] T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, “Communication: Relation of centroid molecular dynamics and ring-polymer molecular dynamics to exact quantum dynamics”, The Journal of Chemical Physics 142, 191101 (2015). [52] J. R. Cendagorta, Z. Bačić, and M. E. Tuckerman, “An open-chain imaginarytime path-integral sampling approach to the calculation of approximate symmetrized quantum time correlation functions”, The Journal of Chemical Physics 148, 102340 (2018). [53] I. R. Craig and D. E. Manolopoulos, “Chemical reaction rates from ring polymer molecular dynamics”, The Journal of Chemical Physics 122, 084106 (2005). [54] I. R. Craig and D. E. Manolopoulos, “A refined ring polymer molecular dynamics theory of chemical reaction rates”, The Journal of Chemical Physics 123, 034102 (2005). 79 [55] A. Witt, S. D. Ivanov, M. Shiga, H. Forbert, and D. Marx, “On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy”, The Journal of Chemical Physics 130, 194510 (2009). [56] J. A. Morrone and R. Car, “Nuclear quantum effects in water”, Physical Review Letters 101, 017801 (2008). [57] B. Cheng, E. A. Engel, J. Behler, C. Dellago, and M. Ceriotti, “Ab initio thermodynamics of liquid and solid water”, Proceedings of the National Academy of Sciences 116, 1110 (2019). [58] T. Zimmermann and J. Vaníček, “Path integral evaluation of equilibrium isotope effects”, The Journal of Chemical Physics 131, 024111 (2009). [59] D. L. Eldridge, R. Korol, M. K. Lloyd, A. C. Turner, M. A. Webb, T. F. Miller III, and D. A. Stolper, “Comparison of experimental vs theoretical abundances of 13 CH3 D and 12 CH2 D2 for isotopically equilibrated systems from 1 to 500 ◦ C”, ACS Earth and Space Chemistry 3, 2747 (2019). [60] G. Strang, “On the construction and comparison of difference schemes”, SIAM Journal on Numerical Analysis 5, 506 (1968). [61] M. P. Calvo and J. M. Sanz-Serna, “Instabilities and inaccuracies in the integration of highly oscillatory problems”, SIAM Journal on Scientific Computing 31, 1653 (2009). [62] V. I. Arnol’d, Mathematical methods of classical mechanics, 2nd edition, Vol. 60, Graduate Texts in Mathematics (Springer Science & Business Media, 2013). [63] B. J. Braams and D. E. Manolopoulos, “On the short-time limit of ring polymer molecular dynamics”, The Journal of Chemical Physics 125, 124105 (2006). [64] N. Bou-Rabee, “Cayley splitting for second-order Langevin stochastic partial differential equations”, Preprint (2017). [65] G. Metafune, D. Pallara, and E. Priola, “Spectrum of Ornstein-Uhlenbeck operators in Lp spaces with respect to invariant measures”, Journal of Functional Analysis 196, 40 (2002). [66] G. A. Pavliotis, Stochastic Processes and Applications, Vol. 60, Texts in Applied Mathematics (Springer, 2014). [67] J. M. Sanz-Serna and M. P. Calvo, Numerical Hamiltonian problems, Vol. 7, Applied mathematics and mathematical computation (Chapman & Hall, 1994). [68] C. Villani, Optimal transport: Old and new, Vol. 338, Comprehensive Studies in Mathematics (Springer Science & Business Media, 2008). [69] C. R. Givens and R. M. Shortt, “A class of Wasserstein metrics for probability distributions”, Michigan Mathematical Journal 31, 231 (1984). 80 [70] C. J. Geyer, “Practical Markov chain Monte Carlo”, Statistical Science 7, 473 (1992). [71] A. Sokal, “Monte Carlo methods in statistical mechanics: Foundations and new algorithms”, in Functional integration: Basics and applications, edited by C. DeWitt-Morette, P. Cartier, and A. Folacci (Springer Science & Business Media, 1997), pp. 131–192. [72] S. Asmussen and P. W. Glynn, Stochastic simulation: Algorithms and analysis, Vol. 57, Stochastic Modelling and Applied Probability (Springer, 2007). [73] R. D. Skeel and Y. Fang, “Comparing Markov chain samplers for molecular simulation”, Entropy 19, 561 (2017). [74] Y. Fang, Y. Cao, and R. D. Skeel, “Quasi-reliable estimates of effective sample size”, preprint (2017). [75] R. Korol, J. L. Rosa-Raíces, N. Bou-Rabee, and T. F. Miller III, “Dimensionfree path-integral molecular dynamics without preconditioning”, The Journal of Chemical Physics 152, 104102 (2020). [76] Y. V. Suleimanov, J. W. Allen, and W. H. Green, “RPMDrate: Bimolecular chemical reaction rates from ring polymer molecular dynamics”, Computer Physics Communications 184, 833 (2013). [77] V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, B. A. Helfrecht, P. Juda, S. P. Bienvenue, W. Fang, J. Kessler, I. Poltavsky, S. Vandenbrande, J. Wieme, C. Corminboeuf, T. D. Kühne, D. E. Manolopoulos, T. E. Markland, J. O. Richardson, A. Tkatchenko, G. A. Tribello, V. Van Speybroeck, and M. Ceriotti, “i-PI 2.0: A universal force engine for advanced molecular simulations”, Computer Physics Communications 236, 214 (2019). [78] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions, with formulas, graphs, and mathematical tables (Dover Publications, Inc., 1965). [79] T. Andreescu, Essential linear algebra with applications (Birkhäuser Basel, 2016). [80] N. Bou-Rabee and J. M. Sanz-Serna, “Randomized Hamiltonian Monte Carlo”, Annals of Applied Probability 27, 2159 (2017). [81] M. B. Priestly, Spectral analysis and time series (Academic Press, 1981). [82] N. Madras and A. Sokal, “The pivot algorithm: A highly efficient Monte Carlo method for the self-avoiding walk”, Journal of Statistical Physics 50, 109 (1988). Part II PARALLEL-IN-TIME STOCHASTIC MOLECULAR DYNAMICS 82 Chapter 4 Path-accelerated stochastic molecular dynamics: Parallel-in-time integration using path integrals Adapted from J. L. Rosa-Raíces, B. Zhang, and T. F. Miller III, “Path-accelerated stochastic molecular dynamics: Parallel-in-time integration using path integrals”, The Journal of Chemical Physics 151, 164120 (2019) DOI: 10.1063/1.5125455. Abstract Massively parallel computer architectures create new opportunities for the performance of long-timescale stochastic molecular dynamics (MD) simulations. Here, we introduce the path-accelerated molecular dynamics (PAMD) method that takes advantage of distributed computing to reduce the wall-clock time of MD simulation via parallelization with respect to stochastic MD time-steps. The marginal distribution for the time evolution of a system is expressed in terms of a path integral, enabling the use of path sampling techniques to numerically integrate MD trajectories. By parallelizing the evaluation of the path action with respect to time and by initializing the path configurations from a non-equilibrium distribution, the algorithm enables significant speedups in terms of the length of MD trajectories that can be integrated in a given amount of wall-clock time. The method is demonstrated for Brownian dynamics, although it is generalizable to other stochastic equations of motion. We apply the method to two simple systems, a harmonic oscillator and a Lennard–Jones liquid, and we show that in comparison to the conventional Euler integration scheme for Brownian dynamics, the new method can reduce the wallclock time for integrating trajectories of a given length by more than three orders of magnitude in the former system and more than two in the latter. This new method for parallelizing MD in the dimension of time can be trivially combined with algorithms for parallelizing the MD force evaluation to achieve further speedup. 83 4.1 Introduction Molecular dynamics (MD) [1, 2] is the central tool for simulating chemical, biological, and materials systems, with new algorithms and hardware expanding the range of accessible scales in space and in time [3–5]. Faster processors have played an important role in this expansion, although the most dramatic improvements in recent years have come from the number of available processors, rather than the clockspeed of the individual cores [6, 7]. In particular, highly multi-threaded computer architectures have been used to parallelize the MD force evaluation, greatly reducing the wall-clock time needed to perform an individual MD step [8–13]. However, despite this progress in the parallelization of MD simulations with respect to the force evaluations (i.e., in space), less attention has been dedicated to the notion of parallelization with respect to the MD time-steps (i.e., in time). The sequential nature of MD (i.e., the need to have access to a given time-step before the next time-step can be computed) would seem to discount the possibility of exploiting parallelization in time; nonetheless, methods for parallel-in-time integration are being developed and applied to MD simulation. Most approaches [14–17] are based on a prediction-correction paradigm that combines fine (i.e., accurate and expensive) and coarse (i.e., inaccurate and inexpensive) solvers to iteratively refine approximations of a trajectory in a convergent and parallel-in-time fashion. A range of coarse solvers and iteration schemes have been employed to evaluate MD trajectories of molecular systems with parallelization in the time domain [18–23], leading to order-of-magnitude reductions in the wall-clock time-to-solution with respect to sequential integration at the fine level of accuracy. Schemes for approximate longtimescale integration via trajectory splicing are an alternative route to parallelization in time, yielding accurate time evolution for systems that exhibit strong timescale separation on well-characterized regions of the potential energy landscape [24]. The current chapter takes a different approach to parallelizing MD in time, focusing on stochastic MD equations of motion for systems in contact with a thermal bath. We demonstrate that by working with ensembles of trajectories in a pathintegral framework, multiple processors can be employed to reduce the wall-clock time needed to evolve a stochastic MD trajectory of arbitrary length, without resorting to parallelization of the MD force evaluation. This method of parallelization for MD trajectories is independent of, and thus entirely complementary to, parallelization of the MD force evaluations, and it creates new opportunities to harness large 84 numbers of available computer processors for the generation of long-timescale MD trajectories. 4.2 Method 4.2.1 MD integration based on path distributions While our approach generalizes to other stochastic equations of motion, we illustrate it in the current chapter for Brownian (i.e., overdamped Langevin) dynamics under potential V at temperature β −1 , √ (4.1) ẋ(t) = −γ −1 V ′ (x(t)) + 2D ẇ(t) , where the diffusion coefficient D and the friction coefficient γ are related by the Einstein relation D = (βγ)−1 , and where w(t) is a standard Wiener process. MD trajectories can be generated by discretizing Eq. (4.1) with various numerical integration schemes [25–29], such as the forward Euler algorithm [2] √ x(t + dt) − x(t) = −γ −1 V ′ (x(t)) dt + 2Ddt ξ , (4.2) where dt is the discretization time-step and ξ is a standard Gaussian random variate. The marginal distribution associated with time evolution of the system by dt according to Eq. (4.2) is [30] K(x(t + dt)|x(t); dt) ∝ (  2 ) dt x(t + dt) − x(t) V ′ (x(t)) + exp − , 4D dt γ (4.3) such that the likelihood of a MD trajectory of length T = N dt that evolves the system along positions X = {x(t0 ), x(t1 ), . . . , x(tN )} at times tn = t + n dt is N −1 Y n=0 K(x(tn+1 )|x(tn ); dt) ≡ e−SN [X] , (4.4) where SN [X] is the action associated with the MD trajectory and, for the case of Eq. (4.1), is equivalent to the Onsager–Machlup action functional in a suitably defined continuous-time limit [31, 32]. From Eq. (4.4), the position of the time-evolved system at time T has a marginal distribution given by the path integral w w K(x(T + t)|x(t); T ) ∝ dx1 · · · dxN −1 e−SN [X] , (4.5) R R 85 where xn = x(tn ). It is clear that this path-integral formulation of the ensemble of MD trajectories provides an equivalent description of the time evolution of the system as Eq. (4.2); for this reason, trajectories generated using the scheme introduced in this section are correctly referred to as MD trajectories. Numerous studies have explored path-based formulations of MD, with variations of the underlying equation of motion and of the discretization of the action [33–45]. Setting aside issues of efficiency until Sec. 4.2.2, we note that the path-integral formulation of the marginal distribution for the time-evolved system offers a simple MD integration scheme, illustrated in Fig. 4.1. First, sampling from the distribution of paths of length T , with likelihood given by Eq. (4.4), is performed using Monte Carlo (MC) or related methods (Fig. 4.1A) [46–48]; by drawing a random realization from this distribution, we obtain a segment of MD trajectory from time 0 to time T (illustrated by the heavy orange path in Fig. 4.1A). Then, by shifting from x(0) to x(T ) along the sampled path, we resolve a trajectory from x(0) to x(T ) (represented by the heavy green path in Fig. 4.1B) that is statistically equivalent to a realization from the Euler algorithm defined in Eq. (4.2). After shifting the left endpoint of the path from x(0) to x(T ), we restart the path sampling to extend the trajectory from time T to time 2T . Iteration of this scheme will lead to the numerical integration of a MD trajectory of arbitrary length in time. Figure 4.2 illustrates a generalization of the integration scheme presented in Fig. 4.1. Figure 4.2A repeats Fig. 4.1A; we first sample a path of length T that is discretized into Npath time-steps (where Npath = T /dt) to obtain a realization of the path that is consistent with the marginal distribution of the time-evolved system for each time ∆t ≤ T . Then, in Fig. 4.2B, we shift the left endpoint of the sampled path (indicated in orange) by Nshift time-steps (where Nshift = ∆t/dt) to the position x(∆t). With the remaining segment of the path now located at positions {x(∆t), x(∆t + dt), . . . , x(T )}, we grow the path out of x(T ) by Nshift time-steps to regenerate the original number of time-steps in the path. The positions of the system at the regenerated time-steps can be drawn from any distribution (and in Fig. 4.2B they are obtained via straight-line extrapolation). Finally, as illustrated in Fig. 4.2C, sampling is again performed to generate a path consistent with evolution from time ∆t to time T + ∆t; this sampling removes any artifacts introduced by the arbitrary distribution used to grow the shifted path. As for the scheme in Fig. 4.1, iterating through the sampling, shifting, and regeneration steps of the scheme in Fig. 4.2 yields a MD trajectory of arbitrary length in time that is statistically equivalent to a realization from 86 A sampled paths x(0) x(T ) B resolved trajectory sampled paths position x(0) 0 dt x(T ) x(2T ) T time Figure 4.1. Illustration of a simple path-based MD integration scheme. A Sampling of the distribution of paths for the period of time from 0 to T . A particular path drawn from this distribution is indicated in orange. B Shifting along a sampled path (orange path in A) from x(0) to x(T ), thereby resolving the segment of MD trajectory indicated in green, and then resuming the sampling of the paths for the period of time from T to 2T . the Euler algorithm. The only difference between these two path-based integration schemes is that Fig. 4.1 involves shifting along the full length of the sampled path, whereas Fig. 4.2 involves shifting only a fraction of the way along the sampled path. Just like the Euler scheme in Eq. (4.2), the schemes illustrated in Figs. 4.1 and 4.2 enable the numerical integration of MD trajectories. Each of these integration schemes consist of sequential iterations of an elementary step that predicts the state of the system at some later time. In the Euler scheme, the prediction can be conducted analytically based on the distribution defined in Eq. (4.3). In path-based integration schemes employing path lengths longer that dt, however, no such analytical expression exists for general systems; path sampling is therefore needed before each shifting event to generate time-evolved system positions consistent with the correct marginal distribution. For the scheme in Fig. 4.2, it is assumed that the path distributions in parts A and C are well sampled. For MC path sampling algorithms, this implies that the number of configurations of the path that are sampled in parts A and C, Nsample , is large 87 sampled paths A x(0) x(T ) resolved trajectory shifted path regenerated path B x(0) x(T ) x( t) resolved trajectory x(T + t) x(T + t) sampled paths C position x(0) x( t) 0 dt time t T (Nshift · dt) (Npath · dt) T+ t Figure 4.2. Illustration of a path-accelerated molecular dynamics (PAMD). A Sampling of the distribution of paths for the period of time from 0 to T . A particular path drawn from this distribution is indicated in orange. B Shifting along a sampled path (orange path in A) from x(0) to x(∆t), thereby resolving the segment of MD trajectory indicated in green, and regenerating the full length of the path by drawing positions for the system from time T + dt to T + ∆t from an arbitrary distribution. C Sampling of the distribution of paths for the period of time from ∆t to ∆t + T . in comparison to the number that is needed to generate uncorrelated realizations of the path. If Nsample is smaller than this decorrelation number, then the distribution of paths that is generated in part C may be biased by the way in which the path was regenerated in part B. However, the only requirement for generating accurate MD trajectories using the scheme in Fig. 4.2 is accurate sampling of paths 88 consistent with the marginal distribution K(x(∆t)|x(0); ∆t); it is not essential that the marginal distribution associated with the full path, K(x(T )|x(0); T ), be sampled without error. Recalling that T = Npath dt and ∆t = Nshift dt, this suggests that for a given path sampling algorithm, there is an interplay between parameters Nsample , Npath and Nshift ; for given values of Npath and Nshift , there is an associated number of path configurations (Nsample ) that must be sampled in order to generate a sufficiently accurate marginal distribution K(x(∆t)|x(0); ∆t). This interplay between Nsample , Npath , and Nshift is illustrated in Figs. 4.3A-C, which plot the error in the marginal distribution generated using the scheme in Fig. 4.2 for the Brownian dynamics of a harmonic oscillator. Full calculation details are provided in Sec. 4.3. The error plotted in Figs. 4.3A-C corresponds to the Kullback– Leibler divergence, w  P (xt |x0 ; t) DKL (t) = dxt P (xt |x0 ; t) log , (4.6) R Q(xt |x0 ; t) where Q(xt |x0 ; t) is the marginal distribution estimated using sampled paths from the scheme in Fig. 4.2, and P (xt |x0 ; t) is the exact marginal distribution. The angled brackets denote averaging with respect to the Boltzmann distribution of positions that is sampled by the exact dynamics, P (x0 ) = Z −1 e−βV (x0 ) , where Z = r dx0 e−βV (x0 ) is the partition function. For a harmonic oscillator with potential R V (x) = 12 kx2 [49], ( ) −1 βk (xt − e−γ kt x0 )2 P (xt |x0 ; t) ∝ exp − , (4.7) 2 (1 − e−2γ −1 kt ) and we employ k = 1 for the oscillator force constant, β −1 = 1 for the temperature and γ = 1 for the friction coefficient. DKL (t) returns non-negative values that approach 0 as Q(xt |x0 ; t) more accurately reproduces P (xt |x0 ; t). As a function of time t along the sampled paths, DKL (t) is plotted in Figs. 4.3A-C for seven simulations that employ the scheme in Fig. 4.2 with different values of Nsample , Npath , and Nshift . The results correspond to sampled paths of length T ≥ 1 that are discretized into time-steps of dt = 1/32, and DKL (t) is evaluated for the numerically generated marginal distributions at times dt ≤ t ≤ 1. Figure 4.3A addresses the case where Nsample and Npath are held fixed and various values of Nshift are used. Comparison of the blue (Nshift = 4), green (Nshift = 8), and red (Nshift = 16) curves shows that for a given value of Nsample , smaller values of Nshift lead to smaller errors in the numerically generated marginal distribution. Using the scheme in Fig. 4.2, a given segment of the path is sampled 89 A B C D E F Figure 4.3. The interplay of the parameters Nsample , Npath , and Nshift in determining the accuracy of the PAMD integration scheme. Panels A, B and C show the Kullback–Leibler divergence DKL (t) defined in Eq. (4.6), from the marginal distribution generated using the scheme in Fig. 4.2 to the exact marginal distribution for the dynamics of an overdamped harmonic oscillator. Panels D, E, and F plot the position distributions P (x) sampled by the MD trajectories integrated using the numerically generated marginal distributions from panels A, B, and C, respectively, in comparison to the exact position distribution shown in black dots. It is seen that the accuracy of the numerically generated marginal distribution dictates that of the integrated MD trajectory, and improves with decreasing Nshift (panels A and D), increasing Nsample (panels B and E), or increasing Npath (panels C and F) for given values of the remaining parameters (indicated at the top of each panel). Npath · Nsample /Nshift times before it is used to generate the marginal distribution for the integration of the MD trajectory; therefore, smaller values of Nshift lead to better sampling of the path distribution and smaller errors in the marginal distribution. Figure 4.3B illustrates a second scenario where Nshift and Npath are held fixed and increasing values of Nsample are used. Comparison of the blue (Nsample = 1), green (Nsample = 2), and red (Nsample = 4) curves shows that for a given value of Nshift , larger values of Nsample (i.e., more sampling per shifting event) lead to smaller errors in the numerically generated marginal distribution for integrating the MD trajectory. This result is intuitive, as more sampling leads to elimination of the bias associated with the arbitrary distribution used in the regeneration of the full length of the path. 90 In Fig. 4.3C, Nshift and Nsample are held fixed as the length of the sampled path (Npath ) is increased while keeping the path discretization time-step unchanged. Comparison of the blue (Npath = 32), green (Npath = 128), and red (Npath = 512) curves demonstrates that increasing the total length Npath ·dt of the sampled paths improves the accuracy of the numerically generated marginal distribution. Like decreasing Nshift for a given Npath (as in Fig. 4.3A), increasing Npath for a given Nshift allows for more sampling of each segment of the path employed to generate the marginal distribution associated with the MD time evolution. While Figs. 4.3A-C illustrate the errors in the marginal distribution generated using the scheme in Fig. 4.2, Figs. 4.3D-F illustrate the corresponding errors in the equilibrium distribution that is sampled by the integrated MD trajectories. For the various employed parameters, the results from path-based MD integration are compared to the exact Boltzmann distribution (dots) and, as expected, the errors in the marginal distribution with given values of Nsample , Npath , and Nshift are reflected in the distribution of positions that are visited in the MD trajectories. Supplementary movies M1 and M2 (see Supplementary Material at the end of this chapter), respectively, illustrate the integration scheme in Fig. 4.2 as it generates the position distributions for the cases Nsample = 1 (blue curve) and Nsample = 4 (red curve) in Fig. 4.3E. In summary, Fig. 4.3 demonstrates that decreasing Nshift , increasing Nsample , or increasing Npath leads to greater accuracy in the integrated MD trajectories; as will be shown in Sec. 4.2.2, the interplay between these three parameters is also critical for determining the computational efficiency of MD integration using the scheme in Fig. 4.2. Before addressing efficiency, however, Fig. 4.4 illustrates that the integration scheme in Fig. 4.2 is a non-equilibrium relaxation process for the segments of the sampled path. For the case of the harmonic oscillator, Figs. 4.3A-C indicate that errors in the numerically generated marginal distributions are typically larger at the right endpoint of the sampled path. This trend emerges because the integration scheme regenerates path segments in configurations that are out of equilibrium (Fig. 4.2B). The light gray paths in Fig. 4.4 correspond to independent realizations of the sampled path (orange) obtained while integrating a harmonic oscillator trajectory. Since the segment at the nose of the path has undergone little sampling after regeneration, it is far from equilibrium with respect to the distribution of segments of an equilibrium harmonic oscillator trajectory (Fig. 4.4A; distributions at right). However, as 91 trajectory sampled paths position x(t) A tail nose P (x(3 t)) trajectory x(3 t) x(0) path time t 0 trajectory 4 t 5 t position distribution P (x(t)) sampled paths position x(t) B 3 t 2 t t nose tail P (x(3 t)) trajectory x(3 t) x( t) path time t t trajectory 4 t 5 t position distribution P (x(t)) sampled paths position x(t) C 3 t 2 t tail nose P (x(3 t)) trajectory x(3 t) x(2 t) path time t 2 t 3 t 4 t 5 t position distribution P (x(t)) Figure 4.4. PAMD integrates equilibrium trajectories by relaxing non-equilibrium path segments. The state of the sampled path is shown for three consecutive iterations of the pathbased integration scheme applied to a Brownian harmonic oscillator on the left of panels A, B, and C. The gray box in each panel highlights the configuration of a particular segment of the path after each iteration of the integration scheme. On the right of each panel are plotted the distributions P (x(3∆t)) of positions x(3∆t) sampled by the boxed path segment at time 3∆t (orange curves), along with the distribution of configurations sampled by the integrated trajectory (green curve). The path segment inside the gray box in each panel reaches equilibrium by undergoing sampling as it shifts from the nose (right endpoint) to the tail (left endpoint) of the sampled path; accordingly, the distribution of sampled positions approaches that sampled by the integrated trajectory. that segment works its way from the nose (right endpoint) to the tail (left endpoint) of the path, it is sampled with increasing accuracy (Figs. 4.4B and C). This relaxation process is illustrated by the distribution of positions, P (x(3∆t)), sampled by 92 the foremost end x(3∆t) of the path segment in the gray box throughout Fig. 4.4; as the orange curves indicate, this distribution approaches that sampled by the harmonic oscillator trajectory (green curve) as the segment relaxes toward equilibrium and simultaneously approaches the tail of the sampled path. 4.2.2 An opportunity for speedup At face value, the path-based integration scheme in Fig. 4.2 may appear to be inefficient, given the difficulties of sampling uncorrelated paths [38, 46, 50, 51]. Yet, it has several potential advantages: Firstly, there is an opportunity for parallelization, given that typical expressions for the path action incur a dominant source of computational cost from the evaluation of the forces in the system along the path (V ′ (x) in Eq. (4.3)). These forces can be evaluated independently, enabling straightforward parallelization of the action with respect to time. Secondly, regeneration of the path to its full length following shifting (Fig. 4.2B) can be performed using an arbitrary distribution to obtain the system positions for the regenerated time-steps; consequently, it is possible to carry out this operation at a cost that is negligible relative to evaluation of the MD forces. Thirdly, MC path sampling provides a numerically more stable way for generating trajectories than integration of the discretized equations of motion [34]; thus, a sufficiently accurate MD trajectory may be obtained with the path-based integration scheme at a larger time-step than a conventional Brownian dynamics integrator would allow. The above considerations suggest that the scheme in Fig. 4.2 could lead to reduction of the wall-clock time associated with MD integration, in comparison with standard methods. To quantify the speedup achieved with the new scheme, we introduce a measure χ, with χ−1 defined as the number of force evaluations per processor per step of time dtE , where dtE is the time-step used by the Euler algorithm to integrate the Brownian dynamics. The wall-clock speedup of the path-based integration scheme is thus χ, assuming that (i) evaluation of the MD forces dominates the cost of the evaluation of the path action, (ii) parallel computer processors are used to independently evaluate the forces along the discretized path, and (iii) regeneration of the full length of the path following shifting (Fig. 4.2B) is performed without evaluating the MD forces. It is clear that for the Euler algorithm, χ = 1, such that this measure provides a simple basis of comparison of the wall-clock time for the proposed path-based integration scheme (which employs parallelization in time) versus the wall-clock time for a conventional stochastic MD integration scheme (which does 93 not). In the current chapter, we set aside the complementary issue of speeding up MD integration via parallelization within the force evaluation. For a general implementation of the integration scheme in Fig. 4.2, the expression for χ is obtained as follows. Recalling previously introduced notation, we employ sampled paths of length T that are discretized with a time-step of dt, which may be different (and is typically larger [34]) than the numerically stable time-step for the Euler algorithm, dtE . Let Nforce be the number of MD force evaluations that are required during path sampling per shifting event, which depends on both Nsample and the details of the path sampling algorithm, and let Nprocs be the number of employed parallel processors. Since the number of force evaluations per processor per shifting event is given by Nforce /Nprocs , the speedup is χ = Nshift · Nprocs dt · . Nforce dtE (4.8) Equation (4.8) shows that the path-based integration scheme in Fig. 4.2 offers the possibility for reduction of the wall-clock time needed to compute MD trajectories, relative to conventional stochastic MD. Factors that enable this speedup include the increase in the path discretization time-step (dt) relative to that possible for conventional stochastic MD (dtE ), maximization of the number of integrated time-steps per shifting event (Nshift ), maximization of the number of parallel processors to perform the independent force evaluations associated with the calculation of the path action (Nprocs ), and minimization of the number of force evaluations needed per shifting event (Nforce ). As will be shown in Sec. 4.4, this approach indeed enables substantial speedups in the integration of MD trajectories while preserving the accuracy of the dynamics, and we henceforth refer to the method as path-accelerated molecular dynamics (PAMD). 4.3 Calculation details We implement the PAMD method with sampling of the path distribution through the multilevel sliding and sampling (S&S) algorithm introduced in Ref. [38]. For a path of Npath time-steps, a total of L = log2 Npath levels are defined (Fig. 4.5A); finer levels (smaller values of the level index 1 ≤ l ≤ L) correspond to partitions of the path into fragments of increasingly smaller length where the local configuration of the path is sampled. In accordance with the S&S algorithm, neighboring path fragments share endpoints that are chosen randomly such that the length of the fragments varies from 1 to 2l time-steps; we call this random fragmentation. For 94 all fragmentations of the path at level l, internal fragments of the path are of length 2l and fragments at the termini of the path have a combined length of 2l . During a MC step for a given fragmentation of the path (Fig. 4.5B), the system positions at shared endpoints of neighboring path fragments are held fixed to permit mutually independent updates of the fragment configurations. Furthermore, the position of the system at the tail endpoint of the path is always fixed throughout the MC step, whereas that at the nose endpoint of the path undergoes sampling together with the nose fragment. Path fragment configurations are updated through application of the standard Metropolis–Hastings criterion [52, 53], with trial configurations drawn from a distribution that satisfies the boundary conditions at the fragment endpoints. Random fragmentation of the path is performed between MC steps, so that fixed system positions at previous fragment endpoints can be sampled during subsequent steps (Fig. 4.5C). Table 4.1. Summary of notation employed to describe the PAMD integration scheme. The last three rows define parameters specific to the path sampling algorithm used. dt Npath Nprocs Nforce Nsample Nshift L lmin lmax time-step for discretization of the sampled path number of time-steps in the sampled path number of processors for parallel-in-time force evaluations number of force evaluations per shifting event number of MC steps per shifting event number of time-steps shifted total number of levels in the sampled path finest sampled level of the path coarsest sampled level of the path At each MC step, a level is randomly selected between lmin and lmax , with 1 ≤ lmin ≤ lmax ≤ L. The calculations reported here employ lmin > 1 and lmax < L, such that not all levels are directly sampled. The choice of lmax < L corresponds to excluding the direct sampling of levels associated with long path fragments, on the basis of negligible acceptance. The choice of lmin > 1 corresponds to excluding the direct sampling of levels associated with short path fragments, as these are trivially updated via the direct sampling of longer fragments at coarser levels. Despite these choices, the sampling remains ergodic due to the random fragmentation of the path that occurs between MC steps [38, 54]. 95 Our choice of the multilevel S&S algorithm over other sampling algorithms offering gradient-based global path updates, based on an effective Hamiltonian or Langevin dynamics on the path space [47, 55–58], is motivated by at least three practical considerations specific to the context where PAMD applications are envisioned. The first pertains to the computational cost of each MC sampling step, which varies drastically across algorithms. While gradient-based algorithms, at least in their position A initial path time l l== 11 l l== 22 l l== 33 l l== 44 B position first MC step time l l== 22 C position second MC step time l=3 Figure 4.5. Illustration of the multilevel sliding and sampling (S&S) algorithm. A Multilevel representation of a path with Npath = 16 time-steps (L = 4 levels), shown for one Cartesian coordinate. An initial configuration of the path is shown in orange. B Update of the path configuration in A via a MC step at level l = 2. Fixed system positions along the path are represented with black dots, configurations of the path before the update in gray, and the configuration of the path after the update in orange. C MC step at level l = 3 following that shown in B. The new fragmentation of the path allows for updates of system positions that were held fixed in previous MC steps. 96 Metropolis-unadjusted variants, can provably exhibit faster mixing than gradientfree counterparts and exhibit robustness to arbitrary refinement of the MD timestep toward the continuous-time limit [55–58], they require the evaluation of higher derivatives of the physical potential that are often unavailable or considered too expensive in applications. On the other hand, gradient-free path sampling algorithms can be combined with a symmetrized choice of the discretized path-action functional [38, 45] to fully bypass all need for higher-order potential derivatives. The second consideration is the requirement that the sampled path ensemble not exhibit spurious artifacts associated with the sampling method employed or the chosen discretization of the Onsager–Machlup functional. Observations have surfaced that gradient-based path sampling algorithms can produced unphysical path ensembles when combined with certain discretizations of the Onsager–Machlup functional [41, 59]; such artifacts are also avoided with gradient-free updates such as the celebrated Crank–Nicolson [60] (or Brownian tube [48]) method. The third and final consideration is to avoid the slow MC mixing associated with global updates in long paths. The above-mentioned gradient-based and gradientfree methods are all formulated to induce updates of the full configuration of the sampled path in each MC step. The acceptance rate of such global updates can scale poorly with sampled path length at fixed time-step size [51, 61], and thus lead to slow equilibration within the long sampled-path lengths envisioned to be necessary for effective application of the PAMD algorithm. This scaling consideration also jeopardizes the performance of semi-local spectral Metropolis-within-Gibbs schemes that sequentially update the components of a normal-mode representation of the sampled path across its full length [62, 63]. The distribution of paths used to generate trial configurations for multilevel S&S in each application is chosen to maximize the statistical efficiency of the sampling (i.e., minimize Nsample ) without requiring evaluation of the MD forces. For the harmonic oscillator, trials are drawn from the distribution of free particle paths. For simulations of the Lennard–Jones liquid, trials are drawn from the path distribution of a fluid of hard spheres with diameter σHS ; this strategy reduces the number of force evaluations needed to obtain likely Lennard–Jones path configurations by excluding those with high interparticle overlap from the ensemble of trial paths. The likelihood of a path at the hard-sphere level is evaluated using an approximation of the pair propagator for diffusive hard spheres [64]. For both the harmonic oscillator and the Lennard–Jones applications, regeneration of the sampled path after shifting 97 is performed with the same distribution used to generate trial configurations for the path sampling. In total, each MC step involves a total number of Npath MD force evaluations. Since Nforce is defined as the number of force evaluations per shifting event, and since Nsample is the number of MC steps per shifting event, we have Nforce = Nsample · Npath . (4.9) As we seek to maximize the wall-clock speedup via parallelization of these independent force evaluations, we employ one processor per force evaluation, and thus Nprocs = Npath . (4.10) Additional parallelization within the force evaluation is of course possible, but is not considered in the current chapter. Thus, we insert Eqs. (4.9) and (4.10) into Eq. (4.8) to arrive at the following expression for the PAMD speedup: χ= Nshift dt · . Nsample dtE (4.11) See Tbl. 4.1 for a summary of terms. 4.4 Results We now apply the PAMD algorithm to two model systems, considering time evolution according to Brownian dynamics (Eq. (4.1)), with β = 1 unless otherwise indicated and with γ = 1 in appropriately reduced units. Like the Euler algorithm (Eq. (4.2)), PAMD is a rigorous and formally exact way to integrate Brownian dynamics, yet the numerical accuracy of the trajectories depends on the parameters employed. In each application, we examine the relationship between the number of parallel processors employed and the speedup in the PAMD algorithm relative to the Euler algorithm (χ in Eq. (4.11)), subject to the requirement that the MD trajectories integrated using both PAMD and the Euler algorithm preserve well-defined measures of accuracy. In the current chapter, we focus exclusively on wall-clock speedups achieved via parallelization of the MD integration in time, setting aside the separate and complementary issue of parallelizing the force evaluation at each time-step. All reported speedups for PAMD in the current study are theoretical; they are obtained from 98 Eq. (4.11) under the stated assumptions. The computational cost associated with each speedup can be obtained by multiplying the reported value with the number of parallel processors, Nprocs , employed in each application. 4.4.1 Harmonic oscillator Here, we consider the example of an overdamped harmonic oscillator, with potential V (x) = 12 x2 . Two measures of the accuracy of the integrated MD trajectories are considered. The first reports on the degree to which the trajectories sample the correct equilibrium distribution, qr dx |P (x) − Ps (x)|2 R Eeq = (4.12) Z where P (x) is the exact Boltzmann distribution, Z is the associated partition function, and Ps (x) is the equilibrium distribution of positions sampled by the numerical integration schemes. The second measure of error reports on the accuracy of the MD time evolution. Specifically, we consider the autocovariance function 1 wT ′ ′ dt x(t )x(t′ + t) , (4.13) C(t) = ⟨x(t′ )x(t′ + t)⟩ = lim T →∞ T 0 which is a simple exponential function for the overdamped harmonic oscillator [49], C(t) = x2 exp(−κt) , (4.14) where the angled brackets indicate Boltzmann averaging. The second measure of error is thus |κ − κs | Edyn = (4.15) k where κ = γ −1 = 1 is the exact decay constant and κs is the decay constant obtained by fitting the exponential decay of the autocovariance from the numerically integrated MD trajectories. Specifically, κs is obtained by averaging over 100 independent trajectories of length 105 time units that are divided into 103 non-overlapping time series, for which the log-autocovariance is linearly fit in the range t ∈ [0, 4]. Simulation parameters for integration of the MD trajectories in this application are chosen to ensure that both measures of error remain below 3%. Table 4.2 indicates that with a time-step of dtE = 0.025, the Euler algorithm yields error values of Eeq = 0.3% and Edyn = 1.0%. Also shown in the table are parameters for three separate PAMD simulations that obtain speedups of χ = 16, 128, and 99 1024 with respect to the Euler algorithm. For both the equilibrium distribution and the autocovariance function, Fig. 4.6 shows the comparison of the exact results and those obtained using PAMD with the aforementioned speedups. The accuracy of the PAMD trajectories is clearly preserved in all simulations, as indicated by the plotted results and the reported values of Eeq and Edyn . A B Figure 4.6. For the harmonic oscillator, comparison of PAMD results (colored lines) with exact results (black dots) for A the Boltzmann distribution P (x) and B the autocovariance function C(t). The PAMD results correspond to Simulation 1 (red), Simulation 2 (green), and Simulation 3 (blue) in Tbl. 4.2, which respectively achieve speedups of χ = 16, 128, and 1024 relative to the Euler algorithm. For all PAMD simulations reported in Tbl. 4.2, a significant fraction of the speedup comes from the 16-fold larger time-step that can be employed in the path-based scheme (dt = 0.4 vs. dtE = 0.025). The larger speedups achieved in Simulation 2 (χ = 128) and Simulation 3 (χ = 1024), in comparison to Simulation 1 (χ = 16), arise from the larger ratios of Nshift to Nsample that are used in these simulations (Nshift /Nsample = 8 for Simulation 2 and Nshift /Nsample = 64 for Simulation 3) in comparison to Simulation 1 (for which Nshift /Nsample = 1). The higher frequency of shifting events associated with larger values of the ratio Nshift /Nsample places greater demand on the efficiency of the path sampling, and a larger number of time-steps in the sampled path (Npath ; hence, a larger number of parallel processors, Nprocs , per Eq. (4.10)) is needed to allow path segments to undergo a sufficient number of MC steps before they are used to generate the marginal distribution for the MD trajectories, as discussed in connection with Figs. 4.3 and 4.4. Accordingly, the value of Npath in each of the simulations in Tbl. 4.2 is set as large as needed to ensure accurate MD integration at the corresponding value of Nsample , which is kept small (Nsample = 2 in Simulation 1, and Nsample = 1 in Simulation 2 and Simulation 3) to enhance the PAMD speedup. 100 Table 4.2. Summary of PAMD simulation parameters used for the harmonic oscillator application. Parameters L, lmin , and lmax are specific to the path sampling algorithm used, described in Sec. 4.3. dtE Eeq Edyn dt Npath Nsample Nshift L lmin lmax χ Eeq Edyn Simulation 1 0.4 16 2 2 4 3 3 16 0.3% 0.3% Euler 0.025 0.3% 1.0% PAMD Simulation 2 0.4 256 1 8 8 4 5 128 0.3% 1.1% Simulation 3 0.4 4096 1 64 12 4 8 1024 0.5% 2.9% 4.4.2 Lennard–Jones liquid Here, we apply PAMD to a model for a molecular liquid. The pairwise interaction between particles is described using the standard energy- and force-shifted Lennard– Jones potential [2], ( u(r) − u(rc ) − (r − rc ) u′ (rc ), r ≤ rc U (r) = , (4.16) 0, r > rc where u(r) = 4ϵ {(σ/r)12 − (σ/r)6 }; throughout, we take ϵ = 1 and σ = 1. The system consists of 27 particles placed in a cubic box at reduced density ρσ 3 = 0.50 and at constant reduced inverse temperature βϵ = 0.74. Simulations are performed with periodic boundary conditions at constant volume, and the cutoff distance rc corresponds to half of the simulation box-length. As described in Sec. 4.3, the reported simulations for the Lennard–Jones liquid employ trial configurations drawn from a distribution of paths for a fluid of hard spheres with diameter σHS . To prevent the path-sampling bias from affecting the accuracy of the integrated trajectories, we employ a hard-sphere schedule that varies as a function of the path-time τ , σHS (τ ); path configurations are then sampled in accordance 101 with the path-time dependent potential (4.17) V (r; σHS (τ )) = U (r) + UHS (r; σHS (τ )) , where U (r) is defined in Eq. (4.16), and UHS (r; σHS ) is the hard-sphere potential ( +∞, r ≤ σHS UHS (r; σHS ) = . (4.18) 0, r > σHS The schedule σHS (τ ) is chosen such that configurational volume is excluded at the nose of the path for enhanced sampling efficiency, and no volume is excluded at the tail of the path where the marginal distribution for MD integration is sampled (Fig. 4.7); in this way, path segments regenerated from the hard-sphere distribution are subsequently relaxed into the Lennard–Jones distribution as they shift from the nose to the tail of the path. A worthy direction for future work is to systematically optimize the schedule σHS (τ ), which in general will depend on the simulation parameters reported in Tbl. 4.1. path sampling path regeneration nose tail scaled path-time Figure 4.7. Hard-sphere schedule σHS (τ ), versus scaled path-time τ , employed in the reported Lennard–Jones simulations. Insets show slices of the path-time dependent potential V (r; σHS (τ )) employed to sample path configurations (Eq. (4.17)) at schedule values corresponding to biased (σHS ≃ 1) and unbiased (σHS = 0) sampling of the equilibrium path distribution, and the hard-sphere potential U (r; σHS ) used to regenerate path segments (Eq. (4.18)). The accuracy of the integrated MD trajectories is evaluated in terms of the radial distribution function g(r) and the self-diffusion coefficient D, using the respective error measures qr rc dr |g(r) − gs (r)|2 0 r rc Eeq = (4.19) dr g(r) 0 102 and |D − Ds | , (4.20) D where g(r) and D are reference quantities obtained using the Euler algorithm with a small time-step (5 × 10−5 Lennard–Jones time units), and gs (r) and Ds are obtained using PAMD and the Euler algorithm with larger time-steps. The diffusion d ⟨R2 (t)⟩, where coefficient is given by D = 16 limt→∞ dt Edyn = 2 R (t) = lim T →∞ wT N 1 X dt |ri (t′ + t) − ri (t′ )|2 0 N i=1 ′ (4.21) is the mean-square displacement, ri the position of the ith particle, and N the number of particles [1]. ⟨R2 (t)⟩ is obtained by averaging over 100 independent trajectories that are divided into 100 non-overlapping time series of length 1, in Lennard– Jones time units, and a linear fit is performed in the range t ∈ [0.2, 1] to evaluate Ds for the PAMD and Euler simulations. Table 4.3 indicates that at a time-step of dtE = 2.5×10−4 , the Euler algorithm yields error values of Eeq = 1.4% and Edyn = 0.1%; larger time-steps were found to lead to unstable Euler trajectories. Also shown in Tbl. 4.3 are two PAMD simulations that lead to 16-fold (χ = 16; Simulation 1) and 128-fold (χ = 128; Simulation 2) reductions of the wall-clock time required to generate equivalently accurate MD trajectories for the Lennard–Jones liquid via the Euler algorithm, using simulation parameters that keep error values below 5%. The radial distribution functions and mean-square displacements obtained from these two simulations are plotted with the corresponding reference quantities in Figs. 4.8A and B. Excellent agreement between the PAMD and reference quantities is evident in the plots and from the values of Eeq and Edyn reported in Tbl. 4.3. As in the harmonic oscillator application, the speedups reported in Tbl. 4.3 for the Lennard–Jones liquid are partially enabled by the use of a larger time-step in PAMD (dt = 5×10−4 ) than is possible for stable numerical integration via the Euler method (dtE = 2.5 × 10−4 ). The remaining speedup in both simulations comes from using shift lengths that integrate Nshift = 8 (Simulation 1) and Nshift = 64 (Simulation 2) time-steps of MD trajectory at a rate of Nsample = 1 MC steps per shifting event. Accurate integration at these speedups requires efficient sampling of path modes that are commensurate with the shifting timescale (8 dt in Simulation 1 and 64 dt in Simulation 2); accordingly, long paths (Npath = 128 in Simulation 1 and Npath = 4096 in Simulation 2) are employed in both simulations so that segments can undergo a 103 A B Figure 4.8. For the Lennard–Jones liquid at βϵ = 0.74 and ρσ 3 = 0.50, comparison of PAMD results (colored lines) with numerically exact results (black dots) for A the radial distribution function and B the mean-square displacement. The PAMD results correspond to Simulation 1 (red) and Simulation 2 (green) in Tbl. 4.3, which respectively achieve speedups of χ = 16 and 128 relative to the Euler algorithm. sufficient number of MC steps before they are used to generate the respective MD trajectories (Fig. 4.4). 4.5 Conclusions The field of MD simulation faces important challenges in harnessing massively parallel computer architectures. Although successful parallelization of the the force evaluation can be expected as the system size grows (i.e., weak scaling), there exists a much more difficult challenge of employing ever-larger numbers of parallel processors to accelerate the simulation of systems of a fixed size (i.e., strong scaling). Remarkable success has been achieved in this vein [8–13], but fundamental limitations are inevitable. The work comprising this chapter suggests that parallelization in the dimension of time via path integrals offers a promising avenue for future progress in terms of strong scaling. We introduce the PAMD approach, which achieves significant speedups over conventional algorithms for stochastic dynamics by parallelizing evaluation of the MD forces along trajectory segments with respect to time. The method strives to be applicable in greater generality than alternative approaches to longtimescale simulation of stochastic trajectories, such as milestoning [65], Markov state modeling [66], and trajectory splicing [24], by enabling the direct simulation of long-timescale processes without assuming separation of timescales, and without re- 104 Table 4.3. Summary of PAMD simulation parameters used for the application to the Lennard–Jones liquid. L, lmin , and lmax are parameters specific to the path sampling algorithm used, described in Sec. 4.3. dtE Eeq Edyn dt Npath Nsample Nshift L lmin lmax χ Eeq Edyn Euler 2.5 × 10−4 1.4% 0.1% PAMD Simulation 1 Simulation 2 −4 5 × 10 5 × 10−4 128 4096 1 1 8 64 7 12 3 4 5 7 16 128 3.5% 3.7% 1.0% 1.9% quiring the specification of a suitable configuration-space tiling or a set of collective variables. Proof-of-principle applications to overdamped systems show that PAMD can accelerate trajectory integration by several orders of magnitude with respect to the conventional Euler scheme for Brownian dynamics, and even greater speedups are possible with the use of larger numbers of parallel processors and enhancement of the MC path sampling efficiency. Looking forward, the PAMD approach will likely require additional methodological developments to become viable for large-scale simulations of complex systems. Central to this effort will be the refinement of path sampling methodologies that lead to the reduction in the number of parallel processors that are needed for a given amount of speedup with the method. Our application to the Lennard–Jones liquid shows that preconditioned sampling of paths from an approximate dynamical model [67], constructed a priori based on known features of the dynamics or generated on-the-fly via adaptive inference techniques for stochastic process approximation [68–70], provides an avenue for efficiency improvements in more general applications. Regardless, we feel that the natural parallelization of path-integral formulations, combined with the increasing availability of massively parallel computer 105 resources, should motivate increased attention to the opportunities of parallelizing molecular dynamics simulation in time. Supplementary Material Supplementary movies M1 and M2, respectively, show execution of the PAMD integration scheme as it accumulates position histograms during simulations of the Brownian harmonic oscillator, for the cases Nsample = 1 (blue curve) and Nsample = 4 (red curve) in Fig. 4.3E. These cases illustrate distinct operational regimes of the PAMD algorithm, respectively corresponding to inaccurate and accurate stochastic MD trajectories as a function number of path sampling steps per PAMD iteration. The movies are available from CaltechDATA at DOI: 10.22002/D1.20111. Acknowledgments. We are grateful for stimulating discussions with Matthew G. Welborn, Eric Vanden-Eijnden, and Gavin E. Crooks. The work comprising this chapter was supported in part by the Department of Energy under Award No. DEFOA-0001912 and the Office of Naval Research under Award No. N00014-10-10884. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the DOE Office of Science under contract DE-AC02-05CH11231. 106 References Cited in Part II [1] D. Frenkel and B. Smit, Understanding molecular simulation: From algorithms to applications (Academic Press, 2002). [2] M. P. Allen and D. J. Tildesley, Computer simulation of liquids (Clarendon Press, 2017). [3] J. D. Durrant and J. A. McCammon, “Molecular dynamics simulations and drug discovery”, BMC Biology 9, 71 (2011). [4] R. O. Dror, R. M. Dirks, J. P. Grossman, H. Xu, and D. E. Shaw, “Biomolecular simulation: A computational microscope for molecular biology”, Annual Review of Biophysics 41, 429 (2012). [5] A. F. Voter, F. Montalenti, and T. C. Germann, “Extending the time scale in atomistic simulation of materials”, Annual Review of Materials Research 32, 321 (2002). [6] D. Keyes, in High performance computing - HiPC 2007, edited by S. Aluru, M. Parashar, R. Badrinath, and V. K. Prasanna (Springer Berlin Heidelberg, 2007), pp. 2–3. [7] J. Dongarra, P. Beckman, T. Moore, P. Aerts, G. Aloisio, J.-C. Andre, D. Barkai, J.-Y. Berthou, T. Boku, B. B, F. Cappello, B. Chapman, X. Chi, A. Choudhary, S. Dosanjh, T. Dunning, S. Fiore, A. Geist, B. Gropp, B. Harrison, M. Hereld, M. Heroux, A. Hoisie, K. Hotta, Z. Jin, Y. Ishikawa, F. Johnson, S. Kale, R. Kenway, D. Keyes, B. Kramer, J. Labarta, A. Lichnewsky, T. Lippert, B. Lucas, B. Maccabe, S. Matsuoka, P. Messina, P. Michielse, M. B., M. S. Mueller, W. E. Nagel, H. Nakashima, M. E. Papka, D. Reed, M. Sato, E. Seidel, J. Shalf, D. Skinner, M. Snir, T. Sterling, R. Stevens, F. Streitz, B. Sugar, S. Sumimoto, W. Tang, J. Taylor, R. Thakur, A. Trefethen, M. Valero, A. van der Steen, J. Vetter, P. Williams, R. Wisniewski, and K. Yelick, “The International Exascale Software Project roadmap”, The International Journal of High Performance Computing Applications 25, 3 (2011). [8] S. Plimpton, “Fast parallel algorithms for short-range molecular-dynamics”, Journal of Computational Physics 117, 1 (1995). [9] J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, and K. Schulten, “Scalable molecular dynamics with NAMD”, Journal of Computational Chemistry 26, 1781 (2005). [10] D. E. Shaw, R. O. Dror, J. K. Salmon, J. P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. Deneroff, B. Batson, K. J. Bowers, E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson, K. LindorffLarsen, P. Maragakis, M. A. Moraes, S. Piana, Y. Shan, and B. Towles, “Millisecond-scale molecular dynamics simulations on Anton”, in Proceedings of 107 the conference on high performance computing networking, storage and analysis, SC ’09 (2009), pp. 1–11. [11] R. Salomon-Ferrer, A. W. Götz, D. Poole, S. Le Grand, and R. C. Walker, “Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald”, Journal of Chemical Theory and Computation 9, 3878 (2013). [12] S. Páll, M. J. Abraham, C. Kutzner, B. Hess, and E. Lindahl, in Solving software challenges for exascale, edited by S. Markidis and E. Laure (Springer International Publishing, 2015), pp. 3–27. [13] J. P. Grossman, B. Towles, B. Greskamp, and D. E. Shaw, “Filtering, reductions and synchronization in the Anton 2 network”, Parallel and Distributed Processing Symposium, 860 (2015). [14] J.-L. Lions, Y. Maday, and G. Turinici, “A parareal discretization in time of PDE’s”, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 332, 661 (2001). [15] C. Farhat and M. Chandesris, “Time-decomposed parallel time-integrators: Theory and feasibility studies for fluid, structure, and fluid-structure applications”, International Journal for Numerical Methods in Engineering 58, 1397 (2003). [16] I. D. Garrido, B. Lee, G. E. Fladmark, and M. S. Espedal, “Convergent iterative schemes for time parallelization”, Mathematics of Computation 75, 1403 (2006). [17] M. Emmett and M. L. Minion, “Toward an efficient parallel in time method for partial differential equations”, Communications in Applied Mathematics and Computational Science 7, 105 (2012). [18] L. Baffico, S. Bernard, Y. Maday, G. Turinici, and G. Zérah, “Parallel-in-time molecular-dynamics simulations”, Physical Review E 66, 057701 (2002). [19] Y. Yanan, A. Srinivasan, and N. Chandra, “Scalable time-parallelization of molecular dynamics simulations in nano mechanics”, in Proceedings of the 2006 international conference on parallel processing, ICPP ’06 (2006), pp. 119– 126. [20] R. Speck, D. Ruprecht, R. Krause, M. Emmett, M. Minion, M. Winkel, and P. Gibbon, “A massively space-time parallel N-body solver”, in Proceedings of the international conference on high performance computing, networking, storage and analysis, SC ’12 (2012), pp. 92-1–92-11. [21] E. J. Bylaska, J. Q. Weare, and J. H. Weare, “Extending molecular simulation time scales: Parallel in time integrations for high-level quantum chemistry and complex force representations”, The Journal of Chemical Physics 139, 074114 (2013). 108 [22] A. L. Blumers, Z. Li, and G. E. Karniadakis, “Supervised parallel-in-time algorithm for long-time Lagrangian simulations of stochastic dynamics: Application to hydrodynamics”, Journal of Computational Physics 393, 214 (2019). [23] F. Legoll, T. Lelièvre, and U. Sharma, “An adaptive Parareal algorithm: Application to the simulation of molecular dynamics trajectories”, SIAM Journal on Scientific Computing 44, B146 (2022). [24] D. Perez, E. D. Cubuk, A. Waterland, E. Kaxiras, and A. F. Voter, “Long-time dynamics through parallel trajectory splicing”, Journal of Chemical Theory and Computation 12, 18 (2016). [25] A. Brünger, C. L. Brooks, and M. Karplus, “Stochastic boundary conditions for molecular dynamics simulations of ST2 water”, Chemical Physics Letters 105, 495 (1984). [26] A. C. Brańka and D. M. Heyes, “Algorithms for Brownian dynamics simulation”, Physical Review E 58, 2611 (1998). [27] G. Bussi and M. Parrinello, “Accurate sampling using Langevin dynamics”, Physical Review E 75, 056707 (2007). [28] N. Bou-Rabee, “Time integrators for molecular dynamics”, Entropy 16, 138 (2014). [29] D. A. Sivak, J. D. Chodera, and G. E. Crooks, “Time step rescaling recovers continuous-time dynamical properties for discrete-time Langevin integration of nonequilibrium systems”, Journal of Physical Chemistry B 118, 6466 (2014). [30] H. Risken and T. Frank, The Fokker–Planck equation: Methods of solution and applications (Springer Berlin Heidelberg, 1996). [31] L. Onsager and S. Machlup, “Fluctuations and irreversible processes”, Physical Review 91, 1505 (1953). [32] K. Yasue, “A simple derivation of the Onsager–Machlup formula for onedimensional nonlinear diffusion process”, Journal of Mathematical Physics 19, 1671 (1978). [33] L. R. Pratt, “A statistical method for identifying transition states in high dimensional problems”, The Journal of Chemical Physics 85, 5045 (1986). [34] R. Olender and R. Elber, “Calculation of classical trajectories with a very large time step: Formalism and numerical examples”, The Journal of Chemical Physics 105, 9299 (1996). [35] D. M. Zuckerman and T. B. Woolf, “Efficient dynamic importance sampling of rare events in one dimension”, Physical Review E 63, 016702 (2000). [36] D. Passerone and M. Parrinello, “Action-derived molecular dynamics in the study of rare events”, Physical Review Letters 87, 108302 (2001). 109 [37] P. G. Bolhuis, D. Chandler, C. Dellago, and P. L. Geissler, “Transition path sampling: Throwing ropes over rough mountain passes, in the dark”, Annual Reviews of Physical Chemistry 53, 291 (2002). [38] T. F. Miller III and C. Predescu, “Sampling diffusive transition paths”, The Journal of Chemical Physics 126, 144102 (2007). [39] J. MacFadyen, J. Wereszczynski, and I. Andricioaei, “Directionally negative friction: A method for enhanced sampling of rare event kinetics”, The Journal of Chemical Physics 128, 114112 (2008). [40] W. E and E. Vanden-Eijnden, “Transition-path theory and path-finding algorithms for the study of rare events”, Annual Review of Physical Chemistry 61, 391 (2010). [41] F. J. Pinski and A. M. Stuart, “Transition paths in molecules at finite temperature”, The Journal of Chemical Physics 132, 184104 (2010). [42] H. Orland, “Generating transition paths by Langevin bridges”, The Journal of Chemical Physics 134, 174114 (2011). [43] S. N. Majumdar and H. Orland, “Effective Langevin equations for constrained stochastic processes”, Journal of Statistical Mechanics: Theory and Experiment 2015, P06039 (2015). [44] P. D. Drummond, “Forward, backward, and weighted stochastic bridges”, Physical Review E 96, 042123 (2017). [45] J. Lee, I.-H. Lee, I. Joung, L. Jooyoung, and B. R. Brooks, “Finding multiple reaction pathways via global optimization of action”, Nature Communications 8, 15443 (2017). [46] D. M. Ceperley, “Path integrals in the theory of condensed helium”, Reviews of Modern Physics 67, 279 (1995). [47] A. M. Stuart, J. Voss, and P. Wilberg, “Conditional path sampling of SDEs and the Langevin MCMC method”, Communications in Mathematical Sciences 2, 685 (2004). [48] G. Stoltz, “Path sampling with stochastic dynamics: Some new algorithms”, Journal of Computational Physics 225, 491 (2007). [49] C. W. Gardiner, Stochastic methods: A handbook for the natural and social sciences, Springer Series in Synergetics (Springer Berlin Heidelberg, 2009). [50] C. Dellago, P. G. Bolhuis, and P. L. Geissler, in Advances in chemical physics (John Wiley & Sons, Inc., 2003), pp. 1–78. [51] T. R. Gingrich and P. L. Geissler, “Preserving correlations between trajectories for efficient path sampling”, The Journal of Chemical Physics 142, 234104 (2015). 110 [52] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of state calculations by fast computing machines”, The Journal of Chemical Physics 21, 1087 (1953). [53] W. K. Hastings, “Monte Carlo sampling methods using Markov chains and their applications”, Biometrika 57, 97 (1970). [54] C. H. Mak and S. Zakharov, “A multigrid algorithm for sampling imaginarytime paths in quantum Monte Carlo simulations”, Journal of Physical Chemistry B 108, 6760 (2004). [55] A. Beskos, G. Roberts, A. Stuart, and J. Voss, “MCMC methods for diffusion bridges”, Stochastics and Dynamics 8, 319 (2008). [56] A. Beskos, F. J. Pinski, J. M. Sanz-Serna, and A. M. Stuart, “Hybrid Monte Carlo on Hilbert spaces”, Stochastic Processes and their Applications 121, 2201 (2011). [57] M. Ottobre, N. S. Pillai, F. J. Pinski, and A. M. Stuart, “A function space HMC algorithm with second order Langevin diffusion limit”, Bernoulli 22, 60 (2016). [58] A. Beskos, M. Girolami, S. Lan, P. E. Farrell, and A. M. Stuart, “Geometric MCMC for infinite-dimensional inverse problems”, Journal of Computational Physics 335, 327 (2017). [59] P. J. Malsom and F. J. Pinski, “Role of Ito’s lemma in sampling pinned diffusion paths in the continuous-time limit”, Physical Review E 94, 042131 (2016). [60] M. Hairer, A. M. Stuart, J. Voss, and P. Wiberg, “Analysis of SPDEs arising in path sampling. Part I: The Gaussian case”, Communications in Mathematical Sciences 3, 587 (2005). [61] A. B. Adib, “Stochastic actions for diffusive dynamics: Reweighting, sampling, and minimization”, The Journal of Physical Chemistry B 112, 5910 (2008). [62] S. L. Cotter, G. O. Roberts, A. M. Stuart, and D. White, “MCMC methods for functions: Modifying old algorithms to make them faster”, Statistical Science 28, 424 (2013). [63] N. Bou-Rabee, “SPECTRWM: Spectral Random Walk Method for the Numerical Solution of Stochastic Partial Differential Equations”, SIAM Review 60, 386 (2018). [64] H. Behringer and R. Eichhorn, “Brownian dynamics simulations with hardbody interactions: Spherical particles”, The Journal of Chemical Physics 137, 164108 (2012). 111 [65] G. Grazioli and I. Andricioaei, “Advances in milestoning. II. Calculating timecorrelation functions from milestoning using stochastic path integrals”, The Journal of Chemical Physics 149, 084104 (2018). [66] L.-T. Da, F. K. Sheong, D.-A. Silva, and X. Huang, “Application of Markov state models to simulate long timescale dynamics of biological macromolecules”, in Protein conformational dynamics, edited by K.-L. Han, X. Zhang, and M.-J. Yang (Springer International Publishing, 2014), pp. 29–66. [67] J. O. B. Tempkin, B. Qi, M. G. Saunders, B. Roux, A. R. Dinner, and J. Weare, “Using multiscale preconditioning to accelerate the convergence of iterative molecular calculations”, The Journal of Chemical Physics 140, 184114 (2014). [68] M. Opper, “Variational inference for stochastic differential equations”, Annalen der Physik 531, 1800233 (2019). [69] V. Harmandaris, E. Kalligiannaki, M. Katsoulakis, and P. Plecháč, “Pathspace variational inference for non-equilibrium coarse-grained systems”, Journal of Computational Physics 314, 355 (2016). [70] S.-T. Tsai, E.-J. Kuo, and P. Tiwary, “Learning molecular dynamics with simple language model built upon long short-term memory neural network”, Nature Communications 11, 5115 (2020).