Archive for the 'Solvation' Category

QM/MM trajectory of an aqueous Diels-Alder reaction

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. 1988, 110, 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. 1994, 59, 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 2002, 106, 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. 2007, 3, 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.

Diels-Alder &Houk &Solvation Steven Bachrach 02 Feb 2016 1 Comment

Microsolvated structure of β-propiolactone

The structure of water about a solute remains of critical importance towards understanding aqueous solvation. Microwave spectroscopy and computations are the best tools we have today to gain insight on this problem. This is nicely demonstrated in the Alonso study of the microsolvated structures of β-propiolactone 1.1 They employed chirped-pulse Fourier transform microwave (CP-FTMW) spectroscopy and MP2(fc)/6-311++G(d,p) computations to examine the structure involving 1-5 water molecules.

The computed structures of these microsolvated species are shown in Figure 1. The deviation of the computed and experimental structures (RMS in the atomic positions) is small, though increasing as the size of the cluster increases. The deviation is 0.014 Å for the 1. H2O cluster and 0.244 Å for the 1.(H2O)5 cluster. They identified two clusters with four water molecules; the lower energy structure, labeled as a, is only 0.2 kJ mol-1 more stable than structure b.




1.(H2O)4 a

1.(H2O)4 b


Figure 1. MP2(fc)/6-311++G(d,p) optimized geometries of the hydrates of 1.

Water rings are found in the clusters having four or five water molecules, while chains are identified in the smaller clusters. One might imagine water cages appearing with even more water molecules in the microsolvated structures.


(1) Pérez, C.; Neill, J. L.; Muckle, M. T.; Zaleski, D. P.; Peña, I.; Lopez, J. C.; Alonso, J. L.; Pate, B. H. Angew. Chem. Int. Ed. 2015, 54, 979-982, DOI: 10.1002/anie.201409057.


1: InChI=1S/C3H4O2/c4-3-1-2-5-3/h1-2H2

MP &Solvation Steven Bachrach 24 Feb 2015 1 Comment

Computationally handling ion pairs

Comparing SN2 and SN1 reactions using computational methods is often quite difficult. The problem is that the heterolytic cleavage in the SN1 reaction leads to an ion pair, and in the gas phase this is a highly endothermic process. Optimization of the ion pair in the gas phase invariably leads to recombination. This is disturbingly the result even when one uses PCM to mimic the solvent, which one might have hoped would be sufficient to stabilize the ions.

The computational study of the glycoside cleavage by Hosoya and colleagues offers some guidance towards dealing with this problem.1 They examined the cleavage of triflate from 2,3,4,6-tetra-O-methyl-α-D-glucopyranosyl triflate 1.

Benchmarking the dissociation energy for the cleavage of 1 and considering computational performance, they settled on M06-2X/BS-III//M06-2X/BS-I, where BS-III is aug-cc-pVTZ basis set for O, F, and Cl and cc-pVTZ for H, C, and S and BS-I is 6-31G(d,p) basis sets were employed for H, C, and S, and 6-31+G(d) for O, F, and Cl. Solvent (dichloromethane) was included through PCM.

Attempted optimization of the contact ion pair formed from cleavage of 1 invariably led back to the covalent bound 1. PCM is not capable of properly stabilizing these types of ions in proximity. To solve this problem, they incorporated four explicit dichloromethane molecules. A minor drawback to their approach is that they did not sample much of configuration space and so their best geometries may not be the lowest energy configurations. Nonetheless, with four solvent molecules, they were able to locate contact ion pairs and solvent-separated ion pairs. Representative examples are shown in Figure 1. This method of explicit incorporation of a few solvent molecules seems to be the direction we must take to treat ions or even highly polar molecules in solution.





Figure 1. Representative examples of microsolvated 1, its contact ion pair (CIP) and solvent separated ion pair (SSIP) computed at M06-2X/BS-III//M06-2X/BS-I, and relative energies (kcal mol-1)


(1) Hosoya, T.; Takano, T.; Kosma, P.; Rosenau, T. "Theoretical Foundation for the Presence of Oxacarbenium Ions in Chemical Glycoside Synthesis," J. Org. Chem. 2014, 79, 7889-7894, DOI: 10.1021/jo501012s.


1: InChI=1S/C11H19F3O8S/c1-17-5-6-7(18-2)8(19-3)9(20-4)10(21-6)22-23(15,16)11(12,13)14/h6-10H,5H2,1-4H3/t6-,7-,8+,9-,10-/m1/s1

Ion Pairs &Solvation &sugars Steven Bachrach 04 Nov 2014 6 Comments

Benchmarked Dispersion corrected DFT and SM12

This is a short post mainly to bring to the reader’s attention a couple of recent JCTC papers.

The first is a benchmark study by Hujo and Grimme of the geometries produced by DFT computations that are corrected for dispersion.1 They use the S22 and S66 test sets that span a range of compounds expressing weak interactions. Of particular note is that the B3LYP-D3 method provided the best geometries, suggesting that this much (and justly) maligned functional can be significantly improved with just the simple D3 fix.

The second paper entails the description of Truhlar and Cramer’s latest iteration on their solvation model, namely SM12.2 The main change here is the use of Hirshfeld-based charges, which comprise their Charge Model 5 (CM5). The training set used to obtain the needed parameters is much larger than with previous versions and allows for treating a very broad set of solvents. Performance of the model is excellent.


(1) Hujo, W.; Grimme, S. "Performance of Non-Local and Atom-Pairwise Dispersion Corrections to DFT for Structural Parameters of Molecules with Noncovalent Interactions," J. Chem. Theor. Comput. 2013, 9, 308-315, DOI: 10.1021/ct300813c

(2) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. "Generalized Born Solvation Model SM12," J. Chem. Theor. Comput. 2013, 9, 609-620, DOI: 10.1021/ct300900e

Cramer &DFT &Grimme &Solvation &Truhlar Steven Bachrach 14 Jan 2013 No Comments

Proper use of computed solvation free energies

I missed this short communication last year, (thanks to the computational chemistry list CCL for bringing this to our attention!) but it is worth commenting on even a year later as this topic is one that frequently confuses users.

Ho, Klamt, and Coote1 note that popular quantum chemistry codes, including the Gaussian series, present the output of continuum solvent models in a way that can be misleading. What is called the free energy is in fact the sum of the electronic energy in solution and the free energy associated with non-electrostatic contributions. What is missing are corrections to the solute to give its free energy. What is assumed (oftentimes without fully recognizing this assumption) is that the thermal corrections for the solute in the gas and solution phase will cancel – but this does not have to be. Let the QM code user (and reader of the literature) beware!


(1) Ho, J.; Klamt, A.; Coote, M. L., "Comment on the Correct Use of Continuum Solvent Models," J. Phys. Chem. A 2010, 114, 13442-13444, DOI: 10.1021/jp107136j

Solvation Steven Bachrach 10 Jan 2012 3 Comments

Thorpe-Ingold Effect

Often gem-dialkyl substitution accelerates a reaction, for example in the formation of an epoxide via reaction 1. Here the relative rates are 1:21:252 in going from 1 to 2 to 3.1 This acceleration is the Thorpe-Ingold effect and had been suggested to arise from a steric reaction: that the methyl groups contract the angle and bring the terminal groups closer together.

1: R1 = R2 = H
2: R1 = Me, R2 = H
3: R1 = R2 = Me

Kostal and Jorgensen2 have examined the reaction of the 2-chloroethoxides 1-3 using computations, especially to look at the effect of solvent. At MP2/6-311+G(d,p) and CBS-Q, the relative rates (based on the activation free energy ΔG) are 1:2.8:17 and 1:0.7:3.7, respectively. Evidently there is no significant rate enhancement afforded by gem-substitution in the gas phase.

However, solution computations give a very different result. Using PCM along with the MP2 method, the computed relative rates are 1:5.8:1100 and with the Monte Carlo-Free Energy Perturbation method, the relative rates for aqueous solution are 1:30:773. Thus, the Thorpe-Ingold acceleration is due to solvent. Analysis of the hydrogen bonded structures and the solute-water pair distributions suggest that increasing alkyl substitution reduces the strength of solvation of the reactant, leading to the lower activation barrier.


(1) Jung, M. E.; Piizzi, G., "gem-Disubstituent Effect:Theoretical Basis and Synthetic Applications," Chem. Rev., 2005, 105, 1735-1766, DOI: 10.1021/cr940337h

(2) Kostal, J.; Jorgensen, W. L., "Thorpe-Ingold Acceleration of Oxirane Formation Is Mostly a Solvent Effect," J. Am. Chem. Soc., 2010, 132, 8766-8773, DOI: 10.1021/ja1023755

Jorgensen &Solvation Steven Bachrach 27 Jul 2010 2 Comments

Mannich reaction

Houk1 examined the Mannich reaction of the enamine formed from acetone and S-proline with N-ethylidine-N-phenylamine (see Chapter 5.3.3 in my book). Parasuk and Parasuk now extend this to the reaction of the enamine of cyclohexanone and S-proline with N-phenylmethanimine (Reaction 1).2 Geometries were optimized at B3LYP/6-31++G(d,p) and single-point energies computed with PCM (for the solvent DMSO) at both B3LYP and MP2.

Reaction 1

First, they examined the formation of the enamine 1, which can be in the syn or anti conformation. The barrier for formation of the syn isomer is 10.2 kcal mol-1. The barrier for the formation of the anti conformer is much higher, 17.9 kcal mol-1, and this is with a single water molecule used to assist the proton migration. However, the rotational barrier between the two conformers is only 4.2 kcal mol-1. So, they conclude that the syn isomer is the only conformer directly formed by the reaction of cyclohexanone and S-proline, and then rotation can produce the anti conformer.

The located the transition state for the reaction of either syn1 or anti1 with phenylmethanimine. The two transition states are shown in Figure 1. The barrier for the reaction of syn1 is 8.5 kcal mol-1, leading to the S product. The other barrier is higher, 13.0 kcal mol-1, and the R product 2R is 6.8 kcal mol-1 higher in energy than the S product 2S. Thus, the reaction to give the S product is both kinetically and thermodynamically favored. This is consistent with experiment3 which gives the S product with 99%ee. Inclusion of solvent makes the S product even more thermodynamically and kinetically favored over the R isomer.



Figure 1. B3LYP/6-311++G(d,p) optimized transition states leading to 2S and 2R.2


(1) Bahmanyar, S.; Houk, K. N., "Origins of Opposite Absolute Stereoselectivities in Proline-Catalyzed Direct Mannich and Aldol Reactions," Org. Lett. 2003, 5, 1249-1251, DOI: 10.1021/ol034198e.

(2) Parasuk, W.; Parasuk, V., "Theoretical Investigations on the Stereoselectivity of the Proline Catalyzed Mannich Reaction in DMSO," J. Org. Chem. 2008, 73, 9388-9392, DOI: 10.1021/jo801872w.

(3) Ibrahem, I.; Zou, W.; Casas, J.; Sundén, H.; Córdova, A., "Direct organocatalytic enantioselective α-aminomethylation of ketones," Tetrahedron 2006, 62, 357-364, DOI: 10.1016/j.tet.2005.08.113.


1: InChI=1/C11H17NO2/c13-11(14)10-7-4-8-12(10)9-5-2-1-3-6-9/h5,10H,1-4,6-8H2,(H,13,14)/t10-/m0/s1/f/h13

2S: InChI=1/C18H24N2O2/c21-18(22)16-10-6-11-17(16)20-12-5-4-9-15(20)13-19-14-7-2-1-3-8-14/h1-3,7-8,15-16,19H,4-6,9-13H2/b20-17+/t15-,16+/m0/s1

2R: InChI=1/C18H24N2O2/c21-18(22)16-10-6-11-17(16)20-12-5-4-9-15(20)13-19-14-7-2-1-3-8-14/h1-3,7-8,15-16,19H,4-6,9-13H2/b20-17-/t15-,16-/m1/s1

Mannich &Solvation Steven Bachrach 04 Jun 2009 No Comments

Protonation of 4-aminobenzoic acid

Molecular structures can differ depending on phase, particularly between the gas and solution phase. Kass has looked at the protonation of 4-aminobenzoic acid. In water, the amino is its most basic site, but what is it in the gas phase? The computed relative energies of the protonation sites are listed in Table 1. If one corrects the B3LYP values for their errors in predicting the proton affinity of aniline and benzoic acid, the carbonyl oxygen is predicted to be the most basic site by 5.0 kcal mol-1, in nice accord with the G3 prediction of 4.1 kcal mol-1. Clearly, the structure depends on the medium.

Table 1. Computed relative proton affinities (kcal mol-1) of 4-aminobenzoic acid.

C=O 0.0 0.0
NH2 7.9 4.1
OH 12.2 9.8

Electrospray of 4-aminobenzoic acid from 3:1 methanol/water and 1:1 acetonitrile/water solutions gave different CID spectra. H/D exchange confirmed that electrospray from the emthanol/water solution gave the oxygen protonated species while that from the acetonitrile/water solution gave the ammonium species.


(1) Tian, Z.; Kass, S. R., “Gas-Phase versus Liquid-Phase Structures by Electrospray Ionization Mass Spectrometry,” Angew. Chem. Int. Ed., 2009, 48, 1321-1323, DOI: 10.1002/anie.200805392.


4-aminobenzoic acid: InChI=1/C7H7NO2/c8-6-3-1-5(2-4-6)7(9)10/h1-4H,8H2,(H,9,10)/f/h9H

Acidity &amino acids &Kass &Solvation Steven Bachrach 30 Mar 2009 No Comments

Solubility in olive oil

Here’s a nice example of the application of computed solvation energies in non-aqueous studies. Cramer and Truhlar have employed their latest SM8 technique, which is parameterized for organic solvents and for water, to estimate solvation energies in olive oil.1 Now you may wonder why solvation in olive oil of all things? But the partitioning of molecules between water and olive oil has been shown to be a good predictor of lipophilicity and therefore bioavailability of drugs! The model works reasonably well in reproducing experimental solvation energies and partition coefficients. They do make the case that fluorine substitution which appears to improve solubility in organics,originates not to more favorable solvation in organic solvents (like olive oil) but rather that fluorine substitution dramatically decreases solubility in water.


(1) Chamberlin, A. C.; Levitt, D. G.; Cramer, C. J.; Truhlar, D. G., "Modeling Free Energies of Solvation in Olive Oil," Mol. Pharmaceutics, 2008, 5, 1064-1079, DOI: 10.1021/mp800059u

Cramer &Solvation &Truhlar Steven Bachrach 17 Feb 2009 1 Comment

Arginine:water cluster

The gas phase structure of the amino acids is in their canonical or neutral form, while their aqueous solution phase structure is zwitterionic. An interesting question is how many water molecules are needed to make the zwitterionic structure more energetically favorable than the neutral form. For glycine, it appears that seven water molecules are needed to make the zwitterion the favorable tautomer.1,2

Arginine, on the other hand, appears to require only one water molecule to make the zwitterion lower in energy than the neutral form.3 The B3LYP/6-311++G** structures of the lowest energy neutral (1N) and zwitterion (1Z) cluster with one water are shown in Figure 1. The zwitterion is 1.68 kcal mol-1 lower in energy. What makes this zwitterion so favorable is that the protonation occurs on the guanidine group, not on the amine group. The guanidine group is more basic than the amine. Further, the water can accept a proton from both nitrogens of the guanidine and donate a proton to the carboxylate group.

1N (1.68)

1Z (0.0)

Figure 1. B3LYP/6-311++G** structures and relative energies (kcal mol-1) of the lowest energy arginine neutral (1N) and zwitterion (1Z) cluster with one water.3


(1) Aikens, C. M.; Gordon, M. S., "Incremental Solvation of Nonionized and Zwitterionic Glycine," J. Am. Chem. Soc., 2006, 128, 12835-12850, DOI: 10.1021/ja062842p.

(2) Bachrach, S. M., "Microsolvation of Glycine: A DFT Study," j. Phys. Chem. A, 2008, 112, 3722-3730, DOI: 10.1021/jp711048c.

(3) Im, S.; Jang, S.-W.; Lee, S.; Lee, Y.; Kim, B., "Arginine Zwitterion is More Stable than the Canonical Form when Solvated by a Water Molecule," J. Phys. Chem. A, 2008, 112, 9767-9770, DOI: 10.1021/jp801933y.


1: InChI=1/C6H14N4O2/c7-4(5(11)12)2-1-3-10-6(8)9/h4H,1-3,7H2,(H,11,12)(H4,8,9,10)/f/h8,10-11H,9H2

amino acids &Solvation Steven Bachrach 15 Dec 2008 1 Comment

Next Page »