Synthesis meets theory: Past, present and future of rational chemistry

3.1 Electronic structure

Hartree–Fock (HF) theory uses a single Slater determinant in its wavefunction definition to comply with the anti-symmetric nature of electron spin. Initial trial molecular spinorbitals are expressed as linear combinations of atomic orbitals and optimized through the Fock operator (Roothaan-Hall equations) and an iterative process known as self-consistent field [43]. The limit of HF resides in its inexistent treatment of electron correlation. Atomic orbitals can be constructed using Slater functions (STO), ϕξ,n,1,m(r,Θ,Φ)=N⋅rn−1⋅e−ξr⋅Y1,m(Θ,Φ), Pople Gaussians (GTO), ϕξ,n,l,m(r,Θ,Φ)=N⋅r(2n−2−1)⋅e−ξr2⋅Yl,m(Θ,Φ), or plane waves, ϕ(r)=A⋅e−iK⋅r. Gaussians are computationally advantageous since they decay with a squared dependence on the radius, but a linear combination of them might be necessary to achieve accurate reproduction of the atomic orbital near and far off the nucleus. Basis sets like 6–311++G(d,p) or aug-ccpvTZ are set of functions using a linear combination of three Gaussian functions to reproduce the valence of the atom; they are complemented by polarization and diffuse functions [44] for accuracy. Post-HF methods are improvement to the HF theory and they usually provide high quality electronic structures, when coupled to saturated basis sets, but they are also time-consuming and resources-intensive [45]. Most of these methods feature multi-determinant expansions of the wavefunction and respect size-consistency [15, 43]. The principal methods [45] are configuration interaction (CI), Moller-Plesset perturbation theory (MP2, MP3, MP4 and MP5, following the expansion of the perturbation term), coupled cluster (CCD, CCSD and CCSD(T)), quadratic configuration interaction (QCIS, QCISD, QCISD(T), QCISD(TQ)), Brueckner Doubles (BD, BD(T), BD(TQ)), compound methods (G2, G3 or CBS).

The early alternative to HF-based methods for large systems was represented by semi-empirical methods: they feature minimal quantum-mechanical treatment of the atomic valences through small basis sets, combined with zero differential overlap (ZDO) approximation [15] and parameterization with terms derived from experimental reference data [15]. Today neglect of diatomic differential overlap (NDDO) methods like PM7 [46] are excellent and inexpensive way to obtain semi-quantitative information on a system, particularly useful when conformational analysis is required. The second alternative is represented by density functional theory (DFT) methods. DFT theory has its roots in the work of W. Kohn, P. Hohenberg and L.J. Sham [47, 48]. Kohn–Hohenberg theorems bind the ground state electron-density of a molecular system, ρ(r) = ∑i⁡ηi|ψi(r)|2, with the ground state energy, E_KS[ρ(r)], through an exact and unique functional F[ρ(r)]. Thus, a simplified version of these theorems could be written as E_KS[ρ(r)] = T_S[ρ(r)] + J[ρ(r)] + E_EN[ρ(r)] + E_XC[ρ(r)], where T_s[ρ(r)] term is the Hamiltonian of the system, J[ρ(r)] term is the Coulomb electron-electron repulsion, E_en[ρ(r)] term is the electrons-nuclei interactions and E_XC[ρ(r)] is the exchange-correlation energy, which includes several corrections including quanto-mechanical treatments of Coulomb and Fermi holes [49]. Once E_XC[ρ(r)] term is approximated, the computation of the other potential energy terms is generally fast, thus DFT methods are usually faster than post-HF counterpart and scale very favorably with number of functions (~10¹ to 10³, versus 10⁴ of MP2 and 10⁷ of CCSD(T) methods). This makes DFT methods very suitable for investigations on large non-truncated systems, where accuracy and speed must go hand in hand [45]. It is precisely the form in which E_x[ρ] and E_c[ρ] are “approximated” and combined together that defines the name of the DFT method (Table 3). J. Perdew envisioned the development of DFT accuracy and performances using the Biblical metaphor of Jacob’s ladder [50]: each rung of the ladder, with its associated strengths and weaknesses, leads to a global improvement of the calculated DFT results over the previous one. The “heaven” is represented in this case by a density functional method using a fully non-local formalism. LDA, or local density approximation, is the first and most rudimentary DFT (Rung 1 on Jacob’s ladder). Proposed by P. Dirac [51], the first local density approximation was used in conjunction with Thomas-Fermi gas model [52]. It is a fully local approach, since the exchange-correlation energy depends upon the local electron density at a certain point, ρ(r). Despite its limitations, it can afford good predictions, mostly due to cancellation of errors when estimating exchange and correlation. GGA, or generalized gradient approximation, introduces a higher degree of space inhomogeneity (rung 2, Table 3). Exchange-correlation term depends upon local electron density and local gradient of the electron density, |∇ρ(r)|, at a certain point. The results obtained by GGA are usually better than those obtained by LDA, thought this is more true for small molecules than for solids. Meta-GGA is an improvement over GGA since the dependencies have been expanded including the Laplacian of the electron density, ∇²ρ(r), and the kinetic energy density, τ (rung 3, Table 3). The next rung in the ladder introduces a degree of non-locality by replacing some local exchange energy density with exact HF exchange energy density; these methods are called hybrid GGAs (rung 4, Table 3). The degree of HF exchange cannot be estimated a priori, but needs to be fitted empirically. The improvement in performances of this class of functionals over the previous ones has been quite drastic, as the hybrid GGAs statistically outperform any previous functional in calculating a wide range of chemical properties. When the kinetic energy density is also taken into account, the result is called hybrid meta-GGA functional (rung 4, Table 3). B3LYP combination (B3 exchange functional by Becke and LYP correlation functional by Lee-Yang-Parr) is a hybrid GGA functional containing 20 % HF-exchange and it is regarded by most as the standard in organic chemistry, although Minnesota functionals like M06 (27 % HF-exchange) and M06-2X (54 % HF-exchange) have received particular attention, since they show superior performances in thermochemistry, non-covalent interactions and kinetics [53, 54]. DFT description of dispersive interactions (i. e., hydrogen bonds, London dispersions), however, is generally regarded as unsatisfactory [49, 55]. Double hybrid functionals are supposed to deliver superior performances for weak interactions by including corrections at perturbation theory level, but they generally show similar computational costs to MP2 [56, 57, 58, 59, 60, 61]. The performances of XYG3 functional in calculating thermochemistry, kinetics, weak interactions are remarkable; it is claimed to fill the 5th position of the ladder [62]. The reproduction of long-range dispersive forces can also be achieved using molecular mechanics corrections to DFT; empirical dispersion formulas have been implemented at no extra computational cost for almost any functional [63, 64, 65, 66]. Functionals like B97D [64, 66] or ωB97XD [67] already contain dispersion corrections. Figure 4 has been imported from R. Peverati and D. Truhlar’s recent work [68] and represents a superb way to capture 30 years of improvements in DFT performances at first glance. Summaries of several aspects of DFT and DFT-related performances and implementations can be found in a few manuscripts that combine in-depth expertise with excellent readability (particularly for non-experts) [68, 69].

