Search tips
Search criteria 


Logo of blackwellopenThis ArticleFor AuthorsLearn MoreSubmit
Journal of Computational Chemistry
J Comput Chem. 2016 January 5; 37(1): 83–92.
Published online 2015 July 31. doi:  10.1002/jcc.24025
PMCID: PMC5324590

Beyond static structures: Putting forth REMD as a tool to solve problems in computational organic chemistry

Section Editor: Charles L. Brooks, III, Masahiro Ehara, Gernot Frenking, and Peter R. Schreiner


Computational studies of organic systems are frequently limited to static pictures that closely align with textbook style presentations of reaction mechanisms and isomerization processes. Of course, in reality chemical systems are dynamic entities where a multitude of molecular conformations exists on incredibly complex potential energy surfaces (PES). Here, we borrow a computational technique originally conceived to be used in the context of biological simulations, together with empirical force fields, and apply it to organic chemical problems. Replica‐exchange molecular dynamics (REMD) permits thorough exploration of the PES. We combined REMD with density functional tight binding (DFTB), thereby establishing the level of accuracy necessary to analyze small molecular systems. Through the study of four prototypical problems: isomer identification, reaction mechanisms, temperature‐dependent rotational processes, and catalysis, we reveal new insights and chemistry that likely would be missed using static electronic structure computations. The REMD‐DFTB methodology at the heart of this study is powered by i‐PI, which efficiently handles the interface between the DFTB and REMD codes. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Keywords: computational organic chemistry, replica exchange molecular dynamics, ab initio molecular dynamics, density functional tight binding


Studies of organic systems frequently utilize computational results as an essential tool to elucidate mechanistic reaction details that are difficult or impossible to access experimentally.1, 2 The literature is rife with examples of reaction pathways in which reactants, intermediates, and products, as well as their associated transition states, are cartooned as a series of static geometries each possessing a specific reaction enthalpy. Often, a picture of this type successfully captures the key aspects of a system allowing, for example, identification of the primary mechanistic pathway3, 4, 5, 6, 7, 8, 9 or rationalization of the presence of a specific intermediate.9, 10, 11, 12, 13 Such descriptions, however, are occasionally insufficient to chemistry occurring in an experimental setting,14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33 in the most extreme cases leading to disastrous failures.34 Inside the computer, a host of factors that govern “real world” chemical reactions must necessarily be approximated or ignored altogether in a static picture. Chiefly among these is a precise description of a reaction's free energy, as opposed the frequently reported enthalpy. While estimates of the free energy within a harmonic approximation are provided by most quantum chemistry codes, techniques based on molecular dynamics (MD) represent a more appealing option to fully access the entropic contribution to the stability of different molecular states, that can also be extended to explicitly include environmental effects (e.g., temperature, pressure, etc.).

From the perspective of the computational organic community, an ideal methodology would directly provide ab initio free energies at accessible computational costs, for example, by combining the physical quantities associated with the free energy obtained from MD simulations with the accuracy of static quantum chemical computations. Roughly speaking, the current state of MD simulations to estimate these quantities can be subdivided into two classes: classical or molecular mechanics (MM) and ab initio (AIMD) methods. On one hand, conventional MM force fields, such as those successfully used in molecular biology and pharmaceutical chemistry,35, 36 rely only on the nuclear coordinates of a system, making them very fast. Unfortunately, in most cases they are incapable of describing chemical processes involving the breaking and formation of bonds (Reactive force fields capable of simulating chemical reactions do exist, but are relatively few in number).37 This makes MD simulations based on MM force fields very useful to examine biological phenomena such as protein folding.38 On the other hand, AIMD is capable of describing chemical reactions involving bond breaking and formation. Indeed, the combination of density functional theory (DFT) with AIMD methods has already been used within the framework of Car‐Parrinello39 and Born–Oppenheimer40, 41 MD to resolve problems associated with reaction pathways, 42, 43, 44 phase transitions,45, 46 and solute/solvent interactions.47, 48 The above mentioned examples utilize GGA functionals and plane‐waves, which benefit from being both relatively fast39, 49 and highly scalable,50, 51, 52, 53 but are, nonetheless not ideal for the organic community as chemically intuitive concepts and properties are lost owing to the delocalized nature of the planewaves and the inevitable use of effective core potentials. In contrast, Born–Oppenheimer MD simulations, using typical localized basis set quantum chemistry codes featuring the more reliable global hybrid functionals and post‐Hartree Fock methods, are rather sparse (see e.g., Refs. 54, 55, 56) (For a seminal implementation on GPU see Ref. [56]). A serious limitation involves the duration for which a simulation can be propagated in time. Naturally, the increased computational expense of AIMD simulations, which require computing potentials from a first principle method such as DFT, limits their applicability to short ~102–103 ps intervals. Depending on the complexity of the potential energy surfaces (PES), chemically relevant interconversions between different possible states may not appear on such short time scales. One feasible solution is to lengthen the simulation time by reducing the computational expense of determining the first principles potential by replacement with a potential derived from semiempirical methods, such as density functional tight binding57, 58 (DFTB). For instance, replacing the DFT by a semiempirical DFTB potential permits access to simulations that are three orders of magnitude longer (on the 1–10 ns time scale) while also accommodating a tenfold increase in the number of atoms.59 Despite these improvements, the complexity of the PES of many organic reactions remains sufficiently large that visiting all chemically relevant regions remains essentially impossible. To ensure that all of the PES is explored, additionally computational tricks must be used.

MD simulations of systems containing a large, complex PES necessitate the use of enhanced sampling techniques that facilitate thorough exploration of the free energy landscape. These techniques are required to overcome problems associated with running insufficiently short simulations by reducing the amount of time that a system spends trapped in a local energy minimum. Such enhanced sampling techniques can be roughly divided into two groups. The first is concerned with the identification of pathways between known initial and final states. This category includes, for instance, transition path sampling60 or constrained dynamics.61 However, these approaches do not facilitate searching for the free‐energy global minima or other important configurations of a system within a complex free energy landscape.62, 63, 64, 65 The second category of enhanced sampling methods is better suited to tackle this problem, as it is aimed at obtaining a thorough exploration of the low‐energy portions of the free energy landscape. Some of these techniques, such as metadynamics62 or accelerated MD66 (aMD), rely on modified potentials. While most of these methodologies were originally envisioned for the study of biological systems using classical MD,36, 67, 68 both metadynamics69, 70 and aMD71 have been used to study chemical reactions requiring quantum mechanical treatments. The principal disadvantage, in this case, is that some insight about the system's reactivity is needed to choose an appropriate reaction coordinate, which could in turn affect the outcome of the simulation.

