ICYMI: Ari and I announced our new company, Rowan! We wrote an article about what we're hoping to build, which you can read here. Also, this blog is now listed on The Rogue Scholar, meaning that posts have DOIs and can be easily cited.
Conventional quantum chemical computations operate on a collection of atoms and create a single wavefunction for the entire system, with an associated energy and possibly other properties. This is great, but sometimes we want to understand things in more detail. For instance, if we have a host A and two guests Bgood and Bbad, a normal calculation would just tell us that E(A•Bgood) is lower than E(A•Bbad), without giving any clue as to why.
Enter EDA. EDA, or “energy decomposition analysis,” is a family of techniques used to dissect interactions in a system with multiple molecules. In this case, an EDA calculation on the AB system would break down the interaction between A and B into various components, which could be used to help scientists understand the origin of the difference, and perhaps used for continued molecular design.
Unfortunately, EDA has always seemed like a pretty troubled technique to me. Wavefunctions are inherently not localized to individual fragments of a multimolecular system—you can’t just slice apart the molecular orbitals or the density matrix and end up with anything that’s physically sane. So you have to do some computational gymnastics to get energetic terms which are at all meaningful. Many such gymnastic workflows have been proposed, leading to a veritable alphabet soup of different EDA methods.
(It’s worth skimming this review on different EDA methods to get a sense for some of the questions the field faces, and also to laugh at how Alston Misquitta & Krzysztof Szalewicz use the review as a chance to relentlessly advertise SAPT and denigrate any and all competing methods.)
I’ll briefly outline how the EDA-NCOV method works for a system AB (following this review), to give a sense for the flavor of the field:
Optimized ground-state fragments A0 and B0 are distorted to the geometries and electronic states (A & B) which they possess in AB, and the energy required for this distortion/excitation is termed Eprep. (The difference between E(AB) and E(A) + E(B) is called Eint, and the total binding energy is equal to Eint + Eprep.)
The distorted fragments A and B are brought together (with frozen charge densities) to form the “promolecule” AB0, and the change in energy is termed Eelstat, the quasiclassical Coulomb interaction energy (typically attractive). The wavefunction for AB0 is ΨAΨB.
The product wavefunction ΨAΨB is antisymmetrized and renormalized to give an “intermediate state” Ψ0 with energy E0, and the change in energy is termed EPauli, originating from Pauli repulsion. This component is always repulsive.
Ψ0 is relaxed to yield the final wavefunction ΨAB. The change in energy is termed Eorb, because it arises from orbital interactions, and is always attractive.
Thus, Eint = Eelstat + EPauli + Eorb. (Dispersion can also be added if an exogenous dispersion correction is employed—that’s pretty trivial.)
The critical reader might observe that the steps taken to obtain these numbers are pretty odd, and that the components of the interaction energy arise from differences in energy between bizarre nonphysical states. Thus, the interpretation of terms like Eelstat in terms of actual physical interactions might not be as easy as it seems. The authors of the above review agree:
It is important to realize that the identification of the three major terms ΔEelstat, ΔEPauli, and ΔEorb with specific interactions is conceptually attractive but must not be taken as genuine expression of the physical forces.
Unfortunately, it seems that imprecise concepts familiar to experimental chemists like “steric repulsion” and “electrostatic attraction” have to be discarded in favor of precise terms like EPauli. Too bad they’re virtually uninterpretable!
And what’s worse is that different EDA-type schemes don’t even give the same results. A paper out today in JACS from Zare/Shaik discusses the use of EDA and related schemes in studying the origin of the hydrogen bond (a pretty fundamental question), motivated by the substantial disagreement between various techniques:
It is important to realize that different methods (e.g., BOVB, ALMO-EDA, NEDA, and BLW) do not fully agree with one another about whether the dominant stabilizing term is ΔEPOL or ΔECT in a particular HB.
While the authors make a good case that the sum of these two terms is relatively conserved across methods, and that it’s this term that we should care about for hydrogen bonds, the conclusions for EDA broadly are not encouraging. (Note, too, that EPOL and ECT don’t even appear in the EDA-NCOV method summarized above—another reason that EDA is a frustrating field!)
And even if the theorists eventually put their heads together and develop a version of EDA that doesn’t have these pitfalls, it’s still not clear that any form of EDA will give the answers that experimental chemists are looking for. Chemistry is complicated, and ground- or transition-state structures arise from a delicate equilibrium between opposing factors: steric repulsion, electrostatic attraction, bond distances, torsional strain, dispersion, &c.
As a result, one can see large changes in the contribution of individual factors even while the overall structure’s stability is minimally perturbed (enthalpy–entropy compensation is a classic example, as is Fig. 2 in this review on distortion–interaction analysis). Looking only at changes in individual factors isn’t always a useful way to gain insight from computation.
For example, imagine a nucleophile adding to two faces of an oxocarbenium, a bulky face and an unhindered face. Based on this description, we might expect to see that TSbulky has higher steric repulsion than TSunhindered (if we’re lucky enough to find a way to extract Esteric out of our EDA method). But it’s also likely that the nucleophile might take a less favorable trajectory towards the oxocarbenium in TSbulky to avoid steric repulsion, thus weakening key orbital interactions. These changes might even end up being larger in magnitude than the destabilization induced by steric repulsion. Is the correct answer, then, that TSbulky is higher in energy because of decreased Eorb, not increased Esteric?
The solution is to recognize that causation is not unique (cf. Aristotle), and so there’s no one right answer here. Within the constraints of the EDA framework, the theorist wouldn’t be incorrect in saying that Eorb is the driving factor—but the experimental chemist might reasonably expect “the bulky TS is destabilized by steric repulsion” as their answer, since this is the root cause of the changes between the two structures. (I side with the experimentalists here.)
And the precisely defined concepts favored by theorists are often hard for experimental scientists to work with. Even if the correct answer in the above scenario were “TSbulky is destabilized by decreased orbital overlap”—what’s an experimentalist supposed to do with this information, add more orbitals? (This is how I feel about Trevor Hamlin’s work on Pauli repulsion.) The steric explanation at least suggests an intuitive solution: make the bulky group or the nucleophile smaller. If the purpose of EDA is to help people to understand intermolecular interactions better on a conceptual level, I’m not sure it’s succeeding in most cases.
(The only use of EDA that led to an actual experimental advance which I’m aware of is Buchwald/Peng Liu’s body of work on ligand–substrate dispersion in hydrocupration: study, new ligand, ligand from Hartwig. I don’t think it’s a coincidence that these papers focus on dispersion, one of the easiest pieces of EDA to decouple and understand.)
I don’t mean to be too critical here. The ability to break intermolecular interactions down into different components is certainly useful, and it seems likely that some version of EDA will eventually achieve consensus and emerge as a useful tool. But I think the utility of EDA even in the best case is pretty limited. Quantum chemistry is complicated, and if we think we can break it down into easy-to-digest components and eliminate the full nonlocal majesty of the Schrodinger equation, we’re lying to ourselves (or our experimental collaborators). Compute with caution!
Recently in off-blog content: I coauthored an article with Eric Gilliam on how LLMs can assist in hypothesis generation and help us all become more like Sharpless. Check it out!
The Pauling model for enzymatic catalysis states that enzymes are “antibodies for the transition state”—in other words, they preferentially bind to the transition state of a given reaction, rather than the reactants or products. This binding interaction stabilizes the TS, thus lowering its energy and accelerating the reaction.
This is a pretty intuitive model, and one which is often employed when thinking about organocatalysis, particularly noncovalent organocatalysis. (It’s a bit harder to use this model when the mechanism changes in a fundamental way, as with many organometallic reactions.) Many transition states have distinctive features, and it’s fun to think about what interactions could be engineered to recognize and stabilize only these features and nothing else in the reaction mixture.
(For instance, chymotrypsin contains a dual hydrogen-bond motif called the “oxyanion hole” which stabilizes developing negative charge in the Burgi–Dunitz TS for alcoholysis of amides. The negative charge is unique to the tetrahedral intermediate and the high-energy TSs to either side, so reactant/product inhibition isn’t a big issue. This motif can be mimicked by dual hydrogen-bond donor organocatalysts like the one my PhD lab specialized in.)
The downside of this approach to catalyst design is that each new sort of reaction mechanism requires a different sort of catalyst. One TS features increasing negative charge at one place, while another features increasing positive charge at another, and a third is practically charge-neutral the whole way through. What if there were some feature that was common to all transition states?
A recent preprint from Diptarka Hait and Martin Head-Gordon suggests an interesting answer to this question: polarizability. (Diptarka, despite just finishing his PhD, is a prolific scientist with a ton of different interests, and definitely someone to keep an eye on.) The authors tackle the question of when precisely a stretching bond can be considered “broken.” An intuitive answer might be “the transition state,” but as the paper points out, plenty of bond-breaking potential energy surfaces lack clear transition states (e.g. H2, pictured below in red).
Hydrogen PES (red), with various polarizability components shown in blue.
Instead, the authors propose that polarizability is a good way to study this question. As seen in the following graph, polarizability (particularly α||, the component parallel to the bond axis) first increases as a bond is stretched, and then decreases, with a sharp and easily identifiable maximum about where a bond might be considered to be broken. This metric tracks with the conventional PES metric in cases where the PES is well-defined (see Fig. 5), which is comforting. Why does this occur?
The evolution of α [polarizability] can be rationalized in the following manner. Upon initially stretching H2 from equilibrium, the bonding electrons fill up the additional accessible volume, resulting in a more diffuse (and thus more polarizable) electron density. Post bond cleavage however, the electrons start localizing on individual atoms, leading to a decrease in polarizability.
In other words, electrons that are “caught in the act” of reorganizing between different atoms are more polarizable, whereas electrons which have settled into their new atomic residences are less polarizable again.
This is cool, but how can we apply this to catalysis? As it happens, there are already a few publications (1, 2) from Dennis Dougherty and coworkers dealing with exactly this question. They show that cyclophanes are potent catalysts for SN2-type reactions that both create and destroy cations, and argue that polarizability, rather than any charge-recognition effect, undergirds the observed catalysis:
Since transition states have long, weak bonds, they are expected to be more polarizable than ground-state substrates or products. The role of the host is to surround the transition state with an electron-rich-system that is polarizable and in a relatively fixed orientation, so that induced dipoles in both the host and the transition state are suitably aligned… Note that in this model, it is not sufficient that polarizability contributes to binding. Polarizability must be more important for binding transition states than for binding ground states. Only in this way can it enhance catalysis.
To defend this argument, the authors prepare a series of variously substituted cyclophanes, and show that while Hammett-type electronic tuning of the aryl rings has relatively small effects, adding more heavy atoms always increases the rate, with the biggest effect observed with bromine (the heaviest element attempted). Heavier atoms are more polarizable, so this supports the argument that polarizability, rather than any specific electrostatic effect, is responsible for catalysis.
The cyclophanes used as catalysts: Y=Br is the best.
The Dougherty work is performed on a very specific model system, and the absolute rate accelerations seen aren’t massive (about twofold increase relative to the protio analog), so it’s not clear that this will actually be a promising avenue for developing mechanism-agnostic catalysts.
But I think this line of research is really interesting, and merits further investigation. Pericyclic reactions, which involve large electron clouds and multiple forming–breaking bonds and often feature minimal development of partial charges, seem promising targets for this sort of catalysis—what about the carbonyl–ene reaction or something similar? The promise of new catalytic strategies that complement existing mechanism-specific interactions is just too powerful to leave unstudied.
Thanks to Joe Gair for originally showing me the Dougherty papers referenced.
Since the ostensible purpose of organic methodology is to develop reactions that are useful in the real world, the utility of a method is in large part dictated by the accessibility of the starting materials. If a compound is difficult to synthesize or hazardous to work with, then it’s difficult to convince people to use it in a reaction (e.g. most diazoalkanes). Organic chemists are pragmatic, and would usually prefer to run a reaction that starts from a commercial and bench-stable starting material.
For instance, this explains the immense popularity of the Suzuki reaction: although the Neigishi reaction (using organozinc nucleophiles) usually works better for the same substrates, you can buy lots of the organoboron nucleophiles needed to run a Suzuki and leave them lying around without taking any precautions. In contrast, organozinc compounds usually have to be made from the corresponding organolithium/Grignard reagent and used freshly, which is considerably more annoying.
The ideal starting material, then, is one which is commercially available and cheap. In recent years, it’s become popular to advertise new synthetic methods by showing that they work on exceptionally cheap and common functional groups, and in particular to compare the abundance of different functional groups to demonstrate that one starting material is more common than another. To pick just one of many examples, Dave MacMillan used this plot to show why cross-coupling reactions of alcohols were important (ref):
This visual works really well.
When I saw MacMillan’s talk at MIT last year, I was curious what it would take to make additional graphics like this. The “number of reactions” plot can be made pretty easily from Reaxys, but I’ve always been uncertain how the “number of commercial sources” plots are made: I haven’t seen references listed for these numbers, nor is anything usually found in the Supporting Information.
I decided to take a swing at getting this data myself by analyzing the Mcule "building blocks" database, which contains about 3.5 million compounds. Although Mcule doesn't define what a building block is (at least, not that I can find), it’s likely that their definition is similar to that of ZINC, which defines building blocks as “those catalogs of compounds available in preparative quantities, typically 250 mg or more” (ref). This seems like a reasonable proxy for the sorts of compounds synthetic chemists might use in reactions. I defined patterns to match a bunch of functional groups using SMARTS/SMILES, and then used RDKit to find matches in the Mcule building blocks database. The code can be found on Github, along with the patterns I used.
The results are shown below. As expected, ethers, amines, amides, and alcohols are quite common. Surprisingly, aryl chlorides aren't that much more common than aryl bromides—and, except for aliphatic fluorides, all aliphatic halides are quite rare. Allenes, carbodiimides, and SF5 groups are virtually unheard of (<100 examples).
Functional Group
Number
Percent
acid chloride
6913
0.19
alcohol
1022229
28.60
aliphatic bromide
42018
1.18
aliphatic chloride
70410
1.97
aliphatic fluoride
650576
18.20
aliphatic iodide
3159
0.09
alkene
176484
4.94
alkyne
35577
1.00
allene
99
0.00
amide
518151
14.50
anhydride
1279
0.04
aryl bromide
451451
12.63
aryl chloride
661591
18.51
aryl fluoride
618620
17.31
aryl iodide
216723
6.06
azide
5164
0.14
aziridine
748
0.02
carbamate
127103
3.56
carbodiimide
28
0.00
carbonate
1231
0.03
carboxylic acid
410860
11.49
chloroformate
250
0.01
cyclobutane
195728
5.48
cyclopropane
349455
9.78
diene
10188
0.29
difluoromethyl
163395
4.57
epoxide
5859
0.16
ester
422715
11.83
ether
1434485
40.13
isocyanate
1440
0.04
isothiocyanate
1389
0.04
nitrile
209183
5.85
nitro
126200
3.53
pentafluorosulfanyl
18
0.00
primary amine
904118
25.29
secondary amine
857290
23.98
tertiary amine
609261
17.04
trifluoromethoxy
18567
0.52
trifluoromethyl
455348
12.74
urea
518151
14.50
Total
3574611
100.00
(Fair warning: I’ve spotchecked a number of the SMILES files generated (also on Github), but I haven’t looked through every molecule, so it’s possible that there are some faulty matches. I wouldn’t consider these publication-quality numbers yet.)
An obvious caveat: there are lots of commercially “rare” functional groups which are easily accessible from more abundant functional groups. For instance, acid chlorides seem uncommon in the above table, but can usually be made from ubiquitous carboxylic acids with e.g. SOCl2. So these data shouldn’t be taken as a proxy for a more holistic measure of synthetic accessibility—they measure commercial availability, that’s all.
What conclusions can we draw from this?
The most common functional groups are the milquetoast ones: alcohols, amines, esters, etc. Perhaps this explains where all the new reactions have gone: unless your new method works on alcohols or amines, it will struggle to get traction in most of chemical space relative to e.g. Williamson ether synthesis or reductive amination. (Kudos to MacMillan for identifying this; vide supra.)
Ureas are much more common than you’d expect from academic methods papers. This I think speaks to the difference between what methodologists want and what medicinal chemists want. Ureas are a bit annoying to work with: they’re pretty polar by the standards of academia, they’re not always soluble in organic solvents, and they have a tendency to stick to transition metal catalysts or get deprotonated by strong bases. But they’re easy to make in libraries, since the isocyanate/amine disconnection is so robust, and they’re excellent hydrogen-bond donors and acceptors.
CORRECTION: There's a SMARTS error, so the match for "ureas" actually matches amides—disregard this section. Thanks to @wmdhn for catching this.
Uncommon functional groups, like SF5 and allenes, are very uncommon. If you want to introduce an SF5 group, you are in for a rough time: there aren’t great ways to add it to molecules (although there have been some steps forward in recent years), and there are only 18 commercial examples. So people can write as many papers as they want about how cool SF5 groups are: I still doubt we’ll see them used very much in the near future.
But also, the abundance of a given functional group is very elastic in the long run. Trifluoromethyl groups used to be extremely rare—they’re not found in nature!—but now 1 in 8 molecules has a CF3 group. CF3 just turns out to be a very good handle for a lot of molecular design tasks, and so people found ways to introduce it all over the place, and now it’s not hard to get molecules that have trifluoromethyl groups. Synthetic chemists should feel good about this.
The functional-group-specific SMILES files are in the previously mentioned Github repo, so anyone who wants to e.g. look through all the commercially available alkenes and perform further cheminformatics analyses can do so. I hope the attached code and data helps other chemists perform similar, and better, studies, and that this sort of thinking can be useful for those who are currently engaged in reaction discovery.
Thanks to Eric Jacobsen for helpful conversations about these data.
In his fantastic essay “The Two Cultures,” C. P. Snow observed that there was (in 1950s England) a growing divide between the academic cultures of science and the humanities:
Literary intellectuals at one pole—at the other scientists, and as the most representative, the physical scientists. Between the two a gulf of mutual incomprehension—sometimes (particularly among the young) hostility and dislike, but most of all lack of understanding. They have a curious distorted image of each other. Their attitudes are so different that, even on the level of emotion, they can't find much common ground.
He reflects on the origins of this phenomenon, which he contends is new to the 20th century, and argues that it ought to be opposed:
This polarisation is sheer loss to us all. To us as people, and to our society. It is at the same time practical and intellectual and creative loss, and I repeat that it is false to imagine that those three considerations are clearly separable. But for a moment I want to concentrate on the intellectual loss.
Snow’s essay is wonderful: his portrait of a vanishing cultural intellectual unity should inspire us all, scientists and otherwise, to improve ourselves, and the elegiac prose reminds the reader that even the best cultural institutions are fragile and fleeting things.
I want to make an analogous—but much less powerful—observation about the two cultures present in atomistic simulation. I’ll call these the “QM tribe” and the “MD tribe” for convenience: crudely, “people who use Gaussian/ORCA/Psi4 for their research” and “people who use Schrodinger/AMBER/OpenMM/LAMMPS for their research,” respectively. Although this dichotomy is crude, I contend there are real differences between these two groups, and that their disunity hurts scientific progress.
The Nature of Energy Surfaces
The most fundamental disagreement between these two cultures is in how they think about energy surfaces, I think. Most QM-tribe people think in terms of optimizing to discrete critical points on the potential energy surface: one can perform some sort of gradient-informed optimization to a ground state, or follow negative eigenvalues to a transition state.
Implicit to this assumption is that there exist well-defined critical points on the PES, and that finding such critical points is meaningful and productive. Conformers exist, and many people now compute properties as Boltzmann-weighted averages over conformational ensembles, but this is usually done for 10–100 conformers, not thousands or millions. Entropy and solvation, if they’re considered at all, are viewed as corrections, not key factors: since QM is so frequently used to study high-barrier bond-breaking processes where enthalpic factors dominate, one can often get reasonable results with cartoonish treatments of entropy.
In contrast, MD-tribe scientists generally don’t think about transition states as specific configurations of atoms—rather, a transition state can emerge from some sort of simulation involving biased sampling, but it’s just a position along some abstract reaction coordinate, rather than a structure which can be visualized in CYLView. Any information gleaned is statistical, rather than concretely visual (e.g. “what is the mean number of hydrogen bonds to this oxygen near this transition state”).
Unlike the QM tribe, MD-tribe scientists generally cannot study bond-breaking processes, and so focus on conformational processes (protein folding, aggregation, nucleation, transport) where entropy and solvation are of critical importance: as such, free energy is almost always taken into consideration by MD-tribe scientists, and the underlying PES itself is rarely (to my knowledge) viewed as a worthy topic of study in and of itself.
Molecular Representations
This divide also affects how the two cultures view the task of molecular representation. QM-tribe scientists generally view a list of coordinates and atomic numbers as the most logical representation of a molecule (perhaps with charge and multiplicity information). To the QM tribe, a minimum on the PES represents a structure, and different minima naturally ought to have different representations. Bonding and bond order are not specified, because QM methods can figure that out without assistance (and it’s not uncommon for bonds to change in a QM simulation anyway).
In contrast, people in the MD tribe generally want a molecular representation that’s independent of conformation, since many different conformations will intrinsically be considered. (See Connor Coley’s presentation from a recent MolSSI workshop for a discussion of this.) Thus, it’s common to represent molecules through their topology, where connectivity and bond order are explicitly specified. This allows for some pretty wild simulations of species that would be reactive in a QM simulation, but also means that e.g. tautomers can be a massive problem in MD (ref), since protons can’t equilibrate freely.
For property prediction, an uneasy compromise can be reached wherein one takes a SMILES string, performs a conformational search, and then Boltzmann-averages properties over all different conformers, to return a set of values which are associated only with the SMILES string and not any individual conformation. (Matt Sigman does this, as does Bobby Paton for NMR prediction.) This is a lot of work, though.
“A Gulf Of Mutual Incomprehension”
These differences also become apparent when comparing the software packages that different tribes use. Take, for instance, the task of predicting partial charges for a given small molecule. A QM-tribe scientist would expect these charges to be a function of the geometry, whereas an MD-tribe scientist would want the results to be explicitly geometry-independent (e.g.) so that they can be used for subsequent MD simulations.
The assumptions implicit to these worldviews mean that it’s often quite difficult to go from QM-tribe software packages to MD-tribe software packages or vice versa. I’ve been stymied before by trying to get OpenMM and openforcefield to work on organic molecules for which I had a list of coordinates and not e.g. a SMILES string—although obviously coordinates will at some point be needed in the MD simulation, most workflows expect you to start from a topology and not an xyz file.
Similarly, it’s very difficult to get the graphics package NGLView to illustrate the process of bonds breaking and forming—NGLView is typically used for MD, and expects that the system’s topology will be defined at the start of the simulation and never changed. (There are kludgy workarounds, like defining a new object for every frame, but it’s nevertheless true that NGLView is not made for QM-tribe people.)
(I’m sure that MD-tribe people are very frustrated by QM software as well, but I don’t have as much experience going in this direction. In general, MD tooling seems quite a bit more advanced than QM-tribe tooling; most MD people I’ve talked to seem to interact with QM software as little as possible, and I can’t say I blame them.)
“A Curious Distorted Image of Each Other”
There are also cultural factors to consider here. The questions that QM-tribe scientists think about are different than those that MD-tribe scientists think about: a somewhat famous QM expert once told me that they were “stuck on an ivory tower where people hold their nose when it comes to DFT, forget anything more approximate,” whereas MD-tribe scientists often seem alarmingly unconcerned about forcefield error but are obsessed with proper sampling and simulation convergence.
It seems that most people have only a vague sense of what their congeners in the other tribe actually work on. I don’t think most QM-tribe scientists I know have ever run or analyzed a regular molecular dynamics simulation using e.g. AMBER or OpenMM, nor do I expect that most MD-tribe scientists have tried to find a transition state in Gaussian or ORCA. In theory, coursework could remedy this, but education for QM alone already seems chaotic and ad hoc—trying to cram in MD, statistical mechanics, etc is probably ill-advised at the present.
Social considerations also play a role. There’s limited crosstalk between the two fields, especially at the trainee level. How many QM people even know who Prayush Tiwary is, or Michael Shirts, or Jay Ponder? How many MD graduate students have heard of Frank Neese or Troy Van Voorhis? As always, generational talent manages to transcend narrow boundaries—but rank-and-file scientists would benefit immensely from increased contact with the other tribe.
“Unite Them”
I’m not an expert on the history of chemistry, but my understanding is that the two fields were not always so different: Martin Karplus, Arieh Warshel, and Bill Jorgensen, key figures in the development of modern MD, were also formidable quantum chemists. (If any famous chemists who read this blog care to share their thoughts on this history, please email me: you know who you are!)
And as the two fields advance, I think they will come closer together once more. As QM becomes capable of tackling larger and larger systems, QM-tribe scientists will be forced to deal with more and more complicated conformational landscapes: modern enantioselective catalysts routinely have hundreds of ground-state complexes to consider (ref), and QSimulate and Amgen recently reported full DFT calculations on protein–ligand complexes (ref).
Similarly, the increase in computing power means that many MD use cases (like FEP) are now limited not by insufficient sampling but by the poor energetics of the forcefields they employ. This is difficult to prove unequivocally, but I’ve heard this in interviews with industry folks, and there are certainly plenty of references complaining about poor forcefield accuracy (1, 2): a Psivant review dryly notes that “historically solvation energy errors on the order of 2–3 kcal/mol have been considered to be accurate,” which is hardly encouraging.
Many QM-tribe professors now work on dynamics: Dean Tantillo and Todd Martinez (who have long been voices “crying out in the wilderness” for dynamics) perhaps most prominently, but also Steven Lopez, Daniel Ess, Fernanda Duarte, Peng Liu, etc. And MD-tribe professors seem more and more interested in using ML mimics of QM to replace forcefields (e.g.), which will inevitably lead them down the speed–accuracy rabbit hole that is quantum chemistry. So it seems likely to me that the two fields will increasingly reunite, and that being a good 21st-century computational chemist will require competency in both areas.
If this is true, the conclusions for individual computational chemists are obvious: learn techniques outside your specialty, before you get forcibly dragged along by the current of scientific progress! There’s plenty to learn from the other culture of people that deals with more-or-less the same scientific problems you do, and no reason to wait.
As a denizen of quantum chemistry myself, I apologize for any misrepresentations or harmful stereotypes about practitioners of molecular dynamics, for whom I have only love and respect. I would be happy to hear any corrections over email.
An important problem with simulating chemical reactions is that reactions generally take place in solvent, but most simulations are run without solvent molecules. This is a big deal, since much of the inaccuracy associated with simulation actually stems from poor treatment of solvation: when gas phase experimental data is compared to computations, the results are often quite good.
Why don’t computational chemists include solvent molecules in their models? It takes a lot of solvent molecules to accurately mimic bulk solvent (enough to cover the system with a few different layers, usually ~103).1 Since most quantum chemical methods scale in practice as O(N2)–O(N3), adding hundreds of additional atoms has a catastrophic effect on the speed of the simulation.
To make matters worse, the additional degrees of freedom introduced by the solvent molecules are very “flat”—solvent molecules don’t usually have well-defined positions about the substrate, meaning that the number of energetically accessible conformations goes to infinity (with attendant consequences for entropy). This necessitates a fundamental change in how calculations are performed: instead of finding well-defined extrema on the electronic potential energy surface (ground states or transition states), molecular dynamics (MD) or Monte Carlo simulations must be used to sample from an underlying distribution of structures and reconstruct the free energy surface. Sufficient sampling usually requires consideration of 104–106 individual structures,2 meaning that each individual computation must be very fast (which is challenging for quantum chemical methods).
The title of this paper makes me so sad, because these techniques are still ignored by most organic chemists.
Given the complexity this introduces, it’s not surprising that most computational organic chemists try to avoid explicit solvent at all costs. The typical workaround is to use “implicit solvent” models, which “reduce the complexity of individual solvent−solute interactions such as hydrogen-bond, dipole−dipole, and van der Waals interactions into a fictitious surface potential... scaled to reproduce the experimental solvation free energies” (Baik). This preserves the well-defined potential energy surfaces that organic chemists are accustomed to, so you can still find transition states by eigenvector following, etc.
Implicit solvent models like PCM, COSMO, or SMD are better than nothing, but are known to struggle for charged species. In particular, they don’t really describe explicit inner-sphere solvent–solute interactions (like hydrogen bonding), meaning that they’ll behave poorly when these interactions are important. Dan Singleton’s paper on the Baylis–Hillman reaction is a nice case study of how badly implicit solvent can fail: even high-level quantum chemical methods are useless when solvation free energies are 10 kcal/mol off from experiment!
An even more important but still open issue is solvation. In the opinion of the authors it is a ‘scandal’ that in 2018 no routine free energy solvation method is available beyond (moderately successful) continuum theories such as COSMO-RS and SMD and classical FF/MD-based explicit treatments.
When computational studies have been performed in explicit solvent, the results have often been promising: Singleton has studied diene hydrochlorination and nitration of toluene, and Peng Liu has recently conducted a nice study of chemical glycosylation. Nevertheless, these studies all require heroic levels of effort: quantum chemistry is still very slow, and so a single free energy surface might take months and months to compute.3
One promising workaround is using machine learning to accelerate quantum chemistry. Since these MD-type studies look at the same exact system over and over again, we could imagine first training some sort of ML model based on high-level quantum chemistry data, and then employing this model over and over again for the actual MD run. As long as (1) the ML model is faster than the QM method used to train it and (2) it takes less data to train the ML model than it would to run the simulation, this will save time: in most cases, a lot of time.
(This is a somewhat different use case than e.g. ANI-type models, which aim to achieve decent accuracy for any organic molecule. Here, we already know what system we want to study, and we’re willing to do some training up front.)
A lot of people are working in this field right now, but today I want to highlight some work that I liked from Fernanda Duarte and co-workers. Last year, they published a paper comparing a few different ML methods for studying quasiclassical dynamics (in the gas phase), and found that atomic cluster expansion (ACE) performed better than Gaussian approximation potentials while training faster than NequIP. They then went on to show that ACE models could be trained automatically through active learning, and used the models to successfully predict product ratios for cycloadditions with post-TS bifurcations.
Their new paper, posted on ChemRxiv yesterday, applies the same ACE/active learning approach to studying reactions in explicit solvent, with the reaction of cyclopentadiene and methyl vinyl ketone chosen as a model system. This is more challenging than their previous work, because the ML model now not only has to recapitulate the solute reactivity but also the solute–solvent and solvent–solvent interactions. To try and capture all the different interactions efficiently, the authors ended up using four different sets of training data: substrates only, substrates with 2 solvent molecules, substrates with 33 solvent molecules, and clusters of solvent only.
Previously, the authors used an energy-based selector to determine if a structure should be added to the training set: they predicted the energy with the model, ran a QM calculation, and selected the structure if the difference between the two values was big enough. This approach makes a lot of sense, but has the unfortunate downside that a lot of QM calculations are needed, which is exactly what this ML-based approach is trying to avoid. Here, the authors found that they could use similarity-based descriptors to select data points to add to the training set: these descriptors are both more efficient (needing fewer structures to converge) and faster to compute, making them overall a much better choice. (This approach is reminiscent of the metadynamics-based approach previously reported by John Parkhill and co-workers.)
With a properly trained model in hand, the authors went on to study the reaction with biased sampling MD. They find that the reaction is indeed accelerated in explicit water, and that the free energy surface begins to look stepwise, as opposed to the concerted mechanism predicted in implicit solvent. (Singleton has observed similar behavior before, and I’ve seen this too.) They do some other interesting studies: they look at the difference between methanol and water as solvents, argue that Houk is wrong about the role of water in the TS,4 and suggest that the hydrophobic effect drives solvent-induced rate acceleration.5
Figure 4B from the paper, showing the change in the PES.
The results they find for this particular system are interesting, but more exciting is the promise that these techniques may soon become accessible to “regular” computational chemists. Duarte and co-workers have shown that ML can be used to solve an age-old problem in chemical simulation; if explicit solvent ML/MD simulations of organic reactions become easy enough for non-experts to run, I have no doubt that they will become a valued and essential part of the physical organic chemistry toolbox. Much work is needed to get to that point—new software packages, further validation on new systems, new ways to assess quality and check robustness of simulations, and much more—but the vision behind this paper is powerful, and I can’t wait until it comes to fruition.
Thanks to Croix Laconsay for reading a draft of this post.
Footnotes
This video from Chris Cramer makes the point nicely.
This obviously depends on the system in question, and what processes are being studied. But in general insufficient sampling is a big issue in molecular dynamics, which I think is underappreciated by organic chemists wading into the area. Jeff Grossman has a nice paper on this.
If you look carefully, many people who claim to be doing big ab initio molecular dynamics studies are actually doing semiempirical molecular dynamics. This isn’t dishonest per se, but it’s a little underwhelming to a computational chemist, especially when it’s only mentioned in the SI. Things get even more confusing when plane wave DFT is employed: in theory, plane wave DFT can be just as accurate as regular DFT, but in practice there are some sneaky approximations that often get introduced.
This argument hinges on whether uphill dynamics (starting from reactants, going to transition state) or downhill dynamics (starting from transition state, going to reactants) are more appropriate. The authors argue that "uphill dynamics allow the solvent sufficient time to reorganise [sic] before the trajectory passes the free energy barrier, providing a more realistic view of solvent behaviour [sic] during the reaction." I'm not fully convinced by this—isn't the idea that the system reorganizes to minimize the energy of the transition state a basic precept of transition state theory? But I'm not convinced I understand these issues deeply enough to have an opinion; I will leave this to the experts.
This argument hearkens back to some old-school computational organic chemistry I love from Bill Jorgensen, studying the hydrophobic effect on conformational preferences of butane. We usually think of the hydrophobic effect as associated with macromolecules (ligands binding to proteins, etc), but it can still matter in tiny systems!
