Accurate ab Initio Spin Densities

Katharina Boguslawski, Konrad H. Marti, Örs Legeza and Markus Reiher Journal of Chemical Theory and Computation 2012, 8, 1970 (Paywall)

The spin density problem
Obtaining accurate spin density distributions and predicting the correct ground state from a number of close lying states of different spin is a challenging problem in quantum chemistry, particularly when transition metal systems are considered. This paper highlights how density-matrix renomarlization group¹ (DMRG) based methods can be used to calculate spin density distributions for molecules that are too large to be treated by complete-active-space self-consistent-field (CASSCF) methods. What I found particularly exciting is the prospect of using DMRG results in the benchmarking and development of new density functionals.

Benchmarking DMRG
The electronic structure of a simple model system of iron nitrosyl is manipulated by surrounding it with point-charges (simulating ligands), adjusting the position of these charges means the system becomes either single- or multi-reference. The results convincingly show how DMRG can converge to a CASSCF(7,7) reference spin density by controlling the number of reference active-system states in the DMRG. The DMRG calculations were carried out using the Reiher group's Qc-Dmrg-ETH code.

Spin densities for large active spaces and comparison to DFT
A problem arrises in CASSCF when the number of active electrons and orbitals required to correctly describe the system becomes too large (this article suggests 18 electrons in 18 orbitals as a practical limit). The paper goes on to show that much larger active spaces are possible with DMRG and, more importantly, that the spin density converges very quickly both with respect to active space and the number of active-system states. The resulting spin densities are then compared to a number of common density functionals, none of which accurately reproduces the DMRG result. Whilst DMRG calculations are not yet commonplace in the computational chemistry community, this paper convinces me that they will have an important role to play in solving the spin density problem.

References
(1) For an introduction to DMRG see: Chan, G. K.-L.; Dorando, J. J.; Ghosh, D.; Hachmann. J.; Neuscamman, E.; Wang, H.;Yanai, T. An Introduction to the Density Matrix Renormalization Group Ansatz in Quantum Chemistry. In Frontiers in Quantum Systems in Chemistry and Physics, 1st ed.; Wilson, S., Grout, P. J.; Maruani, J., Delgado-Barrio, G., Piecuch, P., Eds.; Springer: Dordrecht, The Netherlands, 2008; Vol. 18, pp 49-65. arXiv:0711.1398v1.

Jailbreaking Benzene Dimers

A. Yu. Rogachev, X.-D. Wen, R. Hoffmann Journal of the American Chemical Society 2012, 134, 8062. (Paywall)
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Hoffmann¹ reports on a number of new benzene dimer structures, notably 5-8, whose RIJCOSX-MP2/cc-pVTZ² structures are shown in Figure 1. A few of these new dimers are only somewhat higher in energy than the known dimers 1-4. The energies of these dimers, relative to two isolated benzene molecules, are listed in Table 1.

1	2
3	4
5	6
7	8

Figure 1. RIJCOSX-MP2/cc-pVTZ optimized geometries of 1-8.

Table 1. Energy (kcal mol^-1) of the dimers relative to two benzene molecules and activation energy for reversion to two benzene molecules.

Compound	Erel	Eact
1	50.9	29
2	49.9
3	38.2	9
4	58.7	19
5	71.9	30
6	49.9	36
7	60.8	27
8	98.8	28

The energy for reversion of the isomers 5-8 to two isolated benzene molecules is calculated to be fairly large, and so they should be stable relative to that decomposition mode. They also examined a series of other decomposition modes, including [1,5]-hydrogen migration, all of which had barriers of 21 kcal mol^-1 or greater, retrocyclization ([2+2]), for which they could not locate transition states, electrocyclic ring opening (Cope), with barriers of at least 17 kcal mol^-1 and dimerization – some of which had relatively small enthalpic barriers of 4-5 kcal mol^-1. However, the dimerizations all have very unfavorable entropic activation barriers.

So, the conclusion is that all of the novel dimers (4-8) can be reasonable expected to hang around for some time and therefore are potential synthetic targets.

References

(1) Rogachev, A. Yu.; Wen, X.-D.; Hoffmann, R. "Jailbreaking Benzene Dimers," J.
Am. Chem. Soc., 2012, 134, 8062-8065, DOI:10.1021/ja302597r