Figure 4:

Average of mean unsigned errors (in kcal/mol) for all DFT functionals tested in Truhlar’s database CE345 and grouped according to their year of publication. This database includes 15 subgroups of different chemical properties.

Heavy metals of the third row, especially Au [70], display electrons moving at nearly-relativistic speeds due to a high Z number. The moving mass of these fermions undergoes relativistic increment given by the equation m=m01−(vc)2. This phenomenon leads to important scalar effects, like the contraction of the s and p shells and the expansion of the d and f shells, and vectorial effects, like spin-orbit coupling [71]. These effects, in turn, influence the chemistry of heavy transition metals in terms of reactivity, bonding preferences and strengths, coordination numbers and spectroscopy [71]. Relativistic calculations can be introduced via:

1. Effective core potential, ECP, like LAN [72], SDD [73, 74] and CEP [75]. Basis functions are replaced by potential functions at the core and the extent of the replacement defines large-core (frozen core and valence), medium-core or small-core (frozen core + valence treated by basis set) pseudopotentials. Accurate results are more likely to be obtained using small-core ECPs [76]. The basis sets complementing medium and small-core ECPs are designed to fit the pseudopotential used, generally through appropriate primitive contractions [77].

2. Relativistic Hamiltonians coupled with full-electron basis sets. Relativity can be brought into the Schrodinger equation through the Dirac equation. Douglas–Kroll–Hess (DKH) transformation is commonly used for its accuracy, efficiency and effectiveness [78]. Appropriate basis sets are built by contracting primitives [77].

The goal of electronic structure calculations is to generate potential energy hypersurfaces, alias “mapped” polydimensional surfaces that correlate the potential energy in function of the nuclear coordinates of the atoms. Wells or minima on these surfaces represent stable structures. Calculated metric parameters for these minima provide useful information to synthetic chemists because they can be compared with metric parameters derived from phase scattering of x-ray or neutron single-crystal diffraction experiments [79]. Vice versa, this is also an approach used by theoreticians to validate bona fide the goodness of their calculations. This approach is, however, a risky procedure and must be done carefully for two reasons: the maxima in scattering amplitudes may not coincide with nuclear positions in light atoms and structures in solid state are affected by crystal packing and dielectric fields, not present in the isolated calculated counterparts [80].

Figure 5:

Seidel-Seppelt unique compound. Au•••F contacts with the counterions are shown.

Electronic structures predicted the existence of “impossible” compounds like AuXe⁺ and XeAuXe⁺ in 1995 [81]. The calculations were carried out by P. Pyykkö et al. using methods like MP2, MP4 and CCSD(T): the bond distance Au-Xe calculated at MP2 level was found to be 2.691 Å for AuXe⁺ and 2.66 Å for XeAuXe⁺. The stability of these compounds was attributed to the high electronegativity of gold, due to relativistic effects; the bond dissociation energy for AuXe⁺ decreases from 0.910 to 0.376 eV when relativistic effects are omitted at CCSD(T) level [81]. The existence of AuXe and XeAuXe⁺ was confirmed experimentally in 1998 by Schröder, Pyykkö et al. [82]. Finally, Seidel and Seppelt isolated the first stable gold compound with xenon, [AuXe₄][Sb₂F₁₁]₂ (Figure 5), in 2000[83]. The compound confirmed the stability postulated for this class of compounds being stable up to 40˚C. Metric parameters for Au-Xe bond distance, 2.728(1)-2.750(1) Å, extracted from X-ray diffraction agrees well with distances calculated at MP2 level, 2.787 Å [83]. Electronic structures and enthalpy of formation predicted also the existence of another “impossible” compound in the solid state, Na₂He. The predicting algorithm in the USPEX code [84] suggested the existence of a cubic-phase structure of Na₂He, stable at pressures above 160 GPa. The structure was effectively isolated as a cubic-phase, stable from 113GPa up to 1000 GPa, and characterized in matrix by A.R. Oganov et al. [85]. The electronic structure of the solid was studied by density of states, ELF and Bader’s analysis [85].

