Johannes Kästner

Umbrella sampling: Umbrella sampling

The calculation of free‐energy differences is one of the main challenges in computational biology and biochemistry. Umbrella sampling, biased molecular dynamics (MD), is one of the methods that provide free energy along a reaction coordinate. Here, the method is derived in a historic overview and is compared with related methods like thermodynamic integration, slow growth, steered MD, or the Jarzynski‐based fast‐growth technique. In umbrella sampling, bias potentials along a (one‐ or more‐dimensional) reaction coordinate drive a system from one thermodynamic state to another (e.g., reactant and product). The intermediate steps are covered by a series of windows, at each of which an MD simulation is performed. The bias potentials can have any functional form. Often, harmonic potentials are used for their simplicity. From the sampled distribution of the system along the reaction coordinate, the change in free energy in each window can be calculated. The windows are then combined by methods like the weighted histogram analysis method or umbrella integration. If the bias potential is adapted to result in an even distribution between the end states, then this whole range can be spanned by one window (adaptive‐bias umbrella sampling). In this case, the free‐energy change is directly obtained from the bias. The sampling in each window can be improved by replica exchange methods; either by exchange between successive windows or by running additional simulations at higher temperatures.

DL-FIND: An Open-Source Geometry Optimizer for Atomistic Simulations

Geometry optimization, including searching for transition states, accounts for most of the CPU time spent in quantum chemistry, computational surface science, and solid-state physics, and also plays an important role in simulations employing classical force fields. We have implemented a geometry optimizer, called DL-FIND, to be included in atomistic simulation codes. It can optimize structures in Cartesian coordinates, redundant internal coordinates, hybrid-delocalized internal coordinates, and also functions of more variables independent of atomic structures. The implementation of the optimization algorithms is independent of the coordinate transformation used. Steepest descent, conjugate gradient, quasi-Newton, and L-BFGS algorithms as well as damped molecular dynamics are available as minimization methods. The partitioned rational function optimization algorithm, a modified version of the dimer method and the nudged elastic band approach provide capabilities for transition-state search. Penalty function, gradient projection, and Lagrange-Newton methods are implemented for conical intersection optimizations. Various stochastic search methods, including a genetic algorithm, are available for global or local minimization and can be run as parallel algorithms. The code is released under the open-source GNU LGPL license. Some selected applications of DL-FIND are surveyed.

Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method: "umbrella integration"

We present a method to analyze biased molecular-dynamics and Monte Carlo simulations, also known as umbrella sampling. In the limiting case of a strong bias, this method is equivalent to thermodynamic integration. It employs only quantities with easily controllable equilibration and greatly reduces the statistical errors compared to the standard weighted histogram analysis method. We show the success of our approach for two examples, one analytic function, and one biological system.

Superlinearly converging dimer method for transition state search

Algorithmic improvements of the dimer method [G. Henkelman and H. Jonsson, J. Chem. Phys. 111, 7010 (1999)] are described in this paper. Using the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimizer for the dimer translation greatly improves the convergence compared to the previously used conjugate gradient algorithm. It also saves one energy and gradient calculation per dimer iteration. Extrapolation of the gradient during repeated dimer rotations reduces the computational cost to one gradient calculation per dimer rotation. The L-BFGS algorithm also improves convergence of the rotation. Thus, three to four energy and gradient evaluations are needed per iteration at the beginning of a transition state search, while only two are required close to convergence. Moreover, we apply the dimer method in internal coordinates to reduce coupling between the degrees of freedom. Weighting the coordinates can be used to apply chemical knowledge about the system and restrict the transition state search to only part of the system while minimizing the remainder. These improvements led to an efficient method for the location of transition states without the need to calculate the Hessian. Thus, it is especially useful in large systems with expensive gradient evaluations.

QM/MM Free-Energy Perturbation Compared to Thermodynamic Integration and Umbrella Sampling: Application to an Enzymatic Reaction

We used the free-energy perturbation (FEP) method in quantum mechanics/molecular mechanics (QM/MM) calculations to compute the free-energy profile of the hydroxylation reaction in the enzyme p-hydroxybenzoate hydroxylase (PHBH). k statistics were employed to analyze the FEP sampling including estimation of the sampling error. Various approximations of the free-energy perturbation method were tested. We find that it is adequate not only to freeze the density of the QM part during the dynamics at frozen QM geometry but also to approximate this density by electrostatic-potential-fitted point charges. It is advisable to include all atoms of a QM/MM link in the perturbation. The results of QM/MM-FEP for PHBH are in good agreement with those of thermodynamic integration and umbrella sampling.

ChemShell—a modular software package for QM/MM simulations

ChemShell is a modular computational chemistry package with a particular focus on hybrid quantum mechanical/molecular mechanical (QM/MM) simulations. A core set of chemical data handling modules and scripted interfaces to a large number of quantum chemistry and molecular modeling packages underpin a flexible QM/MM scheme. ChemShell has been used in the study of small molecules, molecular crystals, biological macromolecules such as enzymes, framework materials including zeolites, ionic solids, and surfaces. We outline the range of QM/MM coupling schemes and supporting functions for system setup, geometry optimization, and transition‐state location (including those from the open‐source DL‐FIND optimization library). We discuss recently implemented features allowing a more efficient treatment of long range electrostatic interactions, X‐ray based quantum refinement of crystal structures, free energy methods, and excited‐state calculations. ChemShell has been ported to a range of parallel computers and we describe a number of options including parallel execution based on the message‐passing capabilities of the interfaced packages and task‐farming for applications in which a number of individual QM, MM, or QM/MM calculations can performed simultaneously. We exemplify each of the features by brief reference to published applications.

Analysis of the statistical error in umbrella sampling simulations by umbrella integration

Umbrella sampling simulations, or biased molecular dynamics, can be used to calculate the free-energy change of a chemical reaction. We investigate the sources of different sampling errors and derive approximate expressions for the statistical errors when using harmonic restraints and umbrella integration analysis. This leads to generally applicable rules for the choice of the bias potential and the sampling parameters. Numerical results for simulations on an analytical model potential are presented for validation. While the derivations are based on umbrella integration analysis, the final error estimate is evaluated from the raw simulation data, and it may therefore be generally applicable as indicated by tests using the weighted histogram analysis method.

Locating Instantons in Many Degrees of Freedom.

We implemented and compared four algorithms to locate instantons, i.e., the most likely tunneling paths at a given temperature. These allow to calculate reaction rates, including atom tunneling, down to very low temperature. An instanton is a first-order saddle point of the Euclidean action in the space of closed Feynman paths. We compared the Newton-Raphson method to the partitioned rational function optimization (P-RFO) algorithm, the dimer method, and a newly proposed mode-following algorithm, where the unstable mode is directly estimated from the instanton path. We tested the algorithms on three chemical systems, each including a hydrogen transfer, at different temperatures. Overall, the Newton-Raphson turned out to be the most promising method, with our newly proposed mode following, being the fall-back option.

Hydrogen-Atom Tunneling Could Contribute to H2 Formation in Space†

The interaction of hydrogen atoms with graphite and aromatic hydrocarbons is widely studied because it plays an important role in diverse areas, such as plasma-wall interactions in nuclear fusion reactors, hydrogen-storage materials, petroleum refining, and H2 formation in space. [5–10] H addition to benzene is the rate-determining step in the surface-catalyzed conversion of benzene into cyclohexane. H atoms can also chemisorb on graphite, although the carbon atom to which it binds has to pucker out of the plane, resulting in an activation barrier of about 20 kJmol . Chemisorption on the edge of a polycyclic aromatic hydrocarbon (PAH) is predicted to have a much lower activation barrier and, more strikingly, for some configurations there is no activation barrier for the reaction with a second H atom, yielding either chemisorbed H2 or gaseous H2. [7,15]

Adaptive integration grids in instanton theory improve the numerical accuracy at low temperature

The instanton method allows to accurately calculate tunneling rates down to very low temperature. However, with lowering the temperature, the computational effort steeply increases as many more discretization points are required. This is caused in practical applications by the majority of the discretization points accumulating at a very small region in configuration space. Here, we describe a method to flexibly discretize the instanton path adapted to the temperature. Chosen appropriately, the discretization leads to a much more uniform distribution of the images (control points) along the path which reduces the number of required images by about a factor of two. Combined with a modified Newton-Raphson optimizer and successive updates of the Hessians, the proposed method provides converged reaction rates at computational costs reduced by more than an order of magnitude. We show the success of the method on analytic test potentials and on molecules with energies directly obtained from density functional theory calculations.