Ben Leimkuhler

The adaptive buffered force QM/MM method in the CP2K and AMBER software packages

The implementation and validation of the adaptive buffered force (AdBF) quantum‐mechanics/molecular‐mechanics (QM/MM) method in two popular packages, CP2K and AMBER are presented. The implementations build on the existing QM/MM functionality in each code, extending it to allow for redefinition of the QM and MM regions during the simulation and reducing QM‐MM interface errors by discarding forces near the boundary according to the buffered force‐mixing approach. New adaptive thermostats, needed by force‐mixing methods, are also implemented. Different variants of the method are benchmarked by simulating the structure of bulk water, water autoprotolysis in the presence of zinc and dimethyl‐phosphate hydrolysis using various semiempirical Hamiltonians and density functional theory as the QM model. It is shown that with suitable parameters, based on force convergence tests, the AdBF QM/MM scheme can provide an accurate approximation of the structure in the dynamical QM region matching the corresponding fully QM simulations, as well as reproducing the correct energetics in all cases. Adaptive unbuffered force‐mixing and adaptive conventional QM/MM methods also provide reasonable results for some systems, but are more likely to suffer from instabilities and inaccuracies.

Weakly coupled heat bath models for Gibbs-like invariant states in nonlinear wave equations

Thermal bath coupling mechanisms as utilized in molecular dynamics are applied to partial differential equation models. Working from a semi-discrete (Fourier mode) formulation for the Burgers–Hopf or Korteweg–de Vries equation, we introduce auxiliary variables and stochastic perturbations in order to drive the system to sample a target ensemble which may be a Gibbs state or, more generally, any smooth distribution defined on a constraint manifold. We examine the ergodicity of approaches based on coupling of the heat bath to the high wave numbers, with the goal of controlling the ensemble through the fast modes. We also examine different thermostat methods in the extent to which dynamical properties are corrupted in order to accurately compute the average of a desired observable with respect to the invariant distribution. The principal observation of this paper is that convergence to the invariant distribution can be achieved by thermostatting just the highest wave number, while the evolution of the slowest modes is little affected by such a thermostat.

Numerical Methods for Stochastic Molecular Dynamics

In this chapter, we discuss principles for the design of algorithms for canonical sampling, based on the numerical discretization of stochastic dynamics models (such as Langevin dynamics) introduced in the previous chapter. Before we begin our discussion, let us consider the motivation for computing stationary averages using molecular dynamics.

The Canonical Distribution and Stochastic Differential Equations

Until now, we have considered molecular models for isolated collections of atoms. We have seen how to derive equations of motion, and studied the properties of the relevant dynamical systems. We have also observed that the chaotic nature of a typical system imparts a random aspect which suggests that the dynamics defines, at least in some approximate sense, a limiting probability distribution.

The Stability Threshold

Let us recall the method of studying the asymptotic numerical stability of a linear system \(\displaystyle{ \dot{{\boldsymbol z}} ={\boldsymbol A}{\boldsymbol z}, }\) where \({\boldsymbol z} \in \mathbb{R}^{m}\), \({\boldsymbol A} \in \mathbb{R}^{m\times m}\), when solved by a numerical method.

Phase Space Distributions and Microcanonical Averages

In the previous chapters, we considered the approximation of Hamiltonian trajectories. In this chapter we study the paths emanating from the collection of all initial conditions within a given set. This is the starting point for statistical mechanics which allows the calculation of averages.

Analyzing Geometric Integrators

In the previous chapter, we discussed the growth of error in numerical methods for differential equations. We saw that if the time interval is fixed, the error obeys the power law relationship with stepsize that is predicted by the convergence theory. We also saw that this did not contradict the exponential growth in the error with time (when the stepsize is fixed). The latter issue casts doubt on the reliance on the convergence order as a means for assessing the suitability of an integrator for molecular dynamics.

Extended Variable Methods

In the previous chapter we described numerical methods for solving the equations of Langevin dynamics, a system of stochastic differential equations that sample the canonical ensemble. The methods allow us to compute averages of the form

Optimal control for mathematical models of cancer therapies : an application of geometric methods

\displaystyle{\mathrm{Av}_{\beta }(\phi ({\boldsymbol q},{\boldsymbol p})) = \frac{\int _{\mathcal{D}}\phi ({\boldsymbol q},{\boldsymbol p})\rho _{\beta }({\boldsymbol q},{\boldsymbol p})\mathrm{d}\omega } {\int _{\mathcal{D}}\rho _{\beta }({\boldsymbol q},{\boldsymbol p})\mathrm{d}\omega },}