Molecular orbital methods date back to the origins of quantum mechanics and the Schrödinger equation.
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 BornOppenheimer 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 nonclassical wavefunctions. This means the nuclearnuclear energy contributions are classical mechanical (electrostatic point charge), but the nuclearelectron and electronelectron 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.Ψ
Methods which integrate Hamiltonian (H) and Basis set (Ψ)  



Η: List of Ab initio/Density functional Hamiltonians  Ψ: List of Basis sets  



Methods which extrapolate Hamiltonian and Basis set  

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)
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 BFGSlike algorithms.
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.
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 wavelike 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_{2}^{+} for an electron withdrawing group and CH_{2}^{} for an electron donating group. Its worth considering just for a moment why we are not using more conventional neutral groups (i.e. NO_{2} and NH_{2} 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  mDirecting  

HOMO, R=CH_{2}^{}  HOMO, R=CH_{2}^{+}  HOMO1, R=CH_{2}^{+} 
The Highest occupied π canonical Molecular orbital (HOMO) for R=CH_{2}^{} corresponds to the electronpair 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_{2}^{+}. 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)
The Properties of a Heterocycle: Pyridine.
What about a real system? The differing behaviour of Pyridine and its NOxide 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 3position, whilst the latter nitrates (rather faster) in the 4position. 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 4position even looks larger than that at the 2position.
HOMO, Pyridine, mdirecting  HOMO, Pyridine NOxide, o/p directing 

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 semiempirical 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 Noxide, and also reproduces the preference for 4 rather than 2substitution in the latter (But see DOI: 10.1039/P29750000277 where an argument is presented for 4nitration resulting not by reaction of the free Noxide 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 NOxide, 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 Poxide, 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(PH_{3})_{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.
As a visualisation problem, viewing the Crystal structure of the Pirkle reagent provided a rationalisation of how the intermolecular interactions might proceed via a novel form of hydrogen bonding involving an entire face of one aromatic ring. We saw how Molecular mechanics is not parameterised to deal with such an unusual interaction, and so might not handle it properly. Clearly, the location of the electrons is critical, since the H^{+} binds to regions of high electron density. There are three aromatic rings to choose from in this molecule, and moreover, because of the chiral centre present, two faces to each ring, ie six possibilities in all. Can calculating the wavefunction of this molecule cast light on this? 
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!
Part 2: The integrated σ and πelectrons. In fact the effect we are seeking is actually transmitted via the σ framework (originating in the inductive effect of the CF bonds) and not the π system. The next level of approximation is to recognise that the basicity of a molecule derives from not just a single MO, but a suitably weighted sum of the contributions of all the electrons, and particularly of the σ bonds. This function is called the molecular electrostatic potential or MEP. This is a workfunction, and represents the amount of energy needed to remove a bare proton from any specified position close to a molecule out to an infinite distance away. If this energy is positive, the proton was clearly bound to the molecule, if the energy is negative, it was clearly repelled by the molecule. The function needs to be computed for all points in space around the molecule, and it is conventional to calculate an "isovalue surface", or actually two isosurfaces, one representing a positive value of the MEP, the other a negative value. This is then rendered using computer graphics. In this display the negative potential which in effect attracts a proton (i.e. give rise to a positive work function) is rendered in purple. 
Notice how the two faces of the anthracene ring are quite different, an effect induced by the electron withdrawing effect of the CCF_{3} group which withdraws electrons antiperiplanar to the bond, i.e. on the face opposite to it. But in turn the lone pair of the oxygen atom refocuses 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 stereoelectronic 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:
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 QTAIM Critical point analysis  



The NCI Method: Surface based on (reduced) density gradient ∇ρ(r); surface colour based on ∇^{2}ρ(r)  

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 λ_{2} eigenvalue of the Laplacian of the density ∇^{2}ρ(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 dispersioncorrected 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.
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 (Xray) structure to spot a connection between the (antiperiplanar) orientation of (certain) such idealized lone pairs and the CO 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 KlopmanSalem equation, which uses not so much the twoelectron σbond, as a vacant orbital occupying the (assumed) same orientation known as the σ^{*}bond. This was the basis of the earlier lecture course on stereoelectronics 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:
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,2difluoroethane, 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.
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.


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