Shadow Hamiltonian Research Articles

Unification of Ewald and shifted force methods to calculate Coulomb interactions in molecular simulations.

Three new Ewald series are derived using a new strategy that does not start with a proposed charge spreading function. Of these, the Ewald series produced using shifted potential interactions for the Ewald real space series converges relatively slowly, while the corresponding expression using a shifted force (SF) interaction does not converge. A comparison is made between several approximations of the Ewald method and the SF route to include Coulomb interactions in molecular dynamics (MD) computer simulations. MD simulations of a model bulk molten salt and water were carried out. The recently derived α' variant of Ewald, by K. D. Hammonds and D. M. Heyes [J. Chem. Phys. 157, 074108 (2022)], has been developed analytically and found to be more accurate and computationally efficient than SF in part due to the smaller real space truncation distance that can be used. In addition, with α', the number of reciprocal lattice vectors required is reduced considerably compared with the standard Ewald implementations to give the same accuracy. The invention of the α' method shifts the computational balance back toward using an Ewald construction. The SF method shows greater errors in the Coulomb pressure and time dependent fluctuation properties compared to α'. It does not conserve the shadow Hamiltonian in a microcanonical MD simulation, whereas the α' method does, which facilitates long time stability and insignificant drift of properties over time. The speed of the Ewald computer code is improved by using a new lookup table method.

Shadow Hamiltonians of structure-preserving integrators for Nambu mechanics

Abstract Symplectic integrators are widely implemented numerical integrators for Hamiltonian mechanics, which preserve the Hamiltonian structure (symplecticity) of the system. Although the symplectic integrator does not conserve the energy of the system, it is well known that there exists a conserving modified Hamiltonian, called the shadow Hamiltonian. For the Nambu mechanics, which is a kind of generalized Hamiltonian mechanics, we can also construct structure-preserving integrators by the same procedure used to construct the symplectic integrators. In the structure-preserving integrator, however, the existence of shadow Hamiltonians is nontrivial. This is because the Nambu mechanics is driven by multiple Hamiltonians and it is nontrivial whether the time evolution by the integrator can be cast into the Nambu mechanical time evolution driven by multiple shadow Hamiltonians. In this paper we present a general procedure to calculate the shadow Hamiltonians of structure-preserving integrators for Nambu mechanics, and give an example where the shadow Hamiltonians exist. This is the first attempt to determine the concrete forms of the shadow Hamiltonians for a Nambu mechanical system. We show that the fundamental identity, which corresponds to the Jacobi identity in Hamiltonian mechanics, plays an important role in calculating the shadow Hamiltonians using the Baker–Campbell–Hausdorff formula. It turns out that the resulting shadow Hamiltonians have indefinite forms depending on how the fundamental identities are used. This is not a technical artifact, because the exact shadow Hamiltonians obtained independently have the same indefiniteness.

Shadow Magnetic Hamiltonian Monte Carlo

Incorporating partial momentum refreshment into Magnetic Hamiltonian Monte Carlo (MHMC) to create Magnetic Hamiltonian Monte Carlo with partial momentum refreshment (PMHMC) has been shown to improve the sampling performance of MHMC significantly. At the same time, sampling from an integrator-dependent shadow or modified target density has been utilised to boost the acceptance rates of Hamiltonian Monte Carlo (HMC), which leads to more efficient sampling as the integrator is better conserved by the shadow Hamiltonian than the true Hamiltonian. Sampling from the shadow Hamiltonian associated with the numerical integrator used in MHMC is yet to be explored in the literature. This work aims to address this gap in the literature by combining the benefits of the non-canonical Hamiltonian dynamics of MHMC with those achieved by targeting the modified Hamiltonian. We first determine the modified Hamiltonian associated with the MHMC integrator and use this to construct a novel method, which we refer to as Shadow Magnetic Hamiltonian Monte Carlo (SMHMC), that leads to better sampling behaviour when compared to MHMC while leaving the target distribution invariant. The new SMHMC method is compared to MHMC and PMHMC across various target posterior distributions, including datasets modeled using Bayesian Neural Networks and Bayesian Logistic Regression models.

