Accurate Free Energy Research Articles

Anharmonic thermodynamics of vacancies using a neural network potential

Lattice anharmonicity is thought to strongly affect vacancy concentrations in metals at high temperatures. It is however non-trivial to account for this effect directly using density functional theory (DFT). Here we develop a deep neural network potential for aluminum that overcomes the limitations inherent to DFT, and we use it to obtain accurate anharmonic vacancy formation free energies as a function of temperature. While confirming the important role of anharmonicity at high temperatures, the calculation unveils a markedly nonlinear behavior of the vacancy formation entropy and shows that the vacancy formation free energy only violates Arrhenius law at temperatures above 600 K, in contrast with previous DFT calculations.

Free energy of formation of SPC/E-water and TIP4P-water using BAR and TI in MD and MC

Accurate free energy difference calculation based on the molecular dynamics (MD) and Monte-Carlo (MC) simulations is still a challenge in many different areas of science. In this paper, the free energy difference of the SPC/E-water and the TIP4P-water between two subsequent systems will be calculated using the thermodynamic integration (TI) and the Bennett acceptance ratio (BAR) methods. We compare the accuracy and efficiency of the two methods by calculating the free energy difference at three different temperatures, namely, T = 220, 240, 260, 280 and 300 K. The simulated results show that the free energy difference of small clusters scale nicely to the inverse of TCT−1, with the corresponding right critical temperature, TC, of SPC/E-water of 625 K and of TIP4P-water of 588 K.

Accurate Quantum Chemical Free Energies at Affordable Cost

Free energy sampling methods allow studying the full dynamics of activated processes. Unfortunately, the affordable accuracy of the potential describing the energy and forces of the system is usually rather low. Here we introduce a new method that by combining metadynamics and free energy perturbation allows calculating accurate quantum chemical free energies for chemical reactions. To prove the effectiveness of this new approach we study the SN2 reaction of CH3F + Cl- → CH3Cl + F- in vacuo and solvated by water. Comparisons are made with harmonic transition-state theory to show how this method could provide accurate equilibrium and rate constants for complex systems.

Recent Advances of In-Silico Modeling of Potent Antagonists for the Adenosine Receptors.

The rapid advancement of computer architectures and development of mathematical algorithms offer a unique opportunity to leverage the simulation of macromolecular systems at physiologically relevant timescales. Herein, we discuss the impact of diverse structure-based and ligand-based molecular modeling techniques in designing potent and selective antagonists against each adenosine receptor (AR) subtype that constitutes multitude of drug targets. The efficiency and robustness of high-throughput empirical scoring function-based approaches for hit discovery and lead optimization in the AR family are assessed with the help of illustrative examples that have led to nanomolar to sub-micromolar inhibition activities. Recent progress in computer-aided drug discovery through homology modeling, quantitative structure-activity relation, pharmacophore models, and molecular docking coupled with more accurate free energy calculation methods are reported and critically analyzed within the framework of structure-based virtual screening of AR antagonists. Later, the potency and applicability of integrated molecular dynamics (MD) methods are addressed in the context of diligent inspection of intricated AR-antagonist binding processes. MD simulations are exposed to be competent for studying the role of the membrane as well as the receptor flexibility toward the precise evaluation of the biological activities of antagonistbound AR complexes such as ligand binding modes, inhibition affinity, and associated thermodynamic and kinetic parameters.

Vibrational Spectroscopy in Solution through Perturbative ab Initio Molecular Dynamics Simulations.