Electronic structures can be directly linked to magnetic spectroscopies based on the Zeeman effect [86]. In nuclear magnetic resonance (NMR)[87] and electronic paramagnetic resonance (EPR)[88], low frequencies waves (radiofrequencies for NMR and microwaves for EPR) are pulsed to make net nuclear (NMR) or electronic (EPR) spin magnetization process about an axis in a magnetic splitting field B₀[89]. This precession is sensitive upon the chemical environment of each nuclear or electronic spin. NMR spectra (shift and spin-spin coupling) can be efficiently simulated in gas or in solvent phase by calculating 1D-NMR shielding tensor and susceptibilities using Gauge-Independent Atomic Orbital (GIAO) [90, 91, 92] or Continuous Set Of Gauge Transformations (CSGT)[92, 93, 94] methods. Accurate results versus experiment can be achieved by using high level of theory, usually DFT or MP2[92, 95, 96], in conjunction with large, polarized and diffused basis sets. Similarly, hyperfine coupling constants for EPR spectra can be calculated using EPR-II and EPR-III basis sets [97]. GIAO calculations using a pure GGA functional, BP86 (implemented in G03[98]), in conjunction with SBK basis [99] for metals and 6–311++G(d,p) for carbons and hydrogens have been employed by H.V.R. Dias and T.R. Cundari et al. to provide a comparison to experimental ¹H and ¹³C{¹H} NMR spectra in exceedingly rare adducts of ethane with coinage metals [100, 101], Table 4. A further model [102] was generated using a hybrid GGA, B3PW91, in conjunction with SDD pseudopotential and K.A. Peterson’s correlation consistent basis sets [77]. Theoretical values shown in Table 4 for both models along the homoleptic series [M(η²-C₂H₄)₃][SbF₆] (for M = Cu, Ag, Au) show amazing structural and spectroscopic reproduction of the studied species. Similar computational studies were carried out to explain the NMR and bonding patterns of a peculiar complexation between trans,trans,trans-1,5,9-cyclododecatriene, CuSbF₆ and carbon monoxide [102, 103].

Electronic structures can clarify the dynamic behavior of molecules whose potential energy surfaces change upon the absorption of photons of visible or ultra-violet light. Organic chromophores fall under this category, since light can cause transitions from bonding to antibonding orbitals [104]. Heavy metal complexes show also spin-allowed, Laporte-allowed transitions, ligand-to-metal (LMCT) and metal-to-ligand (MLCT) charge transfers [32]. The transition momentum from ψ₁ to ψ₂ eigenstate, M_2,1 = <ψ₂|μ|ψ₁>, can be calculated for vertical excitations (complying with Born-Oppenheimer approximation) to provide transition energies, multiplicities, symmetry and oscillator strengths, thus simulate UV-Vis and circular dichroism spectra. UV-Vis spectroscopy [89] is used by experimental chemists to identify and quantify an active substance through absorbance and Beer-Lambert law; the description of this approach and its limits are well known [105]. Time dependent density functional theory (TDDFT) coupled with medium-to-large basis sets generally provides reasonably accurate results (i. e., within 0.4 eV from experimental transitions) at a very affordable cost, but it has limitations like the absence of double excitations [106]. Most used DFT functionals are B3LYP and PBE0, although the popularity of functional including long-range corrections to the Coulomb term (e. g., cam-B3LYP, LC-ωPBE and ωB97XD) is rising. More refined and resource-consuming alternatives to TDDFT feature CASSCF and CASPT2, techniques based on the complete active orbital space (including PT2 corrections in the latter) [107]. Modeling excited states encompasses the study of the reorganization energy of the molecule and the solvation layer after the absorption of light; bathochromic shifts can be obtained this way to simulate experimental fluorescence spectra. Study of intersystem crossing phenomena (i. e., conical intersections between different spin potential surfaces) can be valuable to study phosphorescence [108]. Intersystem crossings are more common for heavy metal complexes with strong spin-orbit coupling. Diffuse functions and long-range corrected functionals might be mandatory when studying excited states of higher energy (i. e. other than the first one). Photo-induced single transfer redox reactions via outer sphere mechanism involve bimolecular events between oxidants (electron acceptor) and reducing species (electron donor); Marcus theory was designed to calculate the barrier of activation in terms of free energy difference for such processes, ΔG^‡ = (ΔG+λ)²/4λ (λ is the reorganization energy of all the atoms involved in the reaction)[109, 110]. The thermodynamic function that expresses quantitatively the propensity for the transfer of electron(s) is the standard reduction potential; it can be measured experimentally by cyclic voltammetry and calculated using a variety of methods involving implicit, explicit or mixed solvation schemes [111]. Mean unsigned errors of computed values versus experimental ones are reported in the range of 64 mV (~ 6 kJ·mol⁻¹) and as low as 50 mV for selected classes of compounds [111]. Electrochemical methodologies are fundamental in processes like artificial photosynthesis and water splitting [112, 113].

3.2 Solvation schemes

The solvent influences kinetic and thermodynamic aspects of a chemical reaction, thus the accurate evaluation of solvation effects should be brought into the calculations. Several different solvation schemes are available: implicit, explicit or mixed implicit-explicit. Implicit methods allow for quanto-mechanical treatment of the solute and for bulk treatment of solvents as continuous envelopes around solutes, characterized by macroscopic properties like dielectric constants and microscopic properties like solvent radii. Usually computationally very affordable for spectator solvents and reasonably accurate, implicit solvation models are listed below:

1. Polarizable Continuum Models include D-PCM [114], C-PCM [115] and IEF-PCM [116].

2. Solvation Model Based On Density or SMD [117].

3. Models using the GB approximation for the electrostatics include SM8 [118], SM8AD [119], SM12 [120].

4. Conductor Like Screening Model or (COSMO) [121] and its real solvent variant (COSMO-RS) [122].

5. Others [123].

Explicit methods allow for quantum mechanical treatment of the solute and treat solvents as discrete entities around the solute. They are strongly suggested when interactions between solvent and solute need to be characterized, as in the case of strong hydrogen bonds or coordination to heavy metals, but demanding in terms of resources. The nuclear motion of the solute can be simulated using molecular dynamics or multi-scale models [111]. Mixed methods are probably the best chemical option, although they should be chosen and set-up wisely to provide accurate results at affordable costs.

3.3 Ensemble properties