Bayesian inference of local government audit outcomes.

The scandals in publicly listed companies have highlighted the large losses that can result from financial statement fraud and weak corporate governance. Machine learning techniques have been applied to automatically detect financial statement fraud with great success. This work presents the first application of a Bayesian inference approach to the problem of predicting the audit outcomes of financial statements of local government entities using financial ratios. Bayesian logistic regression (BLR) with automatic relevance determination (BLR-ARD) is applied to predict audit outcomes. The benefit of using BLR-ARD, instead of BLR without ARD, is that it allows one to automatically determine which input features are the most relevant for the task at hand, which is a critical aspect to consider when designing decision support systems. This work presents the first implementation of BLR-ARD trained with Separable Shadow Hamiltonian Hybrid Monte Carlo, No-U-Turn sampler, Metropolis Adjusted Langevin Algorithm and Metropolis-Hasting algorithms. Unlike the Gibbs sampling procedure that is typically employed in sampling from ARD models, in this work we jointly sample the parameters and the hyperparameters by putting a log normal prior on the hyperparameters. The analysis also shows that the repairs and maintenance as a percentage of total assets ratio, current ratio, debt to total operating revenue, net operating surplus margin and capital cost to total operating expenditure ratio are the important features when predicting local government audit outcomes using financial ratios. These results could be of use for auditors as focusing on these ratios could potentially speed up the detection of fraudulent behaviour in municipal entities, and improve the speed and quality of the overall audit.

Separable Shadow Hamiltonian Hybrid Monte Carlo for Bayesian Neural Network Inference in wind speed forecasting

Accurate wind speed and consequently wind power forecasts form a critical enabling tool for large scale wind energy adoption. Probabilistic machine learning models such as Bayesian Neural Network (BNN) models are often preferred in the forecasting task as they facilitate estimates of predictive uncertainty and automatic relevance determination (ARD). Hybrid Monte Carlo (HMC) is widely used to perform asymptotically exact inference of the network parameters. A significant limitation to the increased adoption of HMC in inference for large scale machine learning systems is the exponential degradation of the acceptance rates and the corresponding effective sample sizes with increasing model dimensionality due to numerical integration errors. This paper presents a solution to this problem by sampling from a modified or shadow Hamiltonian that is conserved to a higher-order by the leapfrog integrator. BNNs trained using Separable Shadow Hamiltonian Hybrid Monte Carlo (S2HMC) are used to forecast one hour ahead wind speeds on the Wind Atlas for South Africa (WASA) datasets. Experimental results find that S2HMC yields higher effective sample sizes than the competing HMC. The predictive performance of S2HMC and HMC based BNNs is found to be similar. We further perform hierarchical inference for BNN parameters by combining the S2HMC sampler with Gibbs sampling of hyperparameters for relevance determination. A generalisable ARD committee framework is introduced to synthesise the various sampler ARD outputs into robust feature selections. Experimental results show that this ARD committee approach selects features of high predictive information value. Further, the results show that dimensionality reduction performed through this approach improves the sampling performance of samplers that suffer from random walk behaviour such as Metropolis–Hastings (MH).

Shadow Hamiltonian in classical NVE molecular dynamics simulations involving Coulomb interactions.

Microcanonical ensemble (NVE) Molecular Dynamics (MD) computer simulations are performed with negligible energy drift for systems incorporating Coulomb interactions and complex constraint schemes. In principle, such systems can now be simulated in the NVE ensemble for millisecond time scales, with no requirement for system thermostatting. Numerical tools for assessing drift in MD simulations are outlined, and drift rates of 10-6 K/μs are demonstrated for molten salts, polar liquids, and room temperature ionic liquids. Such drift rates are six orders of magnitude smaller than those typically quoted in the literature. To achieve this, the standard Ewald method is slightly modified so the first four derivatives of the real space terms go smoothly to zero at the truncation distance, rc. New methods for determining standard Ewald errors and the new perturbation errors introduced by the smoothing procedure are developed and applied, these taking charge correlation effects explicitly into account. The shadow Hamiltonian, Es, is shown to be the strictly conserved quantity in these systems, and standard errors in the mean of one part in 1010 are routinely calculated. Expressions for the shadow Hamiltonian are improved over previous work by accounting for O(h4) terms, where h is the MD time step. These improvements are demonstrated by means of extreme out-of-equilibrium simulations. Using the new methodology, the very low diffusion coefficients of room temperature 1-hexyl-3-methyl-imidazolium chloride are determined from long NVE trajectories in which the equations of motion are known to be integrated correctly, with negligible drift.