We have developed a method that allows computing the vibrational spectra at a high quantum mechanical level for molecules in solution or other complex systems. The method is based on the use of configurational samplings from combined QM/MM molecular dynamics simulations and the use of perturbation theory to calculate accurate molecular properties. Such calculations provide in addition accurate free energy gradient vectors and Hessian matrices and thus open the door for the characterization of stationary points in free energy landscapes and the study of chemical reaction mechanisms in large disordered systems. The vibrational spectrum of the water molecule in liquid water has been computed as a test case. It has been obtained using a weighted average of instantaneous signals assuming the instantaneous normal modes approach. Vibrational frequencies are also computed by diagonalizing the Hessian of the free energy surface. Comparison is made with experimental data and with calculations using the Fourier transform of the time autocorrelation function of the dipole moment. The discussion emphasizes the advantages of the developed methodology compared to other techniques in terms of the accuracy/computational cost ratio.

Phq: A Fortran code to compute phonon quasiparticle properties and dispersions

We here introduce a Fortran code that computes anharmonic free energy of solids from first-principles based on our phonon quasiparticle approach. In this code, phonon quasiparticle properties, i.e., renormalized phonon frequencies and lifetimes, are extracted from mode-projected velocity auto-correlation functions (VAF) of modes sampled by molecular dynamics (MD) trajectories. Using renormalized frequencies as input, the code next constructs an effective harmonic force constant matrix to calculate anharmonic phonon dispersions over the whole Brillouin zone and thus the anharmonic free energy in the thermodynamic limit (N→∞). A detailed description of all the input parameters and the subroutines is provided as well. We illustrate the use of this code to compute ab initio temperature-dependent anharmonic phonons of Si in the diamond structure. Program summaryProgram Title:phqProgram Files doi:http://dx.doi.org/10.17632/sk4jsjc6p9.1Licensing provisions: GPLv3Programming language: Fortran 90Nature of problem: Accurate free energy calculations deliver predictive thermodynamic properties of solids. Although the quasi-harmonic approximation (QHA) has been widely used for various materials, its validity at very high temperatures is not guaranteed because of intrinsic anharmonic effects. The QHA cannot be used even at low temperatures for strongly anharmonic systems, e.g., those that are stabilized by anharmonic fluctuations. When anharmonicity is non-negligible, phonon frequencies display a pronounced dependence on temperature, which impacts thermodynamic properties. Therefore, the calculation of anharmonic phonon dispersions is critical for an accurate estimation of the free energy.Solution method: We calculate phonon quasiparticle properties, i.e., renormalized phonon frequencies and lifetimes, by combining molecular dynamics (MD) trajectories and harmonic normal modes (HNM). MD velocities are projected into HNMs and the mode-projected velocity auto-correlation functions (VAF) are then calculated. Quasiparticle properties are then obtained by analyzing the VAFs, if they are well defined. Using the renormalized phonon frequencies and HNMs we build an effective (temperature-dependent) harmonic force constant matrix from which the anharmonic phonon dispersions can be obtained.

Computational Insights into Endo/Exo Selectivity of the Diels-Alder Reaction in Explicit Solvent at Ab Initio Quantum Mechanical/Molecular Mechanical Level.

The microscopic insight into the endo/exo stereoselectivity of the Diels-Alder (DA) reaction between cyclopentadiene and methyl vinyl ketone (MVK) has posed a challenge to the computational chemists, which requires an accurate free-energy (FE) landscape. Although ab initio (ai) quantum mechanical/molecular mechanical (QM/MM) calculations are capable of providing accurate conformation energies, they are far too expensive for free-energy calculations, in which large-scale sampling in the phase space is always indispensable. Recently, on the basis of the idea of reference-potential methods, a new approach termed MBAR+wTP method was proposed by us ( Li ; J. Chem. Theory Comput. 2018 , 14 , 5583 ) for the calculations of ai QM/MM FE profiles with a much less computational expense. In this work, we applied this method to investigate the endo/exo stereoselectivity of the Diels-Alder (DA) reaction at B3LYP/MM level and the results indicate that this method can yield more accurate activation free energies than the semiempirical Hamiltonian PM6 and the "sampling-free" quantum mechanical methods at B3LYP and MP2 levels using a polarizable continuous solvent model (PCM). The stereoselectivity mainly comes from the solvation effect. At the transition state, the first peak of the solvent distribution near the oxygen atom in MVK is slightly closer for the endo pathway than that for the exo pathway, whereas the difference at PM6 level is negligible. Nonetheless, in terms of the stereoselectivity, we obtained an endo/exo ratio of 2 with this method, which is still 1 order of magnitude smaller than the experimental measurement. Although the endo/exo ratio from the MP2/PCM calculation agrees well with the experiment, the absolute reaction barriers deviate by about 5 kcal mol-1. This work issues a warning to the use of the computational methods that are routinely employed for the study of organic reactions in the condensed phase, and it also brings a new method that is affordable and more reliable.