In contrast, modified sampling approaches, such as replica exchange MD72 (REMD), do not require any prior insight. However, they have mostly been used in the context of biological simulations together with empirical force fields, and have only rarely been used in concert with quantum chemical methods. The few existing studies are limited to molecular clusters containing only a handful of atoms.73, 74, 75 The basic idea of parallel tempering replica exchange is to simulate N replicas of a system at a range of different temperatures. Replicas propagated at high temperature freely explore a large amount of the PES in an unencumbered manner by overcoming any barriers present, while low temperature replicas explore local minima regions from which they are unable to escape. The key to the improved sampling seen in REMD involves exchanging complete configurations from replicas at different temperatures via a Metropolis–Teller algorithm,76 thereby enhancing exploration of the entire free energy landscape.

In principle, the coupling of potentials derived from semiempirical methods with REMD would allow access to additional information and larger systems than for the first principles potentials obtained in AIMD. Such a tool would be very useful within the realm of organic chemistry, where the PES landscapes of reactions may be quite complex. Moreover, the importance of directly determining free energies can assist in exploring the chemistry of systems where entropy is known to play an important role. Here, we present results obtained by coupling the i‐PI interface for advanced molecular simulations77 with the DFTB357, 58, 78, 79 semiempirical framework to conduct REMD with the objective of exploring the free energy surfaces of organic systems. Such an approach is appealing as it combines the thorough statistical sampling enabled by enhanced MD with the reliability of approximate quantum chemical techniques, which are capable of accurately describing the energetics of structures of organic systems at an affordable cost. Through coupling with i‐PI, we highlight the abilities of the REMD@DFTB3 method to address prototypical cases relevant to the computational organic community including: (i) exploration of the conformational space of a dithiacyclophane molecule possessing multiple local minima, (ii) estimation of the minimum energy pathway and free energy barrier of the Cope rearrangement (CR) of semibullvalene (SBV), (iii) distinguishing entropically versus enthalpically favored conformational states of a molecular rotor, and (iv) identifying the key conformations of a widely used organocatalyst, cinchona alkaloid.

Computational Details

All forces are computed at the DFTB3/3OB80, 81 using the Universal Force Field parameters to account for dispersion forces as implemented in the DFTB+ program.82 The DFTB+ code was interfaced with the dynamic driver i‐PI after minor modifications. The Hubbard derivatives and the h‐damping factor were chosen according to Refs. [80,81]. A finite electron temperature of 300 K was selected to improve the convergence of the geometries arising from the hottest replicas. While a serious drawback for systems sensitive to solvent effects, the current implementation is restricted to gas phase simulations. Regardless, many valuable aspects of a system's behavior can still be extracted from gas phase data. In turn, this newly revealed information may lead to more informed predictions about how the same system would behave in the condensed phase.

The i‐PI software drives the REMD, which evolves in the NVT ensemble within a cubic box length of 1000 Å. The large box avoids the spurious interaction between replicas within the periodic boundary conditions. For a given system, the simulations were initiated from the same structure and velocities. Gaussian distributions centered at 900 K provide the atomic velocities for each replica. A time step of 0.25 fs was found to be sufficient to integrate the Newton equation without observing any drift on the conserved quantity (see Supporting Information Fig. S6). Exchanges between replicas (Fig. (Fig.1)1) were attempted stochastically every 50 steps83 on average using the following probability (The exchange time was chosen in accordance with the current state‐of‐the‐art, as is frequently done in the literature):


where k B is the Boltzmann constant and Ti and Tj are the temperature of the exchanged replica. In this way, the mapping of the PES is canonical, meaning that the free energy profile obtained at each temperature corresponds to the free energy profile that would be obtained in a normal MD with much improved sampling efficiency.

Figure 1
Illustration of replica exchange: the PES of a hypothetical system is represented in gray. The red lines illustrate the replicas (i.e., trajectories) at different temperatures. Frequent exchanges are attempted between the replicas based on a Metropolis ...

A Langevin thermostat, with a time constant of 100 fs, was used to maintain a constant temperature for each replica. The temperature ranges from 300 K for the coldest replica to 1500 K for the hottest. Sixteen replicas were found to be sufficient (i.e., provided enough round‐trips among the replicas during the simulation time) for the dithiacyclophane molecule and the CR of SBV. Larger molecules—the cinchona alkaloid and the molecular rotor needed 48 replicas to effectively exchange among replicas. A table with all temperatures is provided in the Supporting Information (Table S1). Snapshots of atomic configurations were saved every 50 steps, and used for further processing.

The initial 10,000 steps were used for equilibration, and discarded from subsequent analysis. Data from different replicas were combined by weighting each frame by w=eβrβt and computing the observables separately for each parallel tempering temperature. The different replica were then combined with the weighing factors chosen as w2w2w2 according to the error estimates in reweighed averages given in Ref. 84.

Smooth histograms were constructed using kernel binning with a triangular windows function much smaller than the extent of the main features in the free energy landscape. With respect to timing, for each system considered the REMD@DFTB3 computations presented herein could be performed within 2–4 days (real time) depending on the size of the system on Intel‐based (Xeon E5‐2660) cores. A patch for DFTB+ (version 1.2) is available on demand.

The REMD free energies reported in Table 1 (cinchona alkaloid) and Supporting Information Table S2 (dithiacyclophane) for each relevant region (i.e., a basin i) are computed through evaluation of the following integral:

Fi=kBT logNieF(x)/kBTdx

where each of the nonoverlapping integration regions (Ni) covers a neighborhood of the local free‐energy minima. This equation permits one to account for the thermal fluctuations that distort in a different manner the geometries of the various configurations.