Physical observables represent macroscopic bulk properties of a chemical sample. Simulating those properties at microscopic level requires the average of all the possible states of a chosen ensemble. The taxonomy of statistical thermodynamics encompasses some notable ensembles for N-particle systems: microcanonical or NVE (constant N, volume and energy), canonical or NVT (constant N, volume and temperature), grand canonical or μVT (constant chemical potential, volume and temperature), isoenthalpic-isobaric or NPH (constant number of particles, pressure and enthalpy), isothermal-isobaric or NPT (constant N, pressure and temperature of the system). Calculating all the possible states within an ensemble is impossible, thus few clever simplifications are normally introduced, in order to combine accuracy with computational feasibility. Dynamics methods listed below are excellent examples of methods working on ensemble properties of many-particle systems. Molecular Dynamics, Langevin Dynamics and Quantum Dynamics are particularly suitable for heterogeneous systems, since they are composed of different phases, thus tend to be more affected than homogeneous systems by the multi-scale complexity of the system: the kinetics of diffusion-controlled reaction, for instance, is controlled by particle transport phenomena like diffusion, convection, migration, adsorption and desorption. Chemical reaction happening at the surface of metals, alloys or minerals are many and diverse. The knowledge of the principles ruling transformations in heterogeneous catalysis will allow one to “tailor” a catalyst to carry out the target transformation sought [124]. Figure 6 relates schematically the dimensionality of chemico-physical phenomena to the length/time scale of computational methods [125]

Figure 6:

Multi-scale computational methods most frequently used associated to time-scales of chemico-physical phenomena.

3.3.3 Molecular dynamics

Molecular dynamics is a very popular and largely used technique among researchers in many fields. A many-particle system is left evolved in time through phase-space trajectories that respond to classical laws of motion (Newton): F[x(t)] =  - ∇U[x(t)] = m∙x¨(t), with U being the potential energy of the particle interaction. The equations of motion are usually integrated numerically over the selected interval of time (t+Δt) chosen to describe the system. The initial parameters needed to evolve a system usually are the internal forces among particle, F[x(t)], the composition of the sample (its composing masses, m), the considered ensemble (usually microcanonical) and the time of simulation (t_final – t_initial). The forces between particles are usually derived from parameterized force fields. The standard choice for the integrator falls upon the Verlet algorithm, x(t+Δt) = 2 x(t) + x(t-Δt) + (Δt²/m)•F[x(t)]. This algorithm is particularly useful because it does not allow error accumulation, thus the simulations can be stable for long time [138]. The systems are normally simulated for few femtosecond under periodic boundary conditions to simplify the constraints of a boundary. MD simulations can be efficiently used to obtain static descriptions of the system as well as dynamic evolution of its physical properties. MD simulations can be also run at constant temperature or pressure [139]. One of the major drawback of MD is that the calculation of interatomic forces does not embed any type of quantum mechanical treatment. Without a proper electronic description, MD is not suitable to describe reactions where bonds are formed or broken. The addition of quanto-mechanical electronic calculations extends the use of MD to bond-breaking and forming events (thus scale), but also reduce the size of the sample treatable with this approach by four orders of magnitude (~10² atoms). Quantum dynamics approaches like Born–Oppenheimer molecular dynamics (BOMD) [140] and Car-Parrinello molecular dynamics (CPMD) [140] methods are elegant and rigorous treatment of quantum electronic systems where classical atomic motion is separated, in virtue of Born-Oppenheimer principle, from quanto-mechanic electronic motion. The electronic degrees of freedom are calculated at every step of the optimization in BOMD, but introduced through fictitious variables in CPMD. This allows the simultaneous calculation of both electronic and nuclear degrees of freedom. If from one side the predictive power of these techniques in terms of chemical bonding and its evolution with time is incredible, on the other side, these techniques require great amounts of computational time and resources [140]. This makes them unfeasible for large systems and/or for system where long trajectories and high quality of energy conservation laws are needed. An interesting use of Born-Oppenheimer molecular dynamics (BOMD) to unravel the contribution of entropic disorder on the stabilization of surface polar terminations was published by M. Capdevila-Cortada and N. Lopez [141]. The understanding of surface terminations behavior is very important in the design of new materials. The group performed the study using DFT (PBE+U term) under an NVT ensemble [141]. J.K. Nørskov et al. used the RPBE functional in conjunction with plane-waves, contained in the code DACAPO [142], to model the molecular dynamics of a cheap inorganic surface able to evolve hydrogen in place of the expensive and rare platinum [143]. The evolution of hydrogen is one of the most important electrochemical processes taking place on heterogeneous catalysts for its connection with sustainable energy [144]. Metals like Ni and Mo bind H* too strongly to guarantee a sustainable release of H₂, while Au, for instance, suffers from the opposite problem, being a very poor binder instead (Figure 7). The observation of naturally occurring enzymes able to evolve hydrogen, like hydrogenases and nitrogenases, has sparked the desire to reproduce the catalytic properties of their active centers using mimicking synthetic alternatives [145, 146, 147]. The investigation pointed out that the sites at the edges of MoS₂ are able to bind H* and release H₂, resembling the catalytic structural motif and related activity of the nitrogenase FeMo cofactor (an enzyme that catalyze the evolution of molecular hydrogen, Figure 7). Following the computational lead, the group synthesized nanosized MoS₂ and proved that this material produces molecular hydrogen with an overpotential of 0.1–0.2 V [143].

Figure 7:

Learning from Nature. Mimicking the active site of enzymes like nitrogenase FeMo cofactor (left) leads to rational design of surfaces evolving H₂, like nanosized MoS_2.

Molecular dynamics, coupled with other techniques like umbrella sampling, conformational flooding, metadynamics, and adaptive force bias, plays a pivotal role in drug discovery as well [148]. Drug discovery deals with the energetics, affinity and selectivity of binding/unbinding events between naturally occurring biological receptors and synthetic inhibitors [148]. These events, however, span microseconds to milliseconds, while ab initio molecular dynamics can effectively sample a few femtoseconds at best [149]. G. de Fabritiis et al. published a paper [150] on the quantitative reconstruction of the binding process between β-trypsin and benzamidine using a Markov state model (MSM) [151]. High-throughput classical molecular dynamics parallelized on clusters of graphical processor units (GPUs) [152] allowed 495 simulations of this binding event of 100 ns each. 187 simulations gave productive bindings with root-mean-square deviation of atomic positions within 2Å compared to the crystal structure [153], Figure 8 left. The estimation of the standard free energy of binding is within 1 kcal/mol from the experimental value [150]. Similar routines are employed to screen and/or improve allosteric drugs (drugs that does not bind the orthosteric pocket of the receptor) and numerous methods to predict allosteric sites and their druggability are described and reviewed in literature [154].

