Putting David Craig’s Legacy to Work in Nanotechnology and Biotechnology*

David Craig^[1] pioneered many advances in the field of theoretical chemistry^[2,3] and had a profound impact on the development of chemical and general scientific research in Australia through his roles in the Research School of Chemistry at the Australian National University,^[4] the Australian Academy of Science,^[5] and elsewhere. After his passing in 2015, he will be surely missed. In addition to his own works, he affected Australian and world science indirectly through the careers of his many students. Two of those include Ian Ross^[6] and Gad Fischer^[7] who were to become my Honours supervisors in 1978. His legacy continues through many people such as myself, and on through the careers of our own students and post-doctoral fellows. I was honoured to receive the David Craig Medal of the Australian Academy of Science in 2016.^[8]

Some of the most fundamental papers in the history of the quantum description of chemistry are the ab initio electronic structure method based on Hartree–Fock theory,^[9] the Born–Oppenheimer approximation that separates nuclear and electronic motion to define the concept of molecular potential-energy surfaces,^[10] the Franck–Condon approximation that allows electronic spectra to be understood,^[11] and its extension the Herzberg–Teller approximation that allows electronic and nuclear motions to couple.^[12] However, the idea that it was possible to perform electronic structure calculations using methods with controllable levels of approximation that eventually led to the exact answer, i.e. the concept of ab initio, was first described by Parr, Craig, and Ross in 1950.^[13] In that paper they performed the first ever calculation of electron correlation, the term neglected in the Hartree–Fock approximation to electronic structure, introducing the configuration-interaction (CI) method. Electron correlation has many consequences that are always significant but which sometimes can be ignored in basic thinking, e.g. in describing chemical bonds, and sometimes are critical, e.g. in understanding the London dispersion force that forms the attractive part of the van der Waals interaction, and in understanding exciton transport in light-harvesting proteins and materials. Similarly, while the 1920s Franck–Condon approximation provided the conceptual basis used for the next 20 years in interpreting molecular spectra, it was Craig who in 1950 performed the first ever quantitative calculation for a polyatomic molecule, determining the geometry of benzene in its excited electronic state.^[14] Finally, while weak observed spectroscopic lines were for 20 years qualitatively interpreted in terms of the Herzberg–Teller principle and its associated symmetry selection rules, Craig in 1950 performed the first ever quantitative calculations, considering the excited states of benzene.^[15] It was Craig that first realised that these weak transitions embodied a breakdown of the Born–Oppenheimer approximation^[15] that, in modern language, resulted in entanglement^[16] of the electronic and nuclear wave functions. Very little modern research in molecular spectroscopy, solar energy harvesting, chemical bonding theory, condensed-matter physics, and materials science does not flow from the advances made by Craig.

This review of my own contributions to science starts from Craig’s work and its impact on my career, combining this with classical and quantum simulation techniques^[17–24] learned from my Ph.D. and postdoctoral supervisors Bob Watts,^[25] Kent Wilson,^[26] and Rick Heller,^[27] my long-term collaborators Noel Hush^[28,29] and Max Crossley,^[29] as well as some of the contributions from the next generation. In terms of the development of theoretical chemistry in Australia, Hush’s contributions were also critical, both through the founding of the Theoretical Chemistry Department at The University of Sydney and through his seminal contributions to electron-transfer theory. That a close relationship existed between the works of Craig and of Hush was apparent from the Piepho–Krausz–Schatz quantitative model^[30] of electron-transfer spectroscopy which applied Craig’s methods^[15,31] to Hush’s revolutionary understanding of electron transfer.^[32] Elmars Krausz had been recruited to the Research School of Chemistry by Craig. My work has focussed on the demonstration of this relationship and the power it gives to understanding modern problems in biotechnology and nanotechnology.

Craig’s understanding of the spectrum of benzene,^{[14,15,33–35]} general aromatic molecules,^[36–40] and non-aromatic^[41–43] molecules paved the way for understanding the properties of molecular spectroscopy and electronic structure in not only the gas phase but also condensed phases.^[44–46] Perhaps the most important chromophore, however, is chlorophyll-a (Chl-a), yet its visible absorption spectrum was not authoritatively assigned until 2013.^[47] The absorption, transport, and emission of light are critical to natural photosynthesis, and it is not possible to understand these processes without knowledge of the spectroscopic assignment. Natural photosynthesis has always been an important model system used in designing artificial devices for solar energy harvesting and conversion either into electrical energy or else to perform important functions such as carbon dioxide sequestration and solar hydrogen production.^[48,49] Given that there is nothing fundamentally new in the spectra of Chl-a and modern light harvesting systems^[50] (except for the direct primary charge-separation mechanism discovered by Hush^[32]) that was not already apparent in Craig’s works, such a lengthy delay was not anticipated.

The basic understanding of the spectroscopic properties of the chlorophylls, porphyrins, and related functional molecules was established in the early 1960s by Gouterman^[31,51–54] and is known as the ‘4-orbital’ model.^[51] This model is just an adaption of Craig’s description of benzene: in benzene four orbitals are also critical, the doubly degenerate HOMO and LUMO orbitals, it is just that in porphyrins and chlorophylls the symmetry is lowered and degeneracies often removed. Craig focussed on the role of antisymmetric vibrations in providing intensity to transitions that are ‘forbidden’ by the Franck–Condon approximation^[11] but allowed through Herzberg–Teller (vibronic) coupling.^[12] Gouterman enhanced this by providing a mathematical formalism for the quantitative treatment of these effects, using it to consider three possible options for the interpretation of the available spectroscopic data for Chl-a and many other chlorophyllides, selecting one. A decade later Ross, who had then also moved to the Australian National University, developed the first method for deducing all of the parameters required by Gouterman’s formulae from a priori electronic structure calculations.^[55,56] My honours project (published 1981) involved applying these methods to weak second-order effects on the spectrum of benzene not considered in Craig’s original treatments, showing that indeed every minute detail of the observed spectra could be understood from first principles analyses.^[57] The highly distorted spectra of 1,5-naphthyridine also succumbed to quantitative analysis,^[58] and it seemed that Craig’s ideas were invincible. It was at about the same time that Krausz at the Australian National University applied Craig’s and Gouterman’s ideas to provide^[30] the first quantitative analyses in Hush’s new field of charge-transfer spectroscopy.^[32] Nobody doubted the assignment that Gouterman had made 20 years earlier concerning the interpretation of the spectra of Chl-a, despite the observation of seemingly peripheral results that were yet to be interpreted.

However, this situation changed in 1982 when the first high-resolution spectra were measured for Chl-a by Avarmaa at Tartu University in Estonia.^[59–62] These spectra were inconsistent with Gouterman’s original assignment of the spectrum of Chl-a and led to a new assignment, a different one to that of Gouterman’s original three possibilities. The basic issues are summarized in Fig. 1. The Gouterman 4-orbital model predicts that all chlorophyllides (and porphyrins) show four electronic absorptions in the visible (or near UV) region. These are known as the Q_x and Q_y transitions in the red-orange region and the Soret B_x and B_y transitions in the blue-UV region. A slight shifting of these transitions changes the perceived colour for purple in bacteriochlorophylls to green in chlorophylls to red in the related heme molecules. One or both Q bands can be weak while the Soret bands are always intense, these patterns arising as a result of strong electron correlation coupling individual transitions together. This critical effect was well known to Craig and Ross as their pioneering calculation of electron correlation in 1950 demonstrated this effect for benzene.^[13] Modern nanotechnological applications to which electron correlation is critical will be discussed also in later sections.

Fig. 1.  Observed absorption spectroscopic band spectrum^[63] (A/ν, blue), emission spectrum reflected about the Q_y origin^[64] (E/ν⁵, green), and magnetic circular dichroism spectrum^[63] (ΔA/ν, red) of chlorophyll-a in ether at 298 K, shown as a function of the energy difference to the centre of the Q_y band. The 1960s ‘traditional’, 1980s ‘modern’, and 2010s ‘new’ assignments of the Q_x origin are indicated. Reproduced from ref. [178], with permission of CSIRO Publishing.

Fig. 1 shows the observed^[63] low-resolution Q-band absorption spectrum of Chl-a in the red–orange region and features three resolved bands: an intense band at low frequency followed by two progressively weaker bands shifted to higher frequency by 1400 and 2150 cm^–1. Naively, both the Q_x and Q_y transitions were expected to show origin bands at low energy followed by a much weaker Franck–Condon allowed vibrational sideband in the range of +700 up to +1600 cm^–1, while Q_y was from all estimates expected to be an order of magnitude more intense than Q_x. That the lowest-energy observed band is the Q_y origin was a feature that was never in doubt. Two of Gouterman’s three possible assignments thus were (i) the observed band at +1400 cm^–1 is the Q_y vibrational sideband while that at +2150 cm^–1 is the Q_x origin, the assignment originally preferred,^[10,13–16] and (ii) that +1400 cm^–1 is made up of both the Q_y vibrational sideband and the Q_x origin, the revised assignment of Avarmaa.^[59–61] Known to Gouterman were significant arguments favouring both assignments: the observed emission spectrum,^[64] also shown in Fig. 1 after reflection about the Q_y origin, is highly depleted in the +2150 cm^–1 region, favouring assignment (i) as emission occurs only from Q_y, whereas the observed magnetic circular dichroism spectrum,^[63] again shown in Fig. 1, shows a negative signal indicating the presence of an x-polarized absorption centred at +1000 cm^–1, favouring (ii). If the Q_x origin is at +1000 cm^–1 then this state will contribute significantly to all spectroscopic photosynthetic processes involving Chl-a, but if it is at +2140 cm^–1 then its effect will be minimal. Hence the issue of the assignment is an important one concerning the understanding of natural photosynthesis. Analyses of photosynthesis have always used the traditional assignment,^[65–68] one reason perhaps being that analyses of data based on the modern assignment never yielded basic properties like the Q_x transition moment and its variation with environment. Without knowing these properties, it was very difficult to use the modern assignment to understand photosynthetic mechanisms.

During the 1980s to 1990s, many new experiments^[63,69–71] were performed and interpreted in terms of the ‘modern’ assignment (ii), and this was backed up by density functional theory (DFT) calculations identifying that the +2150 cm^–1 band was an absorption named S₃ not anticipated by the Gouterman 4-orbital model and therefore never previously considered.^[72,73] So in 2000 the situation again seemed clear, with the ‘modern’ assignment replacing the ‘traditional’ Gouterman one. However, neither assignment was actually able to explain all of the data and, as a result, experimentalists desiring a complete answer continued searching for a definitive experiment. In 2009 Freiberg at Tartu University performed new experiments^[64] that were interpreted in terms of the modern assignment while a year later Krausz designed new experiments^[67] interpreted in terms of the traditional assignment. After 50 years of research, the basic challenges to spectroscopic assignment still remained.

I became aware of the issue following the publication of linear dichroism spectra in 1989 that were clearly inconsistent with the traditional assignment.^[69] Based on my previous experience applying the methods of Craig and Ross to benzene and naphthalene, I was confident that appropriate calculations would interpret all spectra and resolve the issues. I wrote a much enhanced version of the complete neglect of differential overlap (CNDO) code of Del Bene and Jaffè^[74] including expansion to include intermediate neglect of differential overlap (INDO)^[75] and multi-reference configuration interaction (MRCI). This failed to deliver the high accuracy required to be helpful in making the assignment, however, but the developed code today remains well used for understanding excited-state solvent shifts,^[76–79] nanoparticle electronic structure,^[80] and conduction between electrodes spanned by single molecules.^[81] It is freely available on the NANOHUB with over 4000 registered users, applied to teach the basics of solar energy capture, charge transport, and photochemical charge separation.^[82]