On the origin of the difference between type A and type B skeletal isomerization of alkenes catalyzed by zeolites: The crucial input of ab initio molecular dynamics

Alkene skeletal isomerization elementary steps catalyzed by acid zeolites are key reactions in refining, petrochemistry and biomass conversion. We unravel the atomic-scale origin of the higher rate constant of type A isomerization (involving a direct alkyl transfer, without any change in the branching degree) than the one of type B isomerization (involving non-classical carbonium ions such as protonated cyclopropane (PCP), inducing a change in the branching degree) of C7 carbenium ions in chabazite. Accurate free energy barriers are calculated at 300 and 500 K for both reactions by means of molecular dynamics in combination with blue moon ensemble approach, whereas the static approach is shown to fail to describe these reactions. The slow transformation between individual rotational isomers, causing non-ergodic sampling of reactant state, largely overlooked in literature, is carefully addressed in the present work. At 500 K (representative of experimental conditions), free energy barriers of 83.4 kJ/mol and 15.0 kJ/mol are determined for type B and type A isomerization respectively. The much lower barrier for type A is thus recovered, and assigned to a loose transition state, with free rotation of the migrating alkyl group, while the transition state of type B isomerization is tighter, with such a rotation blocked, due to the simultaneous hydride shift taking place on the edge of the PCP.

Calculation of Protein-Protein Binding Free Energies using Umbrella Sampling with Dual Resolution Water Models

The binding free energy between two proteins is a measure of how strongly they will interact under physiological conditions. All-atom statistical-mechanics-based free energy calculation methods have the potential to accurately calculate binding free energy but are completely unfeasible for large-scale applications due to their high computational cost. This motivates the development of efficient and accurate free energy calculation approaches. Dual-resolution molecular simulation approach allows one to simulate at the atomic level in the regions of interest and at a coarser level, but much faster, in less relevant parts of the biomolecular systems. This hybrid approach has shown a substantial reduction of computational cost without compromising accuracy in a classical molecular dynamics simulations. Here, we investigate whether combining the statistical-mechanics-based free energy calculation method, umbrella sampling, with dual resolution water models provides a good trade-off between speed and accuracy. We chose Barnase-Barstar protein-protein complex as a test system since it is relatively small to allow convergence of the simulations, has a known structure, and empirical relative binding free energy difference values (ΔΔG) of their observed mutations are reported in the literature. We selected eight different mutations occurring at the different sites with broadly varying empirical ΔΔG values. We then calculated ΔΔG values for each mutation using umbrella sampling simulations, where we treated protein-protein complex and one surrounding water layer at an atomic level and the bulk of the water using WT4 coarse-grained water model. We obtained a high correlation between experimental and calculated ΔΔG values. This dual resolution free energy calculation approach has the potential for future high-throughput studies of protein-protein binding.

Determining Free Energy Differences Through Non-Linear Morphing

