Techniques of Molecular Modelling: Quantum Mechanical Methods

Molecular orbital methods date back to the origins of quantum mechanics and the Schrödinger equation.

1. The first practical chemical solution to this equation was formulated by Erich Hückel in 1931 and dealt exclusively with planar conjugated systems.
2. This molecular orbital theory was thereafter extended by Pauling, Hartree and Fock, Mulliken and many others,
3. and first effectively implemented in the 1960s by Hoffmann and Lipscomb to all valence electron systems (EHT),
4. by Pople in the 1970s who defined a family of semi-empirical methods
5. which were first effectively paramterised by Dewar (MINDO/3, MNDO, AM1), Stewart PM3 and 6)
6. In parallel, effective implementations of ab initio and density functional methods emerged in the 1980s (Gaussian, GAMESS, Turbomole, Orca, etc)

In essence, QM methods differ from the MM methods in providing explicit solutions for the electronic component of a molecule rather than just the nuclear terms.

These methods all adopt the Born-Oppenheimer approximation. Stated approximately, we assume that the motion and (point) position of the nuclei is described by classical equations, whereas the electrons is described by non-classical wavefunctions. This means the nuclear-nuclear energy contributions are classical mechanical (electrostatic point charge), but the nuclear-electron and electron-electron terms are quantum mechanical and the total energy is the sum of these nuclear and electronic components. The starting point for solving the electronic part is the Schrödinger equation:

H.Ψ = E.Ψ

A Quick Summary of the Principle Quantum Mechanical Methods used today

How the various methods scale with size: Let us define N as ~ the number of atoms (slightly more accurately, it represents the number of basis functions associated with all the atoms in the molecule under consideration, and so it also goes up with the quality of the basis set)

• Methods which do not require 2-electron integral evaluations scale as N² (due entirely to diagonalising the Hamiltonian (Fock) matrix, but this only kicks in for very large matrices)
• Semi-empirical methods ~ N² (but can be made to scale as N¹ for larger molecules, say > 500 atoms)
• DFT methods ~N³-N⁴
• MP2 methods ~N⁵-N⁶
• Coupled cluster methods, ~N⁶-N⁸

One also has to consider how the 1st and 2nd derivatives of the energy with respect to the nuclear coordinates behave. The most recent implementations of density functional theory have low scaling 1st and 2nd analytical derivatives, which now routinely enable good geometry optimisations using BFGS-like algorithms.

Case Study 6: Directing effects in the Electrophilic substitution of Aromatic rings

There are two ways of approaching this problem. The simplest is to look at either the properties of the reactant, or the properties of the initial product, the Wheland intermediate. We will look at the first of these initially.