Figure 8:

Left, X-ray crystal structure of the Ca²⁺-enzyme-inhibitor complex between benzamidine and serine protease β-trypsin (PDB ID: 3PTB). Right, X-ray crystal structure of the complex barnase-barstar (PDB ID: 1BRS).

3.4 Transition state theory

H. Eyring, M.G. Evans and M. Polanyi proposed this theory in 1933 [169]. Transition structures are mathematical constructs of dynamical instability, dividing surfaces between short-time intrastate (within basin) from long-time interstate dynamics (basin to basin). They possess 3N-7 vibrational degrees of freedom, where the missing degree represents the reaction coordinate. Thus, transition states (TS) are not real structures and many scientists refers to TS-derived thermodynamics as talking of transition state-derived thermodynamics as quasi-thermodynamics [170]. Fundamental requirement of TST is that reactants (or products) and connected transition structures are thermally equilibrated [170]. Several minimization algorithms are available to locate transition structures: steepest descent, Newton-Raphson, rational function optimization, direct inversion of iterative space (GDIIS and GEDIIS) and synchronous transit-guided quasi-Newton (STQN) [171]. The immediate purpose of Eyring-Polanyi equation is to calculate (or estimate at the very least) rate constants for chemical reactive events. Conventional Transition State Theory (CTST) usually provides an upper-bound estimate of experimental reaction rate constants [170, 172]. Within this approach, the transition surface is placed at the saddle point and 100 % crossing-over efficiency is assumed between reactants and products. The variational optimization of the dividing surface in order to minimize re-crossing leads to the Variational Transition State Theory (VTST) [170, 173, 174, 175, 176]. VTST and its variants usually produce much better agreement with experimental kinetics than CTST [172]. The program POLYRATE [177] implements several variational approaches (e. g., Canonical VTST, Improved Canonical VTST, Microcanonical VTST), as well as several quantum tunneling corrections for light atoms. As a consequence of the dual particle-wave nature of matter, light particles like electrons or protons can “tunnel” through potential energy barriers [178, 179]. De-Broglie equation predicts a wavelength of ~27 Å associated to an electron and a ~0.6 Å wavelength associated to a proton moving with kinetic energy of 20 kJ/mol. Experimental signs of tunneling effect usually entail strong H/D isotopic effect, non-linear behaviors of the Arrhenius equation and temperature independence of the reaction rate. Examples of electron tunneling can be found in semiconductivity, superconductivity, scanning microscopy and biological phenomena related to charge transfer [180]. The first experimental proof of tunneling in a chemical reaction came from F. Williams’s work [181]. There are many experimental evidences claiming tunneling effects of light atoms like protons in elimination reactions [182]. Heavier atoms like carbon or nitrogen or systems with small reduced masses can experience tunneling effect as well, in reactions like conformational changes [183], isomerizations [184] and ring expansions [185]. The computation of tunneling effects is generally accomplished with accurate multidimensional derivations of WKB method [186], like the small-curvature (SCT), the large-curvature (LCT), the optimized multidimensional and the least-action tunneling approximations [170, 187]. The inclusion of a multiplicative factor, α_T, called transmission coefficient, takes care of the corrections to the Eyring equation that relate to effects like tunneling, for instance. The inclusion of the transmission factor into the calculation of rate constants (therefore reaction equations) improves the goodness of calculated reaction rates versus experimental [170, 172, 188]. The following equations show the analogy between transition state theory and experimental kinetic laws in calculating rate constants (for standard state of 1 mol/L): (a) Eyring-Polanyi equation, (b) E-P equation with separated enthalpic and entropic terms and (c) Arrhenius law.

