Anharmonic Vibrational Frequency Calculations Are Not Worthwhile for Small Basis Sets

Jacobsen, R. L.; Johnson, R. D.; Irikura, K. K.; Kacker, R. N.  J. Chem. Theor. Comput. 2013, 9, 951
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Vibrational frequencies are routinely computed within most QM codes assuming the harmonic approximation. To correct for the neglect of higher order terms (anharmonicity), along with correcting for the inherent approximations of whatever quantum mechanical method is used, the harmonic frequencies are typically corrected by using a multiplicative scaling factor. The values of the scaling factor is method-dependent: a different scaling factor is need for every method and basis set combination! Nonetheless, this is a simple approach to computing often quite reasonable vibrational frequencies.

Anharmonic corrections can also be computed, and this is usually done using perturbation theory, which requires computing the third and often fourth derivatives, a mightily expensive proposition for reasonably large molecules even with DFT, let alone with some wavefunction-based post-HF method.

Jacobsen and co-workers¹ examined a set of 176 molecules that include 2738 vibrational modes, using HF, MP2, B3LYP and PBE0 with the 6-31G(d) or 6-31+G(d,p) basis sets. The unscaled anharmonic frequencies are much better than the unscaled harmonic frequencies; for example, using B3LYP/6-31+G(d), the root mean square deviation (RMSD) for the harmonic frequencies is 78 cm^-1 and 36 cm^-1for the anharmonic frequencies. But the scaled harmonic frequencies are just as good as the scaled anharmonic frequencies; using the same QM method, the RMSD for the scaled harmonic frequencies is 38 cm^-1 and 36 cm^-1 for the scaled anharmonic frequencies.

These authors suggest that accurate anharmonic corrections require very accurate potential energy surfaces, and so they recommend that unless you are using a very highly accurate computational model, there is no point in computing anharmonic frequencies; scaled harmonic frequencies will suffice!

References

(1) Jacobsen, R. L.; Johnson, R. D.; Irikura, K. K.; Kacker, R. N. "Anharmonic Vibrational Frequency Calculations Are Not Worthwhile for Small Basis Sets," J. Chem. Theor. Comput. 2013, 9, 951-954, DOI:10.1021/ct300293a.

This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported 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. 

Free Fructose is Conformationally Locked

Cocinero, E. J.; Lesarri, A.; Ecija, P.; Cimas, A.; Davis, B. G.; Basterretxea, F. J.; Fernandez, J. A.; Castano, F. J. Am. Chem. Soc. 2013, 135, 2845
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Sugars comprise a very important class of organic compounds for a variety of reasons, including dietary needs. On the chemical side, their stereochemical variations give rise to interesting
conformational questions. While sugar structures are a well-studied dating back to Fischer, most of these studies are in the solid or solution phase, and these phases can certainly play a role in dictating conformational preferences. This is seen in the differences in conformational distribution with different solvents. Only recently has instrumentation been developed (see these posts for some earlier applications: A, B, C) that can provide structural information of sugars in the gas phase. Cocinero and co-workers describe just such an analysis of fructose.¹

Using a combination of laser ablation Fourier transform microwave spectroscopy and quantum chemical computations, they have examined the gas-phase structure of this ketose. There are
quite a number of important conformational and configurational isomers to consider, as shown in Scheme 1. Fructose can exist in a pyranose form (6-member ring) with the anomeric carbon being α or β. An alternative cyclic form is the 5-member ring furanose form, which again has the two options at the anomeric position. Both the 5- and 6-member rings can ring flip, giving rise to 4 pyranose and 4 furanose forms. Of course there is also the
acyclic form.

Scheme 1. Major Frucotse isomers

α-pyranose

β-pyranose

α-furanose

β-furanose

Open chain

The observed microwave spectrum is quite simple, showing evidence of only a single isomer. In comparing the microwave rotational constants and the quartic centrifugal distortion constants with those obtained from MP2 and M06-2x computations, it is clear that the only observed isomer is the β-pyranose isomer in its ²C₅ conformation.

Both MP2 and M06-2x (with a variety of TZ basis sets) predict this isomer to be the lowest energy form by about 2.5 kcal mol^-1. This structure is shown in Figure 1. Interestingly, B3LYP predicts the open chain configuration as the most stable isomer, with the β-pyranose isomer about 0.5 kcal mol^-1 higher in energy. The authors strongly caution against using B3LYP for any sugars.

Figure 1. MP2/maug-cc-pVTZ optimized structure of β-fructopyranose.

This most stable furanose isomer displays five intramolecular hydrogen bonds that account for its stability over all other possibilities. However, the pyranose form of fructose is very rare in nature, and the Protein Data Bank has only four examples. The furanose form is by far the more commonly found isomer (as in sucrose). Clearly, hydrogen bonding to solvent and other solvent interactions alter the conformational distribution.

References

(1) Cocinero, E. J.; Lesarri, A.; Ecija, P.; Cimas, A.; Davis, B. G.; Basterretxea, F. J.; Fernandez, J. A.; Castano, F. "Free Fructose is Conformationally Locked," J. Am. Chem. Soc. 2013, 135, 2845-2852, DOI:10.1021/ja312393m.

InChIs:

β-fructopyranose:
InChI=1S/C5H10O6/c6-2-1-11-5(9,10)4(8)3(2)7/h2-4,6-10H,1H2
InChIKey=FFDHYUUPNCCTDA-UHFFFAOYSA-N



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

Are Protein Force Fields Getting Better? A Systematic Benchmark on 524 Diverse NMR Measurements

Kyle A. Beauchamp, Yu-Shan Lin, Rhiju Das, and Vijay S. Pande, J. Chem. Theory Comput. 2012 vol. 8 (4) pp. 1409-1414

Contributed by Ric Baron

Beauchamp and co-workers answer: YES! Proteins force fields are indeed getting better. After systematic analysis of 11 recent force field parameter sets, in combination with 5 popular fixed-charge water models and implicit solvent schemes, the authors conclude that two force fields in the pool (ff99sb-ildn- phi, ff99sb-ildn-NMR) that combine recent side chain and backbone torsion modifications achieved high accuracy. A simple metric to compare force field quality was employed and allows easy identificaiton of optimal force field/water model combinations (Figure 1). The benchmark set considered included dipeptides, tripeptides, tetra-alanine, and ubiquitin, and was validated against an extensive set of 524 NMR measurements. Expectedly, explicit solvent models outperform implicit solvent approaches.


Figure 1. Quality of different protein force fields when combined with different treatement of the solvent. This is an unofficial adaptation of Figure 1 that appeared in an ACS publication, J. Chem. Theory Comput. 2012 vol. 8 (4) pp. 1409-1414 . ACS has not endorsed the content of this adaptation or the context of its use

I believe this study is of particular interest for the biomolecular simulation community due to the extensive and rigorous force field comparison presented. It is a well-known challenge to define metrics that fairly capture force field quality. It is also problematic in practice to perform apples-to-apples force field comparison studies. Among the chief obstacles, few software engines enable simulations with alternative force fields using (i) identical treatment of bonded and non-bonded interactions among the force fields compared, and/or (ii) underlying algorithms and energy terms identical to those employed for force field parametrization. Beauchamp et al. employ an extensive benchmark set and reach a compromise following solution (i), i.e. by considering some of the force fields that can be evaluated by simulations in GROMACS.

An interesting implication of Beauchamp’s thorough force-field comparison study is that - for the two optimal force fields in the group considered (ff99sb-ildn- phi, ff99sb-ildn-NMR) – the calculation error is comparable to the uncertainty in the experimental comparison. This observation has remarkable implications for the comparison of molecular simulations with experiments. Therefore, while protein force fields steadily improve, extracting additional force field improvements from NMR data “may require in the future increased accuracy in J coupling and chemical shift prediction.”

At this stage, this work focused on peptide conformational sampling. As the authors note in their concluding remarks, “other experimental data, such as solvation free energies, may be critical for evaluating other aspects of force field performance.”

The Fold-In Approach to Bowl-Shaped Aromatic Compounds: Synthesis of Chrysaoroles

Myśliwiec, D.; Stępień, M. Angew. Chem. Int. Ed. 2013, 52, 1713
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Mysliwie and Stepian report on a new method for creating buckybowls.¹ The usual way had been to build from the inside outward. They opt instead to build from the outside in and have constructed the heterosubstitued bowl chrysaorole 1.

1

B3LYP/6-31G** optimizations reveal two conformers that are very close in energy: one has the butyl chains outstretched (1a) and one has the butyl arms internal or pendant (1b). These structures are shown in Figure 1. The depth of this bowl (1.96 Å) is quite a bit larger than in corranulene (0.87 Å). The agreement between the computed and experimental ¹³C and ¹H chemical shifts are excellent, supporting the notion that this gas phase geometry is similar to the solution phase structure. Though 1 is strained, 53.4 kcal mol^-1 based on B3LYP/6-31G** energies for Reaction 1 (which uses the parent of 1 – replacing the butyl groups with hydrogens), on a per sp² atom basis, it is no more strained than corranulene.

1a

1b

Figure 1. B3LYP/6-31G** optimized geometries of two conformers of 1.

Reaction 1

This new synthetic strategy is likely to afford access to more unusual aromatic structures.

References

(1) Myśliwiec, D.; Stępień, M. "The Fold-In Approach to Bowl-Shaped Aromatic Compounds: Synthesis of Chrysaoroles," Angew. Chem. Int. Ed. 2013, 52, 1713-1717, DOI:10.1002/anie.201208547.

InChI

1: InChI=1S/C54H45N3/c1-4-7-16-55-49-19-31-10-12-33-21-51-45-27-39(33)37(31)25-43(49)44-26-38-32(20-50(44)55)11-13-34-22-52-46(28-40(34)38)48-30-42-36(24-54(48)57(52)18-9-6-3)15-14-35-23-53(47(45)29-41(35)42)56(51)17-8-5-2/h10-15,19-30H,4-9,16-18H2,1-3H3
InChIKey=VUUJVETWVYQACL-UHFFFAOYSA-N

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

Musings about C2

Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

A short note here mainly to call to the reader’s attention a fascinating “trialogue” on the C₂ molecule.¹Shaik, Danovich, Wu, Su, Rzepa, and Hiberty² recently presented a full CI study of C₂ and concluded that the molecule contains a quadruple bond (see my previous post on this paper). This work was inspired in part by a blog post by Henry Rzepa.

The trialogue¹ is a conversation between Sason Shaik, Henry Rzepa and Roald Hoffmann about the nature of C₂, its 4^th bond, its diradical character, and some historical detours to see how some of our theoretical chemistry ancestors came close to proposing a quadruple bond. The discussion weaves together simple MO pictures, simple VB models, and the need for much more sophisticated analysis to ultimately approach the truth. Very much worth pointing out is the careful analysis of trying to tease out bond dissociation energies, especially analyzing the assumptions made here – including the possibility of errors in the experiments and not just errors in the computations! This is a very enjoyable read, following these three theoreticians as they traipse about the complex C₂ landscape!

References

(1) Shaik, S.; Rzepa, H. S.; Hoffmann, R. "One Molecule, Two Atoms, Three Views, Four Bonds?," Angew. Chem. Int. Ed. 2013, 52, 3020-3033, DOI: 10.1002/anie.201208206.

(2) Shaik, S.; Danovich, D.; Wu, W.; Su, P.; Rzepa, H. S.; Hiberty, P. C. "Quadruple bonding in C2 and analogous eight-valence electron species," Nat. Chem. 2012, 4, 195-200, DOI: 10.1038/nchem.1263.

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