Efficient prediction of reaction paths through molecular graph and reaction network analysis

Yeonjoon Kim, Jin Woo Kim, Zeehyo Kim and Woo Youn Kim (2017)

Highlighted by Jan Jensen

Figure 2 from the paper. Reproduced under the CC-BY licence.

The paper presents an interesting approach to generating the most important intermediates given reactants and products. 

1. The reactants and products are mapped to atom connectivity (AC) matrices, where the presence and absence of bonds are represented by 1's and 0's. 

2. Intermediates are then enumerated by changing the matrix elements using simple rules and constraints are employed to keep the number of intermediates manageable: the change in bonding is limited to a maximum of two per steps and only changes to bonds involving "active atoms" (atoms with different bonding in reactants and products) are considered. Identifying active atoms involves mapping reactant atoms to product atoms, using the method proposed by Floudas and co-workers.  The AC matrices are converted to SMILES strings and 3D coordinates to ensure that they correspond to chemically reasonable molecules.

3. A reaction network is then constructed by finding intermediates that best connect reactants and products by computing "chemical distances (CDs)" between AC matrices and identifying intermediates with CDs to reactant and product that are similar to the CD between the reactant and product themselves. The CDs are then used to find the shortest path to reactant and product for each intermediate to yields a minimal reaction network.

4. Finally, TSs connecting the structures in the minimal reaction are found using conventional QM methods.

The method is applied to Claisen ester condensation and Cobalt-catalyzed hydroformylation. The method is implemented in a python program called ACE-Reaction but the availability of this program is unclear.  

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

A Zwitterionic, 10 π Aromatic Hemisphere

Hafezi, N.; Shewa, W. T.; Fettinger, J. C.; Mascal, M., Angew. Chem. Int. Ed. 2017, 56, 14141-14144
Contributed by Steven Bacharach

Reposted from Computational Organic Chemistry with permission

The range of aromatic compounds seems limitless. Mascal and co-workers have prepared the azatriquinacene 1 in a remarkably simple fashion.¹ The molecule is a zwitterion, with the carbon atoms forming a 9-center, but 10 π-electron ring, and the quaternary nitrogen sitting above it. The carbon ring satisfies Hückel’s rule (4n+2) and so should be aromatic. The capping nitrogen should help to keep the carbon ring fixed in a shallow bowl.

As expected, the molecule in fact turns out to possess an aromatic 10 π-electron ring. The B3LYP/6-311++G(d,p) geometry is shown in Figure 1. There is little bond alternation among the C-C distances: the mean deviation is only 0.015 Å with the largest difference only 0.024 &Aring. The x-ray crystal structure shows the same trends. The NICS(1) value is -12.31 ppm, larger even than that of benzene (-10.22 ppm).

Figure 1. B3LYP/6-311++G(d,p) geometry of 1.

References

1) Hafezi, N.; Shewa, W. T.; Fettinger, J. C.; Mascal, M., "A Zwitterionic, 10 π Aromatic Hemisphere." Angew. Chem. Int. Ed. 2017, 56, 14141-14144, DOI: 10.1002/anie.201708521.

InChIs

1: InChI=1S/C10H9N/c1-11-8-2-3-9(11)6-7-10(11)5-4-8/h2-7H,1H3
InChIKey=ZXZPLDVSQUVKTH-UHFFFAOYSA-N

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

Heavy-Atom Tunneling Calculations in Thirteen Organic Reactions: Tunneling Contributions are Substantial, and Bell’s Formula Closely Approximates Multidimensional Tunneling at ≥250 K

Doubleday, C.; Armas, R.; Walker, D.; Cosgriff, C. V.; Greer, E. M., Angew. Chem. Int. Ed. 2017, 56, 13099-13102
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Though recognized to occur in organic systems, the breadth of involvement of heavy-atom tunneling has not been established. Doubleday, Greer and coworkers have examined 13 simple organic reactions sampling pericyclic reactions, radical rearrangements and S_N2 reactions for heavy-atom tunneling.¹ A few of these reactions are shown below.

Reaction rates were obtained using the small curvature tunneling approximation (SCT), computed using Gaussrate. Reaction surfaces were computed at B3LYP/6-31G*. The tunneling correction to the rate was also estimated using the model developed by Bell: k_Bell = (u/2)/sin(u/2) where u = hν^‡/RT and ν^‡ is the imaginary frequency associated with the transition state. The temperature was chosen so as to give a common rate constant of 3 x 10^-5 s^-1. Interestingly, all of the examined reactions exhibited significant tunneling even at temperatures from 270-350 K (See Table 1). The tunneling effect estimated by Bell’s equation is very similar to that of the more computationally demanding SCT computation.

Table 1. Tunneling contribution to the rate constant

Reaction	% tunneling
	35
	17
	28
	95
CN– + CH3Cl → CH3CN + Cl– (aqueous)	45

This study points towards a much broader range of reactions that may be subject to quantum mechanical tunneling than previously considered.

References

1. Doubleday, C.; Armas, R.; Walker, D.; Cosgriff, C. V.; Greer, E. M., "Heavy-Atom Tunneling Calculations in Thirteen Organic Reactions: Tunneling Contributions are Substantial, and Bell’s Formula Closely Approximates Multidimensional Tunneling at ≥250 K." Angew. Chem. Int. Ed. 2017, 56, 13099-13102, DOI: 10.1002/anie.201708489.

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

Understanding and Breaking Scaling Relations in Single-Site Catalysis: Methane-to-methanol Conversion by Fe(IV)=O

