Deflate to Understand Complex Molecular Kinetics

Brooke E. Husic & Frank Noé. arXiv: 1907.04101

Contributed by Jesper Madsen

Dimensionality reduction is at the core of understanding and making intuitive sense of complex dynamic phenomena in chemistry.  It is usually assumed that the slowest mode is the one of primary interest; however, it is critical to realize that this is not always so! A conceptual example hereof is a protein folding simulation (Lindorff-Larsen et al. Science 334, 517-520, 2011) where the slowest dynamical mode is not the folding itself (see Figure). What is the influence, then, of “non-slowest” modes in this process and how can it most appropriately be elucidated?

FIG: Figure 2 from the preprint: "(A) Sampled villin structures from the MD trajectory analyzed. Helical secondary structure is colored and coils are white. Each image represents five structures sampled from similar locations in TIC space as determined by a 250-center k-means model built upon the first three original TICs. The purple structure represents the folded state, and the blue structure represents the denatured state. The green structure is a rare helical misfolded state that we assert is an artifact. (B) Two-dimensional histograms for TICA transformations constructed from villin contact distances. Dashed lines indicate the regions corresponding to the sampled structures of the same color. The first TIC tracks the conversion to and from the rare artifact only. The second TIC tracks the majority of the folding process and correlates well with RMSD to the folded structure."

This work by Husic and Noé show how deflation can provide an answer to these questions. Technically speaking deflation is a collection of methods for how to modify a matrix after the largest eigenvalue is known in order to find the rest. In their provided example of the folding simulation, the dominant Time-lagged Independent Component (TIC) encapsulates the "artifact" variation that we are not really interested in. Thus, a constructed kinetic (Markov-state) model will be contaminated in several undesirable ways as discussed by the authors in great detail.  

In principle, this should be a very common problem since chemical systems have complex Hamiltonians. Perhaps the reason why we don’t see it discussed more is that ultra-rare events – real or artifact – may not usually be sampled during conventional simulations. So, with the increasing computational power available to us, and simulations approaching ever-longer timescales, this is likely something that we need to be able to handle. This preprint describes well how one can think about attacking these potential difficulties.   

Ab Initio Calculation of Rate Constants for Molecule–Surface Reactions with Chemical Accuracy

GiovanniMaria Piccini, Maristella Alessio, and Joachim Sauer (2016)

Contributed by Jan Jensen

Piccini et al. reproduce experimental rate constants for the reactions of methanol with ethene, propene, and trans-2-butene catalyzed by an acidic zeolite (H-MFI), to within one order of magnitude. Key to this is the inclusion of anharmonic effects using the method I highlighted earlier, but it should be noted that the reaction is biomolecular so entropy effects may be larger than for unimolecular reactions such as most enzyme catalysed reactions. However, anharmonic effects also changed the activation enthalpy by as much as 8 kJ/mol.

The PBE/plane wave electronic energy is corrected using MP2/CBS computed for a smaller systems plus a CCSD(T)/TZVP correction computed on an even smaller system.  Such corrections are becoming increasingly feasible for many problems and this study shows that the usual harmonic treatment of the vibrational free energy may become the limiting factor in terms of accuracy. However, anharmonic methods such as the one used here must be implemented, in a black box-fashion, in at least one of the major quantum chemistry packages before we'll see them widely applied.

See also this highlight of the article



This work is licensed under a Creative Commons Attribution 4.0

Enhancing NMR Prediction for Organic Compounds Using Molecular Dynamics

Kwan, E. E.; Liu, R. Y. J. Chem. Theor. Comput. 2015, 11, 5083-5089
Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

The prediction of NMR chemical shifts and coupling constants through ab initio computation is a major development of the past decade in computational organic chemistry. I have written about many developments on this blog. An oft-used method is a linear scaling of the computed chemical shifts to match those of some test set. Kwan and Liu wondered if the dynamics of molecular motions might be why we need this correction.¹

They suggest that the chemical shift can be computed as

<σ> = σ(static molecule using high level computation) + error

where the error is the obtained by using a low level computation taking the difference between the chemical shifts obtained on a dynamic molecule less that obtained with a static molecule. The dynamic system is obtained by performing molecular dynamics of the molecule, following 25 trajectories and sampling every eighth point.

They find outstanding agreement for the proton chemical shift of 12 simple molecules (mean error of 0.02 ppm) and the carbon chemical shift of 19 simple molecules (mean error of 0.5 ppm) without any scaling. Similar excellent agreement is found for a test set of natural products.