1. The Properties of the Reactant: A substituted Benzene.

   Henry Armstrong, working at (the precursor to) Imperial College around 1887, was one of the first chemists to categorise substituents (R, or Radicles as he calls them) on a benzene ring in terms of the effect they had towards electrophilic substitution reactions (of X) of the ring. With the electron yet to be discovered, he attributed the two observed modes of behaviour (i.e. o/p in Table I vs m in Table II of 10.1039/CT8875100258) as being due to an electrically polarizable entity he called an "affinity", of which he suggested a benzene ring had six. He also suggested these affinities acted at a distance over the whole ring, but which differed in their directionality (or Resultant as he puts it, the modern term for which is Vector) according to the nature of R (a difference we nowadays categorize as being the electron donating or withdrawing properties of R). Armstrong's description of the properties of his affinity (See DOI: 10.1039/CT8875100258 where he writes "the introduction of a radical/substituent [onto the benzene ring] doubtless involves an altered distribution of the affinity, much as the distribution of the electric charge in a body is altered by bringing it near to another body". In 10.1039/PL8900600095, p102, he also clearly describes what we now know as a Wheland intermediate) may well constitute the first glimmering of the wave-like properties of the electron!

   

   Erich Hückel some 44 years later was able to derive from Schrödinger's wave equation a simple expression which predicted how these "affinities", now named electrons of course, could be polarized. The resulting (π) molecular orbital functions (eigenvectors) can be used to illustrate graphically how both the directing effects (i.e. the probability of finding electrons in any particular position, i.e. the MOs) and the activating/deactivating effects (i.e. the probability of interacting with the electrons, i.e. the orbital energies) operate.

   To maximise the visual effect, we are going to use two substituents R that were not in fact in Armstrong's original list; CH₂⁺ for an electron withdrawing group and CH₂^- for an electron donating group. Its worth considering just for a moment why we are not using more conventional neutral groups (i.e. NO₂ and NH₂ respectively) for this illustration. Being neutral, the latter groups can only polarize the molecule by separating charge to produce a dipolar ionic species. Such ionic charge separation always takes a fair bit of energy to achieve compared to neutral covalent bonding, and actually requires proper solvation to treat quantitatively. It is much easier if the R group is already ionic, because then moving the charge from one part of the molecule to another takes little energy and no change in ionicity, and hence eliminates the need for any solvation treatment.

   o/p Directing	m-Directing
   HOMO, R=CH2-	HOMO, R=CH2+	HOMO-1, R=CH2+

   The Highest occupied π canonical Molecular orbital (HOMO) for R=CH₂^- corresponds to the electron-pair least bound to the molecule and is hence the most available to react with an Electrophile. It has an energy of -1.7eV, which shows how relatively weakly it is bound (most π electrons have energies around -10ev or even less). This weak binding is also associated with it being a good nucleophile. Note how this wavefunction has density on the ortho and the para positions of the ring, but none at all on the meta position. This of course matches the resonance bond formulation familiar to all organic chemists, and hence the recognition of o/p direction and activation for such a substituent.

   Contrast this with R=CH₂⁺. The HOMO in this case has an energy of -15.0eV, very much more tightly bound, with hence with vastly reduced nucleophilic properties. It has a node in the para position. The next highest orbital is almost degenerate with the HOMO (-15.6 eV) and has a node in the ortho positions. Because these two orbitals are essentially the same in energy, they have to be considered together; the only position where both of them do have density is in the meta position. This substituent therefore is meta directing, but strongly deactivating (due to the very negative orbital energies).

   If we now assume that the reactant has similar properties to the transition state for this reaction (the Hammond principle, in part) we can infer that the transition state will inherit the same properties, and that the electrophilic reactivity of such substituted benzene derivatives can be derived purely from the properties of the reactant (if its an early transition state)

2. The Properties of a Heterocycle: Pyridine.

   What about a real system? The differing behaviour of Pyridine and its N-Oxide towards electrophilic substitution is well known. In particular, these two, apparently almost identical species, have quite different properties when nitrated. The former nitrates (albeit slowly) in the 3-position, whilst the latter nitrates (rather faster) in the 4-position. Does the HOMO for each reflect this? The answer is that the corresponding orbitals look very similar to the ones shown above! Inded, the size of the orbital at the 4-position even looks larger than that at the 2-position.

   
   
   

   HOMO, Pyridine, m-directing	HOMO, Pyridine N-Oxide, o/p directing