The advent of time-dependent density functional theory (TDDFT) a few years later offered the hope of performing calculations accurate enough to assign the spectrum of Chl-a. The first directly relevant publication was on porphyrin and claimed high accuracy,^[83] but the assignment was inconsistent with the Gouterman 4-orbital model and to me seemed unrealistic. The influential subsequent works on Chl-a^[42,43] had similar properties but now the deviation from the 4-orbital model was being used to assign the spectra. Believing something to be fundamentally wrong, I searched for simple ways of isolating and demonstrating the problem,^[84] discovering a significant shortcoming of TDDFT when H₂ was stretched (failure owing to ground-state instability),^[85] that it failed when charge transfer was involved owing to the asymptotic potential error,^[86] and most significantly that it failed for extended conjugation (including polyacetylene, porphyrins, and large acenes) despite many reports to the contrary.^[87] The culmination of this work was the demonstration that all existing TDDFT methods predicted results for plant Photosystem-I that was starkly in contrast with experiment, showing that the non-Gouterman states must actually occur at much higher energies.^[88] Realising the widespread significance of these result to modern technological applications, Kobayashi and Amos from the National Computational Infrastructure (NCI) recognized the potential of a proposed new computational method CAM-B3LYP^[89] and programmed it into GAUSSIAN.^[90] Together we showed that it indeed did work for extended conjugated molecules including chlorophylls.^[91] Today this method is widely applied to chemistry, biotechnology, and nanotechnology applications, the method currently having over 3500 citations.

Over the next few years, we demonstrated that CAM-B3LYP could be used to interpret many sophisticated observations made concerning photosynthetic systems, including a priori prediction^[68] of the visible absorption spectrum of Photosystem-I based on an optimized geometry for the 150 000-atom protein trimer that had previously been fully optimized using a linear-scaling DFT method we had developed.^[92] Also, Freiberg was able, for the first time, to develop methods allowing comparison of both high-resolution absorption and emission spectra for chlorophyllides, obtaining results unprecedented in molecular spectroscopy, that CAM-B3LYP was able to predict to high accuracy.^[93] However, three years’ work by Zhengli Cai attempting to use CAM-B3LYP to assign the spectrum of Chl-a produced no positive conclusions. The approach taken was to calculate the properties of some 32 chlorophyllides and their environments, comparing the predicted Q_x–Q_y spacing to that inferred by the traditional and modern assignments applied to each system. No correlation was found between the calculated data and either assignment,^[47,94] suggesting that either they were both incorrect or else CAM-B3LYP was still not precise enough to be used for the evaluation of this property.

In 2012, Freiberg, Krausz, and myself got together to search for some key element that had possibly been missing from the research of the previous 50 years. Reviewing the extensive literature coming from both Eastern-Block and Western-Block countries, we found no significant error in logic or methods leading to conclusions favouring either the traditional or modern assignments. There was one remaining option, however: Gouterman had suggested three possible assignments but only two had been seriously advanced. The third assignment can be thought of as a quantum-mechanical supposition of the other two assignments, assignments based on application of the Born–Oppenheimer principle that allows Q_x and Q_y to be treated as independent potential-energy surfaces using say classical mechanics. Craig’s original works on aromatic molecules had shown that this was not strictly true and that breakdown of the Born–Oppenheimer approximation had significant repercussions on spectra.^[15] Indeed, my honours project was based on quantifying these effects using CNDO calculations,^[57,58] but the degree of breakdown was small and the effects only perturbative. Gouterman’s enhancements to Craig’s ideas^[31,52] allowed for strong interactions and the complete breakdown of the Born–Oppenheimer equation. The great attraction of this possibility is that it can deliver two bands of x polarization, like what is seen in the MCD spectra in Fig. 1. However, if this effect were operative in Chl-a then it would also have consequences for all other chlorophyllides including highly categorized molecules like pheophytin-a and bacteriochlorophyll-a, and no trace of this effect had ever been reported in the spectra of these and many other chlorophyllides. Reviewing the literature, Gouterman’s argument against this assignment was as valid in 2012 as it was in 1962.

We reasoned that the missing bands must be there in the spectra of all these other chlorophyllides but be obscured by more intense features. Two data analysis methods were developed to reveal them, a computational method for fitting the observed spectra to Gouterman’s model appropriate for very strong coupling,^[47] and an experimental method providing a new way of qualitatively interpreting data coming from MCD experiments.^[95] Both approaches identified the missing bands in all species, with the fitting procedure delivering quantitative information concerning the unperturbed location of the Q_x origin (the property deduced by electronic-structure calculations such as CAM-B3LYP), the Q_x intensity, and the vibronic coupling strength between Q_x and Q_y. These results were in excellent agreement with the previous CAM-B3LYP calculations, and so a new assignment of the spectrum of Chl-a was established.^[47]

The computational fitting procedure embodied a combination of the condensed-phase spectroscopic modelling techniques developed for liquid water and ice during my Ph.D. project with Bob Watts^[17–20,96] with the time-dependent quantum dynamics techniques that I learnt from Eric Heller as a postdoctoral fellow^[21] to aid^[97] the design of NASA’s GALILEO probe to Jupiter.^[98] The core infrastructure for this was built a decade earlier during work with Noel Hush to understand the nature of the special-pair radical cation produced after primary charge separation during bacterial photosynthesis.^[99] This more complex problem involved four electronic states strongly coupled by 70 molecular vibrations,^[100] with the highly quantum non-Born–Oppenheimer properties of the system controlling all observables including the (functionally critical) output voltage, spin distribution, electronic and infrared spectra, and Stark effect.^[99,101,102] Our work on the bacterial reaction centre was also a long quest to understand perplexing experimental data, this time the observation of an apparently intervalence electronic transition peaking at 2700 cm^–1 amongst the normal infrared lines.^[103] It was clear that this spectrum could reveal the critical features controlling primary charge separation during photosynthesis, but its interpretation involved many aspects including reassignment of other observed band systems,^[104] identifying other new band systems,^[105] extensive simulations of mutagenesis effects,^[106] and the identification of long-standing errors with the interpretation of electron-paramagnetic spectroscopy data.^[101] This project was originally conceived as flowing from Hush’s classic theory of adiabatic electron transfer applied to intervalence spectroscopy,^[32] but, as discussed in the next section, its relation to the works of Ross, Krausz, and Craig soon became apparent.^[107,108]

The developed qualitative procedure for the interpretation of MCD spectroscopy involved finding a work around to a fundamental problem that limited the usefulness of the technique. Two independent spectra (absorption and MCD) are measured at n frequencies, while what is desired are component spectra depicting x and y polarization at each frequency, but two susceptibility parameters must also be determined from the spectroscopic data in order to construct an analytical data inversion. The data analysis problem is thus fundamentally underdetermined: there are 2n observables and 2n + 2 unknowns to evaluate. We showed that one of these unknowns can be determined based on other assumptions already used in the analysis, and that it was possible to construct a second unknown whose value acts to rescale component spectra but not significantly change their shape.^[95] The result is a general method that can be used to get critical polarized bandshape information from MCD experiments. Had Gouterman had these simple equations in 1962, he would no doubt have immediately correctly assigned the spectrum of Chl-a.

One last problem remained for the acceptance of the strong vibronic-coupling assignment for the spectrum of Chl-a: high-resolution spectra measured in ether at low temperature, spectra whose qualitative features led to the initial rejection of the traditional assignment,^[59,61] showed quantitative aspects that, from the beginning, were obviously incompatible with a strong-vibronic-coupling assignment. This issue was solved by Krausz and Freiberg who showed that the presence of trace water in the samples meant that the observed spectra were for an ether–Chl-a–water complex rather than for the presumed ether–Chl-a–ether complex.^[109] Extensive free-energy simulations performed by Kobayashi then confirmed that indeed it is feasible for water to displace the solvent when doubly dilute with Chl-a in solution at low temperature.^[110]

What emerges then is a rather unexpected description of photosynthetic machinery as working with an unusually quantum nature. Both the excited states of Chl-a critical for solar-energy harvesting and its transport and the primary charge separation unit from bacterial photosynthesis are subject to extreme failure of the Born–Oppenheimer approximation,^[47,100] the idea that gives us the conceptual basis of modern chemistry. This obviously affects critical features such as the effect of environment and mutagenesis on output voltage,^[100,101] spectra,^[102] and excitonic energy flow.^[47] However, there are many other fascinating consequences and associated features: Chl-a in protein environments is optimized to decohere Q_x excitation as fast as could be possible,^[47] asymmetry unprecedented in molecular spectroscopy is observed between absorption and emission of isolated electronic states,^[93] and the intensity of the Q_x transition varies by a factor of 7 as a function of chlorophyllide and solvent.^[47]

Some time after Craig moved to the Australian National University, Hush moved to Sydney University where he established the Department of Theoretical Chemistry. Their parallel teaching and research platforms in theoretical chemistry were to have a profound impact on chemistry research in Australia. Hush had pioneered many advances in understanding chemical reactivity, paralleling Craig’s work in spectroscopy. Although these research fields both involved breakdown of the Born–Oppenheimer and Franck–Condon approximations as core elements, their differences were sufficient for them to be considered as independent research fields.

To understand why this would be so, we need to trace the origin of Hush’s theories back to the beginnings of quantum mechanics and the understanding of chemical reactivity. The first treatise on reactivity was by London in 1928 using a diabatic description^[111] in which simple potential-energy surfaces described the reactants and the products coupled together to make an adiabatic transition state. The London–Eyring–Polanyi (LEP) potential-energy surface^[112] for general triatomic molecular reactions^[113] was extended from this by Eyring and Polanyi in 1931, Horiuti and Polanyi made a similar extension for proton transfer reactions in 1935,^[114] Wall and Glocker did the same for isomerization reactions in 1937,^[115] and Hush fully generalized it in 1953, considering e.g. redox hydrogen transfer reactions.^[116] Non-adiabatic reactions depicting failure of the Born–Oppenheimer approximation were similarly treated in 1932 by London,^[117] with the more well known but more restricted treatment of Landau also appearing in that year.^[11,12] Today diabatic models are widely used to depict most chemical processes including electron-transfer,^[118–121] natural photosynthesis and organic photovoltaics,^[122–133] chemical quantum computing,^[134,135] racemization processes,^[136,137] non-adiabatic photochemical reactions,^[138–148] proton-transfer, hydrogen bonding, hydrogen-transfer, coupled electron–proton transfer reactions, general reactions,^[149] and aromaticity.^[150–154] An essential aspect of this modern pre-eminence is its connection with Hush’s adiabatic theory of electron transfer.^{[32,155–159]} This approach differed significantly from other approaches to electron transfer such as that developed by Marcus^[160] as Hush focussed on the physical properties of the involved transition states and the use of experimental and ab initio or a priori computational methods to quantify analyses. A centrepiece of Hush’s work was the assignment of the nature of unknown electronic transitions observed in mixed-valence species, leading to the field of intervalence spectroscopy.^[32] This process was identified by Hush as producing primary charge separation, converting solar into electrical energy^[161] and providing a basis for the subsequent understanding of the operation of natural photosynthesis and most processes involved in modern molecular-based artificial solar-energy harvesting systems.