Utilising Partial Momentum Refreshment in Separable Shadow Hamiltonian Hybrid Monte Carlo

Sampling using integrator-dependent shadow Hamiltonian’s has been shown to produce improved sampling properties relative to Hamiltonian Monte Carlo. The shadow Hamiltonian’s are typically non-separable, requiring the expensive generation of momenta, with the recent trend being to utilise partial momentum refreshment. Separable Shadow Hamiltonian Hybrid Monte Carlo (S2HMC) employs a canonical transformation which results in the Hamiltonian being separable and makes use of a processed leapfrog integrator. In this work, we combine the benefit of sampling using S2HMC with partial momentum refreshment to create the Separable Shadow Hamiltonian Hybrid Monte Carlo with Partial Momentum Refreshment (PS2HMC) algorithm which leaves the target distribution invariant. Numerical experiments across various targets show that the proposed algorithm outperforms S2HMC and Shadow Hamiltonian Monte Carlo with partial momentum refreshment. Comprehensive analysis is performed on the Banana shaped distribution, multivariate Gaussian distributions of various dimensions, Bayesian logistic regression and Bayesian neural networks.

Antithetic Magnetic and Shadow Hamiltonian Monte Carlo

Hamiltonian Monte Carlo is a Markov Chain Monte Carlo method that has been widely applied to numerous posterior inference problems within the machine learning literature. Markov Chain Monte Carlo estimators have higher variance than classical Monte Carlo estimators due to autocorrelations present between the generated samples. In this work we present three new methods for tackling the high variance problem in Hamiltonian Monte Carlo based estimators: 1) We combine antithetic and importance sampling techniques where the importance sampler is based on sampling from a modified or shadow Hamiltonian using Separable Shadow Hamiltonian Hybrid Monte Carlo, 2) We present the antithetic Magnetic Hamiltonian Monte Carlo algorithm that is based on performing antithetic sampling on the Magnetic Hamiltonian Monte Carlo algorithm and 3) We propose the antithetic Magnetic Momentum Hamiltonian Monte Carlo algorithm based on performing antithetic sampling on the Magnetic Momentum Hamiltonian Monte Carlo method. We find that the antithetic Separable Shadow Hamiltonian Hybrid Monte Carlo and antithetic Magnetic Momentum Hamiltonian Monte Carlo algorithms produce effective sample sizes that are higher than antithetic Hamiltonian Monte Carlo on all the benchmark datasets. We further find that antithetic Separable Shadow Hamiltonian Hybrid Monte Carlo and antithetic Magnetic Hamiltonian Monte Carlo produce higher effective sample sizes normalised by execution time in higher dimensions than antithetic Hamiltonian Monte Carlo. In addition, the antithetic versions of all the algorithms have higher effective sample sizes than their non-antithetic counterparts, indicating the usefulness of adding antithetic sampling to Markov Chain Monte Carlo algorithms. The methods are assessed on benchmark datasets using Bayesian logistic regression and Bayesian neural network models.

Adaptively Setting the Path Length for Separable Shadow Hamiltonian Hybrid Monte Carlo

Hybrid Monte Carlo (HMC) has been widely applied to numerous posterior inference problems in machine learning and statistics. HMC has two main practical issues, the first is the deterioration in acceptance rates as the system size increases and the second is its sensitivity to two user-specified parameters: the step size and trajectory length. The former issue is addressed by sampling from an integrator-dependent modified or shadow density and compensating for the induced bias via importance sampling. The latter issue is addressed by adaptively setting the HMC parameters, with the state-of-the-art method being the No-U-Turn Sampler (NUTS). We combine the benefits of NUTS with those attained by sampling from the shadow density, by adaptively setting the trajectory length and step size of Separable Shadow Hamiltonian Hybrid Monte Carlo (S2HMC). This leads to a new algorithm, adaptive S2HMC (A-S2HMC), that shows improved performance over S2HMC and NUTS across various targets and leaves the target density invariant.