Transition state theory provides valuable kinetic discrimination of reaction-controlled events, alias events that need to cross an energy barrier. The first predictive use of transition state theory in organic chemistry involved the electro-cyclic ring opening of 3-formylcyclobutene; computational insights predicted for the latter system an unusual reactivity compared to its congener 3-methylcyclobutene [189]. The Hajos-Parrish reaction was the first organocatalysis to be studied from first principles [190, 191, 192]; yet again transition state theory proved to be an invaluable tool in uncovering the energetic of the stereoselection between the more stable s-trans geometry (3.5 kcal/mol) and the less stable s-cis [193]. The first example of computational design was applied to selective Mannich-type reactions. The reaction between propionaldehyde and N-PMP-protected α-iminoethyl glyoxylate is catalyzed by (S)-proline and affords Mannich products of syn configuration (syn:anti 3:1, 99 % enantiomeric excess, entry #1 in Table 5). K.N. Houk, C. Barbas et al. produced a computational model of the reaction at HF/6–31G* level in 2006. Transition state theory provided an excellent estimation of the observed diastereoselectivity, 3.5:1, and enantiomeric excess, ~95%[193]. They produced a second computational model where (S)-pipecolic acid replaced (S)-proline as active catalyst of the reaction (Figure 9 and entry #2 in Table 5) and noticed that a loss of diastereoselectivity in Syn(1) product was taking place in favor of Anti(2) product. Targeted catalytic design, based on the insight built during this strategic approach, lead to the first computational design of an organocatalyst in 2006. The catalyst, an artificial aminoacid (entry #3 in Table 5), was designed at HF/6-31G* level of theory to afford selective anti-Mannich reactions, a completely opposite reactive pathway in this type of systems. The predicted outcome of the reaction, 5:95 syn:anti and 98 % ee, were perfectly matched by the successive experimental findings, 6:94 syn:anti and 99 % ee.

Figure 9:

Schematic representation of Barbas-Houk’s Mannich-type reaction discussed in this section (top). Transition state theory lead to crucial observations (bottom) concerning structures and energetics that sparked computational design of a new catalyst with desired properties.

Table 5:

The strategy to catalytic design. Entry #1 shows the data for the Mannich-type reaction catalyzed by (S)-proline, entry #2 shows those for the reaction catalyzed by (S)-pipecolic acid. Computation designs experiment in entry #3.

While many other excellent examples of synergistic theory-synthesis approaches in chemistry can be found in literature [194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208], we report here a very interesting paper published by A. Kleij, C. Bo et al. on the use of triazabicyclodecene (TBD) as a organocatalyst of chemoselective ring-opening of cyclic carbonates to afford N-aryl carbamates in cheap and mild conditions [209]. Carbamates are products of CO₂-incorporation. Greenhouse gases (GHG) like CO₂ derived from anthropogenic activities are held responsible for the global warming scenario [210, 211, 212, 213, 214]. Sequestrating these gases from the atmosphere and incorporating them in building blocks for fine chemistry is highly sought after [215, 216, 217]. DFT-based (B97-D3/6-311G**/SMD) methods as implemented in G09[98] helped to uncover a favorable proton-relay mechanism as the main source of the fast catalytic effect of TBD at room temperature (Figure 10); two hydrogen bonds imprint a 8-membered ring transition state and lower considerably the barrier of activation at the rate-determining step. Two competitive pathways, namely the direct nucleophilic attack of the amine onto the cyclic carbonate and water-assisted pathway, are either kinetically unfeasible or not as effective as TBD (Figure 10). The crucial contribution of this proton-relay switch has been counter-checked by experiment: replacing TBD by methylated TBD (MTBD) leads to much poorer catalytic effect [209].

Figure 10:

Transition structures for ring-opening of a cyclic carbonate enhanced by n-butylamine: left, non-catalyzed (ΔG^‡ = 33.3 kcal/mol), center, H₂O-catalyzed (ΔG^‡ = 24.2 kcal/mol) and right, TBD-catalyzed (ΔG^‡ = 18.2 kcal/mol).

A priori design of enzymatic pockets relies on transition state theory as well. Theozymes are theoretical catalysts tailored by optimizing the transition state of the reactants in a non-catalyzed reaction [218]. Though natural enzymes are molecular machines of unmatched efficiency in binding and transforming several substrates with high selectivity, synthetic enzymes can be designed to catalyze reactions unknown to natural enzymes. This is the case of the retro-aldol dissociation reaction of 4-hydroxy-4-(6-methoxy-2-naphthyl)-2-butanone [219], the Kemp rearrangement reaction of a benzisoxazole [220] and the Diels-Alder condensation between 4-carboxylbenzyl trans-1,3-butadiene-1-carbamate and N,N-dimethylacrylamide [221]. The group employed a sequential methodology involving data-mining, predictive design (Rosetta hashing algorithm [222]), transition state optimizations, selection of best leads and finally expression and directed evolution to generate synthetic enzymes. The strategic approach and the exceptional results of such computationally-assisted design are reported for the first two reactions in Figure 11 and Table 6, respectively. An experimental technique that could be conceptually linked to transition state theory is femtosecond spectroscopy. This technique allows pulsing wave-packets on a molecule by the means of an ultrafast laser. A vibrational normal mode of interest is excited and the structure rearranges trough an activated complex to the formation of the desired product. This spectroscopy allows the accurate mapping of the entire surface and its dynamics from the reactants to the products, thus providing important experimental details on the characteristics of the “near transition” region [11].

Figure 11:

Nature does it better, but scientists are closing the gap. The route to “artificial” enzymes is an interplay ground for statistical techniques, theoretical computation and targeted synthesis.

3.5 Kinetic models

Complex multi-step chemical reactions are the direct result of parallel, competitive and consecutive events establishing an intricate network of simultaneous equilibria in the medium: within this network, species appear and disappear following (sometimes) complex phenomenological laws [223]. Experimental kinetic laws use rate constants derived from macroscopic quantities like activation energies (e. g., Arrhenius [223] law) to average microscopic quantities deriving from particle transport and collisions (continuum → microscopic, Figure 12).

Figure 12:

The different approach of experimental kinetics versus kinetic model.

A kinetic model uses integrated system concentrations over time of microscopic quantities (i. e. transition state thermodynamics) on elemental events to simulate the macroscopic reaction (microscopic → continuum, Figure 12): it is used to complement and corroborate mechanistic proofs derived from potential and free energy surfaces. In principles, a kinetic model should allow the accurate recovery of information on selectivity, turnover number and kinetics ↔ reactivity ↔ structure relation in conditions where Curtin-Hammett principle [224], Winstein-Holness equation [224] and steady-state approximation [223] are not always so easily applicable.

Kinetic models abide the principle of microscopic reversibility and represent a strong direct link between pure theory and experiment since they can be compared straightforwardly with plots obtained by experimental kinetics [225]. The higher the number of kinetic equations included in the simulation, the better the agreement of the model versus experiment. Kinetic models are particularly useful to model catalytic reactions, where the rate limiting step(s) may change throughout the reaction following the change in concentrations of the reactants/intermediates [223]. Analytical integration is prohibitive on large systems; therefore one should resort to robust numerical algorithms like Euler, 4th-order Runge-Kutta, Runge-Kutta-Fehlberg, Runge-Kutta-Cash-Karp, Runge-Kutta-Dormand-Prince, Bader-Stoer and Adams’ methods [223]. “Stiff” equations can be efficiently treated with Bader-Deuflhard algorithm [226]. An elegant example of kinetic models has been reported by F. Maseras et al. It concerns the computational study at B97-D3/6-311G(d)/SMD level of theory of a known 1,3-dipolar cycloaddition between propargylamine and azidoethylamine [227], Figure 13. The computed thermodynamic parameters reproduced successfully the experimental findings (ΔG^‡ = 27.1 versus 27.3 kcal/mol, for the rate limiting step of the reaction) and the kinetic model helped to re-interpret the experimental mechanism of the reaction in light of the catalytic effect of host-guest encapsulation provided by cucurbit[6]uril (CB6) [227]. Another excellent example has been reported by Maseras et al. The work encompasses the use of non-adiabatic DFT calculations (ωB97X-D) supported by kinetic model in explaining the outcome of photo-activated aromatic perfluoroalkylations [228].

Figure 13:

Schematic representation of the dipolar cycloaddition discussed in this section, top-left. Maseras’s pruned network of calculated competing events proves the “combinatorial” approach behind kinetic model, top right. Molecular structures of the most important intermediates in the reaction (bottom).

3.6 Multi-scale models

The idea of a consistent force-field was born in 1968 [229]. The total potential energy of a molecule can be calculated as a function of simple classical laws involving Hooke’s spring law for bonds, periodic laws for dihedral torsions, Lennard-Johns equation for Van der Walls interactions and Coulomb law for charge interactions [230]. The potential energy can be further related to energy minimization, intrastate and interstate dynamics and random moves (Monte Carlo)[230]. Multi-scale models were originally proposed in 1976. They allow the optimization and calculation of properties of very large systems (encompassing hundreds of thousands of atoms); thus they have been very useful in the calculations of large peptides, proteins and enzymes like lysozyme [231]. The general idea behind this method is that the chemical system under study can be divided into different layers, calculated separately using different methods. Larger areas that do not take direct part in the reaction (e. g., protein side chains) are calculated using inexpensive classical treatments like molecular mechanics (MM), while active site of reaction can be treated with higher-level quanto-mechanical theory (QM). For such reason these methods are called QM/MM methods. The energy of a total system, E_Sys, can be defined as:

E_Sys = E_QMQM + E_MMMM + E_QMinterQM/MM + E_MMinterQM/MM

The steric and polarization effects of the classically-treated low-level layer are contained in the term E_MMMM, while E_QMQM contains the quantum mechanics energy of the high-level layer. E_QMinterQM/MM is called electronic embedding, while E_MMinterQM/MM is called mechanical embedding. Different methods have different ways to calculate the latter two parameters. Two of these QM/MM methods are available and implemented in computational software: IMOMM [232] and ONIOM [233] methods.

There are excellent reviews [234, 235] and highlights [236] and commentaries [194] on QM/MM calculations in asymmetric catalysis via organometallic species; these works include reactions like hydrogenation, hydroboration, epoxidation and hydroformilation. F. Schoenebeck et al. used ONIOM(B3LYP/6-31+G(d)/LANL2DZ:HF/LANL2MB) level of theory to design a new small-bite phosphine ligand able to enhance the reductive elimination of PhCF₃ from Pd(II) complexes. Contrarily to accepted knowledge, the group demonstrated computationally that the reactivity was not correlated with the bite-angle of the phosphine, rather with the repulsive electrostatic interactions between the phosphine and the leaving group (namely CF₃) [237]. Calculations provided a good estimate of the activation enthalpy for elimination of PhCF₃ in [(Xantphos)Pd(CF₃)(Ph)], 25.1 kcal/mol versus the experimental value of 25.9 ± 2.6 kcal/mol [238]. Schoenebeck’s new complex [(dfmpe)Pd(CF₃)(Ph)] (Figure 14), afforded an activation enthalpy range of ΔH^‡ = 20.7–23.5 kcal/mol, depending on various theories [237]. Driven by the successful computational design, the group synthesized [(dfmpe)Pd(CF₃)(Ph)] and demonstrated that the elimination of PhCF₃ is complete after 3 hours at 80 °C after 100 min. Kinetic experiments by NMR measured the activation enthalpy as ΔH^‡ = 27.9 ± 1.6 kcal/mol [237]. Another challenge for QM/MM methods is represented by designing chemotherapeutic agents with in vivo high affinity and selectivity towards target DNA helix. A synergistic approach between molecular dynamics simulations, QM/MM calculation, force fields and experimental resources (i. e., high-resolution structures of nucleosome-drug adduct) has been applied to further the knowledge on binding location, binding modes and selectivity of metallo-drugs versus DNA and proteins during these last years [239].

Figure 14:

Bite-angle is not the issue, says theory. Schoenebeck’s designed palladium complex reductive-eliminates PhCF₃, despite its small bit-angle phosphine. The figure shows the transition structure for such elimination (bond distances in Å).

3.7 Wavefunction analysis

Optimized wavefunctions at the stationary points, Ψi(r)=∑iCl,iΦl(r), reveal a great deal of information concerning chemical bonds structures and properties like nature, directionality, composition, reactivity and selectivity. Among others, some properties derived from wave-function are very useful and routinely used by chemists: molecular orbitals, population analysis, bond order analysis, molecular electrostatic potential, Fukui function [240], charge decomposition analysis [241], natural bond orbital (NBO) analysis [242, 243] and adaptive natural density partitioning analysis [244]. Electron densities obtained from optimized wave-functions, ρ(r)=∑iηi|ψi(r)|2 are also excellent tool of investigation to uncover points of interest within distributions of electronic clouds around atoms and molecules. Unlike wave-functions, electron densities are less sensitive to the choice of theory and basis sets and provide a direct comparison to observable electron density phases obtained from x-ray or neutron single-crystal scattering [79], ρ(xyz)=1V∑hkl⁡F(hkl)e−2πi(hx+ky+lz). The Quantum Theory of Atoms-in-molecule (QTAIM) analysis is a technique based on the topological analysis (critical points) of ρ(r) and ∇²ρ(r)[245]. Figure 15 shows the QTAIM of two isoelectronic closed-shell d¹⁰ complexes, Ni(η²-C₂H₄)₃ and Cu(η²-C₂H₄)₃⁺ generated by MULTIWFN program [246, 247] (other software for QTAIM includes AIM2000 [248] and AIMAll [249]). M-C bond pathway is unique in the case of copper and the bonding pattern in Cu−η²-C₂H₄ is T-shaped, implying a central λ-donation-type interaction between the electron densities of copper and ethylene moieties; in the case of nickel, however, the pathway is doubly degenerated and the bonding pattern assumes a more triangular shape, in line with the increased π-back-bonding character in the Ni-C interaction (increased sp³ character on the carbons). C=C bond critical point (BCP) suggests that ethane loses 10.0 % of C-C internuclear electron density per unit upon complexation with Ni⁰ and only 4.5 % upon complexation with Cu⁺.

Figure 15:

Molecular graphs of Ni(η²-C₂H₄)₃, left, and Cu(η²-C₂H₄)₃⁺, right, overlaid with their respective Laplacian maps, ∇²ρ(r). The structures have been calculated using B3PW91 in conjunction with the non-relativistic full-electron DGDZVP basis set. Blue dots represent bond critical points, or BCP (3,-1), orange dots represent ring critical points, or RCP (3,+1), and brown dots represent nuclear attractors (3,-3). Within green zones kinetic energy density dominates and the electron density is depleted {∇²ρ(r) > 0}. Within red zones potential energy dominates and the electron density is concentrated {∇²ρ(r) < 0}.

Among others, topological analysis functions are very useful and routinely used by chemists: electron localization function (ELF) [250, 251], localized orbital locator (LOL) [252, 253], reduced density gradient (RDG) [254] and source function [255]. ELF was used by M. Kaupp et al. to explain the bond in another uncommon compound, HgF₄ [256]. The bond pattern, consistent with the representation of a low-spin d⁸ complex, defines mercury as real transition metal, rather than a post-transition element. Mercury in high oxidation state was theoretically predicted by Kaupp et al. to be stable as HgF₄ since 1993[257, 258, 259], with weakly-coordinating ligands like SbF₆^- [260] and as HgH₄ by P. Pyykkö et al. [261] and finally isolated and characterized spectroscopically in matrix in 2007 [256].

3.8 Chemoinformatics and machine-learning

Under the name of chemoinformatics falls a large and diverse number of different disciplines spanning from computational chemistry, mathematics, and statistics, data-mining to informatics. The accepted goal of chemoinformatics is to fabricate Hit-to-Lead-to-Candidate [262] predictions on unknown molecular systems using validated and statistically significant data obtained from libraries of known properties and chemical behaviors. Computational chemistry is primarily used in pharmaceutical industries to complement, integrate or prioritize targeted drug synthesis via virtual libraries and virtual screening [263]. Virtual libraries are collections of virtual molecules that can be studied to generate new molecule or to hypothesize chemical properties. Such libraries should be composed of large and fitting populations. These libraries can be increased using “trained” Markov chains like in fragment optimized growth (FOG) algorithm [264]. Virtual screening is the analysis of database of chemical compounds to identify viable candidates. Virtual screening uses in silico techniques like molecular docking [263]. Combinatorial and high-speed syntheses are used as lead or complementary synthetic methods to chemoinformatics, since they are designed to generate or screen large libraries of compounds [263]. High-throughput computational modeling uses languages like the Chemical Markup Language (CML) and the Simplified Molecular-Input Line-Entry System (SMILES). Physico-chemical properties are normally represented by descriptors: a descriptor is a compact transcription of a certain property within a sample or library of species into mathematical language. The goal of a descriptor is to convey the essence of the quantitative structure-activity relationship (QSAR) [265]. The complexity of the property described by the descriptor increases with the dimensionality of the descriptor itself: 1D-descriptors deal with simple macroscopic bulk properties like molecular weight, while 3D-descriptors capture the essence of the donor-acceptor interaction and binding affinity (polar surface area for instance). 1D-descriptors are used as first-approach filters, while 2D-descriptors (also called fingerprints) seem to work well in reproducing many chemical properties [266]. GRID is an interesting program for drug discovery applications and ADME predictions related to pharmacokinetics [43]. QSAR approach is used extensively in pharmaceutical chemistry [267, 268], but its use in catalytic applications involving small molecules results to be still sparse. An excellent example of QSAR, or in this case QSSR (quantitative structure-selectivity relationship), is represented by M.C. Kozlowski et al.’s work [269, 270]. The group aimed to correlate the structure of in situ-generated zinc catalysts (ZnBr₂ + β-aminoacids) with the enantioselectivity of the asymmetric reduction of benzaldehyde to related alcohol. The group generated a grid-based QSAR approach whose descriptors provided a correct replica of the enantioselective footprint of catalysts under study. Since every transition structure have been generated using computationally low-demanding PM3 theory, this model provides important insight at very affordable costs [269].

Machine learning is a branch of computer science that deals with “supervised” or “unsupervised” learning of a machine. The goal of this discipline is the development of artificial intelligences, networks that mimic the human mind’s capability to learn, store and relay data using electrochemical signals [271]. Network analysis, for instance, has been applied in organic synthesis since the Sixties, when E. Corey created software called Logic and Heuristics Applied to Synthetic Analysis (LHASA). The software was meant to suggest the sequences of postulated steps in a de novo synthesis, based on a database of retrosynthetic rules [272]. In 2005, B. A. Grzybowski et al. developed in 2005 a revolutionary concept of a network connecting molecules (“nodes”) by chemical reactions (“arrows”)[273]. The network contains about 6 million organic compounds interrelated by 30,000 retrosynthetic rules. The group focused on the evolution of this network, the relations between connections in the network and reactivity and optimization/pruning of the network itself. The software developed after this network was called Chematica. It can successfully generate and optimize synthetic pathways [274, 275] and detect synthetic pathways leading to dangerous chemicals [276]. Other interesting examples on chemoinformatics and machine learning (i. e., artificial neural networks) can be found in literature [277].