Sunday, February 7, 2016

Molecular Rift: Virtual Reality for Drug Designers

Magnus Norrby, Christoph Grebner, Joakim Eriksson, and Jonas Boström (2015)
Contributed by Jan Jensen

The Oculus Rift will be available soon making it one of the first generally available virtual reality headset. There are some molecular visualization programs, like iView, that has an "oculus rift button", but that will "just" make the image on the screen look more three dimensional. While that's great, "head-tracking, allowing the user to look in all directions, and the wide field of view are key components for providing a real VR experience" (see video below).  "Traditional molecular visualizers show molecules in front of you. With Molecular Rift you can step into ligand-protein complexes and feel like you are really there."

An illustration of head-tracking (not Molecular Rift)

Thus, a virtual reality environment for molecular visualization "objects such as atoms and bonds have a location in 3D-space relative to the user’s position" and this is what Molecular Rift provides.  It is written from scratch using the Unity 3D game engine (which only runs on Windows).

Interacting with the image presents an additional challenge since it's hard to interact with a keyboard of mouse while wearing the headset.  To solve this, the authors also interfaced the program with a Microsoft Kinect.  The code is open source and available on GitHub.

While Molecular Rift is primarily meant for drug design, the Molecular Rift would also clearly be fantastic for teaching, especially if a smartphone version (minus the Kinect interface) could be made to be used in conjunction with Google Cardboard.

Another fantastic extension would be an interface with Reiher's Interactive Chemical Reactivity Exploration tool.

This work is licensed under a Creative Commons Attribution 4.0

Thursday, February 4, 2016

QM/MM Protocol for Direct Molecular Dynamics of Chemical Reactions in Solution: The Water-Accelerated Diels–Alder Reaction

Yang, Z.; Doubleday, C.; Houk, K. N. J. Chem. Theor. Comput. 2015,  5606-5612
Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

I discuss the aqueous Diels-Alder reaction in Chapter 7.1 of my book. A key case is the reaction of methyl vinyl ketone with cyclopentadiene, Reaction 1. The reaction is accelerated by a factor of 740 in water over the rate in isooctane.1 Jorgensen argues that this acceleration is due to stronger hydrogen bonding to the ketone than in the transition state than in the reactants.2-4
Rxn 1
Doubleday and Houk5 report a procedure for calculating trajectories including explicit water as the solvent and apply it to Reaction 1. Their process is as follows:
  1. Compute the endo TS at M06-2X/6-31G(d) with a continuum solvent.
  2. Equilibrate water for 200ps, defined by the TIP3P model, in a periodic box, with the transition state frozen.
  3. Continue the equilibration as in Step 2, and save the coordinates of the water molecules after every addition 5 ps, for a total of typically 25 steps.
  4. For each of these solvent configurations, perform an ONIOM computation, keeping the waters fixed and finding a new optimum TS. Call these solvent-perturbed transition states (SPTS).
  5. Generate about 10 initial conditions using quasiclassical TS mode sampling for each SPTS.
  6. Now for each the initial conditions for each of these SPTSs, run the trajectories in the forward and backward directions, typically about 10 of them, using ONIOM to compute energies and gradients.
  7. A few SPTS are also selected and water molecules that are either directly hydrogen bonded to the ketone, or one neighbor away are also included in the QM portion of the ONIOM, and trajectories computed for these select sets.
The trajectory computations confirm the role of hydrogen bonding in stabilizing the TS preferentially over the reactants. Additionally, the trajectories show an increasing asynchronous reactions as the number of explicit water molecules are included in the QM part of the calculation. Despite an increasing time gap between the formation of the first and second C-C bonds, the overwhelming majority of the trajectories indicate a concerted reaction.


(1) Breslow, R.; Guo, T. "Diels-Alder reactions in nonaqueous polar solvents. Kinetic
effects of chaotropic and antichaotropic agents and of β-cyclodextrin," J. Am. Chem. Soc. 1988110, 5613-5617, DOI: 10.1021/ja00225a003.
(2) Blake, J. F.; Lim, D.; Jorgensen, W. L. "Enhanced Hydrogen Bonding of Water to Diels-Alder Transition States. Ab Initio Evidence," J. Org. Chem. 199459, 803-805, DOI: 10.1021/jo00083a021.
(3) Chandrasekhar, J.; Shariffskul, S.; Jorgensen, W. L. "QM/MM Simulations for Diels-Alder
Reactions in Water: Contribution of Enhanced Hydrogen Bonding at the Transition State to the Solvent Effect," J. Phys. Chem. B 2002106, 8078-8085, DOI: 10.1021/jp020326p.
(4) Acevedo, O.; Jorgensen, W. L. "Understanding Rate Accelerations for Diels−Alder Reactions in Solution Using Enhanced QM/MM Methodology," J. Chem. Theor. Comput. 20073, 1412-1419, DOI:10.1021/ct700078b.
(5) Yang, Z.; Doubleday, C.; Houk, K. N. "QM/MM Protocol for Direct Molecular Dynamics of Chemical Reactions in Solution: The Water-Accelerated Diels–Alder Reaction," J. Chem. Theor. Comput. 2015, 5606-5612, DOI: 10.1021/acs.jctc.5b01029.

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

Sunday, January 31, 2016

Assessment of the accuracy of coupled cluster perturbation theory for open-shell systems. I. Triples expansions

Janus J. Eriksen, Devin A. Matthews, Poul Jørgensen, and Jürgen Gauss (2015)
Contributed by Jan Jensen

Figure 1: Normal distributions of the recoveries of (in percent (%), Figure 1a) and deviations from (in kcal/mol, Figure 1b) CCSDT–CCSD frozen-core/cc-pVTZ correlation energy differences for RHF, UHF, and ROHF references. (Taken from the paper)

This paper argues that we shouldn't expect the kind of accuracy we have come to expect from CCSD(T) for closed shell molecules, for open-shell molecules.  I have a slightly different take on the data.

But first it is fair to ask whether coupled-cluster is an appropriate standard for open shell systems in the first place because of spin-contamination (note there is also spin-contamination in ROHF- coupled-cluster).  If you get significantly different answers for UHF- and ROHF-based calculations it is not clear which one is the most reliable.

The first thing the paper shows - for 18 atoms and small open shell molecules - is that spin-contamination decreases by roughly an order of magnitude for each step in SCF > CCSD > CCSDT > CCSDTQ so that it is negligible for the latter.  For most of the molecules the spin-contamination is also negligible for CCSDT.  This is good to know.

The authors go on to show that CCSD(T) is a worse approximation for CCSDT for open shell systems (see Figure 1) and that CCSD(T) is always a worse approximation to CCSDTQ than CCSDT for open shell systems, in contrast to closed shell systems.  While that's true, the data also shows that the CCSDT correlation energy is closer to the CCSDTQ energy for open shell systems.  In fact the the mean error in CCSD(T) correlation energies relative to CCSDTQ is actually lower for open shell systems (0.80, 0.70, and 0.57 kcal/mol for RHF-, UHF-, and ROHF-based calculations).  So this is good news.

One note of caution in all of this is that the study uses cc-pVTZ (understandable as CCSDTQ calculations are performed). It remains to be seen whether the conclusions are true for the CBS.

On a related note Anacker, Tew, and Friedrich have recently presented an incremental scheme for estimating UHF-based CCSD(T) energies for larger systems.

This work is licensed under a Creative Commons Attribution 4.0

Saturday, January 30, 2016

Seeking Extremes in Molecular Design: To What Extent May Two “Non-Bonded” Hydrogen Atoms be Squeezed in a Hydrocarbon?

Firouzi, R.; Shahbazian, S. ChemPhysChem 2016, 17, 51-54
Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

Setting the record for the shortest non-bonded HH contact has become an active contest. Following on the report of a contact distance of only 1.47 Å that I blogged about here, Firouzi and Shahbazian propose a series of related cage molecules with C-H bonds pointed into their interior.1 The compounds were optimized with a variety of computational methods, and many of them have HH distances well below that of the previous record. The shortest distance is found in 1, shown in Figure 1. The HH distance in 1 is predicted to be less than 1.2 Å with a variety of density functionals and moderate basis sets.

Figure 1. Optimized geometry of 1 at ωB97X-D/cc-pVDZ.


(1) Firouzi, R.; Shahbazian, S. "Seeking Extremes in Molecular Design: To What Extent May Two “Non-Bonded” Hydrogen Atoms be Squeezed in a Hydrocarbon?," ChemPhysChem 201617, 51-54, DOI:10.1002/cphc.201501002.


1: InChI=1S/C41H44/c1-7-13-25-17-9-3-40-5-11-19(35(17)40)27-15-8-2-39(1,33(13)15)34-14(7)26-18-10-4-41-6-12-20(36(18)41)28(16(8)34)32(27)30-23(11)37(40)21(9)29(31(25)26)22(10)38(41)24(12)30/h7-38H,1-6H2

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

Sunday, January 24, 2016

A Practicable Real-Space Measure and Visualization of Static Electron-Correlation Effects

Grimme, S.; Hansen, A. Angew. Chem. Int. Ed. 2015, 54, 12308-12313
Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

Assessing when a molecular system might be subject to sizable static (non-dynamic) electron correlation, necessitating a multi-reference quantum mechanical treatment, is perhaps more art than science. In general one suspects that static correlation will be important when the frontier MO energy gap is small, but is there a way to get more guidance?

Grimme reports the use of fractional occupancy density (FOD) as a visualization tool to identify regions within molecules that demonstrate significant static electron correlation.1 The method is based on the use of finite temperature DFT.2,3

The resulting plots of the FOD for a series of test cases follow our notions of static correlation. Molecules, such as alkanes, simple aromatics, and concerted transition states show essentially no fractional orbital density. On the other hand, the FOD plot for ozone shows significant density spread over the entire molecule; the transition state for the cleavage of the terminal C-C bond in octane shows FOD at C1 and C2but not elsewhere; p-benzyne shows significant FOD at the two radical carbons, while the FOD is much smaller in m-benzyne and is negligible in o-benzyne.
This FOD method looks to be a simple tool for evaluating static correlation and is worth further testing.


(1) Grimme, S.; Hansen, A. "A Practicable Real-Space Measure and Visualization of Static Electron-Correlation Effects," Angew. Chem. Int. Ed. 201554, 12308-12313, DOI: 10.1002/anie.201501887.
(2) Mermin, N. D. "Thermal Properties of the Inhomogeneous Electron Gas," Phys. Rev. 1965137, A1441-A1443, DOI: 10.1103/PhysRev.137.A1441.
(3) Chai, J.-D. "Density functional theory with fractional orbital occupations," J. Chem. Phys 2012136, 154104, DOI: doi: 10.1063/1.3703894.

[Editors note: this paper has also been highlighted by Tobias Schwabe]

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

Wednesday, January 13, 2016

Enhancing NMR Prediction for Organic Compounds Using Molecular Dynamics

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

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

They suggest that the chemical shift can be computed as

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

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

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

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

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



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


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

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

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

Wednesday, January 6, 2016

Attraction or Repulsion? London Dispersion Forces Control Azobenzene Switches

Schweighauser, L.; Strauss, M. A.; Bellotto, S.; Wegner, H. A. Angew. Chem. Int. Ed. 2015, 54, 13436-13439
Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

The role of dispersion in organic chemistry has been slowly recognized as being quite critical in a variety of systems. I have blogged on this subject many times, discussing new methods for properly treating dispersion within quantum computations along with a variety of molecular systems where dispersion plays a critical role. Schreiner1 has recently published a very nice review of molecular systems where dispersion is a key component towards understanding structure and/or properties.

In a similar vein, Wegner and coworkers have examined the Z to E transition of azobenzene systems (1a-g 2a-g) using both experiment and computation.2 They excited the azobenzenes to the Z conformation and then monitored the rate for conversion to the E conformation. In addition they optimized the geometries of the two conformers and the transition state for their interconversion at both B3LYP/6-311G(d,p) and B3LYP-D3/6-311G(d,p). The optimized structure of the t-butyl-substituted system is shown in Figure 1.

a: R=H; b: R=tBu; c: R=Me; d: R=iPr; e: R=Cyclohexyl; f: R=Adamantyl; g: R=Ph



Figure 1. B3LYP-D3/6-311G(d,p) optimized geometries of 1a, 2a, and the TS connecting them.

The experiment finds that the largest activation barriers are for the adamantly  1f  and  t-butyl 1b azobenzenes, while the lowest barriers are for the parent 1a and methylated 1c azobenzenes.

The trends in these barriers are not reproduced at B3LYP but are reproduced at B3LYP-D3. This suggests that dispersion is playing a role. In the Z conformations, the two phenyl groups are close together, and if appropriately substituted with bulky substituents, contrary to what might be traditionally thought, the steric bulk does not destabilize the Z form but actually serves to increase the dispersion stabilization between these groups. This leads to a higher barrier for conversion from the Z conformer to the Econformer with increasing steric bulk.



(1) Wagner, J. P.; Schreiner, P. R. "London Dispersion in Molecular Chemistry—Reconsidering Steric Effects,"Angew. Chem. Int. Ed. 2015, 54, 12274-12296, DOI: 10.1002/anie.201503476.
(2) Schweighauser, L.; Strauss, M. A.; Bellotto, S.; Wegner, H. A. "Attraction or Repulsion? London Dispersion Forces Control Azobenzene Switches," Angew. Chem. Int. Ed. 2015, 54, 13436-13439, DOI:10.1002/anie.201506126.



1b: InChI=1S/C28H42N2/c1-25(2,3)19-13-20(26(4,5)6)16-23(15-19)29-30-24-17-21(27(7,8)9)14-22(18-24)28(10,11)12/h13-18H,1-12H3/b30-29-
2b: InChI=1S/C28H42N2/c1-25(2,3)19-13-20(26(4,5)6)16-23(15-19)29-30-24-17-21(27(7,8)9)14-22(18-24)28(10,11)12/h13-18H,1-12H3/b30-29+

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