Table 1
Relative electronic and free energies (in kcal mol−1) of the cinchona alkaloid conformers determined from gas phase static computations or from the REMD simulations (DFTB3/3OB‐UFF) at 300 K. The geometries for computing the static DFTB3 ...

Static electronic structure computations for the dithiacyclophane and CR included optimizations at the M06‐2X85, 86/def2‐SVP level using the “Ultrafine” grid as implemented in Gaussian09.87 Alternative energy assessments were obtained at the PBE088, 89‐dDSC90, 91, 92, 93/TZ2P (for dithiacyclophane) and CCSD(T)/cc‐pVTZ (for the CR) using ADF94, 95, 96 and Molpro,97, 98 respectively. Reported static free energies include unscaled free energy corrections from M06‐2X/def2‐SVP computations. Reported static DFTB3/3OB‐UFF electronic energies are computed at the DFTB3 optimized geometry.

Results and Discussion


Understanding the conformational analysis of molecules represents a cornerstone of organic chemistry. Free energy mapping provides direct relationships between structure and energy and assists in understanding molecular behavior. In this context, a dynamic exploration of the free energy surface may not only alter, but possibly even reverse pictures provided by static relative energy computations. More importantly, it may also reveal unexpected energetically low‐lying conformations that were not envisioned owing to preconceived user‐based biases about a system. The conformations of dithiacyclophane (Fig. (Fig.2),2), previously investigated by two of the authors,99, 100 perfectly illustrates this aspect. This highly flexible molecule was originally found to possess several low‐lying conformers featuring both π‐stacked (structure 1, meta‐stable) and open conformations (structure 2, lower in energy) using accurate electronic structure methods. Not surprisingly, Born–Oppenheimer MD simulations performed at both the DFT (i.e., PBE) and DFTB3 levels were shown to be highly sensitive to the inclusion of a dispersion correction.99, 100 As might be expected, in the absence of van der Waals corrections, the π‐stacked conformation readily converted into the open conformers (on a scale as fast as 250 fs), whereas the stacked meta‐stable conformer persisted for a few picoseconds in the dispersion‐corrected trajectories. The current REMD@DFTB3 results add a new element to the former picture, revealing a somewhat less intuitive “disarticulated” conformational state (structure 3, Fig. Fig.2).2). More importantly, this new conformational region is thermodynamically favored from the MD free energy and static electronic or free energy perspective (see Supporting Information Fig. S1 and Table S2), and thus would affect any measured properties (e.g., NMR chemical shifts, infrared spectrum, etc.). Note, however, that the conformational entropic contributions are the largest for 2 as indicated by its larger basin (larger number of conformations) and by the small REMD free energy difference between 2 and 3 in comparison to the static picture (see Supporting Information Table S2). The population of 3 remains barely dominant at 300 K, when accounting for the full entropic contribution. The two‐dimensional (2D) plot of Figure Figure22 provides further insights into the conformational dynamics inaccessible from a static picture. For instance, the direct paths connecting the open conformational regions with the two other areas (2  1, 2  3) contrast with the absence of a low‐energy pathway directly connecting the closed and “disarticulated” conformer (1  3). Having access to information of this type could be useful when, for example, aiming to alter dynamic fluctuation through chemical modification. Finally, it is worth noting that “standard” Born–Oppenheimer MD simulations (as opposed to REMD) performed on an even longer timescale (i.e., 1.3 ns see Supporting Information Fig. S2) than those in Refs. [38] and [39] remained, for most of the time, trapped in its original stacked conformational region without identifying the thermodynamically more important region associated with the new conformation.99, 100

Figure 2
Two‐dimensional representation of the free energy landscape obtained from the REMD@DFTB3 of the dithiacyclophane molecule. The relevant collective variables are shown in the plot. 1 kcal mol−1 isocontours are shown in yellow.

Cope rearrangement of semibullvalene

Aside from the identification of chemically important conformers presented in the previous example, REMD can also be used to explore the PES of organic reactions. This second type of case has been explored less often, as most REMD simulations to date have been performed together with a fixed‐bonds, nonreactive force field. The CR of SBV represents a prototypical example for which a great deal of mechanistic detail has already been amassed, including information regarding: the (a)synchronous nature and the role of tunneling in the reaction mechanism,101, 102, 103, 104, 105, 106, 107 the presence and magnitude of homoaromaticity in the transition state,108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121 and the link between molecular properties and structure.122 In its simplest form, the CR transitions between two equivalent structures each characterized by three‐membered ring on one side of the bridging ethyl subunit (Fig. (Fig.3).3). Utilizing a static picture derived from electronic structure theory [CCSD(T)/cc‐pVTZ//M06‐2X/def2‐TZVP], the transition state connecting these structures possesses equivalent C2—C8 and C4—C6 bond distances (2.00 Å) associated with a concerted reaction mechanism with an overall barrier height of 7.3 kcal mol−1.

Figure 3
Cope rearrangement of semibullvalene. a) Two‐dimensional free energy map obtained from REMD@DFTB3 simulations indicating the expected CR (minimum energy pathway, 12, given in yellow) as well as an unexpected region corresponding to ...

From the organic chemist's perspective, this reaction likely appears far too simple to merit study using MD techniques, predominantly because the essential mechanistic components appear well described from static computations. In contrast, our exploration of the PES using REMD@DFTB3 reveals considerably more information than might have been anticipated (Fig. (Fig.3).3). More than 600,000 structures were analyzed to generate a 2D free energy plot, that identifies chemically meaningful regions including the expected SBV minima (1 and 2) and the CR TS (blue central area), as well as the minimum energy reaction pathways connecting these structural regions (yellow). From this data, the relative height of the TS barrier (+8.7 kcal mol−1 is the highest point along the DFTB3 MEP) and any free energy differences between the reactant and products can be estimated (zero in this case) by the upper most reaction coordinate plot of Figure Figure3b,3b, which closely mirror the type of energy plots typically associated with static computations. A quantitative assessment of the reaction rate, however, would require a careful analysis of deviations from transition‐state theory that are to be expected whenever the order parameters chosen to describe the transition are not ideal.