(2) Kossmann, S.; Neese, F. "Efficient Structure Optimization with Second-Order Many-Body Perturbation Theory: The RIJCOSX-MP2 Method," J. Chem. Theory Comput., 2010, 6, 2325-2338, DOI:10.1021/ct100199k

InChIs

1: InChI=1S/C12H12/c1-2-6-10-9(5-1)11-7-3-4-8-12(10)11/h1-12H/t9-,10+,11-,12+
InChIKey=WMPWOGVJEXSFLI-UHFFFAOYSA-N

2: InChI=1S/C12H12/c1-2-6-10-9(5-1)11-7-3-4-8-12(10)11/h1-12H/t9-,10+,11+,12-
InChIKey=WMPWOGVJEXSFLI-IWDIQUIJSA-N

3: InChI=1S/C12H12/c1-2-4-12-10-7-5-9(6-8-10)11(12)3-1/h1-12H/t9?,10?,11-,12+
InChIKey=ONVDJSCNMCYFTI-CAODYFQJSA-N

4: InChI=1S/C12H12/c1-2-10-4-3-9(1)11-5-7-12(10)8-6-11/h1-12H
InChIKey=BCBHEUOKKNYIAT-UHFFFAOYSA-N

5: InChI=1S/C12H12/c1-2-6-10-9(5-1)11-7-3-4-8-12(10)11/h1-12H/t9-,10-,11+,12+/m1/s1
InChIKey=WMPWOGVJEXSFLI-WYUUTHIRSA-N

6: InChI=1S/C12H12/c1-2-4-12-10-7-5-9(6-8-10)11(12)3-1/h1-12H/t9?,10?,11-,12-/m0/s1
InChIKey=ONVDJSCNMCYFTI-QQFIATSDSA-N

7: InChI=1S/C12H12/c1-2-6-10-9(5-1)11-7-3-4-8-12(10)11/h1-12H/t9-,10-,11-,12+/m1/s1
InChIKey=WMPWOGVJEXSFLI-KKOKHZNYSA-N

8: InChI=1S/C12H12/c1-2-6-10-9(5-1)11-7-3-4-8-12(10)11/h1-12H/t9-,10-,11-,12-
InChIKey=WMPWOGVJEXSFLI-NQYKUJLISA-N

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

Interpreting Protein Structural Dynamics from NMR Chemical Shifts

Paul Robustelli, Kate A. Stafford, and Arthur G. Palmer, III Journal of the American Chemical Society 2012, 134, 6365 (Paywall)

The effect of dynamical averaging on chemical shift predictions

In this work snapshots from MD simulations of are used as input to chemical shift predictors Sparta+, ShiftX+, Camshift, and ShiftS to compute average chemical shifts for two homologous ribonuclease H enzymes from Escherichia coli (ecRNH) and Thermus thermophilus (ttRNH).  What I especially liked about this article is not only the thoroughness: all the most popular chemical shift predictors were compared and the averaging is based on 100 ns MD simulations, one of which was tested against a 1000 ns MD simulation, but also how clearly the results were analyzed and presented.

Does averaging improve overall accuracy?

In the case of ecRNH, RMSDs from experiment computed from averaged chemical shifts are always more accurate than those computed from the x-ray structure for all nuclei and methods tested - on average by 23%, while for ttRNH the improvement is much less and for methods like ShiftX+ averaging consistently leads to slightly worse results. The resolution of the x-ray structures of ecRNH and ttRNH are 2.8 Å and 1.5 Å, respectively, so some of the improvement for ecRNH presumably derives from errors in the low resolution x-ray structure that gets fixed by the MD, and the authors describe one concrete example of this.

Identifying and analyzing interesting residues

However, as the authors point out RMSDs only tell part of the story.  One interesting analysis tool employed in the study is a plot of $|\delta_{x-ray}-\delta_{Exp}|-|\delta_{MD}-\delta_{Exp}|$ as a function of residue number to identify where averaging significantly improves the prediction and where it makes it significantly worse compared to using the x-ray structure.  Some of these cases are then chosen for further analysis to highlight the different underlying structural causes.  Here the authors use another powerful tool: a plot of the probability distribution of chemical shifts [$P(\delta)$], which in many of these cases is bi-modal, i.e. it clearly shows that the average chemical shift is not the most likely chemical shift and that it represents at least two distinct conformations.