Shadow Hamiltonian in classical NVE molecular dynamics simulations: A path to long time stability.

The shadow energy, Es, is the conserved quantity in microcanonical ensemble (NVE) molecular dynamics simulations carried out with the position Verlet central-difference algorithm. A new methodology for calculating precise and accurate values of Es is presented. It is shown for the first time that Es rather than E is constant during structural changes occurring within a supercooled liquid. It is also explained how to prepare and conduct microsecond range bulk-phase NVE simulations with essentially zero energy drift without the need for thermostating. The drift is analyzed with block averaging and new drift functions of the shadow energy. With such minimal drift, extremely small and accurate standard errors in the mean for quantities like Es, E, and temperature, T, can be obtained. Values of the standard error for Es of ≈10-10 in molecule-based reduced units can be routinely achieved for simulations of 108 time steps. This corresponds to a simulation temperature drift of ≈10-6 K/μs, six orders of magnitude smaller than generally considered to be acceptable for protein simulations. We also show for the first time how these treatments can be extended with no loss of accuracy to polyatomic systems with both flexible degrees of freedom and arbitrary geometric constraints imposed via the SHAKE algorithm. As a bonus, estimates of simulation-average kinetic and total energies from high order velocity expressions can be obtained to a good approximation from 2nd order velocities and the average mean square force (for polyatomics, this refers to per site, including any constraint forces).

On the accuracy of symplectic integrators for secularly evolving planetary systems

ABSTRACT Symplectic integrators have made it possible to study the long-term evolution of planetary systems with direct N-body simulations. In this paper we reassess the accuracy of such simulations by running a convergence test on 20 Myr integrations of the Solar System using various symplectic integrators. We find that the specific choice of metric for determining a simulation’s accuracy is important. Only looking at metrics related to integrals of motions such as the energy error can overestimate the accuracy of a method. As one specific example, we show that symplectic correctors do not improve the accuracy of secular frequencies compared to the standard Wisdom–Holman method without symplectic correctors, despite the fact that the energy error is three orders of magnitudes smaller. We present a framework to trace the origin of this apparent paradox to one term in the shadow Hamiltonian. Specifically, we find a term that leads to negligible contributions to the energy error but introduces non-oscillatory errors that result in artificial periastron precession. This term is the dominant error when determining secular frequencies of the system. We show that higher order symplectic methods such as the Wisdom–Holman method with a modified kernel or the SABAC family of integrators perform significantly better in secularly evolving systems because they remove this specific term.

Refinement of thermostated molecular dynamics using backward error analysis.

Kinetic energy equipartition is a premise for many deterministic and stochastic molecular dynamics methods that aim at sampling a canonical ensemble. While this is expected for real systems, discretization errors introduced by the numerical integration may lead to deviations from equipartition. Fortunately, backward error analysis allows us to obtain a higher-order estimate of the quantity that is actually subject to equipartition. This is related to a shadow Hamiltonian, which converges to the specified Hamiltonian only when the time-step size approaches zero. This paper deals with discretization effects in a straightforward way. With a small computational overhead, we obtain refined versions of the kinetic and potential energies, whose sum is a suitable estimator of the shadow Hamiltonian. Then, we tune the thermostatting procedure by employing the refined kinetic energy instead of the conventional one. This procedure is shown to reproduce a canonical ensemble compatible with the refined system, as opposed to the original one, but canonical averages regarding the latter can easily be recovered by reweighting. Water, modeled as a rigid body, is an excellent test case for our proposal because its numerical stability extends up to time steps large enough to yield pronounced discretization errors in Verlet-type integrators. By applying our new approach, we were able to mitigate discretization effects in equilibrium properties of liquid water for time-step sizes up to 5 fs.

Stability of velocity-Verlet- and Liouville-operator-derived algorithms to integrate non-Hamiltonian systems