Perhaps, the most relevant feature of Figure Figure3a3a is the presence of three minima, as opposed to the two that we expected. As mentioned earlier, regions 1 and 2 correspond to the two symmetric SBV structures, however, the third black region (3, upper right) represents conformations that possess significant lengthening of both the C2—C8 and C4—C6 interatomic distances. Molecular structures of this type are no longer three‐dimensional cages, instead adopting quasiplanar conformations, as exemplified by dihydropentalene (3, Fig. Fig.3).3). Indeed, such structures are linked to SBV by two distinct, yet directly connected TS structures and a valley‐ridge inflection point indicative of a surface bifurcation.17, 123 Examining the 31 and 32 MEPs reveals this feature: the structure must first pass over a relatively high barrier close to region 3 (Fig. (Fig.3a,3a, also visible 1→3 and 2→3 reaction coordinate plots, Fig. Fig.3b).3b). Overcoming this first, energetically more costly TS barrier leads to the second TS associated with the CR transition between the two symmetric SBV structures. As we observe just a few transitions toward region 3, we cannot deem simulations to be converged and we cannot be confident of the quantitative accuracy of the free energy landscape we computed. More sophisticated REMD implementations or combinations of REMD and metadynamics124 are probably needed to obtain reversible sampling. Similar in nature to the bare SBV picture, placing an electron withdrawing substituents onto selected positions either enhances or suppresses the CR TS barrier.108, 109, 125, 126 For instance, the placement of a CN‐ onto the C8 carbon results in a decrease in the CR transition state barrier as well as creating asymmetrical products and reactants, yet the substituted dihydropentalene is again present (see Supporting Information Fig. S3).

Clearly, the inherent nature of sampling the PES using REMD provides an enhanced view of the ways in which organic systems move between different conformers or along competing energetic pathways. For the example provided here, the three energetic basins corresponding to structural minima, as well as the minimum energy pathways connecting them, were obtained with no prior knowledge of the system. Assessments with no a priori knowledge that are provided by enhanced sampling techniques are likely to become of increasing importance, particularly for systems featuring a more complex PES. Taken to its extreme, REMD can be used to better understand the behavior of systems consisting of thousands of isomers, such as shape‐shifting organic molecules (e.g., bullvalene).127

Molecular rotor

The dynamics of a recently investigated molecular rotor128 (Fig. (Fig.4),4), is a further illustrative case study highlighting the importance of accessing conformational states associated with both entropically and enthalpically favored regions. For this rotor, structural changes occur upon increases in temperature, which are directly visible in the variable temperature 1H NMR spectra. The measured chemical shifts imply that enthalpy favored conformational states, characterized by CH/π interactions between the hydrogen of the central rotating phenylene unit and the aromatic ring of the steroid, are at the origin of the pronounced upfield 1H chemical shift evident at low temperatures. Increasing the temperature results in an increasing population of CH/π unbound rotational states, leading to the displacement of the 1H NMR signal of the rotator C—H group toward lower fields. In such a situation, static computations can help to identify a handful of relevant conformers of both types (bound and unbound), yet the direct one‐to‐one relationship between the most prevalent conformational states at a given temperature and the observed chemical shifts remains undetectable. REMD, conversely, delivers a more comprehensive picture by disclosing key changes in the relative populations associated with temperature‐dependent rotational processes. The REMD@DFTB3 results in Figure Figure5a5a for a temperature of 300 K are analyzed in terms of ring current chemical shielding at the position of the hydrogen atoms H1 and H1 for each structure. The shielding cone created by the aromatic ring of the steroids at proton H1 and H1 were evaluated using Pople's ring current model129 (see Supporting Information Fig. S4). The region above 0 ppm is characterized by CH/π bound states, whereas values around 0 ppm are representative of CH/π unbound states. Accordingly, the results show that the enthalpy driven region corresponds to a set of conformers possessing only one CH/π interaction. In other words, if a CH/π interaction exists with one of the steroid aromatic rings ( δH1orH1> 0 ppm in Fig. Fig.5a),5a), the other side is essentially unbound ( δH1 or δH1 = 0 ppm) at 300 K. Increasing the temperature displaces the population toward the CH/π unbound states (around 0,0) with both hydrogen atoms lying outside the shielding zone (Fig. (Fig.5b).5b). The overall picture is clearly visible in the variable temperature 1H NMR experiment128 but the REMD data provides significantly deeper insight into the dynamic process by revealing the relative conformational population evolution when varying the temperature as well as the direct temperature dependence of the measure properties (Fig. (Fig.5b).5b). For such processes, it is rather unrealistic to neglect the conformational entropic contributions and calculate the “static” population of a few representative conformers randomly extracted from a minimum‐energy structures search, as is often done in computational chemistry.

Figure 4
Schematic representation of the molecular rotor and the key hydrogen atoms.
Figure 5
a) Two‐dimensional map of the shift induced by the current created by the aromatic ring of the steroids on H1 and H1. Pople's ring current model was used to calculate the ring current chemical shifts (see Supporting Information Figure ...

Cinchona alkaloid

The final illustrative application of REMD@DFTB focuses on asymmetric catalysis, more specifically on the use of cinchona‐based primary amines as chiral phase‐transfer catalysts (Fig. 6). This class of compounds, derived from natural sources, has been identified as a promising alternative to other amino‐catalysts such as proline.130, 131, 132, 133 (A) As well as its quinidine diasteteroisomer activate various carbonyl compounds with a consistently high level of stereocontrol. Rather than addressing the origin of the stereoselectivity, which depends heavily on: concentration and nature of the acid co‐catalyst, formation of hydrogen bond motifs, solvent effects, and so forth, here we revisit the gas phase conformational behavior of (A), which is one of the primary elements responsible for the control of stereochemistry.