It is Hush’s work on intervalence spectroscopy^[32] that builds a bridge to Craig’s work. The connection was first realised by Krausz in his 1978 quantitative model for the shapes of intervalence transitions.^[30] My first collaborations with Hush came through this connection when I started as an ARC Research Fellow at Sydney University in the late 1980s. Hush was collaborating with Michael Paddon-Row at the University of New South Wales, synthesizing, measuring (by Jan Verhoeven in The Netherlands) primary charge separation rates, and interpreting the data.^[161] An unexpected result, questioning the conceptual understanding of the process, was that electronic symmetry-forbidden reactions were occurring only slightly slower than Franck–Condon allowed ones. My contribution was to apply the spectroscopic methods of Craig, Ross and Fischer to understand the kinetics process, showing that the observed fast rates were completely consistent with known spectroscopic properties.^[162] This led to predictions of unusual properties that were later observed spectroscopically.^[163] Ideas from spectroscopy were then imbedded into general theories for electron transfer,^[164–166] resulting in the development of a chemical model^[167–169] for interpreting Stark-effect measurements pioneered by Steven Boxer at Stanford for understanding charge transfer processes in general and photosynthetic function in particular,^{[102,170–175]} as well as an understanding of electric-field driven molecular switches.^[176] Parallel work included simulations of the effects of solvent dynamics on molecular spectra^[78,79] and the importance of microheterogenuity in modifying electric fields and optical properties.^[177]

However, the major focus of the decade became interpreting the spectroscopic and chemical properties of the special-pair radical cation from photosynthesis, as mentioned earlier. Hush’s original description of intervalence spectra^[32] and Krausz’ computational model^[30] had focussed on the reaction coordinate as the critical dimension controlling properties. Using the power of unitary transformations in quantum mechanics, I developed a way of transforming the full spectroscopic problem discussed initially by Craig in terms of the effects of symmetric^[14] and antisymmetric^[15] modes into the kinetics representation used by Hush, providing a complete description of the kinetics and spectroscopy of intervalence charge transfer.^[107,108] This allowed the bacterial reaction centre cation to be understood completely,^[99,102] facilitated assignment of the spectrum of Chl-a,^[47] and led to the understanding of a range of phenomena where properties were determined by more than just the reaction coordinate, as has been recently reviewed.^[178] It unified critical aspects of the works of Craig and of Hush, playing on the non-uniqueness of diabatic descriptions for chemical and spectroscopic phenomena.^[179] Unification of basic theories for chemistry and spectroscopy touches the very foundations of chemical understanding. Traditionally, chemistry is taught in separate blocks using different nomenclatures and ideas, e.g. transition state theory, aromaticity, coordination compounds, crystal field theory, Walsh’s rules, valence-shell electron-pair repulsion theory (VSEPR), etc., and the teaching is quite different to approaches taken in condensed-matter physics for describing metals, insulator–conductor transitions, magnetism, etc. Yet inherently the same underlying principles apply.

Recently, we have been able to make progress in finding a unified approach to describing all these phenomena based on an underlying diabatic model, augmented where necessary for complicating effects.^[180] This makes practical the dreams of the original 1920s–1930s developers of diabatic models. This work was done in conjunction with Ross McKenzie (Physics, University of Queensland) who brought his interest in valence-bond processes and quantum coherence, as well as his vast experience in studying coupled quantum nuclear-electronic problems in solid-state physics.^{[135,180–186]} Many key advances in solid-state physics^[187–189] were made at around the same time as Hush’s theory of adiabatic electron transfer was developed, and much was to be gained by merging results from the two disciplines. Also, much of the work was performed by our student Laura McKemmish as undergraduate research projects and during her honours year.^{[135,180–186,190]}

Previously, diabatic models had been very successful only for electron-transfer theory. For it, a simple model with just two or three parameters could qualitatively describe a wide range of spectroscopic and kinetics properties.^[180] Analogous quantitative treatments of chemical properties based on say pseudo Jahn–Teller approaches^[19,94–96] have been very successful but require parameterization, often based on accurate ab initio calculations, for each property of interest. These approaches focus on two diabatic surfaces and the conical intersection that links them, by analogy to electron-transfer theory, but the deduced parameters are not transferrable either for the evaluation of different properties of the same system or else between similar systems. Hence they could describe important fine experimental details well but could not provide a fundamental basis for chemical understanding.

The central feature of electron-transfer theory not replicated in diabatic treatments of other systems was essentially the inability to transfer parameters. Fig. 2 shows the fundamental difference in electronic structure that caused this. Fig. 2a shows the situation for electron transfer, with an analogous one existing for hole transfer. Basically, it involves one electron (or hole) distributed amongst two orbitals, giving rise to two possible electronic states. Antisymmetric vibrations then couple the two resulting diabatic electronic states, producing all the kinetics and spectroscopic properties of electron-transfer theory.^{[30,32,107,108,156,159,167,168,180,191]} General chemical reactions occur between closed-shell species, however, typically involving two electrons being redistributed amongst two orbitals, as sketched in Fig. 2b. This scenario leads naturally to three diabatic electronic states, the ground state G, the singly excited state S, and the doubly excited state D. Antisymmetric vibrations couple G to S in analogy to what happens in electron transfer, but the same vibrations necessarily also couple S to D. Therefore there exists one low-lying conical intersection seam for the G to S interaction and the process can be considered as a pseudo Jahn–Teller effect,^[19,94–96] but there also exists a second seam linking S to D. This seam changes to properties of S significantly, just as the first seam affects G and S, and hence by an indirect mechanism D influences G. If one assumes that the vibronic coupling between S and D is the same as that between G and S (in practice small differences arise owing to the effects of electron correlation), then a renormalization of the full set of coupled equations is possible that allows the properties of the ground-state to be re-represented in terms of a simple two-state model analogous to electron-transfer theory.^[180] This tells why two-state models involving only one conical intersection seam have been so successful in numerical applications, but it also indicates that the parameters appearing in such approaches get renormalized from a more basic set of parameters in different ways. Hence the underlying simplicity is lost and model parameters cannot be transferred.

Fig. 2.  Molecular-orbital diagrams showing representative unique electron occupations for the diabatic ground states G, singly excited states S, doubly excited states D, triply excited states T, and quadruply excited states Q, where appropriate, for (a) electron transfer reactions involving radicals, (b) reactions between closed-shell species without orbital degeneracy, and (c) reactions between closed-shell species with doubly degenerate HOMO and LUMO orbitals. Note that only one spin component of each state is shown, as is one spatial component of the appropriate symmetry; the G, D, and Q states have one symmetry (totally symmetric for closed-shell reactions) while the S and T states have the opposite symmetry, the product of these symmetries being the symmetry of the coupling vibration. Reproduced from ref. [180] with permission of the PCCP Owner Societies.

For aromatic molecules, the HOMO and LUMO orbitals controlling chemistry and spectroscopy become either doubly degenerate (as is the case for benzene) or else nearly so. Fig. 2c sketches the scenario for benzene, it being possible to place the available four electrons in the four frontier orbitals in seven ways with singlet spin. Three of these ways are double excitations D, plus new triplet T and quadruple Q excitations are possible. If the three D states are replaced by one and the vibronic coupling connecting G to S, S to D, D to T, and T to Q are taken to be the same, then renormalizations of the five-state coupled problem in the strong coupling and weak coupling limits can be constructed that again map this problem involving many interacting conical-intersection seams back to a simple two-state one-seam situation.^[180] The notion thus naturally appears that the ground state is ‘twinned’ to an excited state, just as the ground and excited states in electron transfer are intimately connected. This concept in regards to aromaticity has been highlighted for a long time by Shaik.^[150–154] However, quantitative analysis akin to that applied in electron-transfer theory was never attempted. The reason for this is that Shaik mistakenly assumed that the excited state about which much was known, the singly excited state S, was the twin state whereas in fact it is the poorly characterized quadruply excited state Q; S indeed also has a twin state, but it is T not G.^[180]

These results allow the theoretical description of aromaticity, general chemical reactions, and electron transfer to be unified.^[180] Fig. 3 shows the situation for a closed-shell reaction involving isomerization, taking as an example the inversion reaction of ammonia.^[135,180] The ground-state Born–Oppenheimer adiabatic potential-energy surface fitted to experimental data by Swalen and Ibers in 1962^[192] (which is very similar to modern estimates^[193,194]) is shown as green circles and fitted to a diabatic model. This model is described equivalently in the figure in two ways, (i) using the electron-transfer description of Hush as left (L), central (C), and right (R) localized diabatic potentials, and (ii) using the spectroscopic description of Craig, Ross, etc., of ground-state (G), single excited state (S), and double excited state (D) diabatic surfaces. Both diabatic approaches then yield identical descriptions of the ground-state (g), single-excited (s), and double excited (d) Born–Oppenheimer adiabatic surfaces. All potentials are specified by the parameters J (the resonance energy) and λ (the reorganization energy), as shown in the figure. The reorganization energy is related to the angle 2^1/2τ_m of the minimum of the localized diabatic L and R potentials through λ = 2kτ_m², where k is the force constant. As atomic masses are rarely varied, knowing k is the same as knowing the vibration frequency ω for each of the diabatic potentials. The ammonia inversion reaction is a symmetric reaction between equivalent isomers, but symmetry is not an essential property. If the reaction is exothermic then a parameter E₀ must be introduced to indicate the energy difference between the L and R diabatic minima. All chemical and spectroscopic properties are then determined from the four parameters J, λ, ω, and E₀, including the transition state energy shown in Fig. 3. For reactions involving general molecules, the specific torsional reaction coordinate τ used for the ammonia inversion reaction is replaced by a generalized dimensionless normal coordinate Q.

Fig. 3.  Diabatic models for the ammonia inversion reaction. The lower figure shows the Swalen and Ibers adiabatic potential for the torsional motion () of NH₃ (green circles) fitted up to a torsional angle of τ = 44° (NHN bond angle θ = 77°) by a three-state three-parameter diabatic-potential model ( = 4.13 eV, λ = 10.55 eV, τ_m = 25.1°, RMS error 1.2 meV): purple- resulting adiabatic ground-state (g), singly valence-excited state (s), and doubly valence-excited state (d) (the ‘twin state’); red - localized representation of two of the diabatic potentials; blue - delocalized representation of the two of the diabatic potentials (both ). Reproduced from ref. [180] with permission of the PCCP Owner Societies.

Fig. 4 and Table 1 depict 10 quite different applications of these ideas to modern chemical systems, including not just the ammonia inversion reaction^{[135,180,186]} but also aromatic resonance in benzene,^[135] breakdown of aromaticity in the triplet (π,π*) ground-state of pyridine (³PYR),^[135,195] the Creutz–Taube ion in aqueous solution (CT),^[196] possibly including an orthomethoxy substituent on the pyrazine ring (CT-OMe),^[197] a variant of this in which pyrazine is replaced by a dipyridyl polyene (DPP),^[198,199] the bacterial photosynthetic reaction centre from Rhodobacter sphaeroides (PRC),^[99,100] gas-phase BNB,^[200] pertinent to organic light-emitting diodes (OLED)s and artificial photovoltaics the fastest hole-transfer process in the β-phase crystal of the molecular conductor mer-tris(8-hydroxyquinolinato)aluminium(iii) (Alq3),^[201,202] and a ferrocene–porphyrin–fullerene model triad molecule (FcPC₆₀) designed by Crossley that undergoes long-range photochemical charge recombination.^[203,204] Limitations of the application of the one-mode model to these molecular systems can be significant and are discussed elsewhere,^[185] but the model does capture the critical chemical ideas.