We investigate the difference between the velocity Verlet and the Liouville-operator-derived (LOD) algorithms by studying two non-Hamiltonian systems, one dissipative and the other conservative, for which the Jacobian of the transformation can be determined exactly. For the two systems, we demonstrate that (1) the velocity Verlet scheme fails to integrate the former system while the first- and second-order LOD schemes succeed and (2) some first-order LOD fails to integrate the latter system while the velocity Verlet and the other first- and second-order schemes succeed. We have shown that the LOD schemes are stable for the former system by determining the explicit forms of the shadow Hamiltonians which are exactly conserved by the schemes. We have shown that the Jacobian of the velocity Verlet scheme for the former system and that of the first-order LOD scheme for the latter system are always smaller than the exact values, and therefore, the schemes are unstable. The decomposition-order dependence of LOD schemes is also considered.

Tuning the hybrid Monte Carlo algorithm using molecular dynamics forces’ variances

Within the HMC algorithm, we discuss how, by using the shadow Hamiltonian and the Poisson brackets, one can achieve a simple factorization in the dependence of the Hamiltonian violations upon either the algorithmic parameters or the parameters specifying the integrator. We consider the simplest case of a second order (nested) Omelyan integrator and one level of Hasenbusch splitting of the determinant for the simulations of a QCD-like theory (with gauge group SU(2)). Given the specific choice of the integrator, the Poisson brackets reduce to the variances of the molecular dynamics forces. We show how the factorization can be used to optimize in a very economical and simple way both the algorithmic and the integrator parameters with good accuracy.

Time step and shadow Hamiltonian in molecular dynamics simulations

We examine the time step and the shadow Hamiltonian of symplectic algorithms for a bound system of a simple harmonic oscillator as a specific example. The phase space trajectory moves on the hyperplane of a constant shadow Hamiltonian. We find a stationary condition for the time step τ n with which the motion repeats itself on the phase space with a period n. Interestingly, that the time steps satisfying the stationary condition turn out to be independent of the symplectic algorithms chosen. Furthermore, the phase volume enclosed by the phase trajectory is given by n τ n Ẽ n , where Ẽ n is the initial shadow energy of the corresponding symplectic algorithm.

Shadow Hamiltonians, Poisson brackets, and gauge theories

Numerical lattice gauge theory computations to generate gauge field configurations including the effects of dynamical fermions are usually carried out using algorithms that require the molecular dynamics evolution of gauge fields using symplectic integrators. Sophisticated integrators are in common use but are hard to optimise, and force-gradient integrators show promise especially for large lattice volumes. We explain why symplectic integrators lead to very efficient Monte Carlo algorithms because they exactly conserve a shadow Hamiltonian. The shadow Hamiltonian may be expanded in terms of Poisson brackets, and can be used to optimize the integrators. We show how this may be done for gauge theories by extending the formulation of Hamiltonian mechanics on Lie groups to include Poisson brackets and shadows, and by giving a general method for the practical computation of forces, force-gradients, and Poisson brackets for gauge theories.

Combining stochastic and deterministic approaches within high efficiency molecular simulations

Abstract Generalized Shadow Hybrid Monte Carlo (GSHMC) is a method for molecular simulations that rigorously alternates Monte Carlo sampling from a canonical ensemble with integration of trajectories using Molecular Dynamics (MD). While conventional hybrid Monte Carlo methods completely re-sample particle’s velocities between MD trajectories, our method suggests a partial velocity update procedure which keeps a part of the dynamic information throughout the simulation. We use shadow (modified) Hamiltonians, the asymptotic expansions in powers of the discretization parameter corresponding to timestep, which are conserved by symplectic integrators to higher accuracy than true Hamiltonians. We present the implementation of this method into the highly efficient MD code GROMACS and demonstrate its performance and accuracy on computationally expensive systems like proteins in comparison with the molecular dynamics techniques already available in GROMACS. We take advantage of the state-of-the-art algorithms adopted in the code, leading to an optimal implementation of the method. Our implementation introduces virtually no overhead and can accurately recreate complex biological processes, including rare event dynamics, saving much computational time compared with the conventional simulation methods.