According to our simulations, the rich conformational behavior of (A) is not easily rationalized in terms of four relevant (anti/syn and closed/open) conformations.132, 134 The 2D conformational map with 5, 10, 15, and 20 kcal mol−1 isocontours given in Figure Figure66 suggests a somewhat richer conformational picture consisting of four easily accessible conformational regions (14) and one that is less frequently visited (2′). The most populated region, 1, essentially encompasses a 60° range for the C8—C9C4—C4a angle with no significant energy barrier. The most illustrative conformation associated with 1 (i.e., lowest free energy) is characterized by an angle around 90°. Previous NMR measurements135 in apolar solvents along with static computations132 established the anti‐open (1) and syn‐open (4) conformers as being energetically comparable and ~5 kcal mol−1 more stable than the “closed” conformers, with the quinuclidine nitrogen lone pair in gauche arrangement with respect to quinoline group (e.g., 2, see Table 1). The information extracted from REMD follows this interpretation while also revealing an additional accessible conformational space (3) lying between the free energy of 1 (or 4) and 2, and a flatter PES around the lower free energy minimum. The broader conformational space of 1 is particularly illustrative of the importance of capturing the full entropic and anharmonic contributions to the free energy to establish the relative stabilities between 1 and 4. In fact, while region 4 is slightly favored at both static DFTB3 and DFT‐D levels,132 the REMD free‐energy conformational minimum corresponds to 1 (Table 1). Another crucial aspect for catalysis concerns interconversion between the different conformational regions. It is clear from Figure Figure66 that complete rotation around the C8—C9 bond is much more favorable in comparison to rotation around the C9C4 bond, which is hindered by the contact between the two polycyclic groups. Even though regions 1 and 4 lie close in energy the syn‐open conformers (4) are trapped and cannot easily convert into other conformational states. Finally, it is interesting to stress that, in line with previous observation, the diastereomers of quinidine, which carry opposite sterochemistries at carbons 8 and 9, essentially behaves like an enantiomer130 as illustrated by the 2D conformational energy map provided in Supporting Information Figure S5.

Figure 6
Two‐dimensional free energy map of the cinchona alkaloid catalyst a). 5, 10, 15, and 20 kcal mol−1 isocontours are shown in yellow. The representative structural figures of the four relevant minima are shown at the bottom.

Regarding the crucial role played by lowest‐energy conformers and by conformational barriers in the enantiodifferentiation processes, clearly having access to the conformational free‐energy landscape is an undeniable asset. While we stress that the present analysis was performed at a fairly low electronic structure level (DFTB3) and in the gas phase (a limitation of the current implementation), the key findings related to the conformational landscape and population remain valid. Based on these finding, we believe that computational studies in the field of asymmetric catalysis would benefit greatly from accessing more realistic pictures of the conformational and reaction dynamics for each step of the catalytic cycle.


Here, the coupling of REMD with DFTB (REMD@DFTB3) via the i‐PI dynamic driver has been introduced and its utility demonstrated as a tool to solve problems in computational organic chemistry. Through the study of four illustrative examples, we have highlighted the importance of mapping the free energy of the PES to increase chemical understanding. As opposed to static electronic structure computations, exploring chemistry via REMD often reveals unexpected or unimagined features even for seemingly quite simple systems. This fact is exemplified by the search for low energy conformers of dithiacyclophane and analysis of the CR of SBV. Chemically more complex examples further illustrate this point by showing the important role played by dynamics in reproducing experimentally observable properties (such as proton chemical shifts) and the understanding of underlying chemical structure and its role in determining stereochemistry during organocatalysis. In general, we hope to have shown the utility of using enhanced sampling approaches, such as REMD, to better understand and solve problems of interest to computational organic chemists.

Supporting information

Supporting Information


C.C. and M.C. thank the MARVEL platform for initiating the discussion that inspired this work. R. P. thanks Piero Procacci for discussion. This paper is dedicated to the memory of an inspiring mentor Paul v. R. Schleyer, who is deeply missed.


How to cite this article: Petraglia R., Nicolaï A., Wodrich M. D., Ceriotti M., Corminboeuf C.. J. Comput. Chem. 2016, 37, 83–92. DOI: 10.1002/jcc.24025