Fig. 4.  Some sample molecular systems (see text) with electronic states that can be described using two coupled diabatic potential-energy surfaces. OMe is methoxy, PHY is phytyl; ^tBu is tertiary butyl; A is ammonia; FcPC₆₀ is Zinc, [[5,10,16,21-tetrakis[3,5-bis(1,1-dimethylethyl)phenyl]-13-[4-(1′,5′-dihydro-1′-methyl-2′H-[5,6]fullereno-C₆₀-I_h-[1,9-c]pyrrol-2′-yl)phenyl]-1,12-dihydro-23H,25H-diimidazo[4,5-b:4′,5′-l]porphin-2-yl-κN²³,κN²⁴,κN²⁵,κN²⁶]ferrocenato(2-)]-, (SP-4–1); Alq3 is mer-tris(8-hydroxyquinolinato)aluminum(iii); DPP is ruthenium⁵⁺, decaammine[μ-[4,4'-[(1E,3E)-1,3-butadiene-1,4-diyl]bis[pyridine-κN]]]di- (9CI). Reproduced from ref. [185] with permission of the PCCP Owner Societies.

Table 1.  Diabatic-model parameters describing the chemical and spectroscopic properties of the 10 molecular systems shown in Fig. 3

Table 1 also presents five derived characteristic quantities depicting system properties. The ratio 2|J|/λ controls whether the ground-state Born–Oppenheimer surface is double welled (hence allowing reactions from one well to the other over the intervening transition state) or single welled (as in benzene, where the ground-state structure is a quantum resonance between two classical structures). For symmetric systems (E₀ = 0) this parameter alone controls this localization or delocalization^[205–207] and hence the class of a mixed-valence system,^[208] but general expressions are also known.^[184] For the molecules considered it ranges from 0.029 for the non-adiabatic electron transfer charge recombination process in FcPC₆₀ to 0.8 in CT and ammonia to 3.3 in benzene. The parameter

is a global measure of the electronic energy spacing, while is the vibrational energy spacing, so the ratio provides a basic guide to the suitability of the Born–Oppenheimer approximation, this becoming untenable as . In this limit, the Jahn–Teller effect most likely becomes important and so the inclusion of at least two vibrational modes becomes critical to providing a qualitatively useful description of the system. Of all the systems investigated, the PRC is the one most likely to suffer from Jahn–Teller effects, with , though non-Condon effects are also known to be critical for CT.^[196,209] The driving force of a reaction is best thought of in terms of the number of vibrational quanta involved, .

Born and Oppenheimer’s original consideration of ways to separate out nuclear and electronic motion employed an expansion in terms of the quarter root of the ratio of the electronic and atomic masses, a term roughly proportional to . Failure of the Born–Oppenheimer approximation as the Jahn–Teller region is approached is therefore anticipated. It is also known that the Born–Oppenheimer approximation fails locally whenever diabatic surfaces become near degenerate, as is the case when . Hush, however, argued that the critical parameter was the diameter Q_c of the cusp formed at the intersection point of the diabatic surfaces,^[156] and we were recently able to quantify this effect.^[185] For the simple one harmonic mode, two-state problem, the cusp diameter is given by

which unites the two basic ideas into one expression, with the Born–Oppenheimer approximation failing whenever Q_c < 8^–1/2 ≈ 0.35. Table 1 shows that, for classical chemical processes like aromaticity in benzene and ammonia inversion, the Born–Oppenheimer approximation holds well, with Q_c = 12.6 and 6.34, respectively. Perhaps unexpectedly, it even appears useful for electron transport through Alq3 (Q_c = 0.67). This process is, however, the fastest electron process found between molecules on different sites in an Alq3 crystal, with most other processes being much slower and hence smaller J and Q_c, although the failure of non-adiabatic limit theories to describe the fast processes in the crystal has been noted.^[202]

Simple spectroscopic quantities to flow from the equations are the electronic transition energies, dipole moment changes Δμ, and transition dipoles M. For symmetric reactions (E₀ = 0) the electronic absorption vertical excitation energy is given by λ if 2|J|/λ < 1 and by 2|J| otherwise,^[180] while

links the transition moment and dipole change (observable in Stark-effect measurements) to the distance R between centres in a charge-transfer application.^[108,167] In the limit of weak coupling this gives the transition moment as^[32]

The vibrational frequency in the ground state is given by

if the lower Born–Oppenheimer surface is double welled, while either the barrier frequency (if double welled) or the vibration frequency (if single welled) is given by

and the excited-state (the ‘twin state’) frequency is given by

The most important kinetics property is the reaction barrier height which is given in adiabatic electron-transfer theory for symmetric systems as

in the weak coupling non-adiabatic limit as^[156,210]

and in general is close to^[185]

In the adiabatic limit the reaction rate is given by transition-state theory as^[156]

whereas in the non-adiabatic limit it is given by the Golden Rule^[210,211] as

to emphasize the role of the cusp diameter.^[185]

The one-mode two-state model is therefore characteristic of many chemical processes. It yields the above-mentioned analytical results plus many more,^{[101,108,167,168,180,185]} and in addition can easily be solved numerically to high accuracy for any property without use of the Born–Oppenheimer approximation. Hence it can provide a means of assessing the reliability of computational schemes used in non-adiabatic calculations. Application of the Born–Oppenheimer approximation requires the neglect of three terms in the total nuclear-electronic Hamiltonian that scale with either the nuclear momentum or the nuclear kinetic energy. These are usually called the diagonal correction (DC), the first-derivative correction (FD), and the second-derivative correction (SD).^[212,213] The one-mode two-state model leads to analytical expressions for these corrections^[185] that are examples of known multi-mode relationships.^[214,215] Their values depend on nuclear geometry but attain maxima of^[185]

The FD term is thus lowest in order of 1/Q_c and could therefore be expected to be the dominant correction as most of the time the Born–Oppenheimer approximation works well. However, this is not the case as Born–Oppenheimer breakdown is a specific example of a pitchfork bifurcation Cusp Catastrophe,^[216–219] a phenomena known to be not treatable using such perturbation expansions as it embodies a fundamentally unstable mathematical description. These general conclusions have profound consequences when applied to the one-mode model that are easy to demonstrate. The DC term only acts to change the electronic energy levels, allowing wave functions to still be written as simple products of electronic wave functions and nuclear wave functions with retained adiabatic separation of variables. It can therefore be simply added to standard Born–Oppenheimer potential-energy surfaces as calculated by all electronic structure codes. This correction is very spiked at the cusp, which is located at or close to the transition state, and serves to increase the barrier for chemical reactions by the amount ΔΔE^† given above. Values for the 10 chemical systems considered are given in Table 1, expressed relative to the vibrational energy spacing. What is clear is that in non-adiabatic reactions, ΔΔE^† can be larger by far than the original activation energy ΔE^†. Quantum tunnelling through the high but very narrow barrier generated by DC negates this effect to some extent, but we find that over the entire parameter space of the model there is no region in which inclusion of only the FD term leads to a reliable solution of the quantum dynamics equations depicting the reaction process.^[185] Always all three corrections must be used to obtain reliable results. This observation is expected from general considerations of the cusp catastrophe but is contrary to standard practice for the determination of non-adiabatic reaction rates. Success over the last decades by methods including only FD has arisen owing to the accidental cancellation of large errors stemming from the use of an approximated Hamiltonian and its solution by semiclassical rather than full quantum treatments.

The usefulness of this unification of chemical and spectroscopic processes under the one diabatic umbrella is perhaps not so much in the description of important but nevertheless not typical non-adiabatic chemical processes but rather in its ability to unite standard chemical understanding. The difference between an isomerization process like the ammonia inversion reaction and aromaticity in benzene is seen to arise just through a change in the value of 2|J|/λ. This is > 1 in benzene so all C–C bond lengths are equal but < 1 in ammonia so the planar D_3h structure is unstable and pyramidal structures form.^[180] This suggests that the method can be used to provide an explanation of basic chemical structural properties from first principles.

The most fundamental approach to gaining insight into the structure of molecules like XY₃ series such as NH₃ to BiH₃ is by application of Walsh's Rules^[220–227] using understanding of how occupied orbitals change with molecular shape and electronegativity differences. Excited state properties can also be framed in this context^[221,228] and reaction symmetry understood,^[229–233] but the close connectivity^{[139,150–154,234,235]} between ground-state and excited-state properties is lost. To apply the methods, one needs to know beforehand how orbitals behave as a function of geometry, something that may not always be intuitively obvious. A more widely applied but more empirical approach to describe bonding is VSEPR theory^[236–239] and its origins.^[240,241] This approach considers electron pairs associated with the same atom that repel each other, while being attracted to nuclei differently depending on electronegativity. It applies only to ground states and cannot depict excited-state properties. In its original form, VSEPR was applied like this: the XH₃ series has four valence electron pairs and so they are presumed to be arranged tetrahedrally around the central atom making the expected HXH bond angle θ = 109.5° but compressed slightly as lone-pairs take more angular space than do bonding pairs owing to their contraction to the nucleus. Hence the observed HNH angle of 108° in ammonia is explained, but this argument does not explain why the bond angles in PH₃ to BiH₃ range between 93°–90°.^[186] Traditionally, the empirical rule that second-row elements are always octahedral was therefore invented.^[237,239] In the latest version of VSEPR theory, atomic radii are used for every pair of atoms, neglecting the traditional electronegativity arguments used in the earlier approach.^[239] Essentially, VSEPR theory now parameterizes the answer for each problem into the method.

Diabatic descriptions instead offer the description of chemical structures in terms of fundamental properties: the resonance energy and the reorganization energy. It is the factors that control these quantities that control geometry. Considering orbital overlap with the diabatic references orbitals allows the diabatic bond angle in the absence of resonance to be determined. This process is illustrated in Fig. 5. For NH₃⁺, the interaction involves one electron in two orbitals (the N lone pair orbital and the antibonding symmetric NH orbital). The diabatic orbitals are obtained by a 45° rotation of the adiabatic orbitals at the high-symmetry planar geometry, identifying them as sp hybrids (Fig. 5a). In the absence of resonance, the three NH bonding orbitals must be orthogonal to a single one of the two sp diabatic orbitals and so the bond angle becomes θ = acos[3sin²(atan(1/2))/2 – 1/2] = acos[–1/5] =  101.5° (Fig. 5b). Renormalization for the two-electrons in the two-orbitals case of NH₃ increases the associated torsional angle by a factor of 2^1/2 and hence gives a bond angle of θ = acos[3sin²(2^1/2atan(1/2))/2 – 1/2] = 86.7°.^[186] However, resonance drives both of these limiting structures towards the delocalized aromatic-like situation of a planar molecule with 120° angles, achieved whenever 2|J|/λ ≥ 1. Interpretations of the excited-state spectroscopic data of these molecules supported by ab initio calculations indicate that both the resonance energy and the reorganization energy are much larger for N compounds than for P–Bi compounds, but the effect on the resonance energy is second order so that the values of 2|J|/λ are 0.82 for NH₃ but 0.48–0.36 for PH₃–BiH₃, respectively.^[186] It is this discontinuity in 2|J|/λ that causes the geometry of NH₃ to be tetrahedral while the other molecules are octahedral.

Fig. 5.  (a) Modified Walsh diagram showing how sp hybridized diabatic orbitals on X interfere to produce pure p and s adiabatic orbitals at the planar geometry, and how the adiabatic orbitals decouple as XH₃ molecules distort. (b) An X sp diabatic lone-pair orbital and its three orthogonal s^1/2p^5/2 bondable hybrid orbitals at a torsional angle of τ = 26.6° and HXH bond angle θ = 101.5°.). Reproduced from ref. [186] with permission of the PCCP Owner Societies.

Understanding the sources of the large increases in the resonance energy and reorganization energy for N compounds therefore provides the way to understand the bonding from a fundamental perspective. As Fig. 3 shows, these quantities are not just properties of the ground state but are also excited-state properties. The nature of the excited states of molecules is thus intricately related to the properties of their ground states. To provide a good description of the ground-state potential-energy surface, computational methods must therefore also be able to describe the involved excited states. Unoccupied orbitals and their properties are just as important as occupied orbitals when it comes to understanding potential-energy surfaces.

The big difference between the orbital structures of NH₃ and PH₃ is an inverted ordering of the XH symmetric antibonding orbital σ*_A and the lowest Rydberg orbital: in ammonia the Rydberg orbital is lowest in energy whereas in phosphine the valence orbital is lowest.^[186] This change may not seem profound but the consequence is a change in orbital nodal pattern as the two orbitals interact to produce ‘bonding’ and ‘antibonding’ molecular orbitals. If σ*_A is of lower energy then the interaction with the Rydberg orbital stabilizes it, spreading out the electrons over a larger volume, but if it is higher in energy then a nodal surface divides the two, compressing σ*_A on top of itself (Chart 1).

Chart 1. 

This compression results in a large increase of electron–electron repulsion within the same electron pair, forcing up the energy and hence the resonance energy of two interacting sp hybrids. Hence the swap in order of the σ*_A orbital and the lowest Rydberg orbital is responsible for the large change in properties between the first-row element N and the other elements of Group XV. The compression is clearly seen in Fig. 6 where the Hartree–Fock orbitals are shown for NH₃ and AsH₃ obtained using a STO-3G basis, as well as that when augmented by a single Rydberg s-type basis function. Adding the Rydberg orbital compresses the ammonia orbital significantly but has little effect on stilbene.

Fig. 6.  Isodensity surfaces of the unoccupied σ*_A valence orbital of NH₃ and AsH₃ that controls hybridization, at its planar D_3h structure, evaluated using the Hartree–Fock method. Reproduced from ref. [186] with permission of the PCCP Owner Societies.

VSEPR theory provides no fundamental insight into chemical bonding, with its focus on bonding electron pairs causing it to not see the critical features that make first-row compounds different from those of later rows. Walsh’s rules allow ground-state and excited states to be simultaneously interpreted, but the interconnecting properties are omitted and one therefore must determine the key properties individually for each system. Only diabatic analysis can unify a wide range of spectroscopic and chemical properties to provide fundamental physical insight.^[186]

Max Crossley at Sydney University has synthesized many molecules for application in molecular electronics, understanding photosynthesis, and biomedical applications, mostly based on porphyrin chemistry, see Fig. 7 for example. The essential element to the success of his work has been the ability to provide tight synthetic control over structure and function through making synthetically diverse but structurally rigid molecules. Such control allows theories to be developed to describe observed properties, as well as the application of high-level ab initio and density-functional computational methods. An example of this is FcPC₆₀ highlighted in Fig. 4 and Table 1, presenting results obtained by combined synthesis, spectroscopic observation, and fundamental understanding. This and related molecules, at the time of original publication, achieved the longest-lived primary charge separation of any triad, up to 1 ms at room temperature in solution,^[203] following earlier record-breaking performances of related dyad molecules.^[204] Long-lived charge separation is a highly desired property in any artificial device that harvests solar energy and converts it into either electrical or chemical energy.

Fig. 7.  Some of the molecules synthesized by Crossley that led to fundamental understanding in molecular electronics, solar energy harvesting, and related fields. The structures and photochemical properties of the bacterial photosynthetic reaction centre (PRC) and its synthetic analogue (MP)₂Q-Q2HP-C₆₀ are reproduced from ref. [242] with permission from the Royal Society of Chemistry.

Recent advances on this technology include the synthesis of the tetrad molecule (MP)₂Q-Q2HP-C₆₀ (Fig. 7) with structural and functional properties that are the closest-known mimics to the bacterial photosynthetic reaction centre.^[242] In parallel work, simpler compounds such as the tetracyanoporphyrin in Fig. 7 were designed and synthesized for use in organic solar cells, leading to discovery of the dominant mechanism of primary charge separation in poly[2-methoxy-5-(2′-ethylhexyloxy)-1,4-phenylenevinylene] (MEH-PPV) – [6,6]-phenylC₆₁-butyric acid methyl ester (PCBM) blends,^[243] the demonstration of optically controlled switching based on oligoporphyrin oligomers like the tetramer shown in Fig. 7,^[244] the design of molecules like the dichlorotin porphyrin shown in Fig. 7 for reaction with and surface modification of the highly activated silica surfaces found in optical fibres,^[245] the discovery of adatom-bound attachment of nitrogen bases based on 1,10-phenanthroline for the stabilization of gold surfaces and nanoparticles,^[246] solar-energy upconverters 30 times more efficient than had previously been achieved^[247,248] based on the cruciform, T-shaped, and elbow π-expanded porphyrins^[249–252] (Fig. 7), and the development of molecules (see later) and methods for measuring the free energy of polymorphism of self-assembled monolayers on hydrophobic surfaces.^[253]

That we were able to predict and design functional system properties for Crossley’s molecules^{[244,249,250,254–258]} provided an important aspect to their success. These predictions integrated the traditional electron-transfer theories of Hush with the parallel understandings of Craig, as described in the preceding section. Some systems have been designed and shown to be synthetically feasible but not yet developed, however. This includes a molecular n-bit shift register made by the controlled synthesis of large porphyrin assemblies on substrates containing leads fabricated using current silicon device technology.^[259] The established synthetic strategies involve scanning tunnelling microscopy (STM)-induced reactions of molecules^[260–270] on silicon surfaces performed with atomic precision,^[271,272] that are subsequently linkable to porphyrins.

Crossley’s molecules were also important in that they led to fundamental understanding of quantum-chemistry computational procedures now widely used in nanotechnology. The fundamental discoveries leading to the establishment of the CAM-B3LYP method programmed and championed by Kobayashi as part of the chlorophyll-assignment project (see earlier) were demonstrated using Crossley’s porphyrins.^[87,91,244]

The 1950s saw not only the advancement of Craig’s methods for spectroscopy and Hush’s methods for discrete electron transfer events but also the advent of the Landauer method^[273,274] for understanding steady-state conduction through single molecules by coherent quantum transport. The importance of this work was recognized by Mark Ratner^[275–277] (who had first proposed the use of single molecules as a replacement for silicon electronics in 1974),^[278] Supiyro Datta,^[279] and others. The approach was based on Keldysh Green’s function theory^[280] and involved the green’s function G(E) = (H – E1)^–1 where E is the energy of the transmitting electron, and H is the Hamiltonian of the molecule in the field of the semi-infinite external electrodes connecting to it, represented as the unperturbed molecular Hamiltonian plus self-energies Γ^L and Γ^R describing the external contacts to the left (L) and right (R) electrodes, see Fig. 8. Professor Hush and myself showed how this approach could be applied using a priori electronic structure calculations.^[281] Through developments by Hong Guo, DFT implementations of this method combined with non-equilibrium Green’s function techniques^[282] can now model the 14 nm technologies used in current silicon-chip manufacture more accurately than traditional empirical transistor models.^[283] This field is enjoying a renaissance with, for example the recent discovery of a molecular rectifier switchable at 17 GHz.^[284]

Fig. 8.  Left: idealized structure of an isolated benzenedithiol molecule chemisorbed between two Au(111) electrode surfaces (blue – Au, yellow – S, green – C, white – H); and right: conceptualization of the way molecular vibrations scatter transmitting electrons in terms of symmetry properties of the incoming (L) channel, outgoing (R) channel, and vibrational wave functions. The junction channels are described in terms of the molecular conductance point group which is the molecular point group without end-to-end symmetry, while the vibrations are described in terms of the full molecular point group.

While calculation methods for molecular electronics are now well established, this was not always the case and indeed the whole field was once under a cloud in that it had never been shown that single molecules carrying currents of the order of μA remained chemically intact throughout the process. This situation changed in 2004 when two independent measurements^[285,286] reported the vibrational spectrum of a conducting molecule, determined from the second-derivative of the current as a function of voltage in a process discovered at the Ford Motor Co. known as inelastic electron tunnelling spectroscopy (IETS).^[287] Both spectra depicted vibrational levels at about the expected energies for the (very similar) molecules investigated, but the observed intensity patterns were completely different. A new method to understand the intensities in this type of spectra was required.

Our student Gemma Solomon quickly showed that the qualitatively different intensity patterns originated from different binding sites of the conducting molecule on the attached electrodes,^[288] and Mark Ratner showed that the IETS experiment was more likely to emphasize totally symmetric molecular vibrations.^[289,290] Such a discussion of symmetry was unprecedented in the field of molecular electronics, though it was one of the great successes of molecular spectroscopy decades earlier and a central feature of Craig’s works.^{[14,15,33–43]} A reason for this is that it was reasoned that in general the junctions between arbitrarily positioned molecules and electrodes would not show symmetry. However, we demonstrated that the symmetry of the junction has only a minor effect^[291] and can largely be ignored, allowing the full power of group theory to be applied to understand coherent elastic and inelastic electron transport through molecules.^[292,293] Solomon showed that the appropriate point group was not that of the molecule but instead the molecular point group without any end-to-end symmetry as the applied voltage eliminates such symmetry elements. As a result, the symmetry of the electron-transport channels through the junctions can be defined, and the symmetry of the molecular vibrations utilized.

An example of this is shown in Fig. 8 for chemisorbed benzenedithiol spanning two gold electrodes. Current through the junctions and into the molecule is provided mostly by π-molecular orbitals of B₁ symmetry, with coefficients on different atoms weighted as shown in the figure; however, σ channels of A₁ symmetry also provide connection. Typically single channels dominate the junctions so if the same connector material is used on each side then the net product of the two junction symmetries must be totally symmetric and hence symmetric vibrations typically scatter the most.^[289,290] However, knowing the ratios of the junction conductances of secondary channels to the dominant channel, scattering off asymmetric vibrations can occur and be quite intense, as illustrated in the figure.^[292,293] The sensitivity of each particular vibration to scatter the transmitting electron is given by the vibronic coupling constants used by Craig to explain the spectra of benzene and related molecules, the same quantities that we showed also pertained to Hush’s electron-transfer theories. There is thus a very close connection between the understanding of molecular spectroscopy and the understanding of IETS phenomena.

As previously described, central aspects of the works of Craig and of Hush involved the breakdown of the Born–Oppenheimer approximation. This was embodied in the work of Craig through vibronic coupling constants linking delocalized diabatic states or, equivalently, through electronic couplings linking localized diabatic states in the work of Hush (Fig. 3). A critical resulting feature is that the resultant electron-vibration system Eigenfunctions cannot be written as a simple product of electronic and vibrational wave functions. As a result, the electronic and vibrational motions become entangled.^[16] Entanglement is the defining property of quantum systems, allowing them to store more information than could be encoded into an analogous classical system. This feature forms the basis of proposals for quantum information processing and hence chemical qubits based on electron-vibration entanglement that have been proposed for possible technological applications.^[294] A qubit is a quantum memory storage device analogous to bits in modern digital computers.

Our general model for chemical reactions developed in conjunction with Ross McKenzie and Laura McKemmish allows insight to be gained into what types of chemical processes may be useful as the basis of electron-vibration chemical qubits.^[135,184] By varying the parameters ħω/ΔE, 2|J|/λ, and E₀/ħω, it is possible to calculate the resulting entanglement and its stability with respect to external perturbations. The critical properties of a qubit are the amount of entanglement it generates, its ability to preserve that entanglement for long times during which quantum computations can occur through quantum coherence, and its ability to be read and written as classical states through intentional quantum decoherence. In essence, the requirements of quantum coherence during processes, which demands extreme isolation of the qubit from its environment, and easy read/write, which demands strong interactions with some part of the outside world, are seemingly mutually inconsistent, defining a central challenge in this research field. Entanglement can be measured as the von Neumann entropy S of the Eigenfunctions of the qubit, which varies from 0 for a classical system (in this case, a pure diabatic electronic state) and 1 when there is formed a phased 50–50 mixture of two such classical states.^[135,184] At the simplest level of approximation, fluctuations in the environment modulate the energy difference E₀ between the two diabatic states, making the sensitivity of the entanglement to variations in E₀/ħω a quantity of primary interest.

Fig. 9 shows the entanglement as a function of ħω/ΔE and 2|J|/λ for three values of E₀/ħω = 0, 0.01, and 0.1. To help understand these results, the figure also shows the bimodality index of the systems. Bimodality occurs when the ground-state Eigenfunction has a local maximum somewhere in its vibrational density profile. This can only happen when the Born–Oppenheimer adiabatic ground-state potential-energy surface has a double minimum that is sufficiently deep to support zero-point motion in each well. While entanglement occurs to some extent for all values of the parameters, it becomes large only when the ground-state density is bimodal, requiring at least 2|J|/λ < 1. However, as Fig. 9 shows, if the wells become very deep (2|J|/λ << 1) then small levels of asymmetry destroy the quantum interference, making the qubit extremely sensitive to environmental effects. A reasonable level of insensitivity is found in only a small region of the parameter space with 2|J|/λ only slightly less than 1 and ħω/ΔE only slightly less than 1. However, the single vibrational mode model that is used becomes invalid as ħω/ΔE → 1 because then Jahn–Teller effects can dominate. Hence finding a suitable chemical process for application as a quantum qubit is a very difficult challenge.^[135,184]

Fig. 9.  Quantum entanglement (von Neumann entropy) S, colour coded, as a function of the parameter space of general chemical reactions (lower row), emphasizing its relationship to the bimodality of the ground vibronic state vibrational density profile (black is bimodal, white is unimodal, top row) and its sensitivity to environmental effects manifested as small changes in asymmetry E₀. The numbers indicate the location on the parameter space of the 10 model compounds (Fig. 4, Table 1).

Quantum entanglement may be of practical use in understanding quantum systems irrespective of its possible application to qubits in information processing. Relevant to information processing, the entanglement was described above in terms of diabatic states mimicking classical states to which experimental measurement could collapse quantum wave functions. No known experiment mimics collapse of non-adiabatic wave functions to Born–Oppenheimer states, but nevertheless the associated entanglement reveals significant insight into the Born–Oppenheimer approximation, indicating novel ways in which this could become manifest. Such insight is different to the commonly considered manifestations of breakdown in terms of say binding strengths, equilibrium geometries, transition energies, reaction probabilities, etc. What is found^[184] is that the entanglement scales simply with the cusp diameter Q_c and so provides a more intrinsic measure of the breakdown process than is provided by the other measures as these all amplify some particular aspect of the phenomenon. What becomes clear is that ways of calculating Born–Oppenheimer breakdown appropriate to one question of interest may not be appropriate for other properties of the same system. Indeed, the entanglement is the only property we have investigated that is dominated by the FD correction term, the DC term not contributing at all to entanglement with the SD term clearly being of secondary importance.

A dramatic possible manifestation of chemical quantum entanglement appears in the original ‘Orch-OR’ proposal of Penrose and Hameroff^[295–299] that qubits based on this principle form in the tubulin dimers that constitute the microtubular primary structural supports in most cells. Their proposal is that such qubits can remain coherent on the neural timescale, being decohered only by quantum gravitational effects in some way which allows them to modulate neural cellular potentials and hence provide the essential element required for conscious thought. The search for such quantum processing contributions to neural activity stems from Penrose’s Lemma stating that electrochemical processes in neurons, as well as say the logic processes achievable in modern (classical) digital computers, are incapable of manifesting consciousness.^[300,301]

All chemical quantum computing proposals involve coupling of electronic motion to vibrational motions that occur on the fs timescale. These, and all processes occurring on slower timescales, would have to remain quantum coherent on the ms neurological timescale for this process to be feasible in the Penrose–Hameroff ‘Orch-OR’ qubit proposal, something that would be extremely amazing, akin to superconductivity and Bose–Einstein condensation. Penrose and Hameroff evoked the possibility that an unobserved yet widely speculated phenomenon known as Fröhich condensation^[302–304] could provide the required coherence. Fröhlich had considered a non-equilibrium steady-state scenario appropriate to biological environments in which energy flow from a source to a sink through an intervening oscillator system caused a dramatic non-equilibrium energy distribution amongst the oscillators that he likened to a Bose–Einstein condensate. However, we showed that Fröhlich’s equations indicate that in biological systems this effect only occurs after the energy of the non-equilibrium system, expressed as an average ‘temperature’, entered the range of 2000–1 000 000 K. We also showed that the effect was purely classical in nature and hence did not purvey quantum coherence.^[182] Finally, we showed that the coherent vibrational oscillation specified in their proposal was in fact not a vibration at all but instead involved irreversible chemical reactions and the dismantling of the microtubules.^[181] The original proposal was therefore untenable on a variety of biochemical, chemical, and physical principles. It has been withdrawn and the Penrose–Hameroff proposal restated^[299] without specification of its central element, the vibrational qubit.^[183] As a result, we concluded that there exists ‘no basis on which Orch-OR theory can be considered as a proposal worthy of further consideration’.^[183]

On a more positive note, this work on the basics of chemical quantum computing led also to the development of qualitative tools as well as a simple computational model for predicting which gates can be successfully operated in an NMR-based quantum information processor.^[190] This allows synthetic chemists to design new molecules with improved functionality, avoiding the huge supercomputing efforts demanded by the complicated quantum dynamics models used in current qubit simulations. It is based on the observation that the critical qualitative features of the complex experiments performed using NMR pulses shaped in the time domain are provided very easily by transforming the problem into the analogous frequency-domain representation. Some major results are yet to be published.

Craig’s heroic first calculation of electron correlation in 1950^[13] opened the way for many modern research fields, from its dominance of the spectroscopy of aromatic molecules studied by Craig to subtle modifications to classical chemical bonding that facilitate many chemical processes to superconductivity to the London dispersion force that liquefies helium, causes oils to boil at high temperatures, controls biological and polymer structure, and allows sticking to ‘non-stick’ surfaces. While the qualitative effects of electron correlation have been recognized since the inception of quantum mechanics,^[305] how to treat it efficiently in a priori or ab initio calculations has been a dominant theme of research over the past decades. Today there is much excitement concerning the possibility of applying 2D materials in electronics, nanotechnology, and biotechnology applications, with dozens of materials now known and thousands conceived.^[306] These stack together by van der Waals forces, and understanding the ‘van der Waals gap’ that controls their properties is critical.^[307]

Electron correlation is essentially a classical effect rather than a quantum one as it relates to how the instantaneous position of one electron affects the motion of others, although quantum effects do influence its manifestation. To put issues in perspective, Bohr’s semiclassical theory of electronic structure may be applied to molecules and is known to deliver more accurate results than comes from many modern sophisticated DFT approaches.^[308–312] The advantage of Bohr’s approach is that it includes electron correlation simply and easily. Indeed, we have shown that a range of basic chemical phenomena like the nature of the valence states of transition-metal compounds, f-block chemistry, the basic description of NMR spectroscopy, relativistic effects in chemistry, etc., can be simply described at this level, providing an easy way to explain chemistry to freshman students.^[313]

However, the dispersion force is an aspect of electron correlation that is purely quantum in origin. Classical mechanics specifies how permanent electrical charge distributions interact with each other and polarize their surroundings, but quantum mechanics adds something extra: the ability of systems to undergo spontaneous short-lived charge fluctuations as demanded by the uncertainty principle. After these fluctuations occur they interact with the environment in the same way as do permanent classical charges, meaning that the interactions take on a very different form if they occur in metals or in insulators, for example. Such classical aspects of a quantum phenomenon are mostly neglected in standard chemical thinking, despite being critical to understanding.

The reason for this of course is that molecules and proteins are not conductors and so the application of ideas coming from the classical electrostatics of insulators is almost always adequate. This leads to the idea that quantum fluctuations are restricted to occur on atoms, and that these fluctuations interact independently with the atoms in all surrounding pieces of matter. The net dispersion force between two molecules is then written simply as the sum over all inter-atomic interactions of terms scaling as r^–6, where r is the interatomic distance. This generates London’s basic description of the dispersion force.

Consider, however, the situation shown in Fig. 10 in which particles are surrounded by a metallic conductor, and that this system then interacts with an external field. Gauss’ Law of classical electrostatics indicates that the metal conductor acts as a Faraday cage which prevents the electric field from penetrating inside, so the internal particles are screened from the external one and no interaction occurs. Gauss’ Law appears as a consequence of strong electron correlation in classical electrostatics. Another directly relevant feature of classical electron correlation is co-operativity in which a system of interacting polarizable centres interfere with themselves in a constructive manor, leading to much larger dipole moments than would otherwise have been expected. The increase in dipole moment of a water molecule from 1.8 Debye in the gas phase to 2.9 Debye in liquid water is an example of such a cooperative effect.

Fig. 10.  Screening of charges by metallic conductors occurs independent of whether the electric field arises from permanent classical charges, dipole, or electromagnetic radiation, or else from quantum dipole fluctuations as part of the London force. Left: classical fields are screen by metallic ‘Faraday cages’; Right: the van der Waals force is screened between an AFM tip and a substrate by an intervening 2D conductor such as graphene.^[315]

The same principles apply to the dispersion force.^[314] Quantum fluctuations on an external molecule are screened by the Faraday cage and are not perceived by molecules inside the cage. Indeed, this effect has been recently demonstrated by Tsoi et al. for a system in which a silica surface is approached by an atomic force microscope with a silica tip.^[315] Normally, a large dispersion force describable in terms of a sum over atoms as r^–6 London energy contributions is found, and this force is maintained if a single layer of insulating fluorographene is inserted between the tip and surface. However, if conducting graphene is inserted instead, the tip and surface no longer experience van der Waals attraction! Instead, the field lines coming from a spontaneous dipole fluctuation inside the tip are reflected back off the graphene, generating interference patterns and possibly enhanced ‘antiscreening’ interactions with solvent or other matter nearby the tip.^[314]

Expressed in a broad context of modern interests in biotechnology, surface chemistry, and nanotechnology, these issues have been clearly described recently by John Dobson from Griffith University.^[314] His conclusions are presented in Table 2 concerning the dispersion interactions between objects of different shape in the asymptotic region of large separation as a function of the separation distance D. There, results from simple approaches to understanding inter-object dispersion as pairwise additive sums over interatomic London terms scaling as r^–6, where r is the interatomic distance, are compared with the form determined using the sophisticated random-phase approximation (RPA).^{[314,316–322]} Simple London-based approaches work well for intermolecular interactions, explaining their success in chemistry and biochemistry, and these also work for interactions between all 3D materials, but they fail for 1D and 2D conductors such as metallic nanotubes and graphene.

Table 2.  Asymptotic form of the dispersion interaction determined by Dobson between two parallel surfaces separated by distance D as determined from sums of atomic London terms and from RPA calculations^{[314,316–322]}

The screening interactions depicted in Fig. 10 are known as Dobson-B corrections^[314] to pairwise-additive dispersion relationships. In addition, Dobson-C corrections involve quantum fluctuations that are not restricted to atoms but instead occur on length scales limited only by the finite duration of the fluctuations and the effective Fermi velocity of the metal. (In addition, the finite speed of light should also be considered, as in the Casimar effect,^[323] but these influences are usually small.^{[316,318,324]}) Hence quantum dipole fluctuations can spread over distances that are large compared with separations to neighbouring objects, possible up to lengths of 100 nm in order. As sketched in Fig. 11 for interactions between 1D conductors like nanotubes, this converts the usual quantum-fluctuation-dipole to induced dipole interaction described by the London force into something more like interacting quantum-fluctuation-charge to induced charge interactions. Such effects lead to the interactions between conducting nanotubes being very much stronger than expected, changing the asymptotic dependence from D^–5 to D^–2ln(D/D₀)^–3/2 in Table 2. In parallel 1D conductors, such fluctuations are perfectly aligned for maximum effect. Averaging over two dimensions reduces their impact, with fluctuations in the correct direction to have an effect becoming too improbable in three dimensions for there to be any major influence.

Fig. 11.  Top: London forces in insulating materials arise through quantum dipole fluctuations on atoms that induce polarization in neighbouring molecules; Bottom: the analogous situation involving long-range quantum fluctuations in a 1D metal inducing polarization in a neighbouring 1D conductor.

Computational methods that include dispersion by first-principle approaches are computationally demanding. This is highlighted in Table 3 where the scaling of various computational methods as a function of the number of electrons n is shown, as well as the scalings required to qualitatively describe different types of dispersion interactions. If only insulators are involved then n⁵ scaling is adequate whereas n⁶ scaling is required in general.^[314] Ab initio methods like second-order Møller–Plessett perturbation theory (MP2) scale as n⁵ and hence are useful in chemical and biochemical applications but more expensive approaches like CI singles and doubles (CISD), coupled-cluster singles and doubles (CCSD), RPA, or TDDFT are required in general. However, these may not necessarily be applicable to conducting systems as they are single-reference in nature and so their multi-reference analogues are typically required. Traditional DFT methods using the local density approximation (LDA), the generalized-gradient approximation (GGA), or hybrid functionals fail completely, although scaled-down MP2-like dispersion is included in double-hybrid functionals. Note, however, that the scalings shown in the table are formal ones, not ones that may be achieved using some approximate yet very accurate implementation technique. Most significantly, a new TDDFT procedure has just been introduced with ~n⁴ scaling instead of n^[6,325] and this method may prove very useful in the future.

Table 3.  Computational scaling requirements with the number of electrons n requirements for qualitative description of exchange and van der Waals correlation effects in various systems, as well as the treatments afforded in various computational approaches

Notwithstanding this, we have recently reviewed^[326] many available options for accurate calculations, including multi-reference coupled-cluster,^[327–337] embedded approaches,^{[332–334,337–343]} mixed DFT–ab initio approaches,^{[324,337,340,344–349]} and quantum Monte Carlo.^[350,351] For most problems such approaches remain computationally unfeasible, however, and so empirical dispersion correction schemes to low-level GGA or hybrid DFT provide the most commonly used strategies. Of the currently available schemes, we recommend two: the D3 scheme of Grimme et al.^{[266,352–354]} and the many-body dispersion (MBD) scheme of Tkatchenko et al.^{[267,355–360]} The D3 scheme is based on pairwise sums of London terms and so suffers from the problems described in Table 2 for long-distance interactions, although a variant containing three-body interactions, D3(ABC), is available^[352] and underused.^[326] It also has older variants known as D and D2^[266] but is most commonly applied in its form utilising Becke-Johnson damping, D3(BJ),^[354] which herein is simply called just D3. The MBD scheme includes all many-body Dobson-B effects but does not allow for long-range Dobson-C fluctuations and hence also produces the wrong limits for long-range interactions. This scheme scales nominally as n³ in the number of atoms but the prefactor is large and the time required for MBD calculations can exceed that for the base GGA density functional by orders of magnitude. Hence this method can only be applied to small-to-medium-sized systems,^[326] but the very efficient D3 scheme can be applied to any problem treatable by DFT. A promising alternate scheme to D3 and MBD is the exchange-hole dipole moment method (XDM)^[361] that has recently undergone significant improvement.^[362]

Both D3 and MBD have been shown to be very effective for many different types of problems, from ice cluster energies through crystal cohesion energies to basic chemistry and to surface adhesion to interactions with graphite and nanotubes.^{[260–270,341]} This occurs despite fundamental limitations in the description of the interactions (Table 2). Success comes from the fact that these methods are designed to model interaction energies and geometries and are usually only applied for these purposes.

Our particular focus within this topic is self-assembly processes from solution. Such processes are important to, e.g. biology, polymerization, and nanotechnology. Two features of critical relevance are that processes are controlled by free energies not just binding energies, and that solvation effects act to dissolve self-assembled monolayers (SAMs). Our work has focussed on the self-assembly of tetraalkylporphyrin molecules synthesized by Max Crossley on highly ordered pyrolytic graphite (HOPG) surfaces from solution in 1-phenyloctane, octanoic acid, or n-tetradecane.^{[253,363–368]} This involved measurement of STM images and the simulation^{[253,363,369]} of the Gibbs free-energy of formation for a wide range of possible atomic structures for each image. The molecules considered are shown in Chart 2.

Chart 2.  Tetraalkylporphyrins M-C_mP of Crossley used to study self-assembly and polymorphism on HOPG from organic solvents.

An important feature of this system is the observation of various polymorphs for shorter-chained molecules. Key aspects are shown in Fig. 12 and include grain boundary, a rare but informative feature of a polymorphic SAM. This boundary appears between different polymorphs and shows one row of molecular structures in the low density ‘L’ polymorph plus several interface ‘I’ rows that in their middle take on the structure of the medium-density ‘M’ polymorph.^[368,369] The molecular structures determined from modelling of large areas of the L and M polymorphs are shown, as is how they seamlessly meld together to form the interface layers. The L polymorph is observed for all odd-length chain molecules studied (lengths m from 11 to 19)^{[253,363,366,368]} but not for the even chain molecule studied, 2H-C₁₂P,^[364,365] while three M polymorphs^{[253,364,365,368]} of quite different form named M_a, M_b, and M_c are observed for 2H-C₁₁P, 2H-C₁₂P, and 2H-C₁₃P, respectively. The L structures have all four alkyl chains lying flat on the surface while the M polymorphs have two such chains with the other two somehow oriented into the solution. As there are many ways that molecules can orient in solution but only a few ways in which they can kink to lie flat on surface, a much greater variety of possible forms exists for the M polymorph compared with the L one. A high density polymorph with only one chain on the surface is also known for 2H-C₁₁P.^[368]

Fig. 12.  Top: observed STM image^[368] of a region of the SAM of Cu-C₁₁P throughout which isolated lines of molecules assemble as in the L polymorph (all four chains on surface) adjacent to in lines of interface molecules (three chains on surface) that connect with each other using the assembly pattern of the M polymorph (two chains on surface), with overlaid molecules taken from the calculated structure of L (middle) and that as merged with M (outer). Bottom: the construction of an approximate I structure by overlaying fully optimized L and M structures, selecting the required unique chain conformations, and averaging the remaining coordinates. Reproduced from ref. [369] with permission from Taylor and Francis Ltd.

Atomic structures for the SAMs are determined by comparing calculated and observed structures. The calculated structures were obtained by optimizing five possible structures for the L polymorphs^[363,366] and 40 possibilities for the M polymorphs.^[363,369] These possibilities were obtained by considering various conformations of the chain kinks that occur as the alkyl chains come away from the porphyrin at an angle near 60° and must subsequently be bent back in order to lie flat on the surface. The 40 possibilities for the M polymorphs were obtained by varying kink-structure conformers and macrocycle rotation angles, first eliminating possibilities that cannot be shaped into a valid SAM. After optimization, only at most 22 stable SAMs were realised, named Ma to Mv. Structure assignment was made by comparing these structures to experiment, considering:^{[253,363,366,369]}

• Agreement of observed and calculated lattice parameters. To do this, new experimental techniques based on internal calibration and repeated measurements were developed to significantly reduce the experimental uncertainties in the lattice parameters, allowing for discrimination between the calculated structures.^[253]

• Agreement between observed and calculated chain hydrogen atom positions. It is common that images show resolved hydrogen atom locations for the alkyl chains, providing detailed information concerning the structure.

• Agreement of the shape of the porphyrin-ring images. These images always have structure, typically eight resolved units coming from the pyrrole rings and the kinks in the alkyl chains, features that can be related to the calculated structures.

Always only one calculated structure matched all of the observed data, leading to a clear assignment. The observation of interface regions provides strong verification of the assigned structures as the L and M polymorphs must smoothly mesh into each other, a feature which is very difficult to obtain as it requires perfect alignment of many aspects of the structures.^[369]

The calculations also yield the free energy of SAM formation.^[253] This is evaluated using DFT in various forms. In one way, the B3LYP functional is used to describe the porphyrin and the kinks, with molecular mechanics force fields used to describe the porphyrin–porphyrin interactions and a specialized empirical graphite-alkane force field fitted to experimental enthalpy data used for the surface attraction, a QM/MM approach.^[366] After optimization the vibration frequency of the SAM is determined, leading to evaluation of the zero-point energy and entropy of the SAM.^[253] Solvation is then included using a self-consistent reaction-field approach that treats the dispersion interaction between the surface, molecule, and solvent as being simply proportional to the surface area exposed to the solvent.^[253] The other approaches are based on using the Perdew–Burke–Ernzerhof (PBE)^[370] density functional augmented by empirical dispersion corrections to describe all interactions.^[253] The D3 and MBD empirical dispersion correction methods were both used.^[326] Only optimized geometries and interaction energies were evaluated in this way, with the zero-point energy and entropy corrections determined by the QM/MM method simply added to these results.

A new method was also developed for estimating free energies of formation from experimental STM data.^[253] Only a few measurements of free energies of SAM formation have ever been performed,^[371] and while the error bars of our technique are large (of order a few kcal mol^–1), they provide means for analysis of the calculated free energies. Comparison of observed and calculated values is made in Table 4. The observed values are about twice as large as the calculated ones, but the absolute differences in free energy of ~7 kcal mol^–1 are small compared with the likely errors involved. The various computational methods are mostly in good agreement with each other.

Table 4.  Observed and calculated Gibbs free energies of SAM formation, in kcal mol^–1, for L and M polymorphs of H2-C_mP SAMs^[253,326]

The calculations provide insight into the factors controlling SAM polymorphism, as indicated in Fig. 13 for 2H-C₁₃P. The strong (PBE-D3) binding to the surface of –148 and –100 kcal mol^–1 for the L and M polymorphs is opposed by entropy effects of 37 and 25 kcal mol^–1, respectively, combined with desolvation effects of 106 and 70 kcal mol^–1, respectively, yielding free energies of formation of just –8 and –6 kcal mol^–1, respectively. SAM formation is therefore associated with the near-complete cancellation of large SAM forming and SAM dissolving contributions, and subtle variations in this balance controls polymorphism.

Fig. 13.  Components of the Gibbs free energy (a) of SAM formation of the L and M polymorphs of 2H-C₁₃P. The net solvation (b) and entropy-dominated thermal energies are the differences between energies for the solvated SAM and the solvated free molecule and free surface, while the net binding energy for the BPE-D3 method (red) is partitioned into its PBE and D3 components (c); the QM/MM binding energy is shown in blue.

As the L polymorph is observed to form at similar porphyrin concentrations over a wide range of chain lengths m, its free-energy of formation must change only slowly with m. In the asymptotic region of large m (> 15), the calculated chain-length dependences d(ΔG)/dm are shown in Table 5 evaluated using the QM/MM and PBE-D3 methods.^[363] These dependences are also partitioned into the free-energy contributing components. As for the total free energies, these chain-length dependences show large cancelations between the contributing desolvation, entropy, and binding terms. Comparing the QM/MM and PBE-D3 methods, the QM/MM attraction term comes purely from the empirically fitted HOPG-alkane force field,^[372] the AMBER inter-alkane force field,^[373] the B3LYP description^[374] of the intramolecular kinking energy, and the D3 dispersion contribution.^[375,376] The empirical, AMBER, and B3LYP contributions would all be expected to be accurately determined. Hence the similar QM/MM and PBE-D3 results indicate that the D3 correction is also very reliable. The PBE-MBD method also gives similar results.^[326] However, that these approaches predict very small total chain length dependences indicates that the explicit treatments used for the adsorbate–surface interactions compliment the dielectric-continuum calculations used to calculate the solvation of the molecule, bare surface, and SAM. These solvation calculations were based on the model of Floris et al.^[377] and hence indicate that this implicit solvation formalism is highly compatible with D3 (and MBD).

Table 5.  Large m asymptotic contributions to the chain-length dependence d(ΔG)/dm of the free energy of formation of the L polymorphs (in kcal mol^–1)^[363]

In the broader context of Dobson’s description of dispersion forces in nanotechnology^[314] described in Table 2, we note that this system involves large molecules with a delocalized central π system with attached alkane chains interacting with a metallic substrate. The molecules are essentially non-conductors, although when charged the central ring becomes internally conducting. According to Dobson’s analysis, methods like D3 based on summing interatomic London forces should be able to provide realistic descriptions of the chemistry, and the agreement found between D3, the empirical force field used that embodies screening, and the more general MBD approach supports this analysis. There are many applications in self-assembly that provide other examples of systems treatable by methods like D3. We have recently reviewed other applications elsewhere.^[326]

While the dispersion force was always known to be an example of electron correlation with profound consequences for molecular self-assembly, Craig’s first calculation of electron correlation has led in recent times to the discovery of new phenomena. The dispersion force between atoms scales somewhat in proportion to the product of the number of electrons on each atom, making it very important for heavy atoms. While such ideas have always been recognized in Pearson’s concept of hard/soft acids and bases,^[378–381] that the dispersion force can be large enough to compete with covalent and ionic binding is just coming to be understood.^[382,383] Sulfur-bound compounds SR stabilize gold surfaces and nanoparticles, and the chemical form has always been assumed to be Au^I-thiolate, and many hemes are also stabilized in proteins by SR groups (like cysteines) that are assumed to be Fe^III-thiolates, but ionic/covalent bonding is not necessary and Au⁰-thiyl or Fe^II-thiyl species are also possible. These are different chemical species, with different reactivities and spectroscopies. It has recently been shown that cytochrome Fe^II-thiyls can be formed,^[383] and that Au⁰-thiyls are the species that actually protect gold surfaces and nanoparticles while Au^I-thiolate production etches surfaces and inhibits nanoparticle formation.^[382] The competition between covalent/ionic and van der Waals bonding is critical in each case, and Pearson’s ideas predict the main results as S and Au are soft while Fe can be either hard or soft. The Au–Au bond has always been regarded as anomalous, with the concept of ‘aurophilicity’ used to describe it,^[384,385] but this is now known to be just the strong dispersion force acting on top of strong s–s interactions that make the gold atoms in a gold dimer more like those in gold metal than like isolated gold atoms.^[386,387]

Understanding whether or not processes occur on the Au⁰-thiyl or Au^I-thiolate potential-energy surfaces is essential to rationalizing the known synthetic chemistry and spectroscopy of gold surfaces and nanoparticles.^[382] Work with Noel Hush (Sydney University) and Jens Ulstrup and Jingdong Zhang (Danish Technical University) has shown that gold (111) surfaces react with thiol or disulfide reagents to form initially ‘physisorbed’ species and then ‘chemisorbed’ ones (Scheme 1):

Scheme 1. 

The SAM may take on two different types of atomic structure, either bound to gold adatoms above the surface or else bound directly to face centred cubic (FCC) surface sites, with the nature of the R group and the coverage determining which form.^[388] For example, from among the butanethiol reagents, n-butanethiol,^[389] iso-butanethiol,^[390] and enantiomerically pure 2-butanethiol^[391,392] bind to adatoms, but t-butanethiol binds to FCC sites^[393,394] owing to strong steric repulsions in the adatom form and racemic 2-butanethiol shows both motifs in the same SAM^[391,392] so as to maximize coverage; cysteamine in aqueous solution also binds to FCC sites owing to repulsion between the charged tail groups.^[395]

Many properties of the two head-group arrangements are very similar despite the significant change in coordination, suggesting that the binding mechanism is similar in both cases. Also the binding strengths of the chemisorbed forms are not that much greater than those of the physisorbed forms^[382,396] and are consistent with the binding of nitrogen and phosphorous bases to the surface, processes known to be controlled by the dispersion interaction. The postulate that both the ‘physisorption’ and ‘chemisorption’ processes depicted in the above reactions involve primarily dispersion forces and the Au⁰-thiyl valence state provides a consistent analysis of all of these data, including that for other chemical variations such as the replacement of Au with Ag or Cu and the replacement of S with O, Se, or Te.^[382]

However, gold does not react with solutions containing thiolate anions RS^–, but in an electrochemical cell in which Au^I is produced at the anode and RS^– at the cathode, the two species combine in solution to form Au^I-thiolate compounds and thin films that then precipitate onto the surface (Scheme 2).^[397]

Scheme 2. 

Au^I and thiolate ions are present as intermediates and their reaction produces an end product completely different to the protected SAM-covered surface produced by reactions with thiols and disulfides. Instead, this process parallels many known chemical reactions in solution that produce Au^I-thiolate compounds and films and it is reasonable to assume that this is the type of chemistry that occurs whenever the Au^I-thiolate potential-energy surface is engaged.

Sulfur protected gold nanoparticles can be formed using a range of synthetic conditions. We focus on what occurs under mild conditions in aqueous solution, and what occurs in the presence of excess strong reducing agents in non-aqueous media. In aqueous solution,^[398] thiols reduce Au^III halide (X) salts to nanoparticles and cover them according to Scheme 3:

Scheme 3. 

Thiolate anions are not reactants under these conditions, and, if the reaction follows known properties of the two-phase alternate Brust–Schiffrin synthesis^[399] then the Au^I intermediate [AuX₂]^– is implicated.^[400–402] Of concern is how these results can be consistent with the details of the one-phase Brust–Schiffrin synthesis^[403] in a 5 : 1 methanol/water solvent containing a large excess of the strong reducing agent NaBH₄. Au^I-thiolate molecular compounds can be isolated as intermediates that must be subsequently reduced by borohydride to form nanoparticles.^[404] The only interpretation of these data consistent with the known Au(111) surface chemistry is that the one-phase Brust–Schiffrin reaction dominantly proceeds in three steps (Scheme 4):

Scheme 4. 

This mechanism was verified by repeating the synthesis except for the omission of the Au^III salt.^[382] The result was the conversion of the thiol used, hexanethiol, into hexanethiolate, showing that thiolate anions are indeed the active species in the reaction. Suggestions made during the Craig Lecture Tour by Sebastian Marcuccio from La Trobe University and Advanced Molecular Technologies that the produced thiolates complex the boron^[405] as part of a thiolate storage mechanism,^[406,407] even in solution in the reactive methanol solvent,^[408] are currently being investigated. Under the mild conditions used in aqueous solution, radical pathways involving Au⁰-thiyl could dominate, whereas under conditions in which thiolate is produced, nanoparticles are not initially formed as the Au^I-thiolate pathways lead to molecular compounds and films rather than stabilized nanoparticles.

In terms of bonding diagrams, the Au^I-thiolate species involve covalent bonds formed between the s-electrons of isolated gold atoms and sulfur orbitals (Chart 3), whereas the Au⁰-thiyl species, prominent in compounds containing Au–Au bonds, involve strong interactions between filled Au d-orbitals and the sulfur (Chart 4) that affects bonding primarily through a strong Dobson-A effect, enhancing the dispersion interaction. The hybridization of the gold thus critically controls the bonding. When Au–Au bonds form, the resulting s–s interaction produces bonding and antibonding orbital combinations that are well removed from the Fermi energy and hence inaccessible to attacking reagents, making gold surfaces noble. In gold metal, the band broadening is sufficient to spread the occupied (‘valence’) and virtual (‘conduction’) bands so much that they weakly overlap, resulting in a low density of states at the Fermi energy that is responsible for the colour and conductivity of the solid.^[385,409] For understanding chemical reactivity, the low density of states at the Fermi energy is the most critical property, however.

Chart 3. 

Chart 4. 

Sulfur bound to FCC sites clearly attaches to gold atoms having many Au–Au bonds, but the adatom-bound motif has traditionally^[410–413] been drawn showing the adatom bound to only S and not to Au (Chart 5) and called the ‘staple’ motif.^[414] In reality, there are also two direct Au–Au bonds (Chart 6) while superexchange through the S atoms makes for an effective four-bond system (Chart 7).^[382] Hence the gold adatoms are tightly connected to the nanoparticle core and show similar broad s-bands.^[382,415] While differences in the band structures of the adatoms do occur that are significant in terms of the total chemical structure, their similarity is enough for Au^I-thiolates to be generally inhibited, becoming high-energy excited states of the protecting layer.

Chart 5. 

Chart 6. 

Chart 7. 

All spectroscopic studies clearly reflect the Au⁰-thiyl nature of the protecting layers of gold nanoparticles and surfaces. A direct experiment is X-ray photoelectron spectroscopy (XPS) which indicates that all atoms in the bulk and on the surface of gold nanoparticles are Au⁰ and not Au^I.^[399,404] Also the surface dipole moments measured for SAMs of SR groups and ammine groups on Au(111) have similar sign and magnitude, oriented so that the surface appears slightly negative and the ligands slightly positive.^[416] For SR groups, this directionality arises because the only electropositive element in the interface is hydrogen and these atoms are located far from the surface, despite sulfur–gold bonds being polarized with slightly more bond electron density on sulfur than on gold.^[416] This result follows directly from the Au⁰-thiyl nature of the surface and is opposite to that predicted if the surface was Au^I-thiolate. Similarly, near-edge X-ray absorption fine structure (NEXAFS) measurements of SAMs^[417] show electronic transitions to orbitals that are vacant in the Au⁰-thiyl description of the surface but already occupied in the Au^I-thiolate description.^[382] All DFT calculations show electronic structures of Au⁰-thiyl form,^{[415,417–426]} with the sulfur orbitals mixing with gold d-orbitals as predicted only by the Au⁰-thiyl description.^{[415,417,421,427]} Often DFT publications report the structure as being Au^I-thiolate, however, reflecting not the nature of the orbitals generated but rather the most commonly used descriptions in the field at the time.

All chemical species involving Au–S bonds must show some degree of Au⁰-thiyl character and some degree of Au^I-thiolate character as these two limiting valence forms mix via resonance in order to produce the ground and excited states of any system. The Au⁰-thiyl character must always be present as the dispersion interaction cannot be switched off, and the Au^I-thiolate character must always be there as the gold s-band is continuous through the Fermi energy meaning that some covalent interactions with the s-orbitals must always occur. Also, in-equivalence of the gold and sulfur orbital energies means that some degree of ionic Au^I-thiolate character must always be present through polarization effects. A simple way to understand the bonding is through an empirical partitioning scheme based on Hückel theory that partitions the binding energy into dispersion, covalent, and ionic contributions using a fitting parameter η specifying the hybridization of the gold orbital binding to the sulfur as d^1-ηs^η.^[382] This yields values of η in the range of 0.75 to 1 for traditional Au^I-thiolate molecular compounds, near 0.5 for reactive species likely to be seen in the early stages of nanoparticle growth, and below 0.2 for protected gold nanoparticles and surfaces. The elimination of Au^I-thiolate character from the bonding is thus seen as a driving force for nanoparticle growth, this force resulting from the strength of the Au–Au bonds that form as the particles grow and become more noble in character.

Craig’s first calculation of electron correlation has led nearly 70 years later to the understanding that electron correlation does not just act to modify typical covalent and/or ionic bonds describable without it, nor just to provide the strong van der Waals force response for surface adhesion, protein structure, and other intermolecular interactions, but can also form bonds between atoms in its own right.

This article, describing the material presented in 26 lectures during 2016–17 as part of the David Craig Medal Lecture Tour, summarizes many of the important achievements of the later years of my life’s work, expressed in the context of advances made by David Craig, Noel Hush, Maxwell Crossley, Jens Ulstrup, and many others. It highlights that Craig’s Legacy remains alive and well in many current fields of Chemical Research.