The Bennett Acceptance Ratio (BAR) method is widely used to estimate the free energy difference between two given states of a many-body system, such as a macromolecule, where Molecular Dynamics or Monte Carlo simulations are typically used for sampling. Because insufficient phase space overlap of the two states can lead to large biases in the results, the calculation is commonly split up into several smaller steps using intermediate states defined along a so-called morphing path — typically a linear interpolation between the start and the end states' Hamiltonians. Here, using BAR, we develop an optimization scheme to find the path that yields the most accurate free energy values from the much more general class of paths that includes also non-linear interpolations. For a broad range of one-dimensional test systems, we show that the free energy sampling error can be decreased by 25% to 35% compared to linear interpolations, with a theoretical upper limit of 50%. Further, correlations arising from using the same sample in one state to evaluate differences to more than one other state can be taken into account. Finally, we demonstrate that the optimization is still valid for a large number of intermediate states with very sparse sampling for each, enabling further optimization by trading off sampling versus number of intermediate states.

Standardization of virtual-screening and post-processing protocols relevant to in-silico drug discovery.

Structure-based drug discovery has emerged as a powerful tool in computational drug discovery and has gained rapid acceleration due to the development of better algorithms for high-end computation in an affordable time. Molecular docking and virtual screening methods are routinely used for the purpose but computing the ligand binding energies with inbuilt scoring functions accurately is still a limitation. Although, MMPBSA and MMGBSA are routinely employed tools for achieving accurate binding free energies, they are applied on well-equilibrated explicitly solvated systems and are computationally expensive and time-consuming. This study compares different post-processing protocols performed on an in-silico screened benchmarked P. falciparum Dihydrofolate reductase (PfDHFR) dataset with AutoDockVina. The docked and implicitly solvated complexes were subjected to (1) rescoring, (2) energy minimization and (3) Binding Estimation After Refinement (BEAR) algorithm. Subsequently, binding free energies were computed using three different tools-MMPBSA, MMGBSA and "g_mmpbsa". Surprisingly, rescoring alone displays lower accuracy than the inherent scoring function of the AutoDockVina. However, encouraging results were seen after post-processing with the other two protocols. The results suggest that MMPBSA applied on energy minimized conformations is able to achieve 42-fold reductions in computational time as opposed to the BEAR algorithm with comparable accuracy.

Targeted Molecular Dynamics Calculations of Free Energy Profiles Using a Nonequilibrium Friction Correction.

As standard unbiased molecular dynamics (MD) simulations become impractical for sampling rare events, "targeted MD" employs a moving distance constraint to enforce rare transitions along some reaction coordinate x. To calculate free energy profiles from these nonequilibrium simulations via Δ G( x) = W( x) - Wdiss( x), apart from the (readily obtained) work W( x) performed on the system, the dissipated work Wdiss( x) is also required. By employing a second-order cumulant expansion of Jarzynski's equality combined with an analysis within Langevin theory, the dissipated work can be expressed via a nonequilibrium friction coefficient ΓNEQ( x) that may be calculated on-the-fly from constraint force fluctuations. Adopting the ion dissociation of NaCl in water as a test system, this friction correction is shown to result in accurate free energy profiles, even for a modest number of simulations and at high constraint velocities. As a bonus, the analysis of ΓNEQ( x) may yield valuable insight into the microscopic mechanism of friction.

Approaching protein design with multisite λ dynamics: Accurate and scalable mutational folding free energies in T4 lysozyme.

The estimation of changes in free energy upon mutation is central to the problem of protein design. Modern protein design methods have had remarkable success over a wide range of design targets, but are reaching their limits in ligand binding and enzyme design due to insufficient accuracy in mutational free energies. Alchemical free energy calculations have the potential to supplement modern design methods through more accurate molecular dynamics based prediction of free energy changes, but suffer from high computational cost. Multisite λ dynamics (MSλD) is a particularly efficient and scalable free energy method with potential to explore combinatorially large sequence spaces inaccessible with other free energy methods. This work aims to quantify the accuracy of MSλD and demonstrate its scalability. We apply MSλD to the classic problem of calculating folding free energies in T4 lysozyme, a system with a wealth of experimental measurements. Single site mutants considering 32 mutations show remarkable agreement with experiment with a Pearson correlation of 0.914 and mean unsigned error of 1.19 kcal/mol. Multisite mutants in systems with up to five concurrent mutations spanning 240 different sequences show comparable agreement with experiment. These results demonstrate the promise of MSλD in exploring large sequence spaces for protein design.

