Quantum Monte Carlo Benchmark of Exchange-Correlation Functionals for Bulk Water

Miguel A. Morales, John R. Gergely, Jeremy McMinis, Jeffrey M. McMahon, Jeongnim Kim, and David M. Ceperley Journal of Chemical Theory and Computation 2014,  10, 2355–2362

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

I want to start by noting that the modelling of water is enormously complex and controversial, and that I’m far from an expert; in water, hydrogen bonding, dispersion interactions and quantum nuclear effects are all important.  The paper I want to discuss[1] performs detailed comparisons of the energies of different configurations of water molecules (drawn from a variety of dynamic simulations) calculated with several DFT functionals and quantum Monte Carlo (QMC).  While there have been steps in the direction of doing MD with QMC, it is still too expensive to do this routinely; static QMC calculations have been shown to model ice extremely accurately[2,3].

The configurations that are tested are drawn from classical forcefields and path-integral MD (PIMD) with van der Waals corrected DFT (see here for a recent overview of approaches to vdW in DFT).  The only test that is applied is the difference in total energy between the methods.  This is, probably, the only test that could be applied, as force calculations in QMC are very expensive, but it raises an interesting question: is it the convergence of total energy always the best thing to test ? It’s certainly true in normal DFT calculations that forces and energy differences converge faster than absolute energies, and many calculations use this fact to reduce the computational overhead.

The classical forcefields used rigid water molecules (a standard approach to speed up the simulations) and for these configurations, the DFT errors are around 0.1 mHa (3 meV) per water molecule, with only small variations between functionals.  This goes some way to explaining the success that can be found for some areas of water modelling with DFT (it’s important to note that dispersion effects are important and do reduce the errors).

When the configurations from PIMD are tested, however, these errors become much larger – around 0.3 mHa (10 meV) per water molecule for standard and vdW DFT methods.  Only the hybrid methods maintain a low error (but they are still too expensive to use routinely with MD).  Why does this happen ? Well, it appears that the major sources of error between DFT and QMC comes from the structure of the water molecule itself; when it is allowed to be flexible, configurations are generated that DFT describes poorly, but hybrids model better.

This failure can be explored further, using a correction scheme for modelling water in DFT[4] (where the error between DFT and QMC or coupled cluster calculations is decomposed in a many body expansion).  They find that the major term contributing to the error is the one-body term, and show that the correction scheme gives excellent results.  This confirms the difference between the errors for rigid and flexible water molecules.

The paper gives us an excellent insight into the limitations of DFT, and how it can be improved.  It’s important to note that there are different ways of testing: PIMD simulations with vdW-DF2 have been shown to reproduce the experimental water structure almost perfectly[5].  It’s easy to get trapped into thinking that one particular approach or technique is best, or to become obsessed with methods: but experiment is always what must be tested, though it’s often far from trivial to work out how to compare to experiments.

[1] J. Chem. Theory Comput. 10, 2355 (2014) DOI:10.1021/ct500129p

[2] Phys. Rev. Lett. 107, 185701 (2011) DOI:10.1103/PhysRevLett.107.185701

[3] J. Chem. Phys. 138, 221102 (2013) DOI:10.1063/1.4810882

[4] J. Chem. Phys. 139, 244504 (2013) DOI:10.1063/1.4852182

[5] arXiv:1402.2697 http://arxiv.org/abs/1402.2697

INDO/X: A New Semiempirical Method for Excited States of Organic and Biological Molecules

Alexander A. Voityuk Journal of Chemical Theory and Computation 2014, ASAP
Contributed by +Jan Jensen

This paper presents a significantly more accurate re-parameterization of the popular INDO/S method. The re-parameterization is based on the TBE-2 data set: vertical excitation energies and oscillator strengths computed for the valence excited states of 28 medium-sized organic molecules using multistate-CASPT2 and aug-cc-pVTZ basis set.

The author sums it up succinctly: "The mean absolute deviations of the INDO/X excitation energies relative to the TBE-2 data is 0.26 eV for singlet states and 0.33 eV for triplet states. The corresponding values for INDO/S are 0.56 and 0.64 eV. ...  Unlike INDO/S where two different Hamiltonians are employed for singlet and triplet excitations, both types of excited states are calculated within the same INDO/X scheme."  The MADs are thus similar to the MAD for TD-B3LYP calculations (0.27 eV) and significantly better than TD-BP86 (0.53 eV) and TD-BHLYP (0.48 eV) calculations.

Like INDO/S, INDO/X uses the CIS approach and can therefore be applied to quite large molecules. Furthermore, the author notes that the results are quite insensitive to the chosen active space.

Thanks to @CompChemNews for alerting me to this paper.



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

Tunneling Assists the 1,2-Hydrogen Shift in N-Heterocyclic Carbenes

Karmakar, S.; Datta, A. Angew. Chem. Int. Ed. 2014, 53, 9587-9591
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Organic chemists are beginning to recognize that tunneling may be more pervasive than previously thought. This blog has noted a number of interesting occurrences of tunneling, and here’s one more, by Karmakar and Datta.¹

The barrier for the intramolecular earrangement (Reaction 1) taking the carbene 1 into 2 is estimated to be 44.1 kcal mol^-1 at M06-2X/6-31+G(d,p), prohibitively large. However, the intermolecular rearrangement (Reaction 2) has a much smaller barrier of 11.4 kcal mol^-1. The structures of the transition states for these two reactions are shown in Figure 1.

TSintra

TSinter

Figure 1. M06-2X/6-31+G(d,p) optimized transition states for Reactions 1 and 2.

Given that the barrier width is likely to be very small for the intramolecular route, perhaps tunneling may play a role. The rate predicted with canonical variational transition-state theory (CVT) and small curvature tunneling (SCT) at 298K is negligible. However, for the intermolecular process, the rate at 298K including tunneling is 3.56 x 10⁴ s^-1, more than 10 times great than predicted with CVT alone, and tunneling makes a dramatically larger difference at lower temperatures.

The intermolecular barrier for the rearrangement of 3 into 4 is very small, only 1.6 kcal mol^-1.
This manifests in a very interesting rate prediction: the reaction is actually predicted to be slower at temperatures above 150K when tunneling is included than when tunneling is omitted. This is a result of quantum mechanical reflection off of the barrier, and this becomes noticeable with the very small barrier. In addition, the kinetic isotope effects are smaller than expected when D is substituted in for H. These predictions await experimental confirmation.

References

(1) Karmakar, S.; Datta, A. "Tunneling Assists the 1,2-Hydrogen Shift in N-Heterocyclic Carbenes," Angew. Chem. Int. Ed. 2014, 53, 9587-9591, DOI: 10.1002/anie.201404368.

InChIs:

1: InChI=1S/C3H6N2/c1-2-5-3-4-1/h4-5H,1-2H2
InChIKey=JKQUEGZDRZXJNY-UHFFFAOYSA-N

2: InChI=1S/C3H6N2/c1-2-5-3-4-1/h3H,1-2H2,(H,4,5)
InChIKey=MTNDZQHUAFNZQY-UHFFFAOYSA-N

3: InChI=1S/C3H2F2N2/c4-2-3(5)7-1-6-2/h6-7H
InChIKey=LHUPDFSUHVZFPD-UHFFFAOYSA-N

4: InChI=1S/C3H2F2N2/c4-2-3(5)7-1-6-2/h1H,(H,6,7)
InChIKey=KXXZDIFMEWOLPE-UHFFFAOYSA-N

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

The Negative Ion Photoelectron Spectrum of Cyclopropane-1,2,3-Trione Radical Anion, (CO)3•– — A Joint Experimental and Computational Study

Chen, B.; Hrovat, D. A.; West, R.; Deng, S. H. M.; Wang, X.-B.; Borden, W. T. J. Am. Chem. Soc. 2014, 136, 12345-12354
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

As recently explicated by Wang and Borden using NIPE spectroscopy and computations, the potential energy surface of cyclopropyl-1,2,3-trione 1 is remarkably complex.¹ (U)CCSD(T)//aug-cc-pVTZ computations of the D_3h singlet (the ¹A_1’ state shown in Figure 1) is actually a hilltop, possessing two imaginary frequencies. Distorting the structure as indicated by these imaginary frequencies and then optimizing the structure leads directly to dissociation to three CO molecules. Thus, (CO)₃ does not exist as a stable minima on the singlet surface.

The D_3h triplet (the ³E” state shown in Figure 1) is not a critical point on the surface; due to the Jahn-Teller effect is distorts into two different states: the ³B₁ state which is a local energy minimum, and the ³A₂ state which is a transition state between the symmetry-related ³B₁ states.

So, this implies the possibility of a very interesting NIPE experiment. If the radical anion (CO)₃^-.
loses an electron and goes to the singlet surface, it lands at a hilltop(!) and should have a very short lifetime. If it goes to the triplet surface, it lands at either a transition state (³A₂) and again should have a short lifetime, or it can land at the ³B₁ state and perhaps have some lifetime before it dissociates by losing one CO molecule.

1-.
1A1’	3E”
3B1	3A2

Figure 1. (U)CCSD(T)//aug-cc-pVTZ optimized geometries of 1 and its radical anion.

The NIPE spectrum identifies three transitions. By comparing the energies of the electron loss seen in the experiment with the computations, along with calculating the Franck-Condon factors using the computed geometries and vibrational frequencies, the lowest energy transition is to the ¹A_1’ state, and the second transition is part of the vibrational progression also to the ¹A_1’ state. This is the first identification of vibrational frequencies associated with a hilltop structure. The third transition is to the ³A₂ state. No transition to the ³B₁ state is found due to the large geometric difference between the radical anion and the³B₁ state; the Franck-Condon factors are zero due to no overlap of their wavefunctions.

Once again, the power of the symbiotic relationship between experiment and computation is amply demonstrated in this paper.

References

(1) Chen, B.; Hrovat, D. A.; West, R.; Deng, S. H. M.; Wang, X.-B.; Borden, W. T. "The Negative Ion Photoelectron Spectrum of Cyclopropane-1,2,3-Trione Radical Anion, (CO)₃^•– — A Joint Experimental and Computational Study," J. Am. Chem. Soc. 2014, 136, 12345-12354, DOI: 10.1021/ja505582k.

InChIs

1: InChI=1S/C3O3/c4-1-2(5)3(1)6
InChIKey=RONYDRNIQQTADL-UHFFFAOYSA-N

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

Design and synthesis of the first triply twisted Möbius annulene

Schaller, G. R.; Topić, F.; Rissanen, K.; Okamoto, Y.; Shen, J.; Herges, R.  Nat. Chem. 2014, 6, 608-613
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Herges and co-workers have prepared a triply-twisted Möbius molecule. ¹ The key element is recognizing that most of the “twisting” needs to be accomplished through writhe, a twisting that produces figure-8-like crossing, the way an old-school phone cord twists about itself or the way a pretzel is formed. Herges employs three bi-naphthelene subunits to provide the template for the writhe needed. The prepared compound is 1. A clever, relatively straightforword synthesis produces this amazing molecule, along with the single-twisted 2.

The B3LYP/6-31G* optimized geometries of 1 and the single-twisted analogue 2 are shown in Figure 1. Table 1 presents the key topological parameters of 1 and 2, comparing the computed and X-ray crystal structure. The absolute value of the linking number L_k is 3, indicating the three twists, and the reason that this highly twisted molecule can be made is that half of the twist actually results from writhe.

1

2

Figure 1. B3LYP/6-31G* optimized geometries of the two diastereomers if 1. (Be sure to click on these images to launch JMol and interactively manipulate the structures!)

Table 1. Topological parameters of 1. ^a

	|Lk|	|Tw|	|Wr|
X-ray, 1	3.0	1.42	1.58
Comp, 1	3.0	1.33	1.67
Comp, 2	1.0	0.37	0.63

^aL_k is the linking number, T_w is the twist number, and W_r is the writhe number, with the condition that T_k+ W_r = L_k.

References

(1) Schaller, G. R.; Topić, F.; Rissanen, K.; Okamoto, Y.; Shen, J.; Herges, R. "Design and synthesis of the first triply twisted Möbius annulene," Nat. Chem. 2014, 6, 608-613, DOI: 10.1038/nchem.1955.

InChIs

1: InChI=1S/C72H36/c1-2-26-56-44-38-50-20-8-14-32-62(50)68(56)70-58(46-40-52-22-10-16-34-64(52)70)28-5-6-30-60-48-42-54-24-12-18-36-66(54)72(60)71-59(47-41-53-23-11-17-35-65(53)71)29-4-3-27-57-45-39-51-21-9-15-33-63(51)69(57)67-55(25-1)43-37-49-19-7-13-31-61(49)67/h7-24,31-48H/b69-67-,70-68-,72-71-
InChIkey=ICOGWBLQQNBWJA-UWFYZODKSA-N

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

Competition Between Concerted and Stepwise Dynamics in the Triplet Di-π-Methane Rearrangement

Jiménez-Osés, G.; Liu, P.; Matute, R. A.; Houk, K. N. Angew. Chem. Int. Ed. 2014,53, 8664-8667
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

The di-π-methane photorearrangement has been known for many years, first studied by Zimmerman.^1,2The triplet photorearrangement gives an interesting rearranged product; and the mechanism of this photorearrangement of 1 into 2 has now been examined by the Houk group using computational techniques, including trajectory analysis. The proposed mechanism is that excitation to the triplet state1* is followed by rearrangement to the triplet intermediate INT1* which then rearranges to the tripletINT2*. Intersystem crossing then leads to the singlet product 2.

The PES for this rearrangement was explored³ at CASMP2(10,10)/6-31G(d)//CASSCF(10,10)/6-31G(d), with geometries and relative energies shown in Figure 1, as well as at (U)M06-2x/6-31G(d) and (U)B3LYP/6-31G(d); they all give qualitatively the same result. The first TS is the rate limiting step, and the second TS lies only 1-2 kcal mol^-1 above the intermediate INT1. So, the reaction appears to be two steps, but with such a low barrier for the second step, dynamic effects might be important as trajectories might cross INT1* and go over TS2* without residing in the intermediate well for any appreciable time – a seemingly one step reaction. Note than no TS for directly traversing from 1* to INT2* was found.

1* 0.0
TS1* 12.9	INT1* 7.1
TS2* 8.2	INT2* -15.4

Figure 1. CASSCF(10,10)/6-31G(d) geometries and CASMP2 energies in kcal mol^-1.

Now in a follow-up study, Houk and co-workers⁴ performed trajectories analysis on the M06-2x/6-31G(d) PES. A total of 256 trajectories were initiated at TS1* and 241 ended at INT2* within 1500fs. Of these, 24 trajectories resided for less than 60fs within the region of INT1, a time less than a C-C vibration. Furthermore, the lifetime of INT1 that is predicted by RRKM is much longer (about 500fs) than what is observed in the trajectories (about 200 fs). Thus, there is significant dynamic effects in this excited state rearrangement, though INT1 is always sampled.

References

(1) Zimmerman, H. E.; Grunewald, G. L. "The Chemistry of Barrelene. III. A Unique Photoisomerization to Semibullvalene," J. Am. Chem. Soc. 1966, 88, 183-184, DOI: 10.1021/ja00953a045.

(2) Zimmerman, H. E.; Binkley, R. W.; Givens, R. S.; Sherwin, M. A. "Mechanistic organic photochemistry. XXIV. The mechanism of the conversion of barrelene to semibullvalene. A general photochemical process," J. Am. Chem. Soc. 1967, 89, 3932-3933, DOI: 10.1021/ja00991a064.

(3)  Matute, R. A.; Houk, K. N. "The Triplet Surface of the Zimmerman Di-π-Methane Rearrangement of Dibenzobarrelene," Angew. Chem. Int. Ed. 2012, 51, 13097-13100, DOI: 10.1002/anie.201208002.

(4) Jiménez-Osés, G.; Liu, P.; Matute, R. A.; Houk, K. N. "Competition Between Concerted and Stepwise Dynamics in the Triplet Di-π-Methane Rearrangement," Angew. Chem. Int. Ed. 2014,
53, 8664-8667, DOI: 10.1002/anie.201310237.

InChIs

1: InChI=1S/C16H12/c1-2-6-12-11(5-1)15-9-10-16(12)14-8-4-3-7-13(14)15/h1-10,15-16H
InChIKey=VWDKVBGOVYWYFZ-UHFFFAOYSA-N

2: InChI=1S/C16H12/c1-3-7-11-9(5-1)13-10-6-2-4-8-12(10)15-14(11)16(13)15/h1-8,13-16H
InChIKey=RATAQXOLJVRERC-UHFFFAOYSA-N

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