3. Properties of the Wheland/Armstrong Intermediate.

   An prequel to looking at reactivity of aromatics is to consider the relative energy of the possible Wheland intermediates as being better indicators of the outcome of the reaction than simply the wavefunction of the reactant. How can one calculate the relative energies? The mechanics method, since in general the force field constants do not cancel out for the various pararmeters, is going to give a very unreliable relative energy. It will also not reflect on the polarization of the π-electrons as illustrated above. Only proper treatment of the electrons via quantum mechanics can treat this sort of problem. The energies shown below have been obtained from the PM6 semi-empirical method (which calculates the energies and optimizes the geometries of the protonated Wheland intermediate in just a few seconds). The outcome shows clearly the differentiation between Pyridine and its N-oxide, and also reproduces the preference for 4- rather than 2-substitution in the latter (But see DOI: 10.1039/P29750000277 where an argument is presented for 4-nitration resulting not by reaction of the free N-oxide but the protonated species).

   Substrate, X	ortho (Relative to para)	meta (Relative to para)	para
   Pyridine, X=N	-3.9 (234.6)	-15.1 (223.4)	0.0 (238.5)
   Pyridine N-Oxide, X=NO	4.8 (208.9)	25.6 (231.2)	0.0 (205.6)
   X=NS	3.3 (232.4)	23.5 (251.1)	0.0 (227.6)
   Phosphabenzene X=P	0.1 (210.7)	0.0 (210.6)	0.0 (210.6)
   Phosphabenzene P-oxide, X=PO	1.7 (145.4)	45.1 (188.8)	0.0 (143.7)
   X=PS	2.2 (170.7)	45.1 (213.6)	0.0 (168.5)
   Arsabenzene X=As	3.6 (229.6)	19.6 (245.6)	0.0 (226.0)
   Iridiabenzene X=Ir(PH3)3	3.4	39.4	0.0

   
   

   What about other substituents? An almost entirely unstudied problem is the electrophilic substitution of phosphabenzene. It's a great surprise to find that in modelling terms at least, it behaves quite differently to pyridine! Using this approach, one can even predict reactivity for something as unconventional as metallabenzenes (for a review of the chemistry of metallabenzenes, s ee 10.1039/b517928a). For examples of such metallabenzenes with or
   . For details of the bonding in the Ru system, see this blog post.

Case Study 2 Revisited: Where ARE the electrons in the Pirkle Reagent?

Part 1: The π-electrons. One solution of the Schrödinger equation provides a canonical set of energy levels for electron (pairs) in the molecule, that having the highest energy being referred to as the HOMO (Highest occupied molecular orbital). As we saw before, the HOMO corresponds to the electron (pair) which is least strongly bound to the molecule, and hence can be regarded as describing where the most basic (= proton seeking) or nucleophilic (= nucleus seeking) electrons in the molecule are. So a good first start to understand where a hydrogen bond might occur is to plot the HOMO for the Pirkle reagent. Unfortunately, the HOMO is a π type orbital and shows almost no discrimination for binding to a proton!

Notice how the two faces of the anthracene ring are quite different, an effect induced by the electron withdrawing effect of the C-CF₃ group which withdraws electrons anti-periplanar to the bond, i.e. on the face opposite to it. But in turn the lone pair of the oxygen atom re-focuses density on to one ring of this opposite face, and in fact the hydrogen bond forms at almost exactly the centre of this electron density. Thus this hydrogen bond is created by a complex stereo-electronic effect caused by the shaping of the electron density by one withdrawing group and one lone pair.

Part 3: The electron Topology (AIM or Atoms in Molecules). What happens when one places two of these molecules into the close proximity found in the crystal structure? The π-MOs are of little help, and the electrostatic potentials of the two components, as used above, interfere to the point of not really revealing anything. Yet another way of analyzing where the electrons are is needed. This is found in a technique known as Quantum Chemical Topology. It is the study of the topology of electron density in molecules, and uses the language of dynamical systems (critical points, manifolds, gradient vector fields etc, first introduced by the mathematician Henri Poincare) to identify key points in the electron densities.

There are precisely four types of points, differing in the characteristics of the curvature of the electron density ρ(r) in the region of the critical point. This curvature is obtained from the second derivative of the electron density with respect to cartesian space (the density Hessian). The eigenvalues of this 3 by 3 matrix can have exactly four conditions, which are summarised by the signature of this matrix, being the sum of the signs of these eigenvalues:

1. NCP: A critical point (all three eigenvalues -ve) which identifies where an attractor (i.e. a nucleus) is in a molecule
2. BCP: A critical point (two -ve, one +ve) which identifies bonds in a molecule. Normally found between a pair of attractors, but can also involve three attractors (a 3-centre bond).
3. RCP: A critical point (one -ve, two +ve) which identifies rings in a molecule
4. CCP: A critical point (all three eigenvalues -+ve) which identifies cages in molecules
5. The total number of these various points must satisfy a topological relationship known as the Poincaré-Hopf condition, which states that num(NCP-BCP+RCP-CCP) = 1.
6. One property is known as the Laplacian of the electron density, ∇²ρ(r) (being the sum of the diagonal elements of the density Hessian) at a BCP is used to characterise the type of bond in that region
   • Covalent, ∇²ρ(r) negative, ρ(r) large (> 0.2)
   • Ionic, ρ(r) small
   • Charge shift, ∇²ρ(r) positive, ρ(r) large (> 0.2)

The most relevant of the critical points is the second type; a BCP. It tells us in effect if a bond exists between any pair of nuclei. What it does not tell is is what the bond order of that bond is; it could easily for example be less than 1. That is a chemical rather than a topological interpretation.

What happens when such an analysis is conducted on the dimer of the Pirkle reagent? In the resulting analysis (left), the light blue points are the positions of some of the critical bond points identified in the electron density topology (there are many more than shown here). These reveal some fascinating insights into how the Pirkle reagent binds with itself.

The very recent NCI (non covalent interactions) method takes this analysis one step further by developing the concept of surfaces of interaction rather than just critical points. These surfaces are derived from the density gradient ∇ρ(r) (the deviation from a homogenous electron distribution), coloured by the value of the λ₂ eigenvalue of the Laplacian of the density ∇²ρ(r). The NCI method is showing great promise for explaing e.g. stereoselectivity in chemical reactions.

Part 4: Energies again. Our analysis of the Pirkle reagent concludes here with a revisitation of the binding energy. We saw previously how this could be estimated using Molecular Mechanics. This had the advantage of including the so called dispersion, or long range correlation effects, but the disadvantage that there was no ready access to the entropy of the process. A very new procedure which incorporates both advantages is the dispersion-corrected density functional method, one very recent implementation of which is ωB97XD (DOI: 10.1039/b810189b), and this gives the following values:

ΔH dimerisation -25.1; ΔS -49.2, TΔS = -14.7; ΔG dimerisation -10.4 kcal/mol.

Case Study 1 Revisited: How do orbitals interact within a molecule?

We (think we) know where covalent σ-bonds are in molecules, they lie along the straight line connecting the nuclei of two atoms. We can idealize lone pairs on atoms such as oxygen, and in case study 1, all we needed was an (X-ray) structure to spot a connection between the (antiperiplanar) orientation of (certain) such idealized lone pairs and the C-O bonds in the molecule. But are the lone pairs really where we have idealized them? And in dealing the orientation of the lone pair with a σ-bond, we become aware that we are actually dealing with a concept derived from the Klopman-Salem equation, which uses not so much the two-electron σ-bond, as a vacant orbital occupying the (assumed) same orientation known as the σ^*-bond. This was the basis of the earlier lecture course on stereo-electronics and also conformational analysis. It is time to take some of these concepts and show how they can be properly quantified. Let us proceed as follows:

1. We will solve the Schrödinger equation for a molecule (to any required degree of accuracy). This gives us functions (the canonical molecular orbitals) which tell us the electron density probability over the molecule as a whole.
2. We now need to convert this (global) function to some sort of local property, since the concept of a bond is local, and not descriptive of the molecule as a whole. One such procedure is known as the NBO (Natural bond orbitals). A more formal description of this type of orbital is: Natural Bond Orbitals (NBOs) are localized few-center orbitals ("few" meaning typically 1 or 2, but occasionally more) that describe the Lewis-like molecular bonding pattern of electron pairs in optimally compact form. This can be precisely expressed in mathematical form (but we shall not do so here!); it is sufficient to say that the electron density probability predicted by summing all the occupied canonical molecular orbitals and alternatively all the local NBOs comes out the same! (i.e. there is more than one way of constructing a wavefunction for a molecule that predicts its overall electron density distribution, which in fact is the only properly measurable property).
3. Just as with canonical MOs, where we had a (doubly) occupied set and a virtual (unoccupied) set, so the NBOs turn out to geometrically occupy the bonds (BD), the lower energy non-valence cores (CR), and the valence lone pairs (LP) of the molecule, with the unoccupied orbitals being represented by BD* and RY* (don't worry about this last one, it stands for Rydberg).
4. In the NBO formalism, one can calculate an interaction energy between any doubly occupied BD or LP NBO, and any unoccupied BD* NBO. This can be represented in very simple terms by the diagram:
   
5. This is exactly the form of the Klopman-Salem diagram you may well recollect from 2nd year conformational analysis lectures. But we have now made progress. By using the NBO expressions, we can obtain an exact value for the value of the quantity E(2).
6. When this is done for the two dioxepins discussed in case study 1, the following results (at the B3LYP/6-31G(d) level) are obtained (with the keyword pop(nbo)) for the stable (green) and unstable (pink) forms;
   
7. The only two large E(2) terms (both ~15.5 kcal/mol) are between an oxygen LP and a C-O BD* involving the adjacent oxygen. They correspond exactly to the two antiperiplanar interactions we first noticed by visualisation, but now we can precisely quantify these interactions.
8. For the unstable isomer, only one large E(2) term is found (13.4); the other (7.5) is now significantly reduced because the orientation is no longer exactly antiperiplanar. The molecule is less stabilised and the two C-O bonds are no longer equal in length.
   

NBO analysis is very useful in a variety of situations, and is increasingly found in the literature. For example, it can be used to rationalise the gauche preference of 1,2-difluoroethane, and many others. Such E(2) NBO terms are thought to indicate useful chemistry down to energies of ~3 kcal/mol, and ones as high as 130 kcal/mol have been found (~15 is however typical of the anomeric effect). See this blog for another example where the unexpected was highlighted by such techniques.

Case Study 4 Revisited: The Helical structure of DNA

I left this case study dangling, by posing the question: how was it determined that the helix is right handed? Watson and Crick actually answered this by stating (very briefly with no details) that a left handed helix could only be constructed by violating the permissible van der Waals contacts. The vdW contacts are actually very well modelled using molecular mechanics (the 4th term in the force field described here) and so a reasonable model could have been constructed using this technique to test this assertion. It suffers one slight deficiency. Because it uses force constants to define the terms, these are assumed to be the same everywhere throughout the molecule. Thus the rungs of the DNA ladder in the middle are treated the same as at the end of the helix. Is this justified? It is only very recently that vdW effects have been incorporated into quantum mechanical models in a manner which can be applied to large molecules (ωB97XD, see above). So we will kill two birds with one stone by jumping straight to this technique. We will also overcome another defect of mechanics, and that is the general lack of implementation for incorporating entropic effects on the energies. We will pose the new question: what are the relative free energies of and handed DNA helices (which, by the way are diastereomers, not enantiomers!). It turns out the answer depends on the nature of the base pairs.

1. If they are exclusively C-G, as in d(CGCG)₂, ΔG₂₉₈ favours the left handed Z-helix by ~12 kcal/mol.
2. If they are A-T, as in d(ATAT)₂, the right handed B form is favoured by ~4.3 kcal/mol. As it happens, this is exactly what is known for small (tetramers) of DNA!
3. Dispersion terms alone (van der Waals) make Z-d(CGCG)₂ 5.1 kcal/mol LESS stable than B-d(CGCG)₂,
4. Dispersion terms alone (van der Waals) make Z-d(ATAT)₂ 12.5 kcal/mol LESS stable than B-d(ATAT)₂.

We conclude that Watson and Crick's model building, based on van der Waals analysis, did indeed favour the B-DNA helical form, but that there are other aspects that can overcome this for CG-rich DNA strands. See for example this blog for more details.