Protocol for Identifying Accurate Collective Variables in Enhanced Molecular Dynamics Simulations for the Description of Structural Transformations in Flexible Metal-Organic Frameworks.

Various kinds of flexibility have been observed in metal–organic frameworks, which may originate from the topology of the material or the presence of flexible ligands. The construction of free energy profiles describing the full dynamical behavior along the phase transition path is challenging since it is not trivial to identify collective variables able to identify all metastable states along the reaction path. In this work, a systematic three-step protocol to uniquely identify the dominant order parameters for structural transformations in flexible metal–organic frameworks and subsequently construct accurate free energy profiles is presented. Methodologically, this protocol is rooted in the time-structure based independent component analysis (tICA), a well-established statistical modeling technique embedded in the Markov state model methodology and often employed to study protein folding, that allows for the identification of the slowest order parameters characterizing the structural transformation. To ensure an unbiased and systematic identification of these order parameters, the tICA decomposition is performed based on information from a prior replica exchange (RE) simulation, as this technique enhances the sampling along all degrees of freedom of the system simultaneously. From this simulation, the tICA procedure extracts the order parameters—often structural parameters—that characterize the slowest transformations in the material. Subsequently, these order parameters are adopted in traditional enhanced sampling methods such as umbrella sampling, thermodynamic integration, and variationally enhanced sampling to construct accurate free energy profiles capturing the flexibility in these nanoporous materials. In this work, the applicability of this tICA-RE protocol is demonstrated by determining the slowest order parameters in both MIL-53(Al) and CAU-13, which exhibit a strongly different type of flexibility. The obtained free energy profiles as a function of this extracted order parameter are furthermore compared to the profiles obtained when adopting less-suited collective variables, indicating the importance of systematically selecting the relevant order parameters to construct accurate free energy profiles for flexible metal–organic frameworks, which is in correspondence with experimental findings. The method succeeds in mapping the full free energy surface in terms of appropriate collective variables for MOFs exhibiting linker flexibility. For CAU-13, we show the decreased stability of the closed pore phase by systematically adding adsorbed xylene molecules in the framework.

Efficient Computation of Free Energy Surfaces of Diels⁻Alder Reactions in Explicit Solvent at Ab Initio QM/MM Level.

For Diels–Alder (DA) reactions in solution, an accurate and converged free energy (FE) surface at ab initio (ai) quantum mechanical/molecular mechanical (QM/MM) level is imperative for the understanding of reaction mechanism. However, this computation is still far too expensive. In a previous work, we proposed a new method termed MBAR+wTP, with which the computation of the ai FE profile can be accelerated by several orders of magnitude via a three-step procedure: (I) an umbrella sampling (US) using a semi-empirical (SE) QM/MM Hamiltonian is performed; (II) the FE profile is generated using the Multistate Bennett Acceptance Ratio (MBAR) analysis; and (III) a weighted Thermodynamic Perturbation (wTP) from the SE Hamiltonian to the ai Hamiltonian is performed to obtain the ai QM/MM FE profile using weight factors from the MBAR analysis. In this work, this method is extended to the calculations of two-dimensional FE surfaces of two Diels–Alder reactions of cyclopentadiene with either acrylonitrile or 1-4-naphthoquinone at ai QM/MM level. The accurate activation free energies at the ai QM/MM level, which are much closer to the experimental measurements than those calculated by other methods, indicate that this MBAR+wTP method can be applied in the studies of complex reactions in condensed phase with much-enhanced efficiency.

Recent Trends and Applications of Molecular Modeling in GPCR–Ligand Recognition and Structure-Based Drug Design

G protein-coupled receptors represent the largest family of human membrane proteins and are modulated by a variety of drugs and endogenous ligands. Molecular modeling techniques, especially enhanced sampling methods, have provided significant insight into the mechanism of GPCR–ligand recognition. Notably, the crucial role of the membrane in the ligand-receptor association process has earned much attention. Additionally, docking, together with more accurate free energy calculation methods, is playing an important role in the design of novel compounds targeting GPCRs. Here, we summarize the recent progress in the computational studies focusing on the above issues. In the future, with continuous improvement in both computational hardware and algorithms, molecular modeling would serve as an indispensable tool in a wider scope of the research concerning GPCR–ligand recognition as well as drug design targeting GPCRs.

Properties of the free energy density using the principle of maximum conformality**Supported by Natural Science Foundation of China (11625520)

We present a detailed study on the properties of the free energy density at high temperature by applying the principle of maximum conformality (PMC) scale-setting method within effective field theory. The PMC utilizes the renormalization group equation recursively to identify the occurrence and pattern of the non-conformal {βi}-terms, and determines the optimal renormalization scale at each order. Our analysis shows that a more accurate free energy density up to -order level without renormalization scale dependence can be achieved by applying the PMC. We also observe that by using a smaller factorization scale around the effective parameter mE, the PMC prediction is consistent with the lattice QCD prediction derived at low temperature.

Correction to Molecular Dynamics in Mixed Solvents Reveals Protein-Ligand Interactions, Improves Docking and Allows Accurate Binding Free Energy Predictions.

Correction to Molecular Dynamics in Mixed Solvents Reveals Protein-Ligand Interactions, Improves Docking and Allows Accurate Binding Free Energy Predictions.

Protein Folding Free Energy Landscape along the Committor - the Optimal Folding Coordinate.

Recent advances in simulation and experiment have led to dramatic increases in the quantity and complexity of produced data, which makes the development of automated analysis tools very important. A powerful approach to analyze dynamics contained in such data sets is to describe/approximate it by diffusion on a free energy landscape - free energy as a function of reaction coordinates (RC). For the description to be quantitatively accurate, RCs should be chosen in an optimal way. Recent theoretical results show that such an optimal RC exists; however, determining it for practical systems is a very difficult unsolved problem. Here we describe a solution to this problem. We describe an adaptive nonparametric approach to accurately determine the optimal RC (the committor) for an equilibrium trajectory of a realistic system. In contrast to alternative approaches, which require a functional form with many parameters to approximate an RC and thus extensive expertise with the system, the suggested approach is nonparametric and can approximate any RC with high accuracy without system specific information. To avoid overfitting for a realistically sampled system, the approach performs RC optimization in an adaptive manner by focusing optimization on less optimized spatiotemporal regions of the RC. The power of the approach is illustrated on a long equilibrium atomistic folding simulation of HP35 protein. We have determined the optimal folding RC - the committor, which was confirmed by passing a stringent committor validation test. It allowed us to determine a first quantitatively accurate protein folding free energy landscape. We have confirmed the recent theoretical results that diffusion on such a free energy profile can be used to compute exactly the equilibrium flux, the mean first passage times, and the mean transition path times between any two points on the profile. We have shown that the mean squared displacement along the optimal RC grows linear with time as for simple diffusion. The free energy profile allowed us to obtain a direct rigorous estimate of the pre-exponential factor for the folding dynamics.

No system-size anomalies in entropy of bcc iron at Earth\u2019s inner-core conditions

New molecular modeling data show that the entropy of bcc iron exhibits no system-size anomalies, implying that it should be feasible to compute accurate free energies of this system using first-principles methods without requiring a prohibitively large number of atoms. Conclusions are based on rigorous calculations of size-dependent free energies for a Sutton-Chen model of iron previously fit to ab initio calculations, and refute statements recently appearing in the literature indicating that the size of the simulation cell is critical for stabilization of the bcc phase.