They finish up with a discussion of [18]annulene 1. The structure of 1 is controversial. X-ray crystallography indicates a near D_6h geometry, but the computed NMR shifts using a D_6h geometry are in dramatic disagreement with the experimental values, leading Schleyer to suggest a C₂  geometry. Kwan and Liu applied their dynamic NMR method to the D_6h, D_3h, and C₂ structures, and find the best agreement with the experimental chemical shifts are from the dynamic NMR initiated from the D_6h geometry. Dynamic effects thus make up for the gross error found with the static geometry, and now bring the experimental and computational data into accord.

One final note on this paper. The authors indicate that they have filed a provisional patent on their method. I am disturbed by this concept of patenting a computational methodology, especially in light of the fact that many other methods have been made available to the world without any legal restriction. For example, full details including scripts to apply Tantillo’s correction method are available through the Cheshire site and a web app to implement Goodman’s DP4 method are available for free. Provisional patents are not available for review from the US Patent Office so I cannot assess just what is being protected here. However, I believe that this action poses a real concern over the free and ready exchange of computational methodologies and ideas.

 

References

(1) Kwan, E. E.; Liu, R. Y. "Enhancing NMR Prediction for Organic Compounds Using Molecular Dynamics," J. Chem. Theor. Comput. 2015, 11, 5083-5089, DOI: 10.1021/acs.jctc.5b00856.

InChIs

1: InChI=1S/C18H18/c1-2-4-6-8-10-12-14-16-18-17-15-13-11-9-7-5-3-1/h1-18H/b2-1-,3-1+,4-2+,5-3+,6-4+,7-5-,8-6-,9-7+,10-8+,11-9+,12-10+,13-11-,14-12-,15-13+,16-14+,17-15+,18-16+,18-17-
InChIKey=STQWAGYDANTDNA-DWSNDWDZSA-N

[Editors note: this paper has also been highlighted by Jan Jensen]


This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Enhancing NMR Prediction for Organic Compounds Using Molecular Dynamics

Eugene E. Kwan and Richard Y. Liu (2015) DOI 10.1021/acs.jctc.5b00856

Contributed by Jan Jensen

This paper presents a novel way of computing vibrational effects on chemical shifts.  Existing methods generally compute such effects by displacing the coordinates along the normal modes to create PESs, computing the vibrational wavefunctions, computing chemical shifts for the displaced coordinates, and computing the chemical shifts as expectation values. However, it can be challenging to use this approach on low frequency modes since they require relatively large displacements which tend to distort the molecule in unphysical ways. This problem can, in part, be solved by using internal coordinates but results can wary greatly depending on which internal coordinates one chooses.  

Kwan and Liu propose instead to perform a short (125 fs) B3LYP/MIDI! quasi-classical MD simulation "along" each normal mode and then compute an average chemical shift based on the trajectory. The MD simulation is "quasi-classical" in the sense that it is initialized based in the vibrational harmonic oscillator energy levels. The energy level is randomly selected from a Boltzmann distribution and 25 trajectories were found to be sufficient. Since the forces are computed at each point and used to determine the next displacement point unphysical distortions of the molecule are avoided.  Another difference is that the chemical shifts are computed using a higher level of theory (B3LYP/cc-pVDZ) than that used to construct the PES.

The various choices for method, basis set, simulation length, number of trajectories, etc is tested extensively in the supplementary materials, which also contains a more thorough description of steps in the algorithm (page 96).  Here the authors also state that they "are in the process of developing a user-friendly package for carrying out these calculations, which will be reported in due course."

Compared to the previously published methods for computing vibrational corrections this method is significantly more expensive in terms of energy and gradient evaluations.  However, I wonder if the methods can be combined such that this method is used only for lower frequency modes that can problematic for the "displacement" methods.  I also wonder if the method can be adapted to compute the anharmonic corrections to the enthalpy and entropy, which also can prove challenging for displacement-based methods.

This work is licensed under a Creative Commons Attribution 4.0

Effect of Anharmonicity on Adsorption Thermodynamics

GiovanniMaria Piccini and Joachim Sauer Journal of Chemical Theory and Computation 2014, 10, 2476.

Contributed by +Jan Jensen 

Low frequency modes can make a sizable contribution to the binding free energy computed using the harmonic oscillator approximation (HO), but is even appropriate to use the HO-expression for the free energy for these modes.  In a paper I highlighted previously, Kjærgaard and co-workers showed that the HO gave good agreement with the experimental binding free energy for acetonitrile-HCl dimer even though the lowest frequency mode was a mere 38 cm$^{-1}$.