1. Lin Z., Acc. Chem. Res. 2010, 43, 602. [PubMed]
2. Cheng G.‐J., Zhang X., Chung L. W., Xu L., Wu Y.‐D., J. Am. Chem. Soc. 2015, 137, 1706. [PubMed]
3. Bahmanyar S., Houk K. N., J. Am. Chem. Soc. 2001, 123, 11273. [PubMed]
4. Hoang L., Bahmanyar S., Houk K. N., List B., J. Am. Chem. Soc. 2003, 125, 16. [PubMed]
5. Roy D., Najafian K., Schleyer P. v. R., Proc. Natl. Acad. Sci. USA 2007, 104, 17272. [PubMed]
6. Boren B. C., Narayan S., Rasmussen L. K., Zhang L., Zhao H., Lin Z., Jia G., Fokin V. V., J. Am. Chem. Soc. 2008, 130, 8923. [PubMed]
7. Mackie I. D., Johnson R. P., J. Org. Chem. 2009, 74, 499. [PubMed]
8. Yamataka Y., Adv. Phys. Org. Chem. 2010, 44, 173.
9. Zhuo L.‐G., Zhang J.‐J., Yu Z.‐X., J. Org. Chem. 2012, 77, 8527. [PubMed]
10. Himo F., Lovell T., Hilgraf R., Rostovtsev V. V., Noodleman L., Sharpless K. B., Fokin V. V., J. Am. Chem. Soc. 2005, 127, 210. [PubMed]
11. Fristrup P., Ahlquist M., Tanner D., Norrby P.‐O., J. Phys. Chem. A 2008, 112, 12862. [PubMed]
12. Frei R., Wodrich M. D., Hari D. P., Borin P.‐A., Chauvier C., Waser J., J. Am. Chem. Soc. 2014, 136, 16563. [PubMed]
13. Duarte F., Gronert S., Kamerlin S. C. L., J. Org. Chem. 2014, 79, 1280. [PubMed]
14. Li Y., Lin S.‐T., Goddard W. A. III, J. Am. Chem. Soc. 2004, 126, 1872. [PubMed]
15. Kaczor A., Reva I. D., Proniewicz L. M., Fausto R., J. Phys. Chem. A 2006, 110, 2360. [PubMed]
16. Carpenter B. K., Annu. Rev. Phys. Chem. 2005, 56, 57. [PubMed]
17. Ess D. H., Wheeler S. E., Iafe R. G., Xu L., Çelebi‐Ölçüm N., Houk K. N., Angew. Chem. Int. Ed. 2008, 47, 7592. [PMC free article] [PubMed]
18. Schreiner P. R., Reisenauer H. P., Pickard F. C. IV, Simmonett A. C., Allen W. D., Mátyus E., Császár A. G., Nature 2008, 453, 906. [PubMed]
19. Xu L., Doubleday C., Houk K. N., Angew. Chem. Int. Ed. 2009, 48, 2746. [PubMed]
20. Rehbein J., Carpenter B. K., Phys. Chem. Chem. Phys. 2011, 13, 20906. [PubMed]
21. Schreiner P. R., Reisenauer H. P., Ley D., Gerbig D., Wu C.‐H., Allen W. D., Science 2011, 332, 1300. [PubMed]
22. Black K., Liu P., Xu L., Doubleday C., Houk K. N., Proc. Natl. Acad. Sci. USA 2012, 109, 1286.
23. Ley D., Gerbig D., Schreiner P. R., Org. Biomol. Chem. 2012, 10, 3781. [PubMed]
24. Carpenter B. K., Chem. Rev. 2013, 113, 7265. [PubMed]
25. Rossi M., Scheffler M., Blum V., J. Phys. Chem. B 2013, 117, 5574. [PubMed]
26. Sandbeck D. J. S., Kuntz C. M., Luu C., Mondor R. A., Ottaviano J. G., Rayer A. V., Sumon K. Z., East A. L. L., J. Phys. Chem. A 2014, 118, 11768. [PubMed]
27. Samanta D., Rana A., Schmittel M., J. Org. Chem. 2014, 79, 8435. [PubMed]
28. de Souza M. A. F., Ventura E., do Monte S. A., Riveros J. M., Longo R. L., Chem. ‐ Eur. J. 2014, 20, 13742. [PubMed]
29. Jiménez‐Osés G., Liu P., Matute R. A., Houk K. N., Angew. Chem., Int. Ed. 2014, 53, 8664. [PubMed]
30. Karmakar S., Datta A., Angew. Chem. Int. Ed. 2014, 53, 9587. [PubMed]
31. Samanta D., Rana A., Schmittel M., J. Org. Chem. 2015, 80, 2174. [PubMed]
32. Goethe M., Fita I., Rubi J. M., J. Chem. Theory Comput. 2015, 11, 351. [PubMed]
33. Törk L., Jiménez‐Osés G., Doubleday C., Liu F., Houk K. N., J. Am. Chem. Soc. 2015, 137, 4749. [PubMed]
34. Plata R. E., Singleton D. A., J. Am. Chem. Soc. 2015, 137, 3811. [PubMed]
35. Karplus M., McCammon J. A., Nat. Struct. Mol. Biol. 2002, 9, 646. [PubMed]
36. Adcock S. A., McCammon J. A., Chem. Rev. 2006, 106, 1589. [PubMed]
37. van Duin A. C. T., Dasgupta S., Lorant F., Goddard W. A. III, J. Phys. Chem. A 2001, 105, 9396.
38. Freddolino P. L., Harrison C. B., Liu Y., Schulten K., Nat. Phys. 2010, 6, 751. [PubMed]
39. Car R., Parrinello M., Phys. Rev. Lett. 1985, 55, 2471. [PubMed]
40. Barnett R. N., Landman U., Nitzan A., Rajagopal G., J. Chem. Phys. 1991, 94, 608.
41. Wentzcovitch R. M., Solid State Commun. 1991, 78, 831.
42. Hytha M., Štich I., Gale J. D., Terakura K., Payne M. C., Chem. ‐ Eur. J. 2001, 7, 2521. [PubMed]
43. Klein M. L., Ivanov I., J. Am. Chem. Soc. 2005, 127, 4010. [PubMed]
44. Bucko T., Hafner J., J. Phys. Condens. Matter 2010, 22, 384201. [PubMed]
45. Martoňák R., Laio A., Parrinello M., Phys. Rev. Lett. 2003, 90, 075503. [PubMed]
46. Martoňák R., Donadio D., Oganov A. R., Parrinello M., Nat. Mater. 2006, 5, 623. [PubMed]
47. Fois E., Gamba A., Spano E., J. Phys. Chem. B 2004, 108, 154.
48. Trinh T. T., Jansen A. P. J., van Santen R. A., Meijer E. J., Phys. Chem. Chem. Phys. 2009, 11, 5092. [PubMed]
49. Kühne T. D., Krack M., Mohamed F. R., Parrinello M., Phys. Rev. Lett. 2007, 98, 066401. [PubMed]
50. Lippert G., Hutter J., Parrinello M., Mol. Phys. 1997, 92, 477.
51. Lippert G., Hutter J., Parrinello M., Theor. Chem. Acc. 1999, 103, 124.
52. Laino T., Mohamed F., Laio A., Parrinello M., J. Chem. Theory Comput. 2005, 1, 1176. [PubMed]
53. Laino T., Mohamed F., Laio A., Parrinello M., J. Chem. Theory Comput. 2006, 2, 1370. [PubMed]
54. Todorova T., Seitsonen A. P., Hutter J., Kuo I.‐F. W., Mundy C. J., J. Phys. Chem. B 2006, 110, 3685. [PubMed]
55. Wu R., Hu P., Wang S., Cao Z., Zhang Y., J. Chem. Theory Comput. 2010, 6, 337. [PMC free article] [PubMed]
56. Ufimtsev I. S., Martinez T. J., J. Chem. Theory Comput. 2009, 5, 2619. [PubMed]
57. Porezag D., Frauenheim T., Köhler T., Seifert G., Kaschner R., Phys. Rev. B 1995, 51, 12947. [PubMed]
58. Elstner M., Porezag D., Jungnickel G., Elsner J., Haugk M., Frauenheim T., Suhai S., Seifert G., Phys. Rev. B 1998, 58, 7260.
59. Cui Q., Elstner M., Phys. Chem. Chem. Phys. 2014, 16, 14368. [PubMed]
60. Dellago C., Bolhuis P. G., Csajka F. S., Chandler D., J. Chem. Phys. 1998, 108, 1964.
61. Ciccotti G., Ferrario M., Hynes J. T., Kapral R., Chem. Phys. 1989, 129, 241.
62. Laio A., Parrinello M., Proc. Natl. Acad. Sci. USA 2002, 99, 12562. [PubMed]
63. Ciccotti G., Ferrario M., Mol. Simul. 2004, 30, 787.
64. Maragliano L., Vanden‐Eijnden E., Chem. Phys. Lett. 2006, 426, 168.
65. Barducci A., Bonomi M., Parrinello M., WIREs Comput. Mol. Sci. 2011, 1, 826.
66. Hamelberg D., Mongan J., McCammon J. A., J. Chem. Phys. 2004, 120, 11919. [PubMed]
67. Abrams C., Bussi G., Entropy 2014, 16, 163.
68. Bernardi R. C., Melo M. C. R., Schulten K., Biochim. Biophys. Acta 2015, 1850, 872. [PubMed]
69. Zheng S., Pfaendtner J., Mol. Simul. 2015, 41, 55.
70. Ensing B., De Vivo M., Liu Z., Moore P., Klein M. L., Acc. Chem. Res. 2006, 39, 73. [PubMed]
71. Bucher D., Pierce L. C. T., McCammon J. A., Markwick P. R. L., J. Chem. Theory Comput. 2011, 7, 890. [PubMed]
72. Sugita Y., Okamoto Y., Chem. Phys. Lett. 1999, 314, 141.
73. Ishikawa Y., Sugita Y., Nishikawa T., Okamoto Y., Chem. Phys. Lett. 2001, 333, 199.
74. Choi T. H., Chem. Phys. Lett. 2012, 543, 45.
75. Fedorov D. G., Sugita Y., Choi C. H., J. Phys. Chem. B 2013, 117, 7996. [PubMed]
76. Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E., J. Chem. Phys. 1953, 21, 1087.
77. Ceriotti M., More J., Manolopoulos D. E., Comput. Phys. Commun. 2014, 185, 1019.
78. Seifert G., Eschrig H., Bieger W., Z. Phys. Chem. ‐ Leipzig 1986, 267, 529.
79. Yang Y., Yu H., York D., Cui Q., Elstner M., J. Phys. Chem. A 2007, 111, 10861. [PubMed]
80. Gaus M., Lu X., Elstner M., Cui Q., J. Chem. Theory Comput. 2014, 10, 1518. [PubMed]
81. Gaus M., Goez A., Elstner M., J. Chem. Theory Comput. 2013, 9, 338. [PubMed]
82. Dolgonos G., Aradi B., Moreira N. H., Frauenheim T., J. Chem. Theory Comput. 2010, 6, 266. [PubMed]
83. Sindhikara D. J., Emerson D. J., Roitberg A. E., J. Chem. Theory Comput. 2010, 6, 2804. [PubMed]
84. Ceriotti M., Brain G. A. R., Riordan O., Manolopoulos D. E., Proc. R. Soc. A 2012, 468, 2.
85. Zhao Y., Truhlar D. G., Acc. Chem. Res. 2008, 41, 157. [PubMed]
86. Zhao Y., Truhlar D. G., Theor. Chem. Acc. 2008, 120, 215.
87. Frisch M. J., Trucks G. W., Schlegel H. B., Scuseria G. E., Robb M. A., Cheeseman J. R., Scalmani G., Barone V., Mennucci B., Petersson G. A., Nakatsuji H., Caricato M., Li X., Hratchian H. P., Izmaylov A. F., Bloino J., Zheng G., Sonnenberg J. L., Hada M., Ehara M., Toyota K., Fukuda R., Hasegawa J., Ishida M., Nakajima T., Honda Y., Kitao O., Nakai H., Vreven T., Montgomery J. A. Jr., Peralta J. E., Ogliaro F., Bearpark M., Heyd J. J., Brothers E., Kudin K. N., Staroverov V. N., Kobayashi R., Normand J., Raghavachari K., Rendell A., Burant J. C., Iyengar S. S., Tomasi J., Cossi M., Rega N., Millam M. J., Klene M., Knox J. E., Cross J. B., Bakken V., Adamo C., Jaramillo J., Gomperts R., Stratmann R. E., Yazyev O., Austin A. J., Cammi R., Pomelli C., Ochterski J. W., Martin R. L., Morokuma K., Zakrzewski V. G., Voth G. A., Salvador P., Dannenberg J. J., Dapprich S., Daniels A. D., Farkas O., Foresman J. B., Ortiz J. V., Cioslowski J., Fox D. J.. Gaussian 09, Revision D.01; Gaussian: Wallingford, CT, 2009.
88. Perdew J. P., Burke K., Ernzerhof M., Phys. Rev. Lett. 1996, 77, 3865. [PubMed]
89. Adamo C., Barone V., J. Chem. Phys. 1999, 110, 6158.
90. Steinmann S. N., Corminboeuf C., J. Chem. Theory Comput. 2010, 6, 1990. [PubMed]
91. Steinmann S. N., Corminboeuf C., Chimia 2011, 65, 240. [PubMed]
92. Steinmann S. N., Corminboeuf C., J. Chem. Phys. 2011, 134, 044117. [PubMed]
93. Steinmann S. N., Corminboeuf C., J. Chem. Theory Comput. 2011, 7, 3567. [PubMed]
94. te Velde G., Bickelhaupt F. M., van Gisbergen S. J. A., Fonseca Guerra C., Baerends E. J., Snijders J. G., Ziegler T., J. Comput. Chem. 2001, 22, 931.
95. Fonseca Guerra C., Snijders J. G., te Velde G., Baerends E. J., Theor. Chem. Acc. 1998, 99, 391.
96. ADF2013, SCM , Theoretical Chemistry, Vrije Universiteit: Amsterdam, The Netherlandds: Available at:
97. Werner H.‐J., Knowles P. J., Knizia G., Manby F. R., Schütz M., Celani P., Korona T., Lindh R., Mitrushenkov A., Rauhut G., Shamasundar K. R., Adler T. B., Amos R. D., Bernhardsson A., Berning A., Cooper D. L., Deegan M. J. O., Dobbyn A. J., Eckert F., Goll E., Hampel C., Hesselmann A., Hetzer G., Hrenar T., Jansen G., Köppl C., Liu Y., Lloyd A. W., Mata R. A., May A. J., McNicholas S. J., Meyer W., Mura M. E., Nicklaß A., O'Neill D. P., Palmieri P., Peng D., Pflüger K., Pitzer R., Reiher M., Shiozaki T., Stoll H., Stone A. J., Tarroni R., Thorsteinsson T., Wang M.. MOLPRO is a Package of Ab Initio Programs.
98. Werner H.‐J., Knowles P. J., Knizia G., Manby F. R., Schütz M., WIREs Comput. Mol. Sci. 2012, 2, 242.
99. Petraglia R., Steinmann S. N., Corminboeuf C., Int. J. Quantum Chem. (in press). DOI:10.1002/qua.24887.
100. Brémond É., Golubev N., Steinmann S. N., Corminboeuf C., J. Chem. Phys. 2014, 140, 18A516. [PubMed]
101. Zhou C., Birney D. M., Org. Lett. 2002, 4, 3279. [PubMed]
102. Zhang X., Hrovat D. A., Borden W. T., Org. Lett. 2010, 12, 2798. [PubMed]
103. Andrae D., Barth I., Bredtmann T., Hege H.‐C., Manz J., Marquardt F., Paulus B., J. Phys. Chem. B 2011, 115, 5476. [PubMed]
104. Bredtmann T., Manz J., Angew. Chem. Int. Ed. 2011, 50, 12652. [PubMed]
105. González‐Navarrete P., Andrés J., Berski S., J. Phys. Chem. Lett. 2012, 3, 2500. [PubMed]
106. Bredtmann T., Paulus B., J. Chem. Theory Comput. 2013, 9, 3026. [PubMed]
107. Bredtmann T., Ivanov M., Dixit G., Nature Commun. 2014, 5, 5589. [PubMed]
108. Dewar M. J. S., Lo D. H., J. Am. Chem. Soc. 1971, 93, 7201.
109. Hoffmann R., Stohrer W. D., J. Am. Chem. Soc. 1971, 93, 6941.
110. Williams R. V., Kurtz H. A., J. Org. Chem. 1988, 53, 3626.
111. Jiao H., Nagelkerke R., Kurtz H. A., Williams R. V., Borden W. T., Schleyer P. v. R., J. Am. Chem. Soc. 1997, 119, 5921.
112. Williams R. V., Gadgil V. R., Chauhan K., Jackman L. M., Fernandes E., J. Org. Chem. 1998, 63, 3302.
113. Williams R. V., Chem. Rev. 2001, 101, 1185. [PubMed]
114. Williams R. V., Eur. J. Org. Chem. 2001, 227.
115. Brown E. C., Henze D. K., Borden W. T., J. Am. Chem. Soc. 2002, 124, 14977. [PubMed]
116. Goren A. C., Hrovat D. A., Seefelder M., Quast H., Borden W. T., J. Am. Chem. Soc. 2002, 124, 3469. [PubMed]
117. Wu H.‐S., Jiao H., Wang Z.‐X., Schleyer P. v. R., J. Am. Chem. Soc. 2003, 125, 10524. [PubMed]
118. Hrovat D. A., Brown E. C., Williams R. V., Quast H., Borden W. T., J. Org. Chem. 2005, 70, 2627. [PubMed]
119. Greve D. R., J. Phys. Org. Chem. 2011, 24, 222.
120. Griffiths P. R., Pivonka D. E., Williams R. V., Chem. Eur. J. 2011, 17, 9193. [PubMed]
121. Ichikawa Y., Sakai S., J. Phys. Org. Chem. 2012, 25, 409.
122. Jana D. F., Wodrich M. D., Corminboeuf C., J. Org. Chem. 2012, 77, 2548. [PubMed]
123. Castaño O., Frutos L.‐M., Palmeiro R., Notario R., Andrés J.‐L., Gomperts R., Blancafort L., Robb M. A., Angew. Chem. Int. Ed. 2000, 39, 2095. [PubMed]
124. Bussi G., Gervasio F. L., Laio A., Parrinello M., J. Am. Chem. Soc. 2006, 128, 13435. [PubMed]
125. Dewar M. J. S., Jie C., Tetrahedron 1988, 44, 1351.
126. Miller L. S., Grohmann K., Dannenberg J. J., J. Am. Chem. Soc. 1983, 105, 6862.
127. He M., Bode J. W., Proc. Natl. Acad. Sci. USA 2011, 108, 14752. [PubMed]
128. Pérez‐Estrada S., Rodríguez‐Molina B., Xiao L., Santillan R., Jiménez‐Osés G., Houk K. N., Garcia‐Garibay M. A., J. Am. Chem. Soc. 2015, 137, 2175. [PubMed]
129. Pople J. A., J. Chem. Phys. 1956, 24, 1111.
130. Dijkstra G. D. H., Kellogg R. M., Wynberg H., Svendsen J. S., Marko I., Sharpless K. B., J. Am. Chem. Soc. 1989, 111, 8069.
131. Zhou J., Wakchaure V., Kraft P., List B., Angew. Chem. Int. Ed. 2008, 47, 7656. [PubMed]
132. Moran A., Hamilton A., Bo C., Melchiorre P., J. Am. Chem. Soc. 2013, 135, 9091. [PubMed]
133. Prakash G. K. S., Wang F., Ni C., Shen J., Haiges R., Yudin A. K., Mathew T., Olah G. A., J. Am. Chem. Soc. 2011, 133, 9992. [PubMed]
134. Lam Y.‐H., Houk K. N., J. Am. Chem. Soc. 2015, 137, 2116. [PubMed]
135. Bürgi T., Baiker A., J. Am. Chem. Soc. 1998, 120, 12920.

Articles from Wiley-Blackwell Online Open are provided here courtesy of Wiley-Blackwell, John Wiley & Sons