The size of the ensemble: dynamical vs conformational averaging

Most of the averages presented in the paper are computed using 22,200 protein structures, but looking at the data makes me wonder how many conformations it really represents.  For example, in the case of Sparta+ C$_\alpha$ chemical shift predictions for ttRNH, averaging leads to significant improvements over x-ray for only about 3-4 residues out of about 145 (and worse predictions for three other residues!).

Furthermore, the peak or peaks in the $P(\delta)$-plots appear quite symmetric in the examples shown, which means that the most probable chemical shift from that peak is a good approximation to the average value computed using that peak.  I wouldn't be surprised if one could identify 10-20 ttRNH structures that would lead to equally accurate average chemical shift values.

Empirical vs QM-based chemical shift predictors

Finally, the authors note an interesting case (Phe27 C') where ShiftS, which is parameterized based on QM-data, gives a completely different $P(\delta)$-plot compared to the other methods (although an equally accurate average chemical shift value for that particular residue!).  I would love to see a detailed analysis of this issue before completely ruling out the ShiftS results.

Acknowledgements.  I thank Anders Christensen for alerting me to this article and for useful discussions.

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

On the aromatic stabilization of corannulene and coronene

Dobrowolski, M. A.; Ciesielski, A.; Cyranski, M. K. Phys. Chem. Chem. Phys., 2011, 13, 20557-20563 (Paywall)
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission


How should one assess the aromatic stabilization energy of non-planar compounds like corannulene 1? The standard approach might be to employ the homodesmotic reaction, like Reaction 1. The energy of this reaction is however quite different whether one chooses s-cis or s-trans butadiene: 67.5 kcal mol^-1 with the former and 14.8 kcal mol^-1 with the latter. Exactly how one balances the total number of cis/trans relationships is problematic, but worse still is that Reaction 1 does not remove the effect of strain and non-planarity of 1.

1

2

Reaction 1
15 H₂C=CH-CH=CH₂ + 20 H₂C=C(CH₃)₂ → 1 + 10 H₂C=CH₂ + 20 H₂C=CH-CH₃

Dobrowolski, Ciesielski and Cyranski¹ propose a series of reactions that extend the isomerization stabilization energy concept of Schleyer. References are chosen that involve fixed alternate polyenes by appending methylene groups, creating radialene-like compounds. Reaction 2 and 3 are two such reactions that attempt to remove strain and non-planarity effects along with balancing the cis/trans relationships and potential H-H clashes between the pendant methylene groups. They report an additional 18 variations, because there is no unique method for portioning these effects.

Reaction 2

Reaction 3

Using B3LYP/6-311G** energies with zero-point vibrational energy, the reaction energies are 46.7 and 46.3 kcal mol^-1 for Reactions 2 and 3, respectively. Using all of the variations, the mean value is 44.7 kcal mol^-1 with a standard deviation of only 1.2 kcal mol^-1. It is clear that corranulene has a rather substantial artomatic stabilization energy, reflecting its decided aromatic character.
In a similar vein, they have also estimated the aromatic stabilization energy of coronene 2 as 58.4 kcal mol^-1, which, while clearly demonstrating the 2 is aromatic, it does not express any “superaromaticity”.

References

(1) Dobrowolski, M. A.; Ciesielski, A.; Cyranski, M. K. "On the aromatic stabilization of corannulene and coronene," Phys. Chem. Chem. Phys., 2011, 13, 20557-20563, DOI: 10.1039/C1CP21994D

InChIs

1: InChI=1S/C21H14/c1-3-13-6-7-15-10-11-16-9-8-14-5-4-12(2)17-18(13)20(15)21(16)19(14)17/h3-11H,1H2,2H3
InChIKey=ZJQHTVPYWDRMLD-UHFFFAOYSA-N

2: InChI=1S/C24H12/c1-2-14-5-6-16-9-11-18-12-10-17-8-7-15-4-3-13(1)19-20(14)22(16)24(18)23(17)21(15)19/h1-12H
InChIKey=VPUGDVKSAQVFFS-UHFFFAOYSA-N



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

Computational Design and Selection of Optimal Organic Photovoltaic Materials

Noel M. O’Boyle, Casey M. Campbell and Geoffrey R. Hutchison Journal of Physical Chemistry C 2011, 115, 16200 (Paywall)

Quantum chemistry for high throughput screening
In this work over 90,000 pi-conjugated copolymer were computationally screened for new and efficient organic photo-voltaics. Not being an expert in organic photo-voltaics I highlight this paper as a very interesting example of what I believe is an important emerging trend in the use of quantum chemistry: efficient high through-put screening of molecules for desirable properties.  Computers and software have now reached a point where this is computationally feasible to perform computations on thousands of molecules and the challenge is now to make it practically possible, i.e. to find the right combinations of methods and automate their use.

How does one construct 90,000 geometries?  Cheminformatics meets quantum chemistry
OpenBabel is used to construct the 3D structure of each polymer, starting from its SMILES string, and to find the lowest energy conformation using a weighted rotor-search and the MMFF94 force field.  This lowest energy conformation is then minimized with PM6 and used to compute energies and oscillator strengths of the 15 lowest- energy electronic transitions with ZINDO/S, which forms the basis for estimating the energy conversion efficiency (see the paper for more details).  The required CPU time is  8-10 minutes per polymer on a single core.  To put that number in perspective: using 100 cores, 100,000 polymers can be screened in roughly a week.

Starting from 131 different monomers, all possible (19,701) dimers are made and these dimers are then used to construct the corresponding 58,707 tetramers. (A brief description on how exactly this was scripted would have been a welcome addition to the supplementary materials).  The energy conversion efficiency was computed for all these dimer and tetramers.  These results were used to calibrate a genetic search algorithm that was used to identify hexamers and octamers with high energy conversion efficiency without doing an exhaustive search.

Promising candidates and new strategies
The current state-of-the-art is ca 8% energy conversion efficiency.  This study found 621 polymers (mostly hexamers and octamers) with greater than 9%  and 2 with greater than 11 % energy conversion efficiency.  (As the authors point out these polymers still need to be filtered for solubility, crystal packing, and other factors.)  Just as interestingly new design strategies emerged:

"Our analysis of component monomers, dimers, and the copolymer sequence demonstrates important design rules for copolymer photovoltaics. Most importantly, the conventional picture of combining a strong donor and strong acceptor into an alternating copolymer is found to frequently yield poor energy- level alignment. Instead, our top hexamers and octamers reflect a decreased optical band gap due to high coupling between the two-component monomers, not solely due to particular HOMO or LUMO energies of the monomers themselves."

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

Energy Decomposition Analysis in Solution Based on the Fragment Molecular Orbital Method

D. G. Fedorov, K. Kitaura, J. Phys. Chem. A 116 (2012) 704-719.


The solvent screening has long been described relying on the dielectric constants taken arbitrarily based on "experience". In this work the dielectric constants are derived from ab initio calculations using the charges induced by the solute on the solvent surface in the polarizable continuum model. These charges uniquely define the solvent screening energy, and also the dielectric constant.

Induced solvent charges δ surrounding two subsystems I and J.

[Reprinted with permission from the Journal of Physical Chemistry A.  Copyright (2012) American Chemical Society.]

For the rather trivial case of two interacting systems (sodium and chloride ions), the analysis provides rich physical insight into the solvation process, and the induced charge quenching effect is observed, caused by the partial loss of the solvation shell and the potential pressure of the nearby ion.
The real challenge for a physical interpretation is found in systems with multiple fragments, in which case strong many-body electrostatic effects lead to  very interesting and peculiar results, such as a partial solvent screening reduced by the pressure of the solute-induced electrostatic potential of nearby fragments.

The solvent screening of the electrostatic interactions between charged residues in a small protein chignolin (PDB: 1UAO) is discussed based on ab initio (FMO) calculations in solution.

Amino acid-catalyzed aldol and Michael reactions

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


Here are a couple of articles describing computational approaches to catalytic enantioselective reactions using variations upon the classic proline-catalyzed aldol reaction of List and Barbas¹ that started the whole parade. I have discussed the major computational papers on that system in my book (Chapter 5.3).

Yang and Wong² investigated the proline-catalyzed nitro-Michael reaction, looking at four examples, two with aldehydes and two with ketones (Reactions 1-4).

Reaction 1

Reaction 2

Reaction 3

Reaction 4

These four reactions were examined at MP2/311+G**//M06-2x/6-31G**, and PCM was also applied. The key element of this study is that they examined two different types of transition states: (a) based on the Houk-List model involving a hydrogen bond and (b) an electrostatic based model with no hydrogen bond. These are sketched in Scheme 1. For each of the reactions 1-4 there are 8 located transition states differing in the orientation of the attack on to the syn or anti enamine.

Scheme 1. TS models

Model A

Model B

The two lowest energy TS are shown in Figure 1. TS1-β1-RS is the lowest TS and it leads to the major enantiomer. The second lowest TS, TS1-β3-SR, lies 2.9 kJ mol^-1 above the other TS, and it leads to the minor enantiomer. This lowest TS is of the Houk-List type (Model A) while the other TS is of the Model B type. The enthalpies of activation suggest an ee of 54%, in reasonable agreement with experiment.

TS1-β1-RS

TS1-β1-RS

Figure 1. M06-2x/6-31G** optimized geometries of TS1-β1-RS and TS1-β1-RS.

The computations of the other three reactions are equally good in terms of agreement with experiment, and importantly the computations indicate the reversal of stereoselection between the aldehydes and the ketones. These computations clearly implicate both the Houk-List and the non-hydrogen bonding TSs in the catalyzed Michael additions.

Houk in collaboration with Scheffler and Mahrwald investigate the use of histidine as a catalyst for the asymmetricaldol reaction.³ Examples of the histidine-catalyzed aldol are shown in Reactions 5-7.

	Reaction 5
	Reaction 6
	Reaction 7

The interesting twist here is whether the imidazole can also be involved in hydrogen bonding to the acceptor carbonyl group, serving the purpose of the carboxylic acid group in the Houk-List TS when proline is the catalyst (Model A). The transition states for these and a few other reactions were computed at M06-2x/6-31+G(d,p) including with the SMD continuum solvent model for water. The two lowest energy TSs for the reaction of isobutyraldehde and formaldehyde are shown in Figure 2; TS1 has the carboxylic acid group as the hydrogen donor while the imidazole is the donor in TS2. Of note is that these two TS are isoenergetic, indicating that both modes of stabilization are at play with histidine as the catalyst.

TS1

TS2

Figure 2. M06-2x/6-31+G(d,p) geometries of the two lowest energy TSs for the reaction of isobutyraldehde and formaldehyde catalyzed by histidine.

The possible TSs for Reactions 5-7 were also located. For example, with Reaction 5, the lowest energy TS involves the imidazole as the hydrogen donor and it leads to the major product. The lowest energy TS that leads to the minor product involves the carboxylic acid as the donor. The computed ee’s for Reactions 5-7 are in very good, if not excellent, agreement with the experimental values. The study should spur further activity in which one might tune the stereoselectivity by using catalysts with multiple binding opportunities.


References

(1) List, B.; Lerner, R. A.; Barbas, C. F., III; "Proline-Catalyzed Direct Asymmetric Aldol Reactions," J. Am. Chem. Soc., 2000, 122, 2395-2396, DOI: 10.1021/ja994280y.

(2) Yang, H.; Wong, M. W. "(S)-Proline-catalyzed nitro-Michael reactions: towards a better understanding of the catalytic mechanism and enantioselectivity," Org. Biomol. Chem.,
2012, 10, 3229-3235, DOI: 10.1039/C2OB06993H

(3) Lam, Y.-h.; Houk, K. N.; Scheffler, U.; Mahrwald, R. "Stereoselectivities of Histidine-Catalyzed Asymmetric Aldol Additions and Contrasts with Proline Catalysis: A Quantum Mechanical Analysis," J. Am. Chem. Soc. 2012, 134, 6286-6295, DOI: 10.1021/ja2118392




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