However, in this study Piccini and Sauer show that anharmonic effects need to be included to obtain agreement with the experimental binding free energy of methane to H-CHA zeolite.  Specifically, Piccini and Sauer compute the vibrational binding free energy by computing the 1-dimensional (1D) PES for each low frequency mode and compute the vibrational energy levels and corresponding partition function numerically using the method proposed by Beste (as opposed to using the anharmonic fundamental frequency together with the HO partition function).  This decreases the binding free energy by 2.5 kcal/mol compared to the standard HO treatment.  

Actually, Piccini and Sauer show that the usual double difference displacement in Cartesian coordinates can result in errors for low frequency modes and if harmonic frequencies are extracted from the 1D PESs the difference between harmonic and anharmonic binding free energies increase to 3.0 kcal/mol.  (Of course if one has already computed the 1D PESs one might as well compute the anharmonic free energy contribution).

As Piccini and Sauer point out it is crucial to use internal coordinates when displacing the molecule along the normal mode to minimize unphysical distortions in bond lengths and bond angles. The choice of internal coordinates can affect the results and must be chosen carefully.

One note on the general methodology (described in more detail here): I believe the normal mode optimization reduces to a standard RFO quasi Newton-Raphson optimization if all modes are optimized, but using a transformed gradient rather than the Cartesian gradient to check convergence. So if one is using such an optimizer to begin with (rather than e.g. conjugate gradient) implementing the normal mode optimization is not necessary or requires merely a change in convergence criteria.    


This work is licensed under a Creative Commons Attribution 4.0 International License.

Benchmarking Ab Initio Binding Energies of Hydrogen-Bonded Molecular Clusters Based on FTIR Spectroscopy

Nicolai Bork, Lin Du, Heidi Reiman, Theo Kurten, and Henrik G. Kjaergaard Journal of Physical Chemistry A 2014, 118, 5316-5322
Contributed by +Jan Jensen

One of the holy grails of computational chemistry is a method that consistently can compute standard binding free energies in bulk solvent to within 1 kcal/mol.  When using electronic structure theory one of the (many) potential problems is the use of the harmonic approximation because low-frequency modes can make a sizable contribution to the vibrational entropy. This paper provides some valuable benchmarking data.

Kjærgaard and co-workers measure the gas phase binding free energy at 295 K of acetonitrile-HCl to be between 1.2 and 1.9 kcal/mol. CCSD(T)/aug-cc-pV(T+d) predicts 2.3 kcal/mol using the harmonic frequencies and 1.9 kcal/mol using the anharmonic scaling factors suggested by Shields and co-workers.  So, at least in this case, the harmonic approximation works quite well, but can be improved quite easily using anharmonic scaling factors.  If this effect is additive and transferable, then making such a correction will be important for systems that bind using several hydrogen bonds of the 1 kcal/mol accuracy-target is to be met.

The authors write "We find that this energy lowering originates almost entirely from the scaling of the two lowest frequency modes (between 28 and 35 cm$^{-1}$, depending on method), emphasizing the importance of calculating the low frequency modes accurately."

Indeed, such low frequency modes are susceptible to numerical noise (e.g. from coarse DFT grids or sloppy optimization convergence criteria) which can even lead to them turning up as imaginary frequencies in the vibrational analysis. Because imaginary frequencies are not included in the vibrational free energy, this can lead to sizable errors in the binding free energy.  For example, a 28 cm$^{-1}$ frequency contributes 1.8 kcal/mol to the free energy at 295 K.

B3LYP-D3, MP2 and PW91/aug-cc-pV(T+d) calculations also yield binding free energies within the experimentally measured range, but in the last two cases this is due for fortunate cancellation of errors in the electronic and free energy components based on comparison to the CCSD(T) calculations.  Having said that, with the exception of B3LYP, all functionals tested in this study are able to predict the binding energy to within 1 kcal/mol using the aug-cc-pV(T+d) basis set.

It is worth noting that the approach by Shields and co-workers find scaling factors by fitting to harmonic vibrational energy and entropy-expressions evaluated using anharmonic fundamental frequencies.  It is not at all obvious (at least to me) that this approach will result in good agreement with experimental binding free energies involving 28 cm$^{-1}$ frequencies, so this (and a previous) study provides very valuable validation of this general approach.



This work is licensed under a Creative Commons Attribution 4.0 International License.

Calculation of Molecular Entropies Using Temperature Integration

Kim Sharp Journal of Chemical Theory and Computations 2013, 9, 1164
Contributed by +Jan Jensen

Extracting accurate entropy changes from molecular dynamics (MD) simulations is notoriously difficult.  The entropy is usually computed using the Boltzmann-Gibbs-Planck expression$$S=-k \int p(\mathbf{q}) \ln p(\mathbf{q}) d\mathbf{q} $$where $\mathbf{q}$ is the postion vector and $p(\mathbf{q})$ the corresponding probability distribution function, which converges slowly for large systems.  In this paper Kim Sharp makes a convincing argument for using the Clausius-van't Hoff (CvH) expression$$S_{CvH}=\frac{\Delta E}{T_{CvH}}=\frac{\Delta E \left( \ln T_2 - \ln T_1\right)}{T_2-T_1}$$There are two issues which must be considered when using this equation.  One issue is that it assumes the heat capacity is constant between $T_1$ and $T_2$ placing a limit on how large the difference in $T$ can be.  The other issue is that the entropy does not converge to 0 as $T$ approaches 0 for simulations using classical energy functions such as the harmonic oscillator.

To address the latter problem Sharp argues that at sufficiently low $T$ classical systems will behave harmonically so that the entropy can be estimated by either the quasiharmonic method or using normal mode frequencies.  

The approach is tested in two ways. CvH conformational entropy changes on going from 280 K to 320 K for five ligands are shown to be close to corresponding values computed using probability density functions, while converging an order of magnitude faster. Furthermore, water vaporization  entropies at 300 K computed using the CvH approach were within 1 cal/mol/K of experiment and conventional Monte Carlo simulations.  This result can be obtained by cooling in as little as three steps to about 50 K and applying the quasiharmonic method for lower temperatures.

In addition to a more rapid convergence compared to existing methods the method also obviates any need to define suitable coordinates for the probability density functions, which can present a problem for some applications such as QM/MM based MD.  Further, MD simulations are often started by heating from an energy minimized structure and this procedure can be conveniently used to extract the entropy of the system.

Thanks to +Casper Steinmann for alerting me to this paper

 
This work is licensed under a Creative Commons Attribution 3.0 Unported License. 

Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory

Stefan Grimme Chemistry A European Journal 2012, 18, 9955 (Paywall)
Contributed by +Jan Jensen

Computing accurate binding free energies is a fundamental challenge to molecular modeling.  I recently highlighted a study by +Hari Muddana and +Michael Gilson where the binding free energies of 29 ligands to the host CB7 were reproduced with an RMSE of 1.9 kcal/mol using  PM6-DH+, though this involved fitting of the solvation free energy change.

In this study Grimme presents a computational methodology that reproduces binding free energies of 13 chemically diverse host-guest complexes with a MAD of 2.1 kcal/mol.  Though the accuracies of the two studies are similar there are many differences between the two approaches.

The gas phase interaction energies are computed using dispersion corrected DFT calculations.  Practically, this means that, unlike the study by Muddana and Gilson, no conformational averaging is performed.  The dispersion correction includes a three-body term that is shown to contribute up to 4.6 kcal/mol to the gas phase interaction energy.  The similar three-body terms also contribute to the vibrational free energy contributions (see below) and together with the gas phase contributions lower the MAD by 1.3 kcal/mol compared to experimental binding free energies.

Similarly, the solvation free energy is computed using DFT and COSMO and, unlike the Muddana and Gilson paper, is not corrected empirically.  This may be due to the fact that the DFT/COSMO calculations (including the radii used to define the surface) already has been carefully calibrated while the PM6/COSMO calculations have not.

The vibrational free energy corrections are computed in the gas phase using dispersion and hydrogen-bond corrected PM6 and DFTB (PM6-D3H and DFTB-D3H, respectively).  These corrections were developed as part of this study based on the S66 data set and implemented in an external program interfaced with the MOPAC and DFTB codes (Grimme, personal communications).

The vibrational entropy contribution increases exponentially with decreasing vibrational frequency and this can introduce numerical unstable binding free energies for these systems as low vibrational frequencies are common.  Grimme deals with this problem by introducing a hybrid function where the vibrational entropy is smoothly replaced (using a switching function) with a corresponding free-rotor entropy for low frequency modes.  As computed in his study, the free-rotor entropy increases much less steeply as the vibrational frequency decreases and produced binding free energies that are in better agreement with experiment. However, computing numerically stable vibrational frequencies for low frequency modes remains a challenge for these systems and may require manual intervention (Grimme, personal communications).  

Thanks to +Hari Muddana for alerting me to this paper and to Stefan Grimme for personal communications.

 
This work is licensed under a Creative Commons Attribution 3.0 Unported License.