Stability of molecular dynamics simulations of classical systems

The existence of a shadow Hamiltonian H̃ for discrete classical dynamics, obtained by an asymptotic expansion for a discrete symplectic algorithm, is employed to determine the limit of stability for molecular dynamics (MD) simulations with respect to the time-increment h of the discrete dynamics. The investigation is based on the stability of the shadow energy, obtained by including the first term in the asymptotic expansion, and on the exact solution of discrete dynamics for a single harmonic mode. The exact solution of discrete dynamics for a harmonic potential with frequency ω gives a criterion for the limit of stability h ≤ 2/ω. Simulations of the Lennard-Jones system and the viscous Kob-Andersen system show that one can use the limit of stability of the shadow energy or the stability criterion for a harmonic mode on the spectrum of instantaneous frequencies to determine the limit of stability of MD. The method is also used to investigate higher-order central difference algorithms, which are symplectic and also have shadow Hamiltonians, and for which one can also determine the exact criteria for the limit of stability of a single harmonic mode. A fourth-order central difference algorithm gives an improved stability with a factor of 3, but the overhead of computer time is a factor of at least two. The conclusion is that the second-order "Verlet"-algorithm, most commonly used in MD, is superior. It gives the exact dynamics within the limit of the asymptotic expansion and this limit can be estimated either from the conserved shadow energy or from the instantaneous spectrum of harmonic modes.

Energy conservation in molecular dynamics simulations of classical systems

Classical Newtonian dynamics is analytic and the energy of an isolated system is conserved. The energy of such a system, obtained by the discrete "Verlet" algorithm commonly used in molecular dynamics simulations, fluctuates but is conserved in the mean. This is explained by the existence of a "shadow Hamiltonian" H [S. Toxvaerd, Phys. Rev. E 50, 2271 (1994)], i.e., a Hamiltonian close to the original H with the property that the discrete positions of the Verlet algorithm for H lie on the analytic trajectories of H. The shadow Hamiltonian can be obtained from H by an asymptotic expansion in the time step length. Here we use the first non-trivial term in this expansion to obtain an improved estimate of the discrete values of the energy. The investigation is performed for a representative system with Lennard-Jones pair interactions. The simulations show that inclusion of this term reduces the standard deviation of the energy fluctuations by a factor of 100 for typical values of the time step length. Simulations further show that the energy is conserved for at least one hundred million time steps provided the potential and its first four derivatives are continuous at the cutoff. Finally, we show analytically as well as numerically that energy conservation is not sensitive to round-off errors.

Error analysis of molecular dynamics and fractal time approximants from a combinatorial perspective

Trotter's theorem forms the theoretical basis of most modern molecular dynamics. In essence this theorem states that a time displacement operator (a Lie operator) constructed by exponentiating a sum of noncommuting operators can be approximated by a product of single operators provided the time interval is "very small." In theory "very small" implies infinitesimally small (at which point the approximate product becomes exact), while in practical analysis a finite time interval is divided into several small subintervals or steps. It follows, therefore, that the larger the number of steps the better the approximation to the exact time displacement operator. The question therefore arises: How many steps are sufficient? For bounded operators, standard theorems are available to provide the answer. In this paper we show that a very simple combinatorial formula can be derived which allows the computation of the global differences (as a function of the number of steps) between the Taylor coefficients of the exact time displacement operator and an approximate one constructed by using a finite number of steps. The formula holds for both bounded and nonbounded operators and shows, quantitatively, what is qualitatively expected-that the error decreases with increasing number of steps. Furthermore, the formula applies irrespective of the complexity of the system, boundary conditions, or the thermodynamic ensemble employed for averaging the initial conditions. The analysis yields explicit expressions for the Taylor coefficients which are then used to compute the errors. In the case of the algorithmically based practical numerical simulations in which fixed, albeit small, steps are repeatedly applied, the rise in the number of steps does not reduce the size of the steps but increases the total time of interest. The combinatorial formula shows that, here, the errors diverge. Furthermore, this work can be used to supplement other efforts such as the use of shadow Hamiltonians where the truncation of the series expansion of the latter will produce errors in the higher order propagator moments.