Niels Hansen

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Molecular simulation of alkene adsorption in zeolites

The adsorption isotherms of various alkenes and their mixtures in zeolites such as silicalite-1 (MFI-type), theta-1 (TON-type), and deca-dodecasil 3R (DDR-type) were calculated using the grand canonical Monte Carlo (GCMC) approach. Additionally, the adsorption of alkene–alkane mixtures was simulated. The GCMC approach was combined with the configurational-bias Monte Carlo (CBMC) method. Effective Lennard–Jones parameters for the interaction between the oxygen atoms of all-silica zeolites and the sp2-hybridized groups of linear alkenes were determined using a united atom force field. They were adjusted to the experimental adsorption data of silicalite-1 (MFI). The inflection behaviour of the 1-heptene isotherm was investigated in detail. It is shown that, in the inflection region, the 1-heptene molecules alter their end-to-end length depending on their location. The occurrence of a maximum in the mixture adsorption isotherms is attributed to two effects: entropic effects and non-ideality effects. From the mixture simulations some general conclusions concerning the separation of hydrocarbons with silicalite-1 can be drawn. The transferability of the Lennard–Jones parameters to other zeolites was investigated. Simulations of adsorption isotherms in the zeolites theta-1 and DD3R and their comparison with experimental data indicate the possibility of transferring the parameters to other all-silica zeolites.

Combining reactive and configurational-bias Monte Carlo: confinement influence on the propene metathesis reaction system in various zeolites.

In order to efficiently calculate chemical equilibria of large molecules in a confined environment the reactive Monte Carlo technique is combined with the configurational-bias Monte Carlo approach. To prove that detailed balance is fulfilled the acceptance rule for this combination of particular Monte Carlo techniques is derived in detail. Notably, by using this derivation all other acceptance rules of any Monte Carlo trial moves usually carried out in combination with the configurational-bias Monte Carlo approach can be deduced from it. As an application of the combination of reactive and configurational-bias Monte Carlo the influence of different zeolitic confinements (MFI, TON, LTL, and FER) on the reaction equilibrium and the selectivity of the propene metathesis reaction system was investigated. Compared to the bulk phase the conversion is increased significantly. The authors study this reaction system in the temperature range between 300 and 600 K, and the pressure range from 1 to 7 bars. In contrast to the bulk phase, pressure and temperature have a strong influence on the composition of the reaction mixture in confinement. At low pressures and temperatures both conversion and selectivity are highest. Furthermore, the equilibrium composition is strongly dependent on the type of zeolite. This demonstrates the important role of the host structure in catalytic systems.

Reactive Monte Carlo and grand-canonical Monte Carlo simulations of the propene metathesis reaction system

The influence of silicalite-1 pores on the reaction equilibria and the selectivity of the propene metathesis reaction system in the temperature range between 300 and 600 K and the pressure range from 0.5 to 7 bars has been investigated with molecular simulations. The reactive Monte Carlo (RxMC) technique was applied for bulk-phase simulations in the isobaric-isothermal ensemble and for two phase systems in the Gibbs ensemble. Additionally, Monte Carlo simulations in the grand-canonical ensemble (GCMC) have been carried out with and without using the RxMC technique. The various simulation procedures were combined with the configurational-bias Monte Carlo approach. It was found that the GCMC simulations are superior to the Gibbs ensemble simulations for reactions where the bulk-phase equilibrium can be calculated in advance and does not have to be simulated simultaneously with the molecules inside the pore. The confined environment can increase the conversion significantly. A large change in selectivity between the bulk phase and the pore phase is observed. Pressure and temperature have strong influences on both conversion and selectivity. At low pressure and temperature both conversion and selectivity have the highest values. The effect of confinement decreases as the temperature increases.

Comparison of enveloping distribution sampling and thermodynamic integration to calculate binding free energies of phenylethanolamine N-methyltransferase inhibitors.

The relative binding free energy between two ligands to a specific protein can be obtained using various computational methods. The more accurate and also computationally more demanding techniques are the so-called free energy methods which use conformational sampling from molecular dynamics or Monte Carlo simulations to generate thermodynamic averages. Two such widely applied methods are the thermodynamic integration (TI) and the recently introduced enveloping distribution sampling (EDS) methods. In both cases relative binding free energies are obtained through the alchemical perturbations of one ligand into another in water and inside the binding pocket of the protein. TI requires many separate simulations and the specification of a pathway along which the system is perturbed from one ligand to another. Using the EDS approach, only a single automatically derived reference state enveloping both end states needs to be sampled. In addition, the choice of an optimal pathway in TI calculations is not trivial and a poor choice may lead to poor convergence along the pathway. Given this, EDS is expected to be a valuable and computationally efficient alternative to TI. In this study, the performances of these two methods are compared using the binding of ten tetrahydroisoquinoline derivatives to phenylethanolamine N-transferase as an example. The ligands involve a diverse set of functional groups leading to a wide range of free energy differences. In addition, two different schemes to determine automatically the EDS reference state parameters and two different topology approaches are compared.

Current Computer Modeling Cannot Explain Why Two Highly Similar Sequences Fold into Different Structures

The remarkable recent creation of two proteins that fold into two completely different and stable structures, exhibit different functions, yet differ by only a few amino acids poses a conundrum to those hoping to understand how sequence encodes structure. Here, computer modeling uniquely allows the characterization of not only the native structure of each minimally different sequence but also systems in which each sequence was modeled onto the fold of the alternate sequence. The reasons for the different structural preferences of two pairs of highly similar sequences are explored by a combination of structure analyses, comparison of potential energies calculated from energy-minimized single structures and trajectories produced from molecular dynamics simulations, and application of a novel method for calculating free energy differences. The sensitivity of such analyses to the choice of force field is also explored. Many of the hypotheses proposed on the basis of the nuclear magnetic resonance model structures of the proteins with 95% identical sequences are supported. However, each level of analysis provides different predictions regarding which sequence-structure combination should be most favored, highlighting the fact that protein structure and stability result from a complex combination of interdependent factors.

Assessment of enveloping distribution sampling to calculate relative free enthalpies of binding for eight netropsin-DNA duplex complexes in aqueous solution.

The performance of enveloping distribution sampling (EDS) simulations to estimate free enthalpy differences associated with seven alchemical transformations of A‐T into G‐C base pairs at the netropsin binding site in the minor groove of a 13‐base pair DNA duplex in aqueous solution is evaluated. It is demonstrated that sufficient sampling can be achieved with a two‐state EDS Hamiltonian even for large perturbations such as the simultaneous transformation of up to three A‐T into three G‐C base pairs. The two parameters required to define the EDS reference state Hamiltonian are obtained automatically using a modified version of a scheme presented in earlier work. The sensitivity of the configurational sampling to a variation of these parameters is investigated in detail. Although for relatively small perturbations, that is, one base pair, the free enthalpy estimate depends only weakly on the EDS parameters, the sensitivity is stronger for the largest perturbation. Yet, EDS offers various convenient measures to evaluate the degree of sampling and thus the reliability of the free enthalpy estimate and appears to be an efficient alternative to the conventional thermodynamic integration methodology to obtain free energy differences for molecular systems.

Multiscale Modeling of Reaction and Diffusion in Zeolites: From the Molecular Level to the Reactor

A hierarchical multiscale approach is presented to calculate effective rates of reaction for zeolite catalyzed reaction systems. The first step in this approach involves the determination of intrinsic rate coefficients for all elementary reactions by means of quantum chemical calculations combined with transition state theory. The second step is the calculation of adsorption isotherms and diffusion coefficients of all species by means of Monte Carlo and molecular dynamics simulations. The third step comprises a continuum description of a zeolite crystal based on the reaction-diffusion equation. For coupling the intrinsic rate coefficients obtained in the first step to the continuum variables, a model is proposed that calculates the local concentration of reactants at the catalytically active centers from the total concentrations of adsorbed species. The adsorption isotherms and diffusivities are coupled to the continuum variables by analytical theories such as the ideal adsorbed solution theory and the Maxwell-Stefan formulation of multicomponent diffusion. The fourth step involves fixed-bed reactor simulations based on the continuum description of single zeolite crystals. The approach is illustrated by means of the alkylation of benzene with ethene over zeolite H-ZSM-5.

Time-averaged order parameter restraints in molecular dynamics simulations

Abstract A method is described that allows experimental