Terry Z. H. Gani and Heather J. Kulik (2017)

Highlighted by Jan Jensen

This is the first study I have come across that locates TS structures as part of a "high-throughput" single-site catalyst design study. Furthermore, the catalyst contains iron, which is not the easiest of elements to work with computationally. 

The study locates 76 and 43 TSs for the oxo formation (TS1) and hydrogen atom transfer (HAT, TS2) steps of the catalytic cycle. These are relatively small numbers compared to high throughput studies of other properties (hence the quotation marks), but they are roughly an order of magnitude larger than the number of TSs found in typical computational study of catalysts. The number is smaller for HAT due to difficulties in locating TSs for this step.

The TSs were located using either NEB implemented in DL-FIND or Q-CHEM where initial guess structures were generated using a locally modified version of molSimplify.

The studies show that there is a good correlation between reaction energy and barrier for the HAT step (R2 = 0.99) but a poor correlation for the oxo formation (R2 = 0.50 - 0.81). The authors conclude "Overall, our work shows that LFERs can be leveraged in single-site catalyst screening only when the coordination geometry is held fixed. Reliance solely on LFERs for single-site catalysis will thus miss rich areas of chemical space accessible through scaffold distortion."

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

Tunneling Control of Chemical Reactions: The Third Reactivity Paradigm

Schreiner, P. R., J. Am. Chem. Soc. 2017, 139, 15276-15283
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Over the past nine years the Schreiner group, often in collaboration with the Allen group, have produced some remarkable studies demonstrating the role of tunneling control. (I have made quite a number of posts on this topics.) Tunneling control is a third mechanism for dictating product formation, in tandem with kinetic control (the favored product is the one that results from the lowest barrier) and thermodynamic control (the favored product is the one that has the lowest energy). Tunneling control has the favored product resulting from the narrowest mass-considered barrier.

Schreiner has written a very clear perspective on tunneling control. It is framed quite interestingly by some fascinating quotes:

It is probably fair to say that many organic chemists view the concept of tunneling, even of hydrogen atoms, with some skepticism. – Carpenter 1983²

Reaction processes have been considered as taking place according to the laws of classical mechanics, quantum mechanical theory being only employed in calculating interatomic forces. – Bell 1933³

Schreiner’s article makes it very clear how critical it is to really think about reactions from a truly quantum mechanical perspective. He notes the predominance of potential energy diagrams that focus exclusively on the relative energies and omits any serious consideration of the reaction coordinate metrics, like barrier width. When one also considers the rise in our understanding of the role of reaction dynamics in organic chemistry (see, for example, these many posts), just how long will it take for these critical notions to penetrate into standard organic chemical thinking? As Schreiner puts it:

It should begin by including quantum phenomena in introductory textbooks, where they are, at least in organic chemistry, blatantly absent. To put this oversight in words similar to those used much earlier by Frank Weinhold in a different context: “When will chemistry textbooks begin to serve as aids, rather than barriers, to this enriched quantum-mechanical perspective?”⁴

References

1) Schreiner, P. R., "Tunneling Control of Chemical Reactions: The Third Reactivity Paradigm." J. Am. Chem. Soc. 2017, 139, 15276-15283, DOI: 10.1021/jacs.7b06035.

2) Carpenter, B. K., "Heavy-atom tunneling as the dominant pathway in a solution-phase reaction? Bond shift in antiaromatic annulenes." J. Am. Chem. Soc. 1983, 105, 1700-1701, DOI: 10.1021/ja00344a073.

3) Bell, R. P., "The Application of Quantum Mechanics to Chemical Kinetics." Proc. R. Soc. London, Ser. A1933, 139 (838), 466-474, DOI: 10.1098/rspa.1933.0031.

4) Weinhold, F., "Chemistry: A new twist on molecular shape." Nature 2001, 411, 539-541, DOI: 10.1038/35079225.

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

The Cope Rearrangement of 1,5-Dimethylsemibullvalene-2(4)-d1: Experimental Evidence for Heavy-Atom Tunneling

Schleif, T.; Mieres-Perez, J.; Henkel, S.; Ertelt, M.; Borden, W. T.; Sander, W., Angew. Chem. Int. Ed. 2017, 56, 10746-10749
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Another prediction made by quantum chemistry has now been confirmed. In 2010, Zhang, Hrovat, and Borden predicted that the degenerate rearrangement of semibullvalene 1 occurs with heavy atom tunneling.¹ For example, the computed rate of the rearrangement including tunneling correction is 1.43 x 10^-3 s^-1 at 40 K, and this rate does not change with decreasing temperature. The predicted half-life of 485 s is 10¹⁰ shorter than that predicted by transition state theory.

Now a group led by Sander has examined the rearrangement of deuterated 2.² The room temperature equilibrium mixture of d⁴–2 and d²–2 was deposited at 3 K. IR observation showed a decrease in signal intensities associated with d⁴–2 and concomitant growth of signals associated with d²–2. The barrier for this interconversion is about 5 kcal mol^-1, too large to be crossed at this temperature. Instead, the interconversion is happening by tunneling through the barrier (with a rate about 10^-4 s^-1), forming the more stable isomer d²–2 preferentially. This is exactly as predicted by theory!

References

1. Zhang, X.; Hrovat, D. A.; Borden, W. T., "Calculations Predict That Carbon Tunneling Allows the Degenerate Cope Rearrangement of Semibullvalene to Occur Rapidly at Cryogenic Temperatures." Org. Letters 2010, 12, 2798-2801, DOI: 10.1021/ol100879t.

2. Schleif, T.; Mieres-Perez, J.; Henkel, S.; Ertelt, M.; Borden, W. T.; Sander, W., "The Cope Rearrangement of 1,5-Dimethylsemibullvalene-2(4)-d1: Experimental Evidence for Heavy-Atom Tunneling." Angew. Chem. Int. Ed. 2017, 56, 10746-10749, DOI: 10.1002/anie.201704787.

InChIs

1: InChI=1S/C8H8/c1-3-6-7-4-2-5(1)8(6)7/h1-8H
InChIKey=VEAPRCKNPMGWCP-UHFFFAOYSA-N

d⁴–2: InChI=1S/C10H12/c1-9-5-3-7-8(4-6-9)10(7,9)2/h3-8H,1-2H3/i5D
InChIKey=WUJOLJNLXLACNA-UICOGKGYSA-N

d²–2: InChI=1S/C10H12/c1-9-5-3-7-8(4-6-9)10(7,9)2/h3-8H,1-2H3/i7D
InChIKey=WUJOLJNLXLACNA-WHRKIXHSSA-N

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

An automated transition state search and its application to diverse types of organic reactions

Leif D. Jacobson, Art D. Bochevarov, Mark A. Watson, Thomas F. Hughes, David Rinaldo, Stephan Ehrlich, Thomas B. Steinbrecher, S. Vaitheeswaran, Dean M. Philipp, Mathew D. Halls, and Richard A. Friesner (2017)

Highlighted by Jan Jensen

Copyright 2017 American Chemical Society

Finding transition states remains one of the most labor intensive pursuits in computational chemistry.  While interpolation methods are becoming increasingly robust, they usually require that the atom order for reactant and product are identical (atom mapping) and can be sensitive the starting conformations and relative orientation in case of bi-molecular reactions.  Furthermore, one still has to check whether the right TS is found and formulate a strategy if it is not.  All these things to do not immediately lend themselves to automation but this paper proposes solutions for all these problems.

In particular the paper offers a very elegant solution for the atom mapping problem: bonds are broken in both reactants and products until the connectivity of the fragments are identical after which the atoms in the fragments can be easily matched. Both the comparison and atom mapping of fragments can be easily done with modern cheminformatics toolkits such as RDKit using canonical smiles and  maximum common substructure searchers (after atom order and charge has been removed).  Cases where this fails due to equivalent atoms (e.g. the hydrogens in a methylene group) can then be dealt with by searching for the solution with the lowest RSMD between reactant and product.

The study focussed on relatively small and rigid molecules and issues due to multiple conformations is left for a future publication.

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

How To Arrive at Accurate Benchmark Values for Transition Metal Compounds: Computation or Experiment?

Y. A. Aoto, A. P. de Lima Batista, A. Köhn, A. G. S. de Oliveira-Filho, J. Chem. Theor. Comput., 2017
Contributed by Theo Keane

Copyright 2017 American Chemical Society 

When performing calculations of any kind, it is important to establish how accurate a method one intends to use is for a given application. Transition metals (TMs) are often problematic systems for computational chemists, because they exhibit “strong correlation”, i.e. either static or dynamic correlation is significant in systems that contain TMs (usually both). This paper adds to the existing literature of benchmark results for TM compounds by performing some rather high-level calculations on 60 diatomic TM compounds. I believe this is intended to improve upon the recent 3dMLBE20 set of Truhlar and co-workers,[1] which was criticised in a pair of papers published in early 2017.[2,3]

In this new benchmarking set, 43 molecules contain first row TMs, 7 contain Ru, Rh or Ag, and the remaining 10 contain Ir, Pt and Au. The multiplicities range from 1 to 7. For these molecules they have assembled experimental data, including bond length, harmonic vibrational frequency and bond dissociation energies. Single reference benchmark values were obtained via (RO/U)CCSD(T)/aug-cc-pwCVnZ-PP[4d, 5d metals]aug-cc-pwCVnZ[else], with n = T, Q, 5, and were extrapolated to the Complete-Basis-Set (CBS) limit. They also investigated the effect of core-valence correlation on the single-reference values. Furthermore, internally contracted Multi-Reference CCSD(T) (icMRCCSD(T)) calculations were performed in the aug-cc-pwCVTZ(-PP) basis, based on full-valence CASSCF reference wavefunctions to investigate the effect of static correlation – the full details of the chosen active spaces are provided in the SI (Table S6). Finally, relativistic effects were considered: scalar relativistic corrections were obtained by comparing frozen core CCSD(T)/aug-cc-pwCVTZ(-PP) calculations with and without the 2nd-order Douglas-Kroll-Hess (DKH2) Hamiltonian. For the 4d and 5d TM containing molecules, Spin-Orbit corrections were obtained from CASSCF calculations with full valence active spaces and the full, 2 electron Breit-Pauli operator. It is important to note that, with the exception of the SO correction, these corrections were not merely calculated for the equilibrium geometry, rather these were calculated at multiple points along the bond length. Overall, the authors have clearly spent a great deal of care ensuring that their ‘benchmark level’ calculations are truly deserving of the title.

An interesting thing to note is that multi-reference, spin-orbit and core-valence correlation corrections all appear to be very weak and sometimes do not improve the agreement with experiment (Table 3). CBS extrapolation is by far the major way to reduce error. This is very important to bear in mind when looking at previous benchmarking results. The authors also note that the usual ‘multireference’ diagnostics are practically useless: there is weak correlation between diagnostics and, more critically, there is very weak correlation between any of the diagnostics and the magnitude of any MR corrections. The M diagnostic[4] is the best performing one; however, it still fails for approximately 30% of cases and yields both false positives and false negatives. The authors also briefly investigate the effect of including 4f orbitals into the correlation treatment for Ir and Pt and find that this has a very weak effect on their results (SI, Table S5).

Finally, the authors use their new benchmark set to rank some functionals. Overall, at the DFT/aug-cc-pVQZ + DKH2 correction level, it appears that hybrid functionals performs on average the best for bond-dissociation energies and equilibrium distances, when compared to the fully corrected results (Table 7). On the other hand, pure functionals perform better for harmonic frequencies. In agreement with the conclusions of the original 3dMLBE20 paper, it is clear that many functionals beat plain CCSD(T)(FC)/aug-cc-pwCVTZ. This reinforces the critical need for CBS extrapolation when performing CC calculations.

(1) Xu, X.; Zhang, W.; Tang, M.; Truhlar, D. G. Do Practical Standard Coupled Cluster Calculations Agree Better than Kohn–Sham Calculations with Currently Available Functionals When Compared to the Best Available Experimental Data for Dissociation Energies of Bonds to 3 D Transition Metals? J. Chem. Theory Comput. 2015, 11 (5), 2036–2052 DOI: 10.1021/acs.jctc.5b00081.
(2) Cheng, L.; Gauss, J.; Ruscic, B.; Armentrout, P. B.; Stanton, J. F. Bond Dissociation Energies for Diatomic Molecules Containing 3d Transition Metals: Benchmark Scalar-Relativistic Coupled-Cluster Calculations for 20 Molecules. J. Chem. Theory Comput. 2017, 13 (3), 1044–1056 DOI: 10.1021/acs.jctc.6b00970.
(3) Fang, Z.; Vasiliu, M.; Peterson, K. A.; Dixon, D. A. Prediction of Bond Dissociation Energies/Heats of Formation for Diatomic Transition Metal Compounds: CCSD(T) Works. J. Chem. Theory Comput. 2017, 13 (3), 1057–1066 DOI: 10.1021/acs.jctc.6b00971.
(4) Tishchenko, O.; Zheng, J.; Truhlar, D. G. Multireference Model Chemistries for Thermochemical Kinetics. J. Chem. Theory Comput. 2008, 4 (8), 1208–1219 DOI: 10.1021/ct800077r.

An Atomistic Fingerprint Algorithm for Learning Ab Initio Molecular Force Fields

Yu-Hang Tang, Dongkun Zhang, and George Em Karniadakis, arXiv:1709.09235
Contributed by Jesper Madsen

Modeling potential energy landscapes of complex atomic environments is challenging. Conventional interatomic potentials are very useful because the potential energy surface is well approximated by some appropriate smooth function of nuclear coordinates. However, choosing the functional form too simple and closed comes with severe limitations because the true potential energy surface may not be (easily) decomposable. 

Instead of sticking with an explicit functional form, one can use continuous density-fields, formed by superimposition of a smoothing kernel on the atoms of the atomic configuration, in order to represent and compare atomistic neighborhoods. Herein, I highlight a recent example of such a method called Density-Encoded Canonically Aligned Fingerprint (DECAF). 

Figure 1: (A) “Two 1D density profiles, ρ1 and ρ2, are generated from two different atomistic configurations using atom-centered smoothing kernel functions. The ‘distance’ between them is measured as the L2 norm of their difference, which corresponds to the highlighted area in the middle plot.” (B) “Shown here is a 2D density field using smoothing kernels whose widths depend on the distances of the atoms from the origin. Darker shades indicate higher density.” 

The preprint by Tang et al. describes the DECAF algorithm (Fig. 1) and also briefly reviews and critically compares with the recent literature of similar methods [such as Smooth Overlap of Atomic Positions (SOAP), Coulomb Matrix, Graph Approximated Energy (GRAPE), and Atom-Centered Symmetry Functions].

The work rests on the key idea of splitting up conventional functional forms into two separate problems, one of representation and one of interpolation, which appears particularly powerful. Molecular fingerprint algorithms such as DECAF are promising in representing atomic neighborhoods faithfully using kernel regression methods. All the beneficial tools and analyses from modern statistics come into play, but there are still open questions that remain. For instance, it is not clear which smoothing kernel, distance metric (and so on) is superior in relating atomic configurations to one-another -- both in general and in specific situations. It is conceivable that there does not exist a best one-size-fits-all option. Furthermore, there will as always be tradeoffs between resolution and computational costs. For an introductory discussion on these topics, the preprint by Tang et al. (and the references within) is a good place to start.

Analyzing Reaction Rates with the Distortion/Interaction-Activation Strain Model

Bickelhaupt, F. M.; Houk, K. N.,  Angew. Chem. Int. Ed. 2017, 56, 10070-10086
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Bickelhaupt and Houk present a nice review of their separately developed, but conceptually identical model for assessing reactivity.¹ Houk termed this the “distortion/interaction” model,² while Bickelhaupt named it “activation strain”.³ The concept is that the activation barrier can be dissected in a distortion or stain energy associated with bringing the reactants into the geometry of the transition state, and the interaction energy is the stabilization energy afforded by the molecular orbital interactions of the reactant components with each other in the transition state.

The review discusses a broad range of applications, including S_N2 and E₂ reactions, pericyclic reactions (including Diels-Alder reactions of enones and the dehdydro Diels-Alder reaction that I have discussed in this blog), a click reaction, a few examples involving catalysis, and the regioselectivity of indolyne (see this post). They also discuss the role of solvent and the relationship of this model to Marcus Theory.

I also want to mention in passing a somewhat related article by Jorgensen and co-authors published in the same issue of Angewandte Chemie as the above review.⁴ This article discusses the paucity of 10 electron cycloaddition reactions, especially in comparison to the large number of very important cycloaddition reactions involving 6 electrons, such as the Diels-Alder reaction, the Cope rearrangement, and the Claisen rearrangement. While the article does not focus on computational methods, computations have been widely used to discuss 10-electron cycloadditions. The real tie between this paper and the review discussed above is Ken Houk, whose graduate career started with an attempt to perform a [6+4] cycloaddition, and he has revisited the topic multiple times throughout his career.

References

1. Bickelhaupt, F. M.; Houk, K. N., "Analyzing Reaction Rates with the Distortion/Interaction-Activation Strain Model." Angew. Chem. Int. Ed. 2017, 56, 10070-10086, DOI: 10.1002/anie.201701486.

2. Ess, D. H.; Houk, K. N., "Distortion/Interaction Energy Control of 1,3-Dipolar Cycloaddition Reactivity." J. Am. Chem. Soc. 2007, 129, 10646-10647, DOI: 10.1021/ja0734086

3. Bickelhaupt, F. M., "Understanding reactivity with Kohn-Sham molecular orbital theory: E2-S_N2 mechanistic spectrum and other concepts." J. Comput. Chem. 1999, 20, 114-128

4. Palazzo, T. A.; Mose, R.; Jørgensen, K. A., "Cycloaddition Reactions: Why Is It So
Challenging To Move from Six to Ten Electrons?" Angew. Chem. Int. Ed. 2017, 56, 10033-10038, DOI: 10.1002/anie.201701085.

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

Benzophenone Ultrafast Triplet Population: Revisiting the Kinetic Model by Surface-Hopping Dynamics

Marco Marazzi, Sebastian Mai, Daniel Roca-Sanjuán, Mickaël G. Delcey, Roland Lindh, Leticia González, and Antonio Monari (2016)
Highlighted by Ravi Kumar Venkatraman

In 1967, Norrish and Porter were honoured with Nobel prize for their seminal work on understanding the fast chemical reactions using flash photolysis technique.¹ Since then benzophenone (Bzp) serves as an archetypal system for understanding the photochemistry of various aromatic ketones. Aromatic ketones find their use in various technologically significant applications like sunscreen, photocatalysis, etc., apart from their fundamental interest.² Efficacy of aromatic ketones for use in various applications relies upon their photophysics and photochemistry. Therefore, understanding the photophysics and photochemistry of Bzp has attracted several experimental and theoretical investigations.² Despite these myriads of investigations, pathways for populating the lowest triplet state (T₁) after photoexcitation to the S1 state remains still elusive. There are two plausible pathways: i) a direct ISC from S₁(nπ*) to T₁(nπ*); or ii) an indirect process, involving ISC from S₁(nπ*) to T₁(ππ*) with subsequent internal conversion (IC) to T₁(nπ*). The latter pathway is more efficient, according to El-Sayed’s rule, because it entails a change in the orbital character during the spin-orbit coupling mediated process.³

Reproduced with permission from Marco Marazzi, Sebastian Mai, et. al., J. Phys. Chem. Lett., 7, 622 (2016) under Creative Commons Attribution (CC-BY) License

In this work, authors have employed ab initio surface-hopping dynamics simulation for Bzp in gas phase to explore the pathways for the lowest triplet state (T₁) population after photo-excitation to the S₁ state. This study clearly demonstrates that the dominant mechanism for populating the T₁ state is the indirect pathway invoking T₂ state as an intermediate. This study urges reinvestigation of spectroscopic assignment for Bzp in various time-resolved spectroscopic techniques. Furthermore, the mechanism for the photoinduced energy transfer in photocatalysis and DNA damage studies must be revisited as in principle now both channels involving T₂ and T₁ states are available.

References:
1.) F. Ariese, K. Roy, V. R. Kumar, H. C. Sudeeksha, S. Kayal, S. Umapathy, Time-Resolved Spectroscopy: Instrumentation and Applications in Encyclopedia of Analytical Chemistry, edited by R. A. Meyers, 1-55, (2017) John Wiley & Sons, Ltd.
2.) M. C. Cuquerella, V. L-Vallet, J. Cadet, and M. A. Miranda, Benzophenone Photosensitized DNA Damage, Acc. Chem. Res., 45, 1558 (2012).
3.) Elsayed, M. A., Spin-Orbit Coupling and Radiationless Processes in Nitrogen Heterocyclics. J. Chem. Phys., 38, 2834 (1963).

More applications of computed NMR spectra

Grimblat, N.; Kaufman, T. S.; Sarotti, A. M., "Computational Chemistry Driven Solution to Rubriflordilactone B." Org. Letters 2016, 18, 6420-6423
Reddy, D. S.; Kutateladze, A. G., "Structure Revision of an Acorane Sesquiterpene Cordycepol A." Org. Letters 2016, 18, 4860-4863
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

In this post I cover two papers discussing application of computed NMR chemical shifts to structure identification and (yet) another review of computational techniques towards NMR structure prediction.

Grimblat, Kaufman, and Sarotti¹ take up the structure of rubriflordilactone B 1, which was isolated from Schisandra rubriflora. The compound was then synthesized and its x-ray structure reported, however its NMR did not match with the natural extract. It was suggested that there were actually two compounds in the extract, the minor one was less soluble and is the crystallized 1, and a second compound responsible for the NMR signal.

The authors looked at all stereoisomers of this molecule keeping the three left-most rings intact. The low energy rotamers of these 32 stereoisomers were then optimized at B3LYP/6-31G* and the chemical shifts computed at PCM(pyridine)/mPW1PW91/6-31+G**. To benchmark the method, DP4+ was used to identify which stereoisomer best matches with the observed NMR of authentic 1; the top fit (92.6% probability) was the correct structure.

The 32 stereoisomers were then tested against the experimental NMR of the natural extract. DP4+ with just the proton shifts suggested structure 2 (99.8% probability); however, the ¹³C chemical shifts predicted a different structure. Re-examination of the reported chemical shifts identifies some mis-assigned signals, which led to a higher C-DP4+ prediction. When all 128 stereoisomers were tested, structure 2 had the highest DP4+ prediction (99.5%), but the C-DP4+ prediction remained problematic (10.8%). Analyzing the geometries of all reasonable alternative for agreement with the NOESY spectrum confirmed 2. These results underscore the importance of using all data sources.

Reddy and Kutateladze point out the importance of using coupling constants along with chemical shifts in structure identification.² They examined cordycepol A 3, obtained from Cordyceps ophioglossoides. They noted that the computed chemical shifts and coupling constants of originally proposed structure 3adiffered dramatically from the experimental values.

They first proposed that the compound has structure 3b. The computed coupling constants using their relativistic force field.³ The experimental coupling constants for the proton H1 are 13.4 and 7.1 Hz. The computed values for 3a are 8.9 and 1.6 Hz, and this structure is clearly incorrect. The coupling constants are improved with 3b, but the ¹³C chemical shifts are in poor agreement with experiment. So, they proposed structure 3c, the epimer at both C1 and C11 of the original structure.

They optimized four conformations of 3c at B3LYP/6-31G(d) and obtained Boltzmann-weighted chemical shifts at mPW1PW91/6-311+G(d,p). The RMS deviation of the computed ¹³C chemical shifts relative to the experiment is only 1.54 ppm, and more importantly, the computed coupling constants of 13.54 and 6.90 Hz are in excellent agreement with the experiment values.

Lastly, Grimblat and Sarotti present a review of a number of methods for using computed NMR chemical shifts towards structure prediction.⁴ These methods include CP3, DP4, DP4+ (all of which I have posted on in the past) and an artificial neural network approach of their own design. They discuss a number of interesting cases where each of these methods has been crucial in identifying the correct chemical structure.

References

1. Grimblat, N.; Kaufman, T. S.; Sarotti, A. M., "Computational Chemistry Driven Solution to Rubriflordilactone B." Org. Letters 2016, 18, 6420-6423, DOI: 10.1021/acs.orglett.6b03318.

2. Reddy, D. S.; Kutateladze, A. G., "Structure Revision of an Acorane Sesquiterpene Cordycepol A." Org. Letters 2016, 18, 4860-4863, DOI: 10.1021/acs.orglett.6b02341.

3. (a) Kutateladze, A. G.; Mukhina, O. A., "Minimalist Relativistic Force Field: Prediction of Proton–Proton Coupling Constants in 1H NMR Spectra Is Perfected with NBO Hybridization Parameters." J. Org. Chem.2015, 80, 5218-5225, DOI: 10.1021/acs.joc.5b00619; (b) Kutateladze, A. G.; Mukhina, O. A., "Relativistic Force Field: Parametrization of 13C–1H Nuclear Spin–Spin Coupling Constants." J. Org. Chem. 2015, 80, 10838-10848, DOI: 10.1021/acs.joc.5b02001.

4. Grimblat, N.; Sarotti, A. M., "Computational Chemistry to the Rescue: Modern Toolboxes for the Assignment of Complex Molecules by GIAO NMR Calculations." Chem. Eur. J. 2016, 22, 12246-12261, DOI: h10.1002/chem.201601150.

InChIs

1: InChI=1S/C28H30O6/c1-13-9-20(32-26(13)30)25-14(2)24-17-6-5-15-12-28-21(8-7-16(15)18(17)10-19(24)31-25)27(3,4)33-22(28)11-23(29)34-28/h5-9,14,19-22,24-25H,10-12H2,1-4H3/t14-,19+,20-,21-,22+,24-,25-,28+/m0/s1
InChIKey=JGSLSHOXBXVVTQ-NEUKEVNNSA-N

2: InChI=1S/C28H30O6/c1-13-9-20(32-26(13)30)25-14(2)24-17-6-5-15-12-28-21(8-7-16(15)18(17)10-19(24)31-25)27(3,4)33-22(28)11-23(29)34-28/h5-9,14,19-22,24-25H,10-12H2,1-4H3/t14-,19-,20-,21-,22+,24+,25-,28+/m0/s1
InChIKey=JGSLSHOXBXVVTQ-WQIRXNRDSA-N

3c: InChI=1S/C16H28O2/c1-6-11(2)9-14-16(5)12(3)7-8-13(16)15(4,17)10-18-14/h9,12-14,17H,6-8,10H2,1-5H3/b11-9-/t12-,13-,14-,15-,16+/m0/s1
InChIKey=WPQIVUHVYBQTBG-AWEVENECSA-N

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

Efficient DLPNO−CCSD(T)-Based Estimation of Formation Enthalpies for C‐, H‐, O‐, and N‐Containing Closed-Shell Compounds Validated Against Critically Evaluated Experimental Data

Eugene Paulechka and Andrei Kazakov (2017)
Highlighted by Jan Jensen

Copyright 2017 American Chemical Society

A computationally methodology is truly robust when it can be used independently and successfully by other groups.  So Frank Neese was understandably delighted when he saw this paper using his DLPNO-CCSD(T) method, as he mentioned during his talk at WATOC2017.

The paper shows tha that DLPNO-CCSD(T)/quadruple zeta//DFT-D3/triple zeta can be used to predict enthalpies of formation as accurate as you can measure them!  It is actually more accurate than G4, but considerably more computationally efficient. 

DLPNO-CCSD(T) cannot handle open shell systems so the energy of H, C, N, and O are replaced by empirical parameters.  This means that enthalpies of formation for molecules containing other elements cannot be computed without similar parametrisation, but usually atom energies/enthalpies of formation are used "only" to validate the method and are not needed to compute reaction energies.

The largest molecule considered is biphenyl and it is not clear to me that B3LYP-D3 is he optimum choice for more complex molecules requiring a conformer search. But, on the other hand, I also doubt the accuracy is very sensitive to the choice of functional for the small molecules used in the study.  It's easy enough to find out: the most time-consuming calculation (biphenyl) required only 10 hours using 10 CPUs. 


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

The Parameter Uncertainty Inflation Fallacy

Pascal Pernot (2017)
Highlighted by Jonny Proppe

In a recent study on uncertainty quantification, Pernot(1) discussed the effect of model inadequacy on predictions of physical properties. Model inadequacy is a ubiquituous feature of physical models due to various approximations employed in their construction. Along with data inconsistency (e.g., due to incorrect quantification of measurement uncertainty) and parameter uncertainty, model inadequacy only acquires meaning by comparison against reference data. For instance, a model is inadequate if it cannot reproduce reference data within their uncertainty range (cf. Figure 1), given all other sources of error are negligible. While parameter uncertainty is inversely proportional to the size of the reference set, systematic errors based on data inconsistency and model inadequacy remain without explicit identification and elimination.

Figure 1. Illustration of model inadequacy. (a) Reference versus calculated (CCD/6-31G*) harmonic vibrational frequencies reveal a linear trend in the data (red line), which is not the unit line. In this diagram, the uncertainty of the reference data is too small to be visible. (b) Residuals of temperature-dependent viscosity predictions based on a Chapman–Enskog model reveal an oscillating trend, even if the 2 confidence intervals of the reference data are considered. Reproduced from J. Chem. Phys. 147, 104102 (2017), with the permission of AIP Publishing.

In related work, Pernot and Cailliez(2) demonstrated the benefits and drawbacks of several Bayesian calibration algorithms (e.g., Gaussian process regression, hierarchical optimization) in tackling these issues. These algorithms approach model inadequacy either through a posteriori model corrections or by parameter uncertainty inflation (PUI). While a posteriori corrected models cannot be transferred to observables not included in the reference set, PUI ensures that the corresponding covariance matrices are transferable to any model comprising the same parameters. However, the resulting predictions may not reflect the correct dependence on the input variable(s), which is determined by the sensitivity coefficients of the model (the partial derivatives of a model prediction at a certain point in input space with respect to the model parameters at their expected values). Pernot referred to this issue as the “PUI fallacy”1 and illustrated it at three examples: (i) linear scaling of harmonic vibrational frequencies, (ii) calibration of the mBEEF density functional against heats of formation, and (iii) inference of Lennard–Jones parameters for predicting temperature-dependent viscosities based on a Chapman–Enskog model (cf. Figure 2). In these cases, PUI resulted in correct average prediction uncertainties, but uncertainties of individual predictions were systematically under- or overestimated.

Figure 2. Illustration of the PUI fallacy for different algorithms (VarInf_Rb, Margin, ABC) at the example of temperature-dependent viscosity predictions based on a Chapman–Enskog model. In all cases, the centered prediction bands (gray) cannot reproduce the oscillating trend in the residuals. Reproduced from J. Chem. Phys. 147, 104102 (2017), with the permission of AIP Publishing

Pernot’s paper presents a state-of-the-art study for rigorous uncertainty quantification of model predictions in the physical sciences, which only recently started to gain momentum in the computational chemistry community. His study can be seen as an incentive for future benchmark studies to rigorously assess existing and novel models. Noteworthy, Pernot has made available the entire code employed in his study (https://github.com/ppernot/PUIF). 

(1) Pernot, P. The Parameter Uncertainty Inflation Fallacy. J. Chem. Phys. 2017, 147 (10), 104102.
(2) Pernot, P.; Cailliez, F. A Critical Review of Statistical Calibration/Prediction Models Handling Data Inconsistency and Model Inadequacy. AIChE J. 2017, 63 (10), 4642–4665.

Spectroscopic Observation of the Triplet Diradical State of a Cyclobutadiene

Kostenko, A.; Tumanskii, B.; Kobayashi, Y.; Nakamoto, M.; Sekiguchi, A.; Apeloig, Y., Angew. Chem. Int. Ed. 2017, 56, 10183-10187
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Cyclobutadiene has long fascinated organic chemists. It is the 4e analogue of the 6e benzene molecule, yet it could hardly be more different. Despite nearly a century of effort, cyclobutadiene analogues were only first prepared in the 1970s, reflecting its strong antiaromatic character.

Per-trimethylsilylcyclobutadiene 1 offers opportunities to probe the properties of the cyclobutadiene ring as the bulky substituents diminish dimerization and polymerization of the reactive π-bonds. Kostenko and coworkers have now reported on the triplet state of 1.¹ They observe three EPR signals of 1 at temperatures above 350 K, and these signals increase in area with increasing temperature. This is strong evidence for the existence of triplet 1 in equilibrium with the lower energy singlet. Using the variable temperature EPR spectra, the singlet triplet gap is 13.9 ± 0.8 kcal mol^-1.

The structures of singlet and triplet 1 were optimized at B3LYP-D3/6-311+G(d,p) and shown in Figure 1. The singlet is the expected rectangle, with distinctly different C-C distance around the ring. The triplet is a square, with equivalent C-C distances. Since both the singlet and triplet states are likely to have multireference character, the energies of both states were obtained at RI-MRDDCI2-CASSCF(4,4)/def2-SVP//B3LYPD3/6-311+G(d,p) and give a singlet-triplet gap of 11.8 kcal mol^-1, in quite reasonable agreement with experiment.

singlet

triplet

Figure 1. Optimized geometries of singlet and triplet 1.

References

1. Kostenko, A.; Tumanskii, B.; Kobayashi, Y.; Nakamoto, M.; Sekiguchi, A.; Apeloig, Y., "Spectroscopic Observation of the Triplet Diradical State of a Cyclobutadiene." Angew. Chem. Int. Ed. 2017, 56, 10183-10187, DOI: 10.1002/anie.201705228.

InChIs

1: InChI=1S/C16H36Si4/c1-17(2,3)13-14(18(4,5)6)16(20(10,11)12)15(13)19(7,8)9/h1-12H3
InChIkey=AYOHYRSQVCLGKR-UHFFFAOYSA-N

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

How Large is the Elephant in the Density Functional Theory Room?

Frank Jensen (2017)

Highlighted by Jan Jensen

Having highlighted this paper it is only right that I highlight Frank Jensen's response. To recap, the previous study used wavelets to compute benchmark energies for PBE and PBE0 functionals and showed that even aug-cc-pV5Z failed to reach chemical accuracy for some atomization energies.

In this paper Jensen shows that this problem goes away when one uses basis sets specifically designed for DFT calculations. At the pentuple-zeta level the maximum errors are reduced by factors of 5 and 10 for segmented contracted and uncontracted basis sets, respectively and the MAE for atomization energies are well below 1 kcal/mol.


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

Is the Accuracy of Density Functional Theory for Atomization Energies and Densities in Bonding Regions Correlated?

Brorsen, K. R.; Yang, Y.; Pak, M. V.; Hammes-Schiffer, S., J. Phys. Chem. Lett. 2017, 8, 2076-2081
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

I recently blogged about a paper arguing that modern density functional development has strayed from the path of improving density description, in favor of improved energetics. The Medvedev paper¹ was met with a number of criticisms. A potential “out” from the conclusions of the work was that perhaps molecular densities do not fare so poorly with more modern functionals, following the argument that better energies might reflect better densities in bonding regions.

The Hammes-Schiffer group have now examined 14 diatomic molecules with the goal of testing just this hypothesis.² They subjected both homonuclear diatomics, like N₂, Cl₂, and Li₂, and heteronuclear diatomics, like HF, LiF, and SC, to 90 different density functionals using the very large aug-cc-pCVQZ basis set. Using the CCSD density as a reference, they examined the differences in the densities predicted by the various functional both along the internuclear axis and perpendicular to it.

The 20 functionals that do the best job in mimicking the CCSD density are all hybrid GGA functionals, along with the sole double hybrid functional included in the study (B2PLYP). These functionals date from 1993 to 2012. The 20 functionals that do the poorest job include functionals from all rung-types, and date from 1980-2012. A very slight upward trend can be observed in the density error increasing with development year, while the error in the dissociation energy clearly is decreasing over time.

They note that six functionals of the Minnesota-type, those that are highly parameterized and of recent vintage, perform very poorly at predicting atomic densities, but do well with the diatomic densities.

Hammes-Schiffer concludes that their diatomic results support the general trend noted by Medvedev’s atomic results, that density description is lagging in more recently developed functionals. I’d add that this trend is not as dramatic for the diatomics as for atoms.

They pose what is really the key question: “Is the purpose to approximate the exact functional or simply to provide chemists with a useful tool for exploring chemical systems?” Since, as they note, the modern highly parameterized functionals have worked so well for predicting energies and geometries, “the observation that many modern functionals produce incorrect densities could be of no great consequence for many studies”. Nonetheless, “the ultimate goal is still to obtain both accurate densities and accurate energies”.

References

1) Medvedev, M. G.; Bushmarinov, I. S.; Sun, J.; Perdew, J. P.; Lyssenko, K. A., "Density functional theory is straying from the path toward the exact functional." Science 2017, 355, 49-52, DOI: 10.1126/science.aah5975.

2) Brorsen, K. R.; Yang, Y.; Pak, M. V.; Hammes-Schiffer, S., "Is the Accuracy of Density Functional Theory for Atomization Energies and Densities in Bonding Regions Correlated?" J. Phys. Chem. Lett. 2017, 8, 2076-2081, DOI: 10.1021/acs.jpclett.7b00774.

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