Bond! Chemical Bond: Electronic Structure Methods at Work

Chapter 4 Bond! Chemical Bond: Electronic Structure Methods at Work Abstract This chapter plunges into applied quantum chemistry, with various examples, ranging from elementary notions, up to rather advanced tricks of know-how and non-routine procedures of control and analysis. In the first section, the first-principles power of the ab initio techniques is illustrated by a simple example of geometry optimization, starting from random atoms, ending with a structure close to the experimental data, within various computational settings (HF, MP2, CCSD, DFT with different functionals). Besides assessing the performances of the different methods, in mutual respects and facing the experiment, we emphasize the fact that the experimental data are affected themselves by limitations, which should be judged with critical caution. The ab initio outputs offer inner consistency of datasets, sometimes superior to the available experimental information, in areas affected by instrumental margins. In general, the calculations can retrieve the experimental data only with semi-quantitative or qualitative accuracy, but this is yet sufficient for meaningful insight in underlying mechanisms, guidelines to the interpretation of experiment, and even predictive prospection in the quest of properties design. The second section focuses on HF and DFT calculations on the water molecule example, revealing the relationship with ionization potentials, electronegativity, and chemical hardness (electrorigidity) and hinting at non-routine input controls, such as the fractional tuning of populations in DFT (with the ADF code) or orbital reordering trick in HF (with the GAMESS program). Keeping the H2O as play pool, the orbital shapes are discussed, first in the simple conjuncture of the Kohn–Sham outcome, followed by rather advanced technicalities in handling localized orbital bases, in a Valence Bond (VB) calculation, serving to extract a heuristic perspective on the hybridization scheme. In a third section, the H2 example forms the background for discussing the bond as spin-coupling phenomenology, constructing the Heisenberg-Dirac-van Vleck (HDvV) effective spin Hamiltonian. In continuation, other calculation procedures, such as Complete Active Space Self-Consistent Field (CASSCF) versus Broken-Symmetry (BS) approach, are illustrated, in a hands-on style, with specific input examples, interpreting the results in terms of the HDvV model parameters, mining for physical meaning in the depths of methodologies. © Springer International Publishing AG, part of Springer Nature 2018 M.V. Putz et al., Structural Chemistry, https://doi.org/10.1007/978-3-319-55875-2_4 291 292 4 Bond! Chemical Bond: Electronic Structure Methods at Work The final section presents the Valence Bond (VB) as a valuable paradigm, both as a calculation technique and as meaningful phenomenology. It is the right way to guide the calculations along the terms of customary chemical language, retrieving the directed bonds, hybrid orbitals, lone pairs, and Lewis structures, in standalone or resonating status. The VB calculations on the prototypic benzene example are put in clear relation with the larger frame of the CASSCF method, identifying the VB-type states in the full spectrum and equating them in an HDvV modeling. The exposition is closed with a tutorial showing nice graphic rules to write down a phenomenological VB modeling, in a given basis of resonance structures. The recall of VB concepts in the light of the modern computational scene carries both heuristic and methodological virtues, satisfying equally well the goals of didacticism or of exploratory research. A brief excursion is taken into the domain of molecular dynamics problems, emphasizing the virtues of the vibronic coupling paradigm (the account of mutual interaction of vibration modes of the nuclei with electron movement) in describing large classes of phenomena, from stereochemistry to reactivity. Particularly, the instability and metastability triggered in certain circumstances by the vibronic coupling determines phase transitions of technological interest, such as the information processing. The vibronic paradigm is a large frame including effects known as Jahn–Teller and pseudo Jahn–Teller type, determining distortion of molecules from formally higher possible symmetries. We show how the vibronic concepts can be adjusted to the actual computation methods, using the so-called Coupled Perturbed frames designed to perform derivatives of a self-consistent Hamiltonian, with respect to different parametric perturbations. The vibronic coupling can be regarded as interaction between spectral terms, e.g. ground state computed with a given method and excited states taken at the time dependent (TS) version of the chosen procedure. At the same time, the coupling can be equivalently and conveniently formulated as orbital promotions, proposing here the concept of vibronic orbitals, as tools of heuristic meaning and precise technical definition, in the course of a vibronic analysis. The vibronic perspective, performed on ab initio grounds, allows clear insight into hidden dynamic mechanisms. At the same time, the vibronic modeling can be qualitatively used to classify different phenomena, such as mixed valence. It can be proven also as a powerful model Hamiltonian strategy with the aim of accurate fitting of potential energy surfaces of different sorts, showing good interpolation and extrapolation features and a sound phenomenological meaning. Finally, within the symmetry breaking chemical field theory, the intriguing electronegativity and chemical hardness density functional dependencies are here reversely considered by means of the anharmonic chemical field potential, so inducing the manifested density of chemical bond in the correct ontological order: from the quantum field/operators to observable/measurable chemical field.


Keywords Computational chemistry ab initio methods Hartree–Fock (HF) Density functional theory (DFT) Complete active space self consistent field (CASSCF) Valence bond (VB) Configuration interaction (CI) Electron correlation Quantum chemistry codes Gaussian General atomic and molecular






4 Bond! Chemical Bond: Electronic Structure Methods at Work
293
electronic structure system (GAMESS) Amsterdam density functional (ADF) Input files Keyword controls Fractional occupation numbers Electronic density Electronegativity Chemical hardness (electrorigidity) Ionization potentials Spin coupling Heisenberg–Dirac–van Vleck (HDvV) spin hamiltonian Resonance structures Graphic rules for the VB phenomenological hamiltonian Vibronic coupling Jahn–Teller effect Pseudo Jahn–Teller effect Vibronic orbitals Coupled perturbed techniques Molecular dynamics Vibronic phenomenological models Mixed valence Potential energy surfaces Symmetry breaking Chemical field Electronic potential




4.1
















Molecular Structure by Computational Chemistry: A Brief Synopsis Up to now, we have presented the principles of basic electronic structure models: Hartree–Fock (HF) based on single determinant wave function, Density Functional Theory (DFT), Multi-Configuration Self-Consistent Field (MCSCF) with its most representative version—Complete Active Space Self-Consistent Field (CASSCF), and Valence Bond (VB), in some modern avatars. At the same time, not denying their importance, details on some other used techniques were not provided. Thus, we confine ourselves only to mention procedures devised as perturbation corrections to the self-consistent methods, such as the second-order Møller–Plesset (MP2), with respect to orbital energy from HF (sum of occupied MO eigenvalues), or perturbation varieties devised as post-CASSCF treatments (CASPT2, NEVPT2). Top quality in accuracy of computed molecular energy (with consequences in good prediction of thermochemical parameters), is assigned to the class of Coupled Cluster (CC) methods, which are currently available as improvements starting from an HF singledeterminant reference (e.g. CCSD stands for Coupled Cluster with single and double excitations with respect of HF Slater determinant). There are also prospects for using CC in multi-configuration frames. For this chapter, we will use the MP2 or CC methods without explaining their principles, just for provisional comparative purpose. In principle, the existing methods, applied in accordance to their technical prescriptions, can really reinvent the molecules from first principles, as states the etymology of ab initio (in free translation, meaning something like “from the very beginning”). Pedantically speaking, the DFT methods are sometimes denied the ab initio quality, because the actual functionals are based on empirical ingredients. However, since the existence a “pure” form of these fundamental objects is ensured at the level of basic theorems (and because the practice with DFT cannot be ranked lower than the ab initio HF), there is no problem in placing the spirit of DFT into the first-principle approaches. As a quick illustration of the power of today’s computational tools, we will “synthesize” the methanol molecule by the so-called geometry optimization procedures, starting from a chaotic bunch of atoms. The movie of the action, performed 294 4 Bond! Chemical Bond: Electronic Structure Methods at Work 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 Fig. 4.1 Frames from the steps of the structure optimization of methanol molecule by B3LYP density functional, with 6-311+G* basis set, starting with a strongly distorted arbitrary geometry. The bar on the left side (visible only in the first frames) shows the energy gradient, in relative scale in DFT frame (with B3LYP/6-311+G* setting) is shown in Fig. 4.1. It looks rather similar when other methods (HF, MP2, CCSD) are used, all being able to reconstitute a realistic equilibrium geometry, starting from a very arbitrary one, as seen in the first frame (left upper corner) in Fig. 4.1. The left-side bar of each frame shows, in relative scale, the evolution of the energy gradient, the magnitude that drives the geometry optimization (considered ended when a value close enough to zero is obtained). Thus, with the universal condition of reaching a minimum of the total energy, obeying the principles of quantum mechanics, the molecules, existent or 4.1 Molecular Structure by Computational Chemistry: A Brief Synopsis 295 imaginary, can be obtained inside the memory of today’s computers, rendered as a list of coordinates of the atoms, together with several basic properties (dipole moment, vibration or electronic spectra, if instructed by specific keywords). With non-demanding methods, such as HF or DFT, even the small computers, designed for personal use, can tackle a basic calculation, such as the computation of the molecular orbitals at a defined geometry, or undertaking the structure optimization. It is true that we have chosen a simple molecule, without a large palette of isomers, but this does not change the point about the power of first-principle methods. Table 4.1 shows results of different technical settings of the geometry optimization, the outcome being quite illustrative for the accepted ranking of several consecrated methods. The atom labeling is derived from Fig. 4.2, which compares directly the experimental and the B3LYP geometries. Visually, the experimental vs. computation match is excellent, the small mutual differences being practically imperceptible. The conformation of the molecule is with the in-plane trans configuration of a H–O–C–H (1) sequence (180° dihedral angle). Focusing concretely on the computed C–O bond length, compared with the value from a crystallographic characterization (Kirchner et al. 2008), one finds (see Table 4.1) a very good agreement from the side of the celebrated B3LYP method, the most popular option in the DFT approaches. This B3LYP result is almost coincident with those of more costly methods, such as the CCSD, suggesting the appeal for such a functional, belonging to the so-called hybrid class. The other functionals, illustrating different classes (BP86 with gradient corrections, while LDA confined to the simple local density approach) perform less well, in line with the level of their designed sophistication. Here, it happens that LDA and HF perform almost the same, a fortuitous result that may allow us to recall that at its very origin, before the DFT itself was stated, the local density approximation of Slater was intended just to find a cheaper calculation recipe, instead of ab initio HF, by an Table 4.1 The structural parameters (bond lengths, in Å, and angles, in degrees) for CH3OH molecule from experimental crystallography data (first numeric column) versus several calculation methods Param. Exp. DFT B3LYP BP86 LDA WFT HF MP2 CCSD CO 1.415 1.418 1.426 1.399 1.398 1.418 1.417 OH 0.840 0.963 0.974 0.971 0.939 0.959 0.959 CH(1) 0.980 1.091 1.100 1.100 1.081 1.089 1.092 CH(2) 0.980 1.099 1.109 1.109 1.087 1.096 1.100 CH(3) 0.980 1.099 1.109 1.109 1.087 1.096 1.100 COH 109.47 108.40 107.58 108.51 110.04 107.97 108.13 OCH(1) 109.48 106.79 106.62 107.02 107.21 106.56 106.68 OCH(2) 109.47 112.52 112.81 113.25 111.95 112.28 112.20 OCH(3) 109.47 112.52 112.81 113.25 111.95 112.28 112.20 The labeling of the atoms is shown in Fig. 4.2. The first half of the geometry optimization results shows outcomes from selected functionals, in the frame of density functional theory (DFT), while the other part is devoted to methods from the wave functional theory (WFT) frame 296 4 Bond! Chemical Bond: Electronic Structure Methods at Work (a) CH(2) CH(3) 0.980 Å 0.980 Å OCH(2) OH 0.840Å OCH(3) 109.5o COH 109.5 OCH(1) 109.5o CH(1) 0.980 Å CO 1.415 Å CH(3) 1.099 Å (b) CH(2) 1.099 Å OCH(2) OH 0.963 Å OCH(3) 112.5o COH 108.4 OCH(1) 106.8o CO 1.418 Å CH(1) 1.091 Å Fig. 4.2 Comparison of a experimental (crystallographic) and b computed (B3LYP level) geometries of the methanol molecule. The notation of the parameters corresponds to those in Table 4.1. empirical approximation of the exchange. In general, the LDA and HF perform differently, possibly with a better score for the DFT option. The HF computed geometry can be judged as the less accurate one, but, as a matter of fact, the relative deviation of the computed Hartree–Fock C–O bond length with respect of the experiment is about 1%, ultimately not so bad. Going now to other parameters, the OH and CH bond lengths, or the COH and OCH angles, something intriguing is observed. Apparently, all the methods are deviating rather much from experiment, while they are quite consistent with each other, in the range of values and trends. Although the comparison with experiment is considered supreme, and should be kept this way, we must however note that the limits of interpretation can affect the experiment itself, not the theoretical approaches only. Namely, the hydrogen atoms are barely visible in the X-ray diffraction methods, although not completely impossible to detect, with dedicated effort. Then, the X-ray analyses, which are themselves a fit of a model to detected diffraction peaks, are tributary, sometimes in a hidden manner, on presets in the variation of certain parameters. Therefore, in certain circumstances, the theoretical estimates can be more realistic than the experiment itself. As a tiny detail, note that the hydrogen atoms of the methyl group are slightly different, because the H(1) is in staggered conformation with the H from the hydroxyl group, while the others, H(2) and H(3), are at dihedral angles of about ±60°. Although this situation is expected to bring only small differences in bond lengths or angles, via inductive effects, or by trough space interaction, a tiny alteration is yet expected to be sensed. The experiment seems insensitive in finding differentiations in C–H parameters, yielding the same bond length, 0.98 Å. In turn, the calculations are detecting the infinitesimally different identities. The experiment finds the variation of O–C–H in a very narrow range, around the ideal tetrahedral angle, 109.47°, while the computation plays more realistically, retrieving a certain variation, as intuitively expected. At this series of parameters, once again, one notes a close parallelism between the computationally cheap B3LYP and the expensive CCSD. 4.1 Molecular Structure by Computational Chemistry: A Brief Synopsis 297 Although not suggesting dethroning the experiment from its deserved position, we emphasize with this occasion the need to view critically this side too. At the same time, one points to the potentially increasing capability of the perfectible ab initio methods to reproduce the molecular properties, keeping in mind also the remarks made in the introduction to this book, that the aim to photographically reproduce and predict the world should not be set as the supreme goal, heuristic and qualitative insight being sometimes more valuable. The capability of the methods to account basic molecular features is seen also for more complicated systems. We must accept, however, that the claim to ab initio construction of molecules (with the complete predictive account of their properties), even outside the field of experiment, is limited to some borders (e.g. the size of molecules, up to several tens of atoms). Other limitations may come not from molecular size and the related demands in computational effort (running job time and allocated memory), but from the intricacy of the electronic structure problems (e.g. needing many states, with delicate balance in their ordering, as may happen in magnetism or electronic spectroscopy). In molecules with firm covalent bonds, as for most of the organic species, the electrons are placed in well rationalized bonding schemes, in rather inert closed shells. In inorganic or coordination systems based on transition metals or lanthanide ions, the electrons are more “nervous”, moving easily between closely spaced states. In such situations, the calculation approach is more problematic. However, the obstacles in reaching good computational interpretation or prediction, either determined by technical limits of hardware, or because of incomplete precision or poor leverages, are not insurmountable, the theoretical approach being fully acknowledged nowadays as a counterpart to experiment. 4.2 Hartree–Fock Versus Density Functional Theory Computation Simple Samples This section will exemplify the simplest calculations, HF and DFT, with one of the most important small molecules of our world: water. As is known, or can be guessed, there are specific programs that can do this, when supplied with appropriate instructions. The input, basically a text file, provides the molecular geometry (experimental or optimized, or a certain guess, if the optimization is included as a step) and keywords choosing the method and specifying the desired state of molecule (charge and spin multiplicity, at the most elementary level). Since quantum chemical programs can be served by graphic interfaces, the producing of a text input can be replaced by mouse clicks, in picking fragments from databases, or constructing in 2D or 3D sketches, selecting options by pushing buttons on menus or windows. With graphic interfaces, driving a quantum chemical calculation, at least at routine levels, is not much different than driving a physical measurement, e.g. with spectrometers or diffractometers. 298 4 Bond! Chemical Bond: Electronic Structure Methods at Work User friendly interfaces are making basic quantum chemistry available for chemists from synthetic or instrumental branches, who can use the computational software as extension of their experimental hardware (e.g. as helper to assign a general pattern or specific peaks of the vibration spectra). Going back to the text-type input, Table 4.2 illustrates the calculation of the water molecule at fixed geometry (single point) with three representative codes: Gaussian (Frisch et al. 2009, Gaussian 09), GAMESS (General Atomic and Molecular Electronic Structure System) (Schmidt et al. 1993), and ADF (Amsterdam Density Functional) (ADF2012.01, SCM; te Velde et al. 2001; Fonseca et al. 1998). Table 4.2 Examples of input files for Gaussian, general atomic and molecular electronic structure system (GAMESS), and Amsterdam density functional (ADF) codes, taking the example of H2O molecule, computed by DFT using the BP86 functional and basis set at the level of triple-zeta with polarization Gaussian and GAMESS %chk=h2o.chk #P BP86/6-311+G* title water molecule in Gaussian ADF $ADFBIN/adf << eor XC GGA BP86 END 01 O 0.000000 0.000000 -0.011508 H 0.764982 0.000000 0.591320 H -0.764982 0.000000 0.591320 Basis Type TZP Core none End ############################## Atoms O 0.000000 0.000000 -0.011508 H 0.764982 0.000000 0.591320 H -0.764982 0.000000 0.591320 End $contrl scftyp=rhf units=angs icharg=0 mult=1 maxit=150 dfttyp=bp86 $end $BASIS gbasis=n311 ngauss=6 ndfunc=1 npfunc=1 $END $DATA title water molecule in GAMESS CNV 2 O 8 0.000000 0.000000 -0.011508 H 1 0.764982 0.000000 0.591320 $END Occupations A1 2 2 2 A2 0 B1 2 B2 2 End scf iterations 150 end End Input eor cp TAPE21 h2o.t21 4.2 Hartree–Fock Versus Density Functional Theory Computation Simple Samples 299 The presented input examples are arranged for DFT calculations, with the BP86 functional, in restricted Slater determinant wave function. If we want to do HF calculations in the Gaussian case, we must replace the BP86 keyword with HF or RHF. In GAMESS, we must erase the dfttyp specification (keeping the already introduced scftyp=rhf). In ADF, although not very natural in the philosophy of the code to work HF, one can replace the GGA BP86 line by the “Hartree–Fock” keyword. The geometry is kept the same in all the codes, being actually obtained by optimization with the ADF setting. Revisited in the other codes, the optimization may yield slightly different parameters, but the shifts are not really important. The Gaussian and GAMESS use the 6-311+G* GTO-type basis set, while ADF, the TZP (Triple Zeta with polarization) STO, a level somewhat comparable with the above GTO setting. The Gaussian and ADF have the possibility to produce binary files (defined here h2o.chk, respectively h2o.t21), useful in post-computational analyses or as predefined databasis for starting other calculations. In GAMESS, there is an implicit text file accompanying the output, usually named with DAT extension, containing the black box info of the calculations (e.g. orbitals coded as a stack of coefficients, usable for restarts). We will illustrate now the calculations done by the ADF, exploiting the occupation number controls. Here the symmetry is used, benefiting from classification of states according to the irreducible representations, specific for the point group of the molecule. We did not provide till here an introduction in symmetry groups. Succinct expositions, which can be visited, if necessary, are given in Sect. 4.3, and at the beginning of Chap. 8. Otherwise, we will assume the reader is accustomed with the basics of symmetry representations, or, at least provisionally will accept that a useful dichotomization of interaction types is realized with the help of such labels, specific for different patterns of molecular geometry. Thus, in the ADF input, the occupation numbers can be specified by lines labeled by the irreducible representations of the orbitals. In the C2v point group of the water molecule, we have a total symmetric representation, a1, for the 1 s core of oxygen, the a1 + b1 representation for the couple of O–H bonds, and a1 + b2 for the two lone pairs on the oxygen atom. If the orientation of the molecule is changed from the actual xz plane to yz, the labeling of bonds and lone pairs is interchanged. In the actual orientation, the occupied orbitals of the water molecule are labeled by their irreducible representation and a prefix of repetition, having 1a1, 2a1, and 3a1 for respective core, bond, and lone-pair types of the total symmetric components. The line “A1 2 2 2” visible in the ADF input example in the “Occupations” block enlists the double occupation of each of the abovementioned three a1 functions. The apparitions of the b-type representations are: 1b1 for the bond components and 1b2 for the lone pairs (non-bonding elements), corresponding to the “B1 2” and “B2 2” occupations in the input. The representation a2 does not have occupied elements, being specified as “A2 0” in input. The computed ordering of occupied orbitals is 1a1, 2a1, 1b1, 3a1, 1b2. As pointed out, in ADF it is possible to provoke fractional variations of occupations. Thus, running calculations with the “A1 2 2 1.99” and “A1 2 2 1.98” occupations, aside the regular “A1 2 2 2” closed shell occupation of the a1 subspace, one may trace an interpolation of the total energy regarding the last occupied 300 4 Bond! Chemical Bond: Electronic Structure Methods at Work total symmetric orbital, 3a1. Doing a similar change on the second position, “A1 2 1.99 2” and “A1 2 1.98 2”, data regarding the 2a1 orbital are gathered. In a different sort of numeric experiment, namely placing infinitesimal extra-electron population, e.g. “A1 2 2 2 0.01” and “A1 2 2 2 0.02”, one probes the first a1-type virtual, namely the fourth total symmetric component, 4a1. Similar numeric experiments can be done on the other components, e.g. “B1 2”, “B1 1.99” and “B1 1.98”, testing the occupied 1b1 orbital, while the “B1 2 0.01” and “B1 2 0.02” input lines are experimenting the 2b1 virtual, starting to occupy it incrementally. The fractional occupations can be tackled either with frozen orbitals, kept from the regular closed shell occupations, or accepting the iteration over the altered occupation numbers. Given the small variation in the occupation numbers, the frozen and variational results are numerically close. In Table 4.1, we confine ourselves to the frozen version. Any quantum chemistry program prints the total energy. In ADF, usually, this is not the absolute value, but relative to predefined fragments, which are the neutral atoms, if no preamble action is taken for their definition. It is possible to demand the absolute total energy, but this is not essential, keeping here the basic setting. Thus, picking the bonding energies from the suggested numeric experiments with occupation numbers in the orbitals “i”, one may fit the variation of the energy, relative to the integer occupation scheme: EðDni Þ  E0  vi
Dni þ 1 g
Dn2i : 2 i ð4:1Þ where E(Dni = 0) = E0 is the reference energy, the coefficients of quadratic fit taking the meaning of orbital electronegativity (vi) and hardness (ηi), as defined in Eqs. (3.40) and (3.41), as energy derivatives. To obtain the coefficients, the Dni = −0.01 and Dni = −0.02 differences are provoked in occupied orbitals (corresponding to the above discussed ni = 1.99 and ni = 1.98 fractional populations), while Dni = +0.01 and Dni = +0.02 for virtuals. The results of the quadratic interpolation are shown in Table 4.3. The first numeric column gives the energies of the computed Kohn–Sham orbitals, at the integer configuration. One observes that these show values very close to the first-order coefficients taken with reverted sign, ei = −vi, having a check of the Janak theorem from Eq. (3.40), with numeric derivatives. The second-order coefficients, assignable to the chemical hardness (electrorigidity) ,should be taken under the provision of somewhat unstable values, since, at the performed small displacement of population, the graphic representation of energy variation looks rather linear, without a curvature perceptible to the eye. This means that their role in the fit is not essential and the values are affected by a certain indetermination. To fix this issue, we must induce larger population shifts, which, on the other hand, start to affect the discussed match of genuine orbital energies with the numerical estimation of orbital electronegativities. Now, we can turn the computational play to the estimation of the ionization potentials. In DFT by ADF we can tackle this by handling occupation lists, this time with integers, e.g. subtracting the energy of the “A1 2 2 2” ground state from the 4.2 Hartree–Fock Versus Density Functional Theory Computation Simple Samples 301 Table 4.3 Orbital energy parameters: ei Kohn–Sham eigenvalues computed at regular occupation scheme and coefficients of quadratic interpolation, at infinitesimal population changes, assignable to orbital electronegativity (vi) and hardness (ηi) i ei vi ηi 2a1 −25.354 25.346 13.366 1b1 −13.174 13.166 12.570 −9.478 9.471 12.429 3a1 −7.308 7.300 12.327 1b2 4a1 −0.577 −0.619 3.907 2b1 1.434 1.300 0.343 10.054 10.069 6.822 2b2 14.829 14.832 7.476 1a2 The zero-order coefficient is E0 = −14.050 eV, for all the orbitals. The approximate ei = −vi relationship is a numeric test of the Janak theorem. The calculations are done with BP86/TZP setting in ADF. All values are in eV result of the run with the “A1 2 2 1” occupation, to get the ionization potential for the 3a1 orbital. The difference between “A1 2 2 2 1” energy and the “A1 2 2 2” one, will correspond to the electron affinity in the 4a1 function. However, this is not completely accurate in the physical sense, since in a restricted frame, the ADF treats the occupation number 1 as a half of a doubly occupied function (0.5a and 0.5b), a situation which is not really a restricted open shell canonicalization. More concretely, if we extrapolate back to HF energy formulas, in this way we will have (ai|ai)–(1/2)(ai|ia) couples of two-electron integrals (where i is the single occupied function and “a” a doubly occupied one), namely the scaling by (1/2) of closed shell scheme, instead of the genuine ROHF type term, (ai|ai)–(ai|ia). Then, we must switch to the option “unrestricted”, writing this keyword in the body of the input and specify the occupation in lists providing first the a occupations, then the b, separated by a slash symbol, e.g. “A1 1 1 1/1 1 0” for positive ion and “A1 1 1 1 1/1 1 1” for an extra-electron in the a1 set. Thus, using the unrestricted frame, the ionization potentials taken from the sequence of last occupied orbitals of the water molecule are I(1b1) = 18.94 eV, I(3a1) = 15.06 eV, and I(1b2) = 12.85 eV, matching well the respective experimental values: 18.51, 14.74, and 12.62 (all in eV). (Dutuit et al. 1985) Some other reports may have interchanged the b1 and b2 labels, because of a different convention in the orientation of the Cartesian coordinates. Observe that, doing finite differences to get the ionization potentials, we refrain from using DFT orbital energies as surrogates of these amounts. Although a certain correlation exists, occasionally a possible good match, the belief that Kohn– Sham orbital energies may render ionization potentials is not well substantiated. One may verify that there is a certain departure, indeed (compare the above I values with orbital energies from Table 4.3). At the same time, there is no need to charge the Kohn–Sham orbitals with this demand, while the obtaining of ionization potentials as proper energy differences is possible, with no big costs, and seems to work well. 302 4 Bond! Chemical Bond: Electronic Structure Methods at Work Table 4.5 VB2000 calculations for methanol in the GAMESS environment GAMESS preamble ! VB calculations on CH3OH $contrl scftyp=rhf units=angs icharg=0 mult=1 maxit=150 nosym=1 vbtyp=vb2000 $end $basis gbasis=n311 ngauss=6 ndfunc=1 npfunc=1 $end ! Defined the 6-311G* basis set $data title Methanol CN 1 C 6 -0.019804 0.000000 0.663187 O 8 0.121916 0.000000 -0.748146 H 1 0.987440 0.000000 1.080980 H 1 -0.543945 0.891877 1.033440 H 1 -0.543945 -0.891877 1.033440 H 1 -0.756997 0.000000 -1.141933 $end ! In the following VBGA and VBORB ! groups, the numeric indices 1-6 ! correspond to the above defined atoms. VB2000 calculations $VB2000 #! VB (4).VB(10) PRINTALL $VBGA 2:(1) => 2 2:(2) => 2 1-2 => 3 1-3 => 3 1-4 => 3 1-5 => 3 2-6 => 3 $02VBORB 0.70711 (2:(1)) + (2:(2)) 0.70711 (2:(1)) (2:(2)) $03VBORB 1->2 2->1 1->3 3->1 1->4 4->1 1->5 5->1 2->6 6->2 $02VBSTR 1 1122 $03VBSTR 1 1 2 3 4 5 6 7 8 9 10 $END 0.70711 0.70711 4.2 Hartree–Fock Versus Density Functional Theory Computation Simple Samples 303 Because of different numerical integration schemes in various codes, one may obtain results slightly dependent on the used program, even within the same setting. In comparing ADF with the other mentioned codes, there is also the reason of different type of basis sets, STOs, versus GTOs. However, the departure is small. Thus, the GAMESS input, shown in Table 4.2, yields for the sequence of the last four occupied KS orbitals the {−24.789, −12.658, −8.751, −6.577} values (in eV), close to the ei from the left column of Table 4.3, with ADF results. In actual implementation of GAMESS there are no controls for handling fractional or nonaufbau occupation schemes. However, aiming ionization potentials from specific orbitals, one may use the trick of reordering the orbital list, loaded from a previous calculation (a $VEC numeric block, introduced with $guess input instruction). Thus, adding in $guess the reordering option and the permuted list, e.g. iorder(2) = 3, 4, 2 saying that, starting from the second position, the new ordering is the mentioned series of indices, one places the 1b1 orbital as the last on the list, instead of the 1b2 HOMO. Then, with the option in rstrct=.true. in the $scf group, one enforces the run to hold this orbital list (instead of the energy ordered one), so that, putting the positive charge (and the spin multiplicity to doublet) in the header, the restricted open shell calculation will extract the electron from the desired 1b1 orbital. Here we can keep the restricted-open frame, since the occupation number 1 is treated in GAMESS as described in a definite ROHF format. Thus, repeating the ionization potentials with BP86/6-311+G* setting in GAMESS, one obtains for the {1b1, 3a1, 1b2} sequence the respective {18.69, 14.75, 12.60} estimation, extremely close to the experimental data, {18.51, 14.74, 12.62} (all in eV) (Dutuit et al. 1985). Once arrived at the GAMESS code, we can do a test on the Hartree–Fock side, checking the meaning of orbital energies as ionization potentials taken with reversed sign. Thus in the {1b1, 3a1, 1b2} occupied sequence, the HF orbital energies are −18.93, −15.53, and −13.58 (in eV), relatively well matching the range of experimental ionization potentials. Then, using the abovementioned handling in setting ionization potentials by customizing the order of MOs in GAMESS (and keeping the orbitals frozen, i.e. by taking one iteration only), we can check now the Koopmans theorem. By the corresponding energy differences between the positive molecules, with one electron extracted from the respective 1b1, 3a1, and 1b2 orbitals, and the neutral ground state, the 18.93, 15.53, and 13.58 values are obtained, retrieving exactly the statement of Koopmans theorem (with frozen orbitals). Thus, in this section, we partly emulated a practical hands-on session, showing the specifics of several codes and the routes of possible non-standard handling of calculations, beyond the simplest level of loading the molecule and choosing the method. In this way, basic theorems of Koopmans (for HF) and Janak (for DFT) were illustrated, finding also the route to chemical meaningful parameters, such as electronegativity and hardness. Recall, from the previous chapter, that we proposed a new term for chemical hardness, namely “electrorigidity”, to underline the status of companionship with electronegativity, and the meaning of resistance of the molecule against the deformation of the electronic cloud. 304 4.3 4 Bond! Chemical Bond: Electronic Structure Methods at Work Orbitals, the Building Blocks of Electronic Structure The molecular orbitals (MOs) are the icons of theoretical chemistry, venerated in various degrees, from chemists who seek heuristic interpretations in qualitative sketches, or simply admire the round and shiny shapes drawn by modern graphic interfaces, to the computational chemists who see them from a rather “desecrated”, utilitarian, and technical perspective. In either case, the orbitals are indispensable objects of the whole molecular science. A useful key to read MO diagrams is the symmetry, whenever the system has one (e.g. axes of rotation, reflection planes, the inversion in the sign of Cartesian coordinates, or a combination of such processes, that render the molecule identical to itself). The symmetry aspects are very powerful, as suggested also in the chapters dedicated to the atom, allowing partial predictions, even in advance of quantum calculation, or illuminating the post-computation scene. In certain cases the symmetry acts by producing sequences of equal orbital energies, such as the doubly degenerate representations labeled e, in corresponding point groups, which are reminiscent of ±m couples met in the discussion of the axial patterns in atomic orbitals. The symmetries associated with the Platonic solids (tetrahedron, octahedron, cube, icosahedron, dodecahedron) allow also triple degenerate series labeled by t. The exceptional cases of icosahedron and dodecahedron (actually, the same point group, Ih) can keep degeneracies of fourth and fifth order, labeled by h and g. The degeneracies play an important role when placed in the frontier zone, where the occupation with electrons ends. An incompletely filled degenerate orbital generates degenerate many-electronic ground states, which are not stable situations. In such cases, a break of symmetry occurs (Jahn–Teller effects), the system ending in non-degenerate states, produced along distorted molecular symmetries. Then, the trend to avoid such situations may generate magic numbers of preferred electron counts, that stabilize the occupation schemes having degenerate or quasi-degenerate frontier orbitals. Even when the symmetry does not imply degeneracy, the classification of the orbitals by irreducible representations of the point group is useful. Thus, the a and b labels behave as even and odd with respect of a rotation, the sign (coloring) of orbital lobes on equivalent molecular moieties being, respectively, the same or reverted. A similar situation occurs with the g and u subscripts, which stand for even or odd at inversion (changing position to antipodes). The prime and double prime superscripts mark the odd–even conjugation with respect of the reflection through a symmetry plane. After this short guided tour on symmetry, it is worth finally mentioning the creator of group theory, Évariste Galois, and the tragic sacrifice made at its very inception. Galois died as a result of a duel (in 1832, at age 21) and wrote down, desperately, the bases of the group theory on the night before he died, instead of resting or trying to fight for his life. The theory survived. In the following, a contemplation of Kohn–Sham molecular orbitals of the water molecule is proposed (Fig. 4.3). It resulted from a DFT calculation with BP86 4.3 Orbitals, the Building Blocks of Electronic Structure 305 functional and the 6-311+g* basis set, but the visual aspect remains the same over other functionals and basis sets, and even at the Hartree–Fock level.The graphics of orbital shapes or the electronic density maps was realized, throughout this book, by the Molekel (Flükiger et al. 2000-2002) or WxMacMolPlt (Bode and Gordon 1998) codes. When the use of ADF is mentioned, the related plots are realized with the routine implemented in the package. (ADF2012.01, SCM). The orbitals are labeled by their representation in the C2v point group and a prefix giving the ordering number in its class. Skipping the 1a1 deepest level, which is just the 1s core function of the oxygen, observe (Fig. 4.3) that 2a1 has the same sign (coloring) over the whole surface, its density wrapping all the atoms of the molecule. Due to this feature, the 2a1 is the most bonding function. With respect of a given pair of proximal atoms, an orbital acts as bonding when the function builds a continuous bridging electronic density, due to atomic components overlapping with the same sign of the lobes. Sed contra, the overlap of local lobes, colored in different signs, destroys the accumulation of the electronic density in the inter-nuclear zones and determines antibonding (destabilizing) effects. This is because the depletion of electronic density between nuclei let their repeal contribute more to the total energy. The more local bonding areas appear in an MO, the lowest is its energy, as is the case of 2a1. The 1b1 function has bonding characteristics on each OH axis, but has a nodal plane. As noticed in the atomic case, the apparition of nodal planes increases the energy of the orbitals. The 1b2 is a bonding function, but slightly less cohesive than 2a1 (see also the energy costs for ionization, discussed in the previous section). The 2a1 and 1b1 look like in-phase (with the same sign) and out-of-phase (opposite signed) combinations of the bonds located on the two OH lines. The former is totally symmetric (remains unchanged at rotation around the axis bisecting the HOH angle, or at reflection through the molecular plane and its Occupied 1b 1 2a1 Bonding Virtuals 1b 2 3a1 Nonbonding 2b 1 4a1 Nonbonding Fig. 4.3 The orbitals of H2O molecule, from a BP86/6-311+G* DFT calculation. The qualitative pattern is the same in other single-determinant settings, including the Hartree–Fock case 306 4 Bond! Chemical Bond: Electronic Structure Methods at Work perpendicular). The 1b1 is asymmetric at rotation and at the perpendicular plane, having, in the circumstances of the C2v point group, the same symmetry as the px orbital (when the molecular plane is xz). The 3a1 and 1b2 orbitals can be characterized as non-bonding, carrying electron pairs pointing outside the inter-atomic axes. The last one is almost a pure py orbital, spanning a symmetry channel distinct from other occupied orbitals. With this occasion we note the role of symmetry as “separator” of different interaction paths. The situation 3a1 and 1b2 orbitals as lodges for the two lone pairs of the oxygen seems quite different from the customary idea of sp3 hybridization of the water molecule, with tetrahedral orientation of four electron couples, two from the bonds and two from equivalent lone pairs, with angular pattern of their prominent lobes. Actually, the MO picture looks merely like an sp2 hybrid. From the symmetry point of view the situation should be sp2-alike, without meaning however that this is the true hybridization of the oxygen in the water molecule. As will be discussed later, chemists’ beloved concepts can be tricky, or superfluous, from the point of view of computational people. A single sp3 hybrid, from the two lone pairs, does not obey alone the twofold symmetry axis. However, taking the two sp3-type lone pair hybrids, one can obtain symmetry adapted combinations, a sum (which points like the sp2 function), and a difference (looking like a p orbital). Such an argument may be consolatory, not dismissing the sp3 idea, but not yet comforting, once we see the lability of the hybridization assignment, with the “double personality”, sp2 and sp3, of the water molecule. But, after all, the indiscernibility between different definitions of the reality is not so strange in the quantum world. Finally, one may quickly note the 4a1 and 2b1 couple of empty orbitals, which, located mostly on the oxygen, are of non-bonding type. These can be accessed by physical processes, or by a computation setting, like CASSCF (including them in the active space). Otherwise, in a single-determinant frame, the virtuals play no direct role in the bonding. In the first instance, the virtuals are left over from the optimization processes that regard the occupied or active orbitals, but play a subtle role in the properties of the molecule, since no variation in density can occur (e.g. those modulated by vibrations) without the occupied-virtual MOs remixing. Also, of course, the electron excitations (the processes giving the color) will not be possible without the empty orbital rooms. In the following, the orbital choice in the frame of Valence Bond (VB) calculations will be illustrated. In this way, we move closer to a chemically intuitive picture, with the cost of a somewhat more intricate input control. Thus, the orbitals will be prepared to resemble two lone pairs, plus two hybrids aimed for two unpaired electrons on oxygen, which will be paired with a companion couple of functions resembling the hydrogen atoms. The leverages of current VB programs, e.g. VB2000 (Li et al. 2007; Li and McWeeny 2002) or XMVB (Song et al. 2005, 2012) codes, allow producing such functions in a user-friendly manner, starting with localized orbitals from a Hartree–Fock preamble. However, for a first-principles aspect, we will initiate the input with a set of merged atomic orbitals, mixed “by hand”, to tune the s–p hybridization for the VB initial functions dedicated to the lone pairs and active SOMOs on the oxygen. 4.3 Orbitals, the Building Blocks of Electronic Structure 307 Placing the hydrogen atoms in the xz plane (with z as symmetry axis), the oxygen hybrids prepared for bonds will be made with px, pz, and s functions. In turn, the lone pairs are constructed from py, pz, and s. Observe that, while px and py functions are distinctly dedicated to bonds or lone pairs, respectively, the s and pz are shared among all the hybrids. In the circumstances, there is a simple recipe for such a mixture. Thus, say that r is the coefficient of s in the hybrids prepared for bonds. pffiffiffiffiffiffiffiffiffiffiffiffiffi Then, for normalization reasons, the coefficient of pz should be 1  r 2 . At the same pffiffiffiffiffiffiffiffiffiffiffiffiffi time, the portion of s remaining for the lone pair hybrids is 1  r 2 . The amount of p in lone pairs z is the normalization counterpart with respect of s coefficient, namely −r, the minus being needed for the orthogonality to the previous functions, and also for determining the orientation of the lobes opposite to hydrogens. Finally, the r and pffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi 1  r 2 coefficients discussed for the two groups of hybrids must be divided by 2, because is shared by two elements in each set. This construction is outlined again in the next chapter, in the section devoted to hybridization. Figure 4.4 shows three points from the r-dependent tuning of the basis. The case (a) occurs at r = 0, when the bond has no s content, being made of two p-type orbitals with lobes at 90° (the in-phase and out-of-phase combinations of px and pz), while the lone pairs are enforced as sp-alike hybrids along the y axes, with 180° between their main lobes (in-phase and out-of-phase mixing of s and py functions). The (b) case, in the middle panel of Fig. 4.4, stands for equal sharing of s perpffiffiffi centage among the orbitals at r ¼ 1= 2, which corresponds to the sp3 hybridization. The (c) is the reversion of case (a), at r = 1, when the s content is entirely 1 2 3 4 5 6 (a) (b) (c) Lone Pairs O-side bonding H-side bonding Fig. 4.4 Different choices of VB starting orbitals, tuning the mutual variation of s:p ratio in oxygen hybrids dedicated lone pairs (components 1 and 2) versus those dedicated to spin coupling (components 3 and 4) against the hydrogen type functions (components 5 and 6). The Rumer diagram corresponds to the following pairing: {1,1}, {2,2}, {3,5}, {4,6} 308 4 Bond! Chemical Bond: Electronic Structure Methods at Work invested into OH bonds, letting the lone pairs stand as non-hybridized py ± pz functions. This is quite inconvenient for making bonds, since the prepared functions are pointing at a 180° angle, not toward the hydrogen positions. In a complete active space option, like CASSCF or CASVB, the transformation of the basis is superfluous, leading to the same poly-electronic states. However, here the complete space is not taken, preferring a manner driven by the following Rumer coupling schemes: 1 1 2 2 3 5 4 6 1 1 2 2 3 3 4 6 1 1 2 2 3 5 4 4 1 1 2 2 3 3 4 4 The proposed Rumer configurations are interpreted as follows. The function represented by the table from the left contains in the first two couples, {1,1}, {2,2} the lone pairs, where the opposed spins are occupying the same orbitals, while the {3,5} and {4,6} columns are describing the two O–H bonds as spin-paired factors. The other configurations can described as ionic, transforming in lone pairs the hybrids oriented toward the hydrogen atoms, {3,3} or {4,4}, as in the second and third Rumer diagrams, or both of them, as in the fourth configuration. As primary interpretation, the all-lone pairs configuration, {1,1}, {2,2},{3,3}, {4,4}, (the right-side table) reduces the role of hydrogens to two point charges. However, in a self-consistent regime, the hydrogen AOs are mixing into the result, this single determinant VB function becoming entirely equivalent to the Hartree–Fock case. In variational VB, the change in the orbital set is reflected in the coefficients of resonance structures, all the options leading to the same self-consistent result in the total energy, represented by the horizontal line in the lower part of Fig. 4.5. However, hindering the variational step, we designed numeric experiments probing the idea of directed valence. These are illustrated by the curves in Fig. 4.5, as a function of the previously described mixing parameter, r. The line (a) shows the non-iterative approach, just merging the hybridized atomic orbitals of oxygen with Fig. 4.5 Different VB calculation of electron pairing for the water molecule in active space options considered in the text, i.e. (a)– (d) cases -75.5 E(Hartree) -75.6 VB(8) CI only VB(4) freeze HF VB(4).VB(4) VB$(8) -75.7 -75.8 (a) -75.9 (b) -76 (c) -76.1 (d) r -76.2 0 0.2 0.4 0.6 0.8 1 4.3 Orbitals, the Building Blocks of Electronic Structure 309 the 1s AOs of the hydrogens and doing the configuration interaction (CI) in the described basis of the four Rumer functions. In this case, the enforced electronic density shows a minimum nearby the r = 0 edge that places the two lone pairs as distantly as possible. It is not exactly r = 0 because the electrostatic attraction from the hydrogen nuclei still determines a trend for the aligned bond hybrids, without encompassing the repeal between the clouds of the filled lone pairs. Another numeric experiment keeps the doubly occupied components (core and lone pairs) frozen, while the {3,5}{4,6} part of the Rumer wave function is allowed to relax. Shown in the (c) curve, this situation produces a minimum at the middle of the r-scale, reflecting in better conditions the needs of the established OH bonds. The (b) curve, which shows a definite minimum nearby the ideal tetrahedral angle (r = 0.707), corresponds to the situation where the lone pairs are allowed to do orbital optimization, but the intermix with the subspace of bonds is yet restricted, treating the {1,1}, {2,2} and {3,5}, {4,6} sequences as different VB blocks. The series of different controls of the hindered variational experiments suggests the hidden role of the hybridization, as effective engine, fueled by electron correlation, if a phenomenology based on atomic orbital configurations is adopted. The discussion made here in terms of VB, a rather niche technique, can be transferred to other many-configuration methods, of more frequent use (e.g. using the same type of atomic and hybrid bases in a CASSCF). As a general conclusion: the orbitals are palpable results of the calculations. The interpretation of their shape and energy ordering can be relevant in qualitative respects, serving as ideograms in the messages exchanged with experimental chemists, who expect heuristic clues, as quintessential distillation of computation and theoretical outcome. The interpretation can be done on canonical orbitals, automatically produced by a given procedure (HF, DFT, CASSCF, CASVB), but a bit of advanced handling (e.g. numeric experiments with transformed sets) can bring a supplement of insightfulness. 4.4 4.4.1 The H2 Molecule: The Simplest Bond Prototype. Phenomenological Models and Calculation Methods The Spin-Coupling Phenomenology of the Chemical Bond The simplest model of a bond can be constructed taking two electrons in two orbitals. Even one electron will be enough to bind the system, but we choose to start with something representative for the common heuristics of the chemical bond, seen as a result of the electron pairing. Having in mind the H2 example, the model is drawn, however, a bit more generally, considering different atomic functions, a and b, on the sites of the AB generic diatomic molecule. The problem is, in certain parts, similar to the two-electrons two-orbitals case of the helium atom (see Chap. 2), treated as configuration interaction over the 1s and 2s orbitals, except that now the functions 310 4 Bond! Chemical Bond: Electronic Structure Methods at Work are located on different centers and considered non-orthogonal. There are six Slater determinant configurations to be formed in this frame, four keeping one particle per site: |aaba|, |aabb|, |abba|, and |abbb|, while two are implying charge-transfer configurations: |aaab| and |babb|. The Sz spin-projections of the first four elements are respectively 1, 0, 0, and −1, while the ionic configurations are both with Sz = 0. As learnt from the situation of the helium atom, the former four components are determining two states, a singlet and a triplet. The Sz = ± 1 determinants are each isolated from the others, in the configuration interaction matrix, because of the spin projection selection rules (all integrals will imply products of a a electron with a b one, vanishing). At the same time, aside the Sz = ± 1 basis components, the triplet must collect a third one, with Sz = 0, resulting from the combination of |aabb| and | abba| elements. For tractability, one assumes now that the ionic configurations can be ruled out, in the first approximation, setting then the problem in the basis of the |aabb| and |abba| Slater determinants. Considering then non-orthogonal orbitals, with sab ¼ hajbi, but normalized, so that saa = 1 and sab = 1, one forms the Slater determinant overlap matrix: S¼  1 s2ab  s2ab : 1 ð4:2Þ This is done with the rules mentioned in Sect. 2.5.2, dedicated to the non-orthogonal wave functions. Namely, the poly-electronic overlap between two Slater determinants is the determinant of the matrix running on rows the orbitals from “bra” and on columns those from “ket” (or vice versa). Then, making the matrix on the indices due to |aabb| and |abba|, one observes that the diagonal elements are null, having opposed spins, while the non-diagonal positions are entries for the same spin, but different sites, containing then, each, the sab elements. The expression of the 2  2 determinant with null elements on the main axis and sab as the two non-diagonal elements is −s2ab, of course. The Hamiltonian matrix is: H¼  haa þ hbb þ Qab 2hab sab  Jab  2hab sab  Jab : haa þ hbb þ Qab ð4:3Þ Here the h’s are one-electron matrix elements, the two-electron integrals being the Coulomb and exchange types, Qab = (ab|ab) and Jab = (ab|ba), respectively. The eigenvalue problem, equivalent to the solving of det|H − ES| = 0 with respect of E, leads to: E ¼ ðhaa þ hbb þ Qab Þ  ð2hab sab þ Jab Þ : 1  s2ab ð4:4Þ The E− solution corresponds to the in-phase combination of the defined Slater  pffiffiffi    basis determinants,ðaa bb  þ ab ba Þ= 2, pertaining to the triplet, because it can be produced from the jaa ba j higher projection applying the decrement ladder operator: 4.4 The H2 Molecule … 311  pffiffiffi    ^ S jaa ba j. It remains that the orthogonal complement, ðaa bb   ab ba Þ= 2, is the singlet, with the E+ energy. We must point out that it is not an error in assigning the minus subscript [inspired from the sign of the second parenthesis in the numerator of Eq. (4.4)] to the energy of the wave function, having the plus sign between the determinants. Actually, the sign of wave function components can be changed by convening a permutation within the second determinant: |baab| = −|abba|. This changes the sign of the non-diagonal elements in the previously expressed overlap and Hamiltonian matrices, keeping the same solutions. The sab, Qab, and Jab are positive quantities, while haa, hbb, and hab are negative. Then, it appears that a vanishing or small overlap integral determines the triplet as ground state, while the singlet becomes lower when the negative 2habsab quantity predominates over the positive Jab. The singlet ground state is the situation describing the bonded state, if we invoke, for concreteness, the hydrogen molecule as prototype, observing then the role of overlap in establishing a chemical bond. In the basic chemistry curriculum it is taught that the chemical bond is made by the spin pairing. Conversely, the parallel spins may be regarded as the source of anti-bonding effects, a rule kept in the circumstance of the overlap terms predominating over the pure exchange. The bonding as electron pairing is then a phenomenological statement, in the sense of taking a thing as it appears to be. The term follows, etymologically, the Greek word phainómenon, meaning “the thing that appears”. In the actual conjuncture, it can be taken as interpretation putting the emphasis on certain visible characteristics that may represent, in simplified manner, a hidden deeper complexity. Characterizing the chemical bond as spin coupling is the visible end of a more complicated causality. In fact, the spin pairing is itself a consequence of the mechanisms determining the bond. The electron pair as signature of the bond is, somewhat, a figure of speech. Nominally, the spin coupling can refer to dipolar interaction of the magnetic moments carried by individual electrons. But, the dipolar coupling is, by orders of magnitude, lower than the energy invested in chemical bond, so that it cannot really determine the cohesion. However, one may take a hint from the dipolar coupling formula, that contains the scalar product of the involved magnetic moments, ~ la
~ lb , ^ which should lead, in the quantum mechanics of spin-based moments, to the Sa
^Sb operator (discarding the corresponding gyromagnetic factors). In the classical acceptation, the sign of the spin product is positive when these are aligned parallel and negative for the anti-parallel case. Noticing that the above formulas for the singlet and triplet are twinned by a sign switch relationship, the idea would be to use the opposite signs generated by a scalar product, to comprise in a single formulation both cases from (4.4). Thus, we want something like:   ES ¼ p þ ^ Sa
^ Sb S q; ð4:5Þ which should account simultaneously for both ES=1=E− and ES=0=E+ solutions. The scalar product of two spins is obtained by an artifice involving the square of their vector sum: 312 4 Bond! Chemical Bond: Electronic Structure Methods at Work ^ S2 ¼ ^ S2a þ ^ S2b þ 2^ Sa
^ Sb : ð4:6Þ The square of a momentum operator can be transformed into a scalar yielding the corresponding eigenvalue, having then a closed expression for the spin product:    1 ^ Sa
^ Sb S ¼ SðS þ 1Þ  Sa ðSa þ 1Þ  Sb ðSb þ 1Þ ; 2 ð4:7Þ as function of the total spin, S. Replacing the spin quantum numbers, Sa = 1/2 and Sb = 1/2, the above expression yields 1/4 for S = 1 and −3/4 for S = 0. The signs correspond to the qualitative expectations for parallel and anti-parallel vectors, but the concrete values are due to the quantum specifics. Expanding in series the Eq. (4.4) and identifying the results as a function of the spin operator from (4.6), the following factors are identified: 1 1 p ¼ haa þ hbb þ Qab  Jab  hab sab þ s2ab ðhaa þ hbb þ Qab Þ þ


2 2 q ¼ 2Jab  4hab sab þ 2s2ab ðhaa þ hbb þ Qab Þ þ


ð4:8Þ ð4:9Þ The p term is the spinless part, depending only on the orbital placement of the electrons (one on the site a, the other on the site b), while the second term, q, senses the spin configuration. Actually, the expression from (4.5) is an effective Hamiltonian, which can be used to describe the spin states in the model system. Equation (4.5) is reformulated as ^ ¼ H   1 eff eff ^ haa þ hbb þ Qab  Jab Sa
^Sb ;  2Jab 2 ð4:10Þ where the effective exchange coupling cumulates other terms, e.g. eff ¼ Jab þ 2hab sab þ


. The ignored configuration interaction with higher terms Jab can be also thought as tacitly incorporated in the effective coupling parameter. The above approach is in the style of the Valence Bond (VB) method, working with orbitals localized on the sites, allowed to be non-orthogonal, without the help of an intermediate level of construction by molecular orbitals (MO) methods, taken as linear combination of atomic orbitals (AOs). The MO theory is a generic name for a large variety of variational methodologies, as discussed in previous chapters and sections. When constructed for the H2 alike molecule, the spins on the a and b sites are both 1/2, but the formalism carried by Eq. (4.10) is more general, known as Heisenberg–Dirac–van Vleck (HDvV) spin Hamiltonian (Heisenberg 1928; van Vleck and Sherman 1935; Anderson 1959). In the purposes useful for molecular magnetism, where this effective model is used, the local spins implied in the 2J ^ Sa
^ Sb coupling can take any non-null value (positive integers or half integers). With the help of Eq. (4.10), and discarding the constant terms, one finds that, in 4.4 The H2 Molecule … 313 general binuclears, the Hamiltonian becomes directly resolved into the ES = −JS (S + 1) spin-dependent energy levels, with the spin quantum number varying between the S = |Sa − Sb| and Sa + Sb limits. 4.4.2 Model Calculations on H2 In the following, the dihydrogen molecule is proposed as an example for different calculations, taking the simplest setting: the STO-3G basis set. According to all the current standards, this is the worst level of numerical precision, but it is not accuracy that is searched for now. It is a compromise for keeping the situation close to the above discussed conceptual representation, in the two-orbitals with two-electrons problem, while using standard codes and procedures representing the current level of quantum chemistry. The energy curves for different calculations are shown in Figs. 4.6 and 4.7, the graphic of related species of orbitals being illustrated in Fig. 4.8. For the H2 case, in the setting of one orbital per center, the canonical molecular orbitals are simply decided by symmetry, as in-phase and out-of-phase combinations of equivalent atomic orbitals, with factors dependent on their overlap integrals, sab: 1 rg ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ða þ bÞ; 2ð1 þ sab Þ ð4:11Þ 1 ru ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ða  bÞ: 2ð1  sab Þ ð4:12Þ The notation of the orbitals recalls that the overlap of the s-type orbitals has the r pattern (circular perimeter, with unique sign on the surface of a section 1.5 1Σ g(1) E(Hartree) Fig. 4.6 The CASSCF  CI energies, as function of inter-atomic distance, for the simplified model of H2 molecule as the two-orbitals two-electron problem, taken with the crude STO-3G basis set CASSCF 1 1Σ u 0.5 3Σ u 1Σ g(0) 0 0.5 1 R(Å) 1.5 2 2.5 Fig. 4.7 First and second states of H2 molecule, as function of inter-atomic distance, accounted with different methods, in the simplified test of STO-3G basis set 4 Bond! Chemical Bond: Electronic Structure Methods at Work 1 ROHF, UHF(Sz=1) E(Hartree) 314 0.5 RHF UHF(Sz=0) RHF, UHF(Sz=0) CASSCF R(Å) 0 0.5 1 1.5 2 2.5 perpendicular to the bond). The g and u are traditional labels for symmetric or asymmetric behavior at inversion through the origin of axes. Namely, the g function (from German gerade) remains the same when all the Cartesian coordinates are changing the sign (from x, y, z to −x, −y, −z), while in the u case, the inversion leads to a function with changed sign (otherwise, identical as module with the initial one). It is clear that the plus and minus signs forming the in-phase and out-of-phase combinations are leading to the g and u behaviors. The even g- type orbital is lower in energy, having bonding character (a map with contours enveloping the nuclei as seen in the left–bottom corner of Fig. 4.8). In this particular situation, of minimal basis set, with one rg and one ru, the MOs are independent from the Self-Consistent Field (SCF) method: restricted or unrestricted Hartree–Fock (RHF or UHF), Density Functional Theory (DFT), or Complete Active Space (CAS). Later on, a methodology intentionally enforcing a non-equivalence of the sites, called Broken Symmetry (BS), workable in unrestricted HF or DFT frames, is briefly illustrated. With richer basis sets this simple relationship of symmetry uniquely determined MOs is no longer kept, the SCF procedures entering in action.         In terms of the molecular orbitals, there are six configurations:rag rbg , rag rau ,        a b  b a  b b rg ru , rg rg , rg ru , and jrau rbu j. Although the symmetry issues are not detailed here, one may suggest that the r nature of all the orbitals is reflected in the same symmetry label for the poly-electronic wave functions, ascribed as capital (R). The g and u labels are obeying, in the orbital products, formed as lines of Slater determinants, the same multiplication rules as the + and − signs. Thus, the first and last Slater determinants in the above enumeration show the g symmetry, because gg = uu = g, while the components occupying orbitals with opposed symmetries are of u type, since gu = u. As in the previous discussion, with localized atomic orbitals, a set of four spin-flip combinations—aa, ab, ba and bb—gives rise 315 Broken Symmetry Orbitals 4.4 The H2 Molecule … µ2 µ1 Localized Orbitals χb Canonical Orbitals χa σu σg R=0.75 Å R=1.5 Å R=2.5 Å Fig. 4.8 Different sets of orbitals for H2 molecule, taken at selected inter-nuclear distances. The bottom line shows the canonical orbitals from a CASSCF calculation. Because of the minimal basis set and symmetry, the MOs are the same as those resulting from RHF. The middle line shows functions localized on sites a and b, obtained by the transformation of canonical MOs with the inverse LCAO matrix. The upper line shows the orbitals resulted in the Broken Symmetry Unrestricted Hartree–Fock (BS-UHF) attempt. At small distances, implying large overlap between AOs, the BS-MOs resemble the canonical MOs, while in the weakly interacting situations, they are becoming localized 316 4 Bond! Chemical Bond: Electronic Structure Methods at Work to a singlet and a triplet. Therefore, combining the above reasons, one finds two symmetric spin singlet functions and two asymmetric (one singlet and one triplet). The lowest function is those mainly based on the double occupation of the lowest molecular orbital, having a small contribution from the configuration interaction (CI) with the Slater determinant with the same g-symmetry, produced by the double occupation of the anti-bonding function:  1 W1 Rgð0Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi jrag rbg j  kjrau rbu j ; 1 þ k2 ð4:13Þ with k as mixing coefficient.DThe k got a negative E factor because the non-diagonal a b a b element of the CI, namely jrg rg jjH jjru ru j ¼ ðrag rbg jrau rbu Þ, is positive, being dominated by the large one-center Coulomb integral, (aa|aa) = (bb|bb) > 0. In such a case, the negative mixing coefficient gives the lowest expectation value of the energy. As a parenthesis, this situation is in reverse analogy with the MO diatomic schemes, where the negative non-diagonal Fock elements are determining the combination with positive mixing, a + b, to form the lowest eigenvalue. As energy ordering, the states following upwards the ground state are of odd symmetry, because these are generated from configurations with one electron in rg and one ru. In the series of excited states, the triplet is lower, because it benefits from the stabilization resulted from a negative sign in the front of a positive exchange integral: ðrg ru jrg ru Þðrg ru jru rg Þ. The spin triplet is represented by a set of three functions: W3 Ru 8 a a jr r j >  < g u 1ffiffi a b b a p p1ffiffi jra rb j  jra rb j ; jr r j þ jr r j ¼ ¼ g g u u g u g u 2 2 > : b b jrg ru j ð4:14Þ The asymmetric spin singlet is equated in this simple model by the ðrg ru jrg ru Þ þ ðrg ru jru rg Þ energy term, having the function: 1  1  W1 Ru ¼ pffiffiffi jrag rbu j  jrbg rau j ¼ pffiffiffi jrag rbu j þ jrau rbg j : 2 2 ð4:15Þ The highest spectral term is those produced having as dominant term the double occupation of the anti-bonding function, linked by CI to the ground state:  1 W1 Rgð2Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi kjrag rbg j þ jrau rbu j : 1 þ k2 ð4:16Þ In this case, the CAS (implying orbital optimization) is equivalent to the CI, because the orbitals are already prepared, by symmetry. Moreover, the full-CI limit is already fulfilled, because of the imposed small basis set. Inside a CAS space, or 4.4 The H2 Molecule … 317 in the full-CI, the remixing of orbitals by unitary transformations produces the same eigenvalues, while, obviously, the eigenvectors are changed to represent the new Slater determinant basis. In this case, one may propose the unitary transformation that leads to functions localized on the two sites:  1  va ¼ pffiffiffi rg þ ru ; 2  1  vb ¼ pffiffiffi rg  ru : 2 ð4:17Þ ð4:18Þ Note that the va and vb orbitals are not the a and b pure AOs. The v orbital dedicated to one site contains a tail on the other center, which accomplishes the orthogonalization relationship. The genuine AOs are allowed to overlap, as described previously in the Valence Bond type of modeling. Rewriting the CI (or CAS) in the new basis, the ground state is: 1k  W1 Rgð0Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi jvaa vba j þ jvab vbb j 2 1 þ k2 1þk  þ pffiffiffiffiffiffiffiffiffiffiffiffiffi jvaa vbb j þ jvab vba j : 2 1 þ k2 ð4:19Þ As the inter-center separation grows, the energy gap between occupied and virtual MOs is lowered. At infinite separation, the rg and ru orbital energies coincide, so that the implication of the two configurations from (4.13) becomes equivalent, reaching k = 1. One observes that at k = 1 the (4.19) form of the ground state tends to keep the Slater determinants with one electron per site. At large separation, the orthogonal localized orbitals are also becoming closer to the AOs themselves. At the same time, as seen from Fig. 4.6, the energy of the triplet is coalescing with the singlet ground state. Then, after the bond breaks, all the states made of one particle per site tend toward equal energies, meaning that these eigenvectors can be mixed arbitrarily, to obtain equivalent descriptions of the system. In this circumstance, one may interpret that, when the interaction is switched off, the most convenient situation is reached with atomic configurations, any of the four spin-flips being equally probable. The ionic situations with two electrons at one site (hydride anion), while the other is the naked proton, are higher in energy (because of large Coulomb repulsion in the doubly occupied AOs). At large distances, the one-site inter-electron repulsion is not counterbalanced by the inter-center electrons–nuclei Coulomb attraction, while at smaller separations this effect can bring the charge-transfer configurations into the composition of the ground state. The Hartree–Fock is the single determinant approximation, based here on the double occupation of the rg orbital. In Fig. 4.7 one observes that the RHF curve is 318 4 Bond! Chemical Bond: Electronic Structure Methods at Work above the CAS line, along the whole distance dependence, the departure being an illustration of the correlation energy, as it is retrieved in the crude computational experiment. The RHF curve has an abnormal trend after the bond breaks: it is not meeting asymptotically the triplet level. In the circumstances of intentionally small active space, the triplet appears only once in the spectrum, the CASSCF result being therefore the same with those produced by Hartree–Fock in restricted-open-shell (ROHF) or unrestricted (UHF) settings. In general, with larger basis sets and for more complex systems, this coincidence does not happen. To understand the drawbacks of the RHF solution in a long range regime, we must take a look at the expansion of MO into AOs-type basis: WRHF ¼ jrag rbg j ¼  a b  1 ja a j þ jaa bb j þ jba ab j þ jba bb j : 2ð1 þ sab Þ ð4:20Þ The four AO configurations have equal weights, meaning equal probability for having zero, one or two particles per site. For two hydrogen atoms placed at one kilometer distance it will, of course, not be reasonable to think of the electron transfer as permitted so easily, at equal footing with the spin flipping. The triplet in localized basis is: W 3 Ru 8 a a jv v j >  < ab b b ¼ p1ffiffi2 jvaa vb j þ jvba vab j ¼ p1ffiffi2 jvaa vb j  jvab vba j ; > : b b jva vb j ð4:21Þ Permuting the a and b indices (which represents the inversion operation) the sign of function is changed, justifying the u label. It appears then that, in the weakly interacting regime, which occurs not only at bond breaking, but is the very nature of other bonding regimes (such as in the coordination compounds), the single determinant HF is not a good model. Several problems, such as magnetism, demand multi-configuration approaches to reflect the physics of the spin-driven exchange interactions, or the strong CI regime imposed by quasi-degenerate states (i.e. series with small energy spacing). Although the single determinant methods are, in principle, unsatisfactory, the unrestricted frame allows numeric experiments which can recover a part of the handicap. Namely, the Broken Symmetry (BS) treatment, as its name suggests, attempts to induce non-equivalence of the orbitals, so that, at large distances, the wave function tends (in the given case) to the |aabb| or |abba| Slater determinants. These solutions are able to mimic the preference of the dissociated H2 molecule for neutral atoms with one unpaired electron per site, as seen in upper-right corner of Fig. 4.8. The BS configurations are not real states. However, the BS determinants are good surrogates for the weakly interacting systems, since the |aabb| ± |abba| symmetry combinations will have practically the same energy as the expectation values of the |aabb| or |abba| unrestricted BS single determinants. From another perspective, the BS is an experiment enabling parameters which can lead to the 4.4 The H2 Molecule … 319 complete description of the system, when replaced in corresponding models. More precisely, the J exchange coupling value obtained in the BS frame can be used as parameter in the HDvV Hamiltonian, which is a VB-type phenomenology. In other words, one may emulate from properly designed single determinant calculations the parameters needed for the multi-configuration description. Broken Symmetry is mostly used in conjunction with DFT, as a very practical tool nowadays in molecular magnetism, being applied on transition metal complexes with coupled paramagnetic sites. However, for the first illustrative purpose, we used it on an HF wave function and on a system a bit non-orthodox, since the covalent systems like the H2 molecule are not in the usual domain of the BS applications. However, at large separations, when the interaction between spins is weak, the BS treatment falls into its domain of validity. Even for the simplest species, the binuclears, the existent literature presents several technical options, e.g. by Noodleman (1981), Noodleman and Davidson (1986), Ruiz et al. (1999), and Bencini et al. (1997). The most general is, however, the formula laid down by Yamaguchi and Onishi (Mitani et al. 2000; Onishi et al. 2001): J¼ EHS  ELS ; 2 ^ S HS  ^ S2 LS ð4:22Þ where the HS and BS denote, correspondingly, the High Spin (the state with parallel spins) and Broken Symmetry Spin configurations. The BS state is the unrestricted Slater determinant with localized anti-parallel spins. The formula is very convenient, because the ^ S2 expectation values are printed out by the regular output of the ab initio packages, which give also the E energies. On the other hand, this form is limited to cases with two spin carriers in the molecule, but yet general for any Sa and Sb couple of spins. The HS configuration is close to the state with the S = Sa + Sb spin, although the unrestricted computational realization makes a slight difference, because the ^ S2 HS can differ, usually not much, from the (Sa + Sb)
(Sa + Sb + 1) value. The BS configuration corresponds to a Slater determinant with the Sz = |Sa − Sb| spin projection, but is not the S = |Sa − Sb| spin state, the ^ S2 BS being much different from the |Sa − Sb|
(|Sa − Sb| + 1) quantity. The true S = |Sa − Sb| state has a multi-configuration character, while the BS is a single determinant and EBS is the energy expectation value, as ^S2 BS is the expectation value of the spin squared operator, not its eigenvalue. At the same time, one may observe that Eq. (4.22) is quite general and can be used in circumstances where a pair of true spin states are involved. As sketched previously, the spin states in a binuclear have the ES = –JS(S + 1) energies and the ^S2 S ¼ SðS þ 1Þ expectation values of the spin square operator. It turns out, then, that the Yamaguchi-type formula is valid also applied for the states with maximal and minimal spin quantum number, S = Sa + Sb and |Sa − Sb| (Fig. 4.9). 320 4 Bond! Chemical Bond: Electronic Structure Methods at Work J (Hartree) 0 BS CAS -0.25 R(Å) -0.5 0.5 1 〈S2〉 2 1.5 2 2.5 Sz=1 (HS) 1.5 1 Sz=0 (BS) 0.5 R(Å) 0 0.5 1 1.5 2 2.5 1 σH 0.5 a α b β 0 -0.5 R(Å) -1 0.5 1 1.5 2 2.5 Fig. 4.9 Synopsis of the quantities related with the broken symmetry (BS) treatment for the simplified model of H2 molecule with STO-3G basis set. The bottom panel shows Mulliken spin populations (rH) on the two atomic sites for UHF calculations with Sz = 0, attempting the artificial asymmetry by reading localized orbitals as initial guess. The proper Broken Symmetry (BS) regime starts when the self-consistent orbitals disproportionate into a and b spin polarizations, here, at distances larger than 1 Å. The middle panel shows the S2 expectation values for the unrestricted calculations, performed at higher spin projection (HS) with Sz = 1 and low spin projection, Sz = 0, conducted in the BS regime. One notes that, in the BS, the expectation values depart from the S(S + 1) estimation. The upper panel shows the exchange coupling obtained properly (in the limited model) from CASSCF solutions as J = (ES=1 − ES=0)/2 (continuous line) and estimated from BS treatment (dashed line) with Yamaguchi-type formula. Due to the simplicity of the case, the match between the two types of J modeling is good in the entire domain, even at smaller distance, when the BS regime is not fulfilled. One notes a discontinuity between the two branches of unrestricted J estimation. To legitimate the comparison with the CASSCF, the BS is used in the HF form (not in the DFT often practiced version) For a touch of concreteness, Table 4.5 shows the samples of several input files, one by Gaussian (Frisch et al. 2009), and the others by GAMESS code (Schmidt et al. 1993). Intentionally, the simplest basis set, the STO-3G, was used, to safeguard the simplicity. The primary information contained in the input file concerns the used method and basis set, the geometry of the molecule, its charge and spin multiplicity. As seen on the top of left-side column from Table 4.5, the Gaussian input is quite minimalist, announcing the method (RHF) and the basis set in initial command lines (defined by the preceding #P or # keys). Although optional, the definition of a so-called checkpoint file, the black box of the calculation, is useful for post-computational deals (visualizing the orbitals or restarting the process with 4.4 The H2 Molecule … 321 Table 4.4 Illustration of inputs for the calculation of H2 molecule in Gaussian and GAMESS codes, with different methods (restricted and unrestricted Hartree–Fock) Restricted Ha rtree Fock and Unrestricted Broken -Symmetry tests. A CASSCF simple sample of calculation %chk=H2 $CONTRL scftyp=mcscf #P RHF/STO-3G icharg=0 units=angs mult=1 $END ! choose the algorithm for Title: H2 restricted HF ! generating the configurations $MCSCF cistep=aldet 0 1 H 1 0.000000 0.000000 0.325 H 1 0.000000 0.000000 -0.325 canonc=.t. $END ! define the active space ================================= $CONTRL scftyp=rhf units=angs icharg=0 ! number of core and active ! orbitals, number of electrons ! and the set of states mult=1 $END $BASIS gbasis=STO ngauss=3 $END !(by spin projection sz) $DATA $DET group=c1 ncore=0 Title: H2 restricted HF CN nact=2 nels=2 1 nstate=4 sz=0 H 1 0.000000 0.000000 0.325 H 1 0.000000 0.000000 -0.325 $END $END $BASIS gbasis=STO ngauss=3 $END ================================= $CONTRL scftyp=uhf icharg=0 units=angs $DATA title H2 mult=1 $END $BASIS gbasis=STO ngauss=3 $END CN 1 $DATA title H2 CN 1 H 1 0.000000 0.000000 H 1 0.000000 0.000000 -0.325 0.325 $END H 1 0.000000 0.000000 1.25 H 1 0.000000 0.000000 -1.25 ! asymmetric Broken-Symmetry guess $guess guess=moread norb=2 ! CASSCF needs usually a guess ! Here we provide the non-normalized $END ! a+b and a-b types of combinations $end $VEC $guess guess=moread norb=2 $end 1 1 1.00000000E+00 0.00000000E+00 $VEC 2 1 0.00000000E+00 1.00000000E+00 1 1 1.00000000E-00 1.00000000E-00 1 1 0.00000000E+00 1.00000000E+00 2 1 1.00000000E-00-1.00000000E-00 2 1 1.00000000E+00 0.00000000E+00 $END $END 322 4 Bond! Chemical Bond: Electronic Structure Methods at Work other methods or basis sets). Here, the %chk = H2 line produces an H2.chk checkpoint file. In more complex systems, the definition of demanded memory, the number of parallel processes may imply other specific initialization lines. The input continues with a title line placed between two blank lines. The charge and spin multiplicity (here 0 and 1) are declared before the coordinate lists. Most often, the Cartesian coordinates are used (a list of atom symbols and the x, y, z values), the alternate option being the so-called Z-matrix, a protocol specifying bond lengths, bond angles, and dihedrals. Here, we exemplify the most basic calculation level, the Hartree–Fock on the closed shell ground state configuration. The next section on the left column from Table 4.5 shows the same calculation in GAMESS format. The text is more verbose, apparently more complicated, in comparison to the Gaussian input. The simplicity of Gaussian is one of the points of appeal that led to the popularity of this program. At the same time, in the lengthier format of GAMESS, it is easier to program supplementary or advanced options, an action which becomes more cryptic in Gaussian case, using the keywords called IOP, comprising numeric coding, for various operations. In turn, the GAMESS further input options are text keys, resembling acronyms or abbreviations of the desired procedure. The text of GAMESS input for the RHF on H2 molecule is quite self-explanatory, noting here only the less visible detail, that this calculation is made without symmetry, this being coded in the CN 1 line of the $DATA section (meaning a Cn point group with trivial n = 1 first-order axis, equivalent to a NoSymm keyword in Gaussian case). When symmetry exists, it is sufficient to put only the unique atomic species, the other ones being retrieved according to the specified point group. The GAMESS has a text-type checkpoint file, usually named with the “dat” extension, where details of the calculation are punched out, and can be reused in restarts and analyses. The most important content is the $VEC section, which stores the formatted block of LCAO (linear combination of atomic orbitals) matrix. As is illustrated in the other examples, this block can be inserted in the input file, for a restart or for educated guess, after appropriate changes. The unrestricted HF example (UHF), in the last text module from the left column of the discussed table, has the somewhat more complicated situation of Broken Symmetry (BS) calculation. First, to have clear success with this type of calculation, a large separation of atoms is induced, as specific to the BS regime. Then, an orbital guess pointing toward the desired spin asymmetry should be introduced, with the $VEC group requested by the $GUESS line. The $VEC group is formatted, with a certain pattern of indices and fixed positions for digit sequences. In our simplified cases, we have two MOs and two AOs in the $VEC group, namely two lines with two coefficients, the blocks being doubled in the case of unrestricted calculation, with the first half referring to the a orbitals, while the second refers to the b ones. In the BS calculation, we aim to put one a electron at one H atom and one b electron on the other site. Having one orbital per site, the first a MO is formally introduced as {1.0, 0.0} while the first b MO is proposed like {0.0, 1.0}. The second MOs, in both sets, are actually not important, being virtuals. In the exemplified input, such coefficients are ignoring the normalization, which involves also the role of overlap between AOs. However, given the large separation between 4.4 The H2 Molecule … 323 sites, the overlap will be small and the automatic normalization, performed at the guess reading, is sufficient for the correct course of the convergence. Otherwise, at smaller distances, or in more complicated systems, a carefully managed set of orbitals should be introduced as guess, e.g. prepared by the localization, or the back-symmetry transformation, of the MOs from a restricted calculation. The right side of Table 4.5 illustrates the CASSCF input, worked with the more general multi-configuration self-consistent field keyword of GAMESS, MCSCF, which, in combination with the given “cistep” option (to be read like: configuration interaction step) is equivalent to complete active space procedure. The active space is defined in the “$DET” block (dedicated to the generation of determinants): number of electrons, orbitals, and states. The CASSCF should be ideally started with proper tailored orbitals, the operations possibly implying reordering of the orbital set (so that the active set is picked from the frontier sequence), or other previous orbital transformations (localization or “handmade” multiplication with unitary matrices). In the exemplified case, the guess can be the unaltered output of a previous restricted HF calculation. We even exaggerated, putting in the format of $VEC not-normalized orbitals, namely the {1.0, 1.0} and {1.0, −1.0} pairs of coefficients, instead of proper LCAOs, since, in the relatively simple situation, the automatic post-guess normalization is restoring well the matrix. 4.5 4.5.1 Computational and Conceptual Valence Bond: The Spin Coupling Paradigm at Work Overture to the Valence Bond Calculations The Valence Bond (VB) theory, crowned with the renown of the first quantum theory of molecules, exists nowadays in revisited forms, conceived in various technical ways. Modern VB avatars can be devised borrowing the technical body from the Complete Active Space (CAS) self-consistent field (SCF) procedures, introducing the orbital optimization, together with the use of many-electron bases, formatted in the style of resonance structures, arriving at the Complete Active Space Self-Consistent Valance Bond (CASVB) methods (Raimondi and Cooper 1999; Cooper et al. 1999; Thorsteinsson and Cooper 1998; Hirao et al. 1996; Nakano et al. 2002). There are several VB programs, standalone or devised as routines inside larger computation packages, confining ourselves here to the mention of the VB2000 (Li et al. 2007; Li and McWeeny 2002) and XMVB (Song et al. 2005, 2012) codes. The VB is less used, in comparison to other methods of choice, but it cannot be judged as less deserving, in conceptual or technical respects. Sed contra, it deserves enhanced attention, because of its virtues for communication in the words of chemical language. 324 4 Bond! Chemical Bond: Electronic Structure Methods at Work 11 12 12 12 11 1 11 2 2 1 2 4 4 3 10 4 3 3 10 9 5 9 7 8 5 6 7 6 8 6 5 7 8 Fig. 4.10 The superposed self-consistent VB orbitals for the methanol molecule, shown from different orientations and labeled in line with the {1,1}, {2,2}, {3,4}, {5,6}, {7,8}, {9,10}, {11,12} Rumer diagram For a brief entry into the specifics of the Valence Bond method, namely the use of localized orbitals, the methanol molecule is taken as first example. Figure 4.10 shows the superposed self-consistent VB orbitals, their labeling being related with the following imposed single Rumer configuration: 1 2 3 5 7 9 11 1 2 4 6 8 10 12 which can be ascribed in condensed form like {1,1}, {2,2}, {3,4}, {5,6}, {7,8}, {9,10}, {11,12}. The input of this case is shown in Table 4.4, handled with the VB2000 module, as implemented in the GAMESS code. The problem would imply 14 active electrons, with a VB(14) command. However, practically equivalently, it can be made as product of two VB subspaces, one with four electrons for the two lone pairs, and another with ten electrons in the five bonds, produced by the VB(4). The VB(10) keyword in the input is illustrated in the following. These two orbital sets are labeled as 2 and 3 (this is why the 02 and 03 prefixes in the VBSTR keywords and the =>2 or =>3 suffixes in the valence bond group assignment are delineated by the VBGA key), while the group no. 1 is tacitly assigned to the core electrons. Then, the $02VBSTR and $03VBSTR keys are corresponding to the previous Rumer structure, rewritten as the merging of the two VB subspaces with the electrons renumbered within each group, {1,1}, {2,2}, for the lone pairs and {3,4}, {5,6}, {7,8}, {9,10}, for the bonds. One may note the two lines of the $02VBORB block as the in-phase and out-of-phase combinations (normalized sum and differences) of two primary lone pairs on oxygen (atom no. 2), which are produced primarily in accordance with the symmetry of the molecule with respect to the C–O–H reflection plane, one function 4.5 Computational and Conceptual Valence … 325 being placed in this plane and another perpendicular to it, in a pattern resembling the sp2 situation. As discussed in the section dedicated to hybrid orbitals convention in the case of the water molecule, the recombination of symmetry adapted localized doubly occupied orbitals, identified in the input by the colon symbols:(1) and:(2), retrieves the “rabbit ears” sp3-type equivalent lone pairs. The 1 1 2 2 sequence in $02VBSTR marks the situation of pairing electrons in these lodges. The $03VBORB block defines the spin-coupled orbitals, the 1 -> 2 and 2 -> 1 denoting, for instance, the hybrids pointing from atom 1 to atom 2 and vice versa, their overlap making the C–O bond, based on the assignment established in the first sequence, 1 2, of the $03VBSTR Rumer table. In principle, everything interacts with everything, but it is reasonable to admit that the overlapping orbital pairs contain the maximum of bonding effects. Thus, the VB is the way for giving a quantum account for ideas of additive bond energies. 4.5.2 Benzene: Valence Bond Versus Complete Active Space Entering deeper into the Valence Bond craft, the benzene molecule is taken, which, with its iconic connotation to the aromaticity and resonance effects, is the perfect illustration for the VB paradigm. In Table 4.6, the GAMESS-VB2000 input for the Spin Coupled Valence Bond (SCVB) treatment of the C6H6 is shown, the results being illustrated in Fig. 4.11. In the input part dedicated to the VB commands, one may recognize the lines corresponding to the five resonance structures described by Rumer graphic diagrams, or by two-row notations resembling Young tableaus. To be distinguished from the rules of the Young representations, here the tables are allowed to correspond directly to the Rumer graphical representations. Having the double bond symbol for spin coupled terms, the Rumer diagrams stand for the two celebrated Kekulé and three Dewar structures.. The VBORB section of the input specifies, by the 1^–6^ entries, the intention to obtain p-type orbitals, the VB2000 code having this facility, specially prepared for the case of planar conjugated systems. The pictures added to the right side of Fig. 4.11 are illustrating the pz-alike nature of the converged VB orbitals. The intersecting contours of the superposed localized functions are suggesting their non-orthogonality, a specific of the VB procedures, as previously noticed. The values of overlap integrals for the respective ortho, meta, and para mutual positions are: 0.525, 0.028, and −0.159. The negative value for the para couple is due to the fact that the main lobes show diffuse outskirts, with opposed sign. This distant overlap substantiates the non-trivial role of Dewar resonances. The poly-electronic basis functions of Rumer type are non-orthogonal, too. In the discussed computational setting, the overlap between the two Kekulé structures has the absolute value of 0.370, the sign depending on the order taken in the columns of Young-type tableaus for the coupled spins (e.g. the swap of 1 and 2 in the first column of first structure switches the sign). The magnitude of the overlap of one Kekulé and one Dewar structure is about 0.616, larger than the inter-Kekulé 326 4 Bond! Chemical Bond: Electronic Structure Methods at Work Table 4.6 VB calculation of benzene by the VB2000 routine implemented in GAMESS code GAMESS preamble ! VB calculation on benzene $contrl scftyp=rhf units=angs icharg=0 mult=1 ispher=1 VBtyp=VB2000 $end VB2000 calculations $VB2000 #! VB(4.86) PRINTALL NoSym=1 $basis gbasis=N311 ngauss=6 diffsp=.t. ndfunc=1 $end $DATA title DNH 6 C 6 1.394783 0.000000 0.000000 H 1 2.480313 0.000000 0.000000 $END ! Using symmetry (D6h point group, defined in the second $DATA line), the geometry input contains only the unique species of atoms. However, the symmetry cannot be used in VB, being deactivated in the MO part by NoSym keyword. $02VBORB 1^ 2^ 3^ 4^ 5^ 6^ $02VBSTR 5 1 2 3 4 5 2 3 4 5 6 1 4 2 3 6 1 6 2 5 3 2 1 3 6 4 6 1 5 4 5 $02ROOT 1 $END The Rumer resonance structures resonance structures in VBSTR block correspond to those represented in Fig. 4.11 amount, apparently counter-intuitively, because of topological reasons debated later on. The overlap between two Dewar-type basis components is 0.419. Because of the large absolute values in the poly-electronic overlap matrix, the coefficients displayed in the table contained in the synopsis from Fig. 4.11 are deviating from the sub-unitary values, customarily appearing for the CI in orthogonal bases. In other words, the sum of square of coefficients is not retrieving the unity, the normalization implying also sums over triple products of two coefficients and one overlap, each. We must recall that the taking of resonance structures follows an algorithm producing the needed number of basis components for a given spin state. The five resonances of the benzene are then covering the five dimension of the singlet states for the case of a system with six electrons. Thus, the Kekulé and Dewar patterns are selected on algebraic grounds, not on the inspiration of chemical intuition, though the procedure is in line with this quality too. The philosophy of Young tables was explained previously. The Rumer graphical algorithm consists in putting the spin carriers on a circle and drawing pairs of lines among sites, in a way avoiding their 4.5 Computational and Conceptual Valence … E (cm-1) 0 327 42587.8 74235.0 74235.0 120381.1 1 6 2 5 3 4 1 2 6 3 5 4 1 6 2 5 3 4 1 2 3 4 5 6 0.442 1.239 0.000 0.000 4.311 2 3 4 5 6 1 0.442 -1.239 0.000 0.000 4.311 1 4 2 3 6 5 -0.070 0.000 0.000 1.658 2.886 1 6 2 5 3 4 -0.070 0.000 1.436 -0.829 2.886 2 1 3 6 4 5 -0.070 0.000 -1.436 -0.829 2.886 1 2 6 3 5 4 1 6 2 5 3 4 Fig. 4.11 Synopsis of the VB calculations on benzene. Left side: graphical and table-type representation of the Rumer resonance structures. Middle panel: coefficients of resonance structures in the five singlet states, given as column eigenvectors, with corresponding eigenvalues on the top line. Right side: superposed contours of the self-consistent localized p-type orbitals, taken at an elevation of 0.5 Å above the molecule plane (upper-right corner) and (lower-right corner) the same thing, from a different viewpoint, adding a 3D representation of one VB orbital (at atom 1), drawn at 0.1 electrons/Å3 isosurface crossing. The non-crossing rule is another strategy (somewhat similar to ordering criteria for standard Young tableaus) to accomplish the necessary count of basis components suited for a given spin. Structures with intersecting lines are linear combinations of non-crossing ones. In the case of non-singlet states, the resonance diagrams contain “free radical” sites, namely products of a spins, represented as unconnected dots. For benzene and other small members of the CnHn series, the Rumer ring coincides with the real structure, but the circular topology should also be used for the complete construction of the spin bases in general situations, irrespective of the molecular geometry. Besides, the sites may represent different orbitals, not atoms only, in the more complex cases the graphic representations going beyond the immediate chemical intuitive grasp. Thus, while for the real skeleton of naphthalene only three Kekulé structures can be drawn, the system of ten p electron forming singlet states demands a basis of 42 resonances. The graphical approach is reasonable for few spins, while in larger systems other equivalent ways of coding should be considered in practical instances. The non-crossing graphic rule is limited to singlet states. For instance, in the case of VB triplets of benzene should be nine states. From each of the singlet Rumer structures there is, in principle, an offspring of three triplet-type resonances, decoupling the three double bonds in bi-radicals, successively, all these respecting the non-crossing pattern. However, in this way one obtains a count of 15 independent graphs, more than the 9 dimension of the triplet set. Therefore, in non-null spin cases we must rely on Young tableaus. On the other hand, one may still use the graphical approach limiting the number of non-crossing structures to the necessary dimension, with the cost of inducing a 328 4 Bond! Chemical Bond: Electronic Structure Methods at Work certain non-symmetric aspect of the basis. As another example, consider the quintet states on benzene, for which one may imagine the traveling of a “double bond” over the six edges, the other sites being radicals, while the needed dimension of the set is five. Then, one may arbitrarily eliminate one of the six resonances suggested by the primary chemical intuition. An alternative procedure is to consider imaginary dots to which the unpaired radicals are making fake links like coupled spins. In this case, one can apply again the non-crossing rules taking out, at the end, the resonances with links between the virtual points. For instance, returning to triplet benzene, with two virtual aid points, it can be treated as an eight particle system, for whose singlets there are 8!/(4!4!)–8!/(5!3!) = 14 states. However, the resonances with an 8–9 link are non-physical and should be silenced. There are five such structures (like for singlet benzene itself) which, subtracted from the set with the 14 dimension, one arrives at the correct 9-dimension of the triplet in a six-electron system with unique orbital configuration. The columns from the table displayed in the middle of the synopsis from Fig. 4.11 give the coefficients of resonance structures in the different singlet states. The first column describes the ground state, as the superposition with equal coefficients of the two Kekulé components and a small portion of Dewar ones, each with the same contribution. The first VB excited state is made from the out-of-phase combination of Kekulé structures (coefficients with opposed sign and equal magnitudes), without Dewar components, ruled out for symmetry reasons. The next two columns correspond to a doubly degenerate state, their coefficients being mutually transformable by an arbitrary rotation. Only the Dewar structures are participating to this orbital doublet. The last state is totally symmetric, like the ground state, with the same resonance actors, but with higher weight of Dewar components. As pointed out previously, the coefficients are supra-unitary, because of the overlap integrals between the resonance structures. It is more intuitive to judge the corresponding weights, which are, in the last state, 0.95% for the Kekulé pieces and 27% for each Dewar-type one. The following discussion will put in parallel the Complete Active Space Self-Consistent Field (CASSCF) calculations with the Valence Bond (VB) paradigm and the phenomenological spin Hamiltonian called Heisenberg– Dirac–van Vleck (HDvV) (Heisenberg 1928; van Vleck and Sherman 1935; Anderson 1959). The HDvV model is often used to fit the magnetic susceptibilities in transition metal complexes. It can produce a spectrum of spin states, object to the Boltzmann statistics over the spin magnetic moments, in simulating properties like magnetization or susceptibility, fed with parameters assigned to couples of interacting centers. However, the HDvV Hamiltonian, having a phenomenology isomorphic to the Spin Coupled Valence Bond method (SCVB, when a single orbital configuration is defined as a factor of spin functions), is a much more powerful tool than the currently perceived credits (diverted, because of long routine use in the molecular magnetism). It can be applied to several prototype cases of chemical bonding, such as the states resulting from the p system of benzene and, furthermore, to the problem of aromaticity (see Chap. 5). While, in magnetism, the spin 4.5 Computational and Conceptual Valence … 329 Hamiltonian may concern any spin state of paramagnetic sites, in VB one deals only with the 1/2 quantum numbers of electrons. Taking the CASSCF as calculation scheme, in a CAS(6,6) setting, aiming to account six electrons in six p orbitals, we must perform a preliminary selection of the relevant orbitals. Guided by symmetry reasons, one knows that the following set of representations in the D6h point group is needed: a2u, e1g, e2u, b2g (namely six functions, considering that the e labels represent double degeneracy). The intuition developed alongside the Hückel simple model seeds the expectation of finding the p-type orbitals at the frontier (three occupied MOs and three virtuals). Such regularity is not retrieved at the use of rich atomic bases, augmented with many diffuse and polarization parts. The non-p frontier components may give rise to states with physical relevance, such as Rydberg levels, which however are ignored in the defined quest. Using the 6-311+G* basis, the Hartree–Fock step needed in the preamble of CASSCF, as primary source of orbitals, produces the desired p-type functions at the following positions: {17, 20, 21, 27, 28, 53}, not quite in line with respect of the occupied–unoccupied frontier, located between the 21–22 couple. Collecting these orbitals and reordering them on the 19–24 successive positions, the CASSCF calculations can be started, the results being displayed in Tables 4.7, 4.8 and 4.9. Table 4.7 Singlet terms from CASSCF(6,6)/6-311+G* calculation on C6H6 and the identification of the HDvV-type states Order No. E CASSCF (cm−1) 1 3 0 43,426.53 E HDVV fit (cm−1) 0.00 41815.25 9 71,202.43 73912.29 10 71,202.43 73912.29 11 13 14 19 22 73,167.69 83,765.98 83,765.98 102,301.1 106,839.1 115,727.54 HDVV 0 pffiffiffiffiffi   13  1 J pffiffiffiffiffi   13 þ 1 J pffiffiffiffiffi   13 þ 1 J pffiffiffiffiffi 2 13J 25 111,496.2 26 111,496.2 27 113,132.9 30 124,246.6 31 124,246.6 40 138,251.2 41 138,251.2 42 139,509.4 The order number is related to the placement of the given state in the whole energy spectrum (from singlets to septets) 330 4 Bond! Chemical Bond: Electronic Structure Methods at Work Table 4.8 Triplet terms from CASSCF(6,6)/6-311+G* calculation on C6H6 and the identification of the HDvV-type states Order no. E CASSCF (cm−1) E HDVV fit (cm−1) 2 33,620.12 21,978.18 4 43,876.19 48,851.68 5 43,876.19 48,851.68 6 62,092.95 57,863.77 7 62,092.95 57,863.77 8 17 66,380.17 94,549.94 93,749.36 18 96,821.38 106,009.34 20 21 23 103,944.3 103,944.3 111,018.6 115,021.43 24 111,018.6 115,021.43 HDVV pffiffiffiffiffi pffiffiffi  13  5 J pffiffiffiffiffi  pffiffiffiffiffi  13  3=2 þ 17=2 J pffiffiffiffiffi  pffiffiffiffiffi  13  3=2 þ 17=2 J pffiffiffiffiffi  13J pffiffiffiffiffi  13J pffiffiffiffiffi pffiffiffi  13 þ 5 J pffiffiffiffiffi   13 þ 3 J pffiffiffiffiffi  pffiffiffiffiffi 13 þ 3=2 þ 17=2 J pffiffiffiffiffi  pffiffiffiffiffi  13 þ 3=2 þ 17=2 J  28 118,611.7 29 118,611.7 34 128,157.6 35 128,157.6 36 128,715.1 37 129,304.4 38 129,304.4 39 130,839.0 43 142,299.5 44 142,299.5 47 144,782.6 The order number is related to the placement of the given state in the whole energy spectrum (from singlets to septets) The CAS(6, 6) calculations involve a large number of states: 175 spin singlets, 189 triplets, 35 quintets, and one septet (in total, 400 levels, including the orbital degeneracies). Then, a non-trivial problem lies in identifying the states following an HDvV phenomenology, solving implicitly the question whether there is such an effective behavior. In the HDvV and SCVB frame, we should identify five singlets, nine triplets, five quintets, and one septet. Without entering into details, the clue to selecting the HDvV-type levels from the many CASSCF ones lies in repeating the CAS calculation, after an orbital transformation on canonical MOs, bringing them to localized functions, looking like pz components on carbon atoms. To be distinguished from the above discussed genuine VB calculations, these pz alike functions are orthogonal each to other, though well localized. This is because they are obtained by the unitary transformation from canonical MOs, keeping then the orthogonality feature. 4.5 Computational and Conceptual Valence … 331 Table 4.9 Quintet and septet terms from CASSCF(6,6)/6-311+G* calculation on C6H6 and the identification of the HDvV-type states. The order number is related to the placement of the given state in the whole energy spectrum (from singlets to septets) Order no. 12 E CASSCF (cm−1) 75,609.63 E HDVV fit (cm−1) 73,912.29 15 93,321.19 89,960.81 16 93,321.19 89,960.81 32 125,457.7 122,057.86 33 125,457.7 122,057.86 45 46 48 143,604.7 143,604.7 147,971.5 138,106.38 HDVV pffiffiffiffiffi  13 þ 1 J pffiffiffiffiffi   13 þ 2 J pffiffiffiffiffi   13 þ 2 J pffiffiffiffiffi   13 þ 4 J pffiffiffiffiffi   13 þ 4 J   pffiffiffiffiffi  13 þ 5 J A unitary transformation of the canonical MOs does not change the computed energy levels, but is based on other pattern of many-electron Slater determinants, which are relevant for the sake of comparison to the considered phenomenological model. Then, one may assign to the HDvV type the states based, as much as possible, on Slater determinants with one electron per local orbital, running only in their spin flips. The states with spin coupling pattern are found among the first 50 levels of the whole spectrum, comprising all the spin multiplicities. Focusing on singlet energies (see Table 4.7), one notes the comparability, though not exact match, of the values resulted from the VB calculations, namely the series {0, 42587.80, 74235.02, 74235.02, 120381.11} (all values in cm−1) against the CASSCF output for the VB-type levels: {0, 43426.53, 71202.43, 71202.43, 106839.10} (in cm−1). One may observe, under this note, that the VB, in its spin coupled version (SCVB), limited to a single orbital scheme, is different from the CAS conceptual and practical approach. However, a full CASVB, allowing all the orbital configurations, would be, at the end, the same with CASSCF, in energy spectrum, although it may imply different orbitals and therefore another system of eigenvectors. Another observation is the quite good closeness of the CASSCF states identified as VB type to the HDvV phenomenology. Here, the simplest parameterization is by a unique J, effectively describing the interaction between two linked carbon atoms. It is rather remarkable that a single parameter, estimated at about J = −16048.5 cm−1, fits reasonably the whole series of CASSCF states having HDvV nature (compare the values marked in bold, nearby the corresponding analytic expressions, in Tables 4.7, 4.8 and 4.9). The fact that not all the states are accounted by the spin coupling phenomenology illustrates the limitations of a standard Valence Bond perspective, while the relatively good match of the corresponding states to the simple model proves the reliability of the VB concepts, in their domain of validity. Another global test on states identified as HDvV is the relative ratio of the barycenters, in the series from singlets to septet, found 1: 3.01: 6.09 for the selected CASSCF states, very close to ideal HDvV pattern: 1: 3: 6. 332 4 Bond! Chemical Bond: Electronic Structure Methods at Work 4.5.3 Playing with Graphic Rules for Setting a VB Modeling Combining the VB and HDvV perspectives, advanced detailing of the benzene example will help the deeper insight into the phenomenology of spin coupling, including the issue of non-orthogonality of the resonance structures. Of particular beauty is the fact that, after the Rumer graphical construction of the basis, there are simple rules to find the overlap integrals and coefficients building the Hamiltonian from exchange integrals, by interpreting the patterns resulting from the superposition of resonance structures. This may become clear inspecting Fig. 4.12. Observe first the resonance structures, this time developed with arrows instead of the double bounds depicted in the previous Fig. 4.11. This serves to define a sign of the wave function, the tail and tip of the arrow between the i and j sites representing, respectively, the first and second member in the (a(i)b(j) − b(i)a(j)) parenthesis, acting as factor in a given resonance wave function. Or, in other words, the i ! j arrow corresponds to the i site standing in the first row of a Young-type representation, while j is placed, below it, on the second line. Writing on horizontal and vertical entries the Rumer structures, signed with arbitrary defined arrow directions, one produces a table with superposition patterns (see left panel) prefiguring the overlap and Hamiltonian matrices. The right panel illustrates the further handling, which brings the superposition figures to the so-called matching situation, with the arrows meeting either tail-to-tail or head-to-head. This condition ensures that there exist common products of spin states in the “bra” and in the “ket” sides. The lines confluent at a given point describe different partners, one from rows and one from columns, because in a 1 1 6 2 6 2 5 3 5 3 4 4 Fig. 4.12 The superposition of resonance structures for the singlet states of the benzene. The left side contains the direct patterns, the right panel showing the modification to the matching form, with all arrows head-to-head or tail-to-tail. The reverted arrows are marked in dark ink, their count providing the sign of overlap integrals: positive or negative, for even or odd number of switches, respectively 4.5 Computational and Conceptual Valence … 333 given resonance structure a site, i, is used only once, namely in a single arrow. The number of arrow reversals, mPQ, bringing the superposition pattern of P and Q resonance structures to the standard match, will decide the sign of the matrix components, by a factor ð1ÞmPQ . If the orbitals making the space-dependence background of the spin coupled wave functions are orthogonal, then a very simple rule for absolute value of the overlap SPQ can be proposed. In spite of the fact that the VB is based on the idea of overlap as driving force, the orthogonality condition is not so restrictive, since, as discussed previously in comparing the genuine VB and the CASSCF results, one may emulate a VB regime using orthogonal basis too. The key feature determining the magnitude of the overlap integrals is the number of so-called islands, namely closed circuits made at the superposition of P and Q graphs, the smallest element of this sort being the stacking of one arrow with itself. Thus, if nPQ is the number of islands formed at the superposition of the P and Q Rumer diagrams, and g is the count of electron pairs of the system (contained in each resonance structure), the absolute value of their overlap is 2nPQ g . A hint for this formula is first caught in the fact that the expansion of a resonance structure with g electron pairs contains 2g micro-configurations, consisting in spin flips. Each product of elementary spin functions is orthogonal to the other primitive elements of this type. The overlap of the resonance function with itself brings 2g situations where each micro-configuration from “bra” meets its own image in “ket”, so that pffiffiffiffiffi the normalization factor of the expanded product of paired functions is 1= 2g . Thus, the product of normalization factors from each resonance structure brings the 2−g factor in a matrix element of overlap or Hamiltonian. A resonance structure may contain an amount of unpaired electrons, acting as supplementary factor to the unfolded product of (a(i)b(j) − b(i)a(j)) parentheses, this composition not affecting the above described count of expanded wave function terms and the normalization factor. Now, we shall understand the formation of an island as a consequence of finding “bra” and “ket” sequences with the same content of primitive micro-configurations. A closed polygon, with h edges and vertices in the superposition figure, means that the encircled collection of h centers is found in a given factorization of the matrix elements in both “bra” and “ket” sides. Besides, the h should be even, because only in this way may one achieve the perfect dichotomy of lines concurring to one point: one from the left matrix side, one from the right matrix component. Considering, for instance, a i–j–k odd-numbered sequence of a superposition island then it is clear that if i–j belonged to “bra”, the j–k line must be related to the “ket”, since the electrons of j were already used in the left component, so that this cannot take the j–k line. The matching condition imposes that the spin on j coincides in “bra” and “ket”. For each island there are only two spin products found both in “bra” and in “ket”, namely those ensuring the complete spin alternation along the contour, a1b2a3b4


ah−1bh and b1a2b3a4


bh−1ah. Then, each island contributes with a factor of 2 to the total overlap, arriving to the 2nPQ g absolute value. In the case of resonances containing r unpaired electrons, prepared for states with r + 1 spin multiplicity, the superposition patterns may form open chains, 334 4 Bond! Chemical Bond: Electronic Structure Methods at Work incorporating the radical sites. These can comprise an even or odd number of sites, situations labeled E or O. The radicals are forming the ends of the superposition strings. An open spin in a resonance structure is linked to nothing. Then, a link touching this site in the superposition pattern must come from the partner resonance structure. Then, from the side of this companion, the concerned site is not involved in the other link, as it was not from the side of the former, the radicals acting as stoppers in the growth of superposition patterns. In the superposition pattern, a radical site, conventionally fixed to a spin can receive only the tail of the arrow. Then, in an O-chain, having an even number of bonds, the number of spin products common to “bra” and “ket” is restricted to one, it being—by convention—not possible to switch the spin on the radicals. For an E-chain, there is a topological impossibility to form matching arrows. As an example, take the r–i–j–s sequence, with r and s open radicals. Let us start with the tail from r, respecting its a nature, having the following tentative matching sequence: r ! i ← j ! s. The last arrow offends the assumed a spin projection on the site s, representing a non-physical situation. It is suggested, in this way, that in E-chains it is impossible to find identical micro-configurations in the corresponding sequences of “bra” and “ket”, leading to the vanishing of the overlap. The O-chains are allowing the formation of overlaps, but are not visible in the final formula, contributing only with unity factors. Ascribing by dE = 0 the quenching due to apparition of at least one E-chain and by dE = 1 the absence of these patterns, the overlap of the P and Q resonance structures is SPQ ¼ dE ð1ÞmPQ 2nPQ g , with mPQ the number of arrow reversals needed to bring the picture to matching form, nPQ the number of island, and g the total number of coupled electron pairs in both diagrams (excluding the r unpaired electrons). As illustration of the described rules, one may inspect the first line in Fig. 4.12, namely the overlap of one Kekulé structure with all the other resonances. The left side panel shows the first instance superposition patterns, while the other is arranged to the matching form, marking the reverted arrows by dark coloring. Reading the {mPQ,nPQ} indices along the upper line of the right side panel, one finds: {0,3}, {1,0}, {2,2}, {1,2},{2,2}. One notes the first couple of indices, with three minimal islands created from the superposition of each arrow on itself. Applying the above formula, with g = 3, one obtains the respective series of overlap integrals: 1, 1/4, 1/2, −1/2, 1/2. This is an important step toward the phenomenological modeling of benzene by VB and a valuable example for other cases. The reach of Hamiltonian matrix elements is based also on graph-type rules. In the spirit of HDvV modeling, the 2Jij ^ Si
^ Sj term of the Hamiltonian is equivalent ^ to the Jij Pij permutation operator. Then, it acts on a resonance structure produced permuting the tail and apex points of the arrow. Thus, applied to a i ! j arrow, this operation simply changes the sign of the resonance structure. On k ← i!j or k ← i!j ← l superposition sequences (with i and j from “ket” while k and l from “bra”), it leads to k ← j←i and k ← j←i ← l patterns, which are taking a skewed aspect on the graphical representation, if we keep fixed the geometrical positions of 4.5 Computational and Conceptual Valence … 335 the i and j points. It results then that the parametric expression of the matrix elements can be expanded estimating the overlap integrals taking the P resonances ^ ij XQ element from “bra”, while the “ket” Q is mutated by permutations, the XP jP representing the coefficient of the Jij effective exchange parameter in the HPQ Hamiltonian matrix element: HPQ ¼ hXP j X j\i   X  1 ^ ij XQ : Si
^ Jij XP jP Jij   2^ Sj  X Q ¼  2 j\i ð4:23Þ These elements, altogether with the overlap: SPQ ¼ XP jXQ ; ð4:24Þ are making the object of a generalized eigenvalue problem: HC = SCE. Figure 4.13 shows the superposition patterns, the primary graphs and their adjusted matching form, between benzene Rumer diagrams and the series of mutated resonances, with the connections of 1 and 2 sites swapped. Applying to these diagrams the rules described for the overlap integrals, one obtains the coefficient of J12 integral in the corresponding Hamiltonian matrix elements, similar procedures going on for the handling the factors of all the Jij parameters. For instance, the intersection of the first row with the first column has an element which is, practically, the superposition of a Kekulé structure with itself, except a sign ^ 12 permutation. Then, change, due to reversal of the arrow induced in “ket” by the P 1 1 6 2 6 5 3 5 4 2 3 4 Fig. 4.13 The superposition of resonance structures with the elements from “ket” (upper row ^ 12 permutation. The overlaps resultingd from matching structures (right entries) operated by the P side panel) are yielding the coefficients of the J12 effective exchange parameter in the HDvV ^ ij modeling of SCVB. The coefficients of other Jij parameters are obtained in similar manner, by P mutations 336 4 Bond! Chemical Bond: Electronic Structure Methods at Work the H11 matrix element has a J12 term, after the sign reversal factorizing the operator (aside those provoked by the operator itself to the superposition pattern). The second element from the first line of superpositions has one island (although skewed, topologically there is a single contour), with two arrows switched to the matching form (positive overlap, therefore), leading to the 1/4 overlap factor and finally to the −J12/4 contribution into the H12 entry. The H13 matrix element implies a pattern with two islands (one with skewed aspect) and three arrow reversals, leading to the −1/2 overlap and the J12/2 content. The following H14 and H15 elements show also two islands, the former with no need for arrow operations, the other with one switch, so that these receive the −J12/2 and J12/2 quantities, respectively. The algorithm can be continued on the other matrix elements and for the other Jij factors. It is also general for the situation of resonances incorporating radical sites, needed for non-singlet spin states, where the factors are also dictated by the number of the formed islands, censored to extinction if any E-type chain occurs. The graph-type rationalization can be distilled to even subtler regularities (McWeeny 2001). Thus, the coefficient of Jij integral in the HPQ matrix element in the HDvV or, equivalently, SCVB, modeling is: (a) SPQ if the i and j sites belong to the same substructure element (island or O chain) separated by odd numbers of links in the superposition pattern; (b) −2SPQ if, in similar placement, there is an even number of links intercalated; (c) −SPQ/2 if the sites belong to different islands or chains; (d) −SPQ if i and j are disjoint radical sites. In a shortcut of reasoning, the factor for case (a) results because the action of the operator permuting two sites in “ket” resonance does not change the number of islands. It provokes only arrow reversal, whose negative signature is counterbalance by the negative factor in front of the exchange parameter, as illustrated before. Therefore, the coefficient of the given exchange coupling remains the overlap itself. In case (b), the permutation causes the split of the island containing the coupling, in two contours, increasing the initial nPQ index to nPQ + 1 and therefore the mutated overlap by a factor of 2, involving also a sign change in the rearrangement. The one half factor in case (c) is because, by contrary, the permutation merges two cells, decreasing the number of islands to nPQ − 1. The (d) situation can be understood by the fact that the independent radical sites, not integrated in chains, in the overlap topology, can be taken as factorized dimer couples that show the −Jij energy, as companion of the Jij singlet couple. In the situation of overlap vanishing because of E-chain formations, there are still contributions to Hamiltonian matrix elements, because permutations between two E-chains may convert them to an O-chain. So, if there are no other Echains able to impose a quenching, the overlap of permuted pattern determines the contributions of all the coupling parameters between the sites linking the two Echains. The absolute value of the coefficient is, similar to previous definition, 2nPQ g for all the elements. The sign is negative if imply head-to-head or tail-to-tail couples and positive, oppositely. Finally, recall that in these formalisms one deals only with exchange integrals, since all the other quantities (Coulomb and kinetic terms) are the same, under all the spin-flip configurations, therefore an equal shift to all the HDvV-type levels. 4.5 Computational and Conceptual Valence … 337 For a consolidation of the enounced rules, the elements of the 55 matrix for singlet states of benzene are given in the following, assuming as non-vanishing only the exchange parameters along the molecular polygon (keeping their indices distinct) and evidencing the corresponding overlap as the factors in the front of each right-side member (for non-diagonal components): 1 1 1 H11 ¼ J12  J23 þ J34  J45 þ J56  J16 ; 2 2 2 1 H12 ¼ ðJ12 þ J23 þ J34 þ J45 þ J56 þ J16 Þ; 4   1 1 1 J12 þ J23 þ J34  J45 þ J56  J16 ; H13 ¼ 2 2 2   1 1 1 J12  J23 þ J34  J45 þ J56 þ J16 ; H14 ¼  2 2 2   1 1 1 J12  J23 þ J34 þ J45 þ J56  J16 ; H15 ¼ 2 2 2 1 1 1 H22 ¼  J12 þ J23  J34 þ J45  J56 þ J16 ; 2 2 2   1 1 1  J12 þ J23  J34 þ J45 þ J56 þ J16 ; H23 ¼ 2 2 2   1 1 1 H24 ¼ J12  J23  J34  J45 þ J56  J16 ; 2 2 2   1 1 1 J12 þ J23  J34 þ J45  J56 þ J16 ; H25 ¼ 2 2 2 ð4:25Þ ð4:26Þ ð4:27Þ ð4:28Þ ð4:29Þ ð4:30Þ ð4:31Þ ð4:32Þ ð4:33Þ 1 1 1 1 H33 ¼  J12 þ J23  J34  J45 þ J56  J16 ; 2 2 2 2 ð4:34Þ 1 H34 ¼  ðJ12 þ J23 þ J34 þ J45 þ J56 þ J16 Þ; 4 ð4:35Þ 1 H35 ¼ ðJ12 þ J23 þ J34 þ J45 þ J56 þ J16 Þ; 4 ð4:36Þ 1 1 1 1 H44 ¼  J12  J23 þ J34  J45  J56 þ J16 ; 2 2 2 2 ð4:37Þ 338 4 Bond! Chemical Bond: Electronic Structure Methods at Work 1 H45 ¼  ðJ12 þ J23 þ J34 þ J45 þ J56 þ J16 Þ; 4 ð4:38Þ 1 1 1 1 H55 ¼ J12  J23  J34 þ J45  J56  J16 : 2 2 2 2 ð4:39Þ By setting all the coupling parameters equal to J, one obtains a matrix which can be written in condensed format, compatible with MathematicaTM (Wolfram 2003, 2014) computer algebra input: H = {{1, 1, 1, -1, 1}, {1, 1, 1, -1, 1}, {1, 1, 0, -1, 1}, {-1, -1, -1, 0, -1}, {1, 1, 1, -1, 0}} Taking, in the same way, the overlap matrix: S = {{1, 1/4, 1/2, -1/2, 1/2}, {1/4, 1, 1/2, -1/2, 1/2}, {1/2, 1/2, 1, -1/4, 1/4}, {-1/2, -1/2, -1/4, 1, -1/4}, {1/2, 1/2, 1/4, -1/4, 1}} one can resolve them analytically with a Mathematica procedure for the secular equation, Solve[Det[H-x*S]==0,x], where x are giving the eigenvalues, as the following set: n pffiffiffiffiffi pffiffiffiffiffi o E ¼ J þ 13J; 0; 2J; 2J; J  13J : ð4:40Þ Here, the level displayed with zero energy is not the ground state. Considering that J is a negative effective parameter, the above list is ordered with respect of energy. The energy of a single Kekulé structure can be taken as the diagonal elements H11 or H22, namely 3 J/2. Then, one may formally estimate the resonance energy, as the difference between nominal ground state and those of a single Kekulé element, Eres ¼   5 pffiffiffiffiffi  þ 13 J  1:1055J; 2 ð4:41Þ equal with the amount obtained in an early work (Pauling and Wheland 1933). With the previously estimated J parameter (see the discussion after Tables 4.7, 4.8, 4.9), equivalent to about –2 eV, the resonance energy results 2.2 eV, or 50.7 kcal/mol, in absolute value. This is higher than the experimental (also empirical) amount, 36 kcal/mol, deduced as the difference between hydrogenation energies of the benzene and those of an hypothetical cyclohexatriene, with three “non-resonating” double bonds, taking as surrogate for the last amount three times the hydrogenation energy of the cyclohexene. However, in these chemical processes there are several other mechanisms involved (e.g. the carbon–carbon bond lengths variation), so that cannot be taken as the undisputed measure for the decoupling of p-type resonance. After the shift of the ground state to zero: 4.5 Computational and Conceptual Valence … n pffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffi o E ¼ 0; J  13J; J  13J; J  13J; 2 13J ; 339 ð4:42Þ one obtains a list that matches the formulas given previously in Table 4.7. The first and second eigenvalues are those assignable to the in-phase and out-of-phase combinations of Kekulé structures, as discussed previously about the direct VB results. As a hint of how the HDVvV modeling, as conceptual VB, can be used for approaching more complex problems, Fig. 4.14 shows an excerpt from the simplified treatment of polyacenes, molecules consisting in concatenated benzene rings. The number of resonance structures needed for description of such systems increases tremendously with the number of rings. However, in first crude instance, one may propose to take only structures with Kekulé appearance. Then, the situation becomes, for linear polyacenes, rather simple. Namely, drawing the double bonds, one may see that only one can be placed “vertically”, while the others form conjugated sequences, in upper and lower parts of the zig-zag carbon chains. For polyacenes with N rings there are N + 1 “vertical” connections and therefore N + 1 Kekulé resonances. Applying the previously discussed topological rules, closed Fig. 4.14 Illustration of the superposition patterns of resonance structures in linear polyacenes with N rings, labeled by the position of the “vertical” link, here taking: a n1 = 4 and b n2 = 7 for the N = 8 molecule. c The primary superposition pattern, marking by circular spots the positions where non-matching arrows appear and by dashed encircling the sequence that must be re-matched. d The matched superposition pattern, showing in dark color the reverted arrows 340 4 Bond! Chemical Bond: Electronic Structure Methods at Work formulas can be obtained for the overlap and Hamiltonian matrix elements of this minimal set of structures. Thus, labeling the resonances by the position of the itinerant double bond (arrow), take the n1 and n2 resonance structures (n1 < n2) out of a N + 1 set, according to the exemplification from Fig. 4.14. One notes that an island is formed between the n1 and n2 positions, with the arrows of the upper part in regular matching (head to head and tail to tail), while the lower n1 and n2 junctions are non-matching. To regularize this sequence, we must return a number of 2(n1 − n2) arrows on the lower margin, the even count ensuring the positive signature of the pattern. Then, one notes that all the arrows at the left of n1 wall, a number 2(n1 − 1) number, and all those at the right side of n2, namely 2(N + 1 − n2) pairs of arrows, are coincident in both resonance structures. Then, counting the one central island, plus the number of superimposed pairs of arrows, one finds a total of 2(N + n1 − n2) + 1 islands. The number of electron pairs is 2N + 1, so that the overlap integral is: S½n1 ½n2  ¼ ð1=4Þjn2 n1 j : ð4:43Þ Assuming a single type of coupling parameter, Ja, for all the links placed in the zig-zag lines, on the horizontal of the drawn polyacene skeleton, while Jb for all the vertical links, the Hamiltonian matrix elements are constructed as described in the following. In the left side of the n1 position there are 2(n1 − 1) intra-island couplings of +Ja magnitude (the superposed arrows), 2(n1 − 1) inter-island contacts weighted by −Ja/2 (the lines along the horizontal, not occupied by arrows), as well as (n1 − 1) inter-island interactions on vertical lines, receiving the −Jb/2. The part contributed by the left side of the central island is summed to (n1 − 1)(Ja − Jb/2). In similar manner, the moiety at the right side gives (N − n2 + 1)(Ja − Jb/2). Counting now the central island, there are 4(n2 − n1) contacts with the Ja coupling strength and (n2 − n1 + 1) with Jb contribution. Summing the interactions, and installing the proportionality to the overlap integrals, the Hamiltonian matrix element is: H½n1 ½n2   jn2 n1 j   1 1 ¼ ðN þ 3jn2  n1 jÞJa1  ðN  3jn2  n1 j  2ÞJb : ð4:44Þ 4 2 This example suggests that the spin-coupling concept can be a rather elegant source of modeling, at the confluence of theory and chemical meaning. Though of limited practical importance, in comparison with the existing palette of computation methods, the HDvV approach to SCVB particular treatment, in principle phenomenological, with empirically adjustable parameters, it keeps the nostalgic savor of the “good old times”. It can be said that this is, at the level of poly-electronic multi-configuration modeling, what is the Hückel model in counterpoint with the ab initio Hartee–Fock or DFT treatment of molecular orbitals: oversimplified, but heuristically relevant. This effective sort of modeling never reached the Hückel level of popularity, in parallel to the fact that methods like CASSCF are not as often used 4.5 Computational and Conceptual Valence … 341 as DFT, for instance, the discrimination being deepened by the cone of shadow fall on the VB theory itself. However, as the reappraisal of VB is contoured in the modern modules of computational chemistry, the perspective of HDvV, as valuable modeling of chemical bonding, may take in some future its quantum of solace. 4.6 4.6.1 Mobilis in Mobile: Electrons Moving Around Mobile Nuclei. Floppy Molecules, Unstable Systems, and Chemical Reactions Jahn–Teller and Related Effects. Vibronic Coupling Although this section is only the lite version of ample and rather difficult problematics, the reader may perceive its importance, particularly considering the growing incidence of molecular dynamics simulations, ported on the power of actual generation of computers. Like every time when computers take over a domain, it is important to urge the users to keep the conceptual compass, holding the track of principles, being aware of possible methodological limitations, digging more for simple clues, extracted from the sea of big data afferent to such calculations. Most often, quantum chemists are practicing the idea of fixed nuclei (Born– Oppenheimer approximation). Of course, when we want to model chemical reactions, we must willingly go beyond this hypothesis. In other circumstances, nature may force us to renounce this convenience, as in the case of so-called Jahn–Teller effect (Bersuker 1984, 2001), when nuclei can become surprisingly strongly coupled with the electron movements, in spite of the large mass difference. The nominal Jahn–Teller theorem refers to poly-atomic non-linear molecules, stating that these are unstable in orbitally degenerate states (particularly ground states), tending to distort in a way that removes the degeneracy (ending with a lower symmetry and a non-degenerate situation). Of course, the Jahn–Teller effect can occur only at molecular geometries having point groups that admit degenerate representations. A related manifestation is the so-called pseudo Jahn–Teller effect, (Bersuker 2013) formally interpreted as a ground state in relative closeness to excited states, to be regarded, together, as quasi-degenerate. Then, one may conceive it as a sort of residual Jahn–Teller manifestation. However, the pseudo Jahn– Teller case is more general than some almost Jahn–Teller occurrence, being driven by couples of states that are not supposed to be tied in a degenerate form, in a higher symmetry, to which a given molecular geometry can be formally idealized. The Jahn–Teller and pseudo Jahn–Teller manifestations are generally called vibronic instability. The keyword is a composite term for vibration-electronic coupling. The vibronic coupling is basically always manifesting in molecules, even when its strength is not sufficient to cause extreme floppiness of a given geometry, driving some less spectacular phenomena, such as the broadening of electronic 342 4 Bond! Chemical Bond: Electronic Structure Methods at Work transitions in bands, that are the convolution of so-called vibronic spectral progressions (in certain cases, as fine structure of visible or ultraviolet spectra). In solid-state physics, the equivalent of vibronic interactions is known as electron– phonon coupling (phonons being collective modes of vibration and deformation, with infinite repetitive patterning). Figure 4.15 shows the patterns of vibronic instability of a nuclear configuration (molecular geometry) conventionally represented as the zero point of the abscissa, which is evasively conceived as a deformation that destroys a given high symmetry. Since it implies degenerate states, the Jahn–Teller energy landscape is also multi-dimensional, here representing a one-dimensional section in such an object, as seen in panel (a) from Fig. 4.15. The curves crossing at the center of the abscissa represent the degeneracy point [here being represented as a double degeneracy, as the left and right side parabolas of the panel (a) are showing two states, whose energies coincide at crossing (in a very simplistic interpretation)]. The reaction coordinate from abscissa is a generic molecular deformation (qualitatively decided by the symmetry of the involved states) which, for infinitesimal displacements is a vibration mode (or a combination of vibration coordinates). The Jahn–Teller instability point can also be named conical intersection, considering a generalization of the linear course of the crossing energy curves, for the case of two or many dimensions. On the other hand, the conical intersection is also known to computational chemists for other situations, not emerging from a Jahn– Teller effect, but yet pertaining to a vibronic-type of problematics. Namely, it is about states keeping different symmetries along a certain molecular deformation (which has a symmetry not allowing their coupling). In this situation, the states look as non-interacting along the given reaction coordinate, having a linear or conical pattern at interaction. However, in such cases are ignored the deformations with symmetries that may couple the states with crossing energy profiles, which will bring the pattern merely resembling the pseudo Jahn–Teller case (avoided crossing), or a general version of it, with asymmetric potential wells. The Jahn–Teller instability is driven by non-diagonal elements, whose magnitude is at least linearly proportional with the amplitude of the symmetry-determined deformation modes. In the degeneracy circumstances, the eigenvalues are also evolving with linear pattern, around the instability point. The higher order Fig. 4.15 Generic patterns of vibronic instability, representing the potential energy as a function of a reaction coordinate: a Jahn–Teller effect (conical intersection), b pseudo Jahn–Teller effect, c Renner–Teller effect 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 343 polynomial dependence of the non-diagonal coupling elements may decide a further way of warping the potential energy surface, but the intersection is dominated by the linear terms. In the pseudo Jahn–Teller case, the non-diagonal elements are at least linearly depending on the deformation amplitude, but the repercussion in eigenvalues will be, in the spirit of perturbation dependence, at least second order, having then the instability point as a smooth maximum on the non-degenerate ground state energy profile (instead of a sharp encounter of degenerate states, like in the Jahn–Teller form). Many stereochemical or reactivity problems can be formulated as a pseudo Jahn–Teller case. For instance, the hill of transition state (activated complex) conceived (qualitatively, or in corresponding computation experiments) as sitting between energy valleys representing the reactants and products, can be considered a general case of pseudo Jahn–Teller type, because it implies critical interaction between the spectral states of the system, coupled by a distortion coordinate, even when the higher states are not explicitly accounted, and even when no symmetry factor is involved. The Renner–Teller occurs only in linear molecules. It is basically conditioned like Jahn–Teller one, taking place on degenerate ground states, but, for symmetry reasons, the non-diagonal element cannot be lower than the second rank, with respect of distortion coordinates. Then, the pattern is like a mixing of Jahn–Teller and pseudo Jahn–Teller features; having two states meeting at the zero of the abscissa, but the profile near this vicinity is rounded, due to the second-order dependence (see Fig. 4.15c). In spite of their somewhat exotic incidence in the landscape of quantum chemistry, yet dominated by the Born–Oppenheimer static approximation, or their presence obscured by the heavy machinery of molecular dynamics calculation, the vibronic effects are very important drivers of many special properties of high practical importance. The bistabilities of all sorts implied in devices used in processing or storing information, practically all the phase transitions, superconductivity (Ceulemans et al. 1997), stereochemistry (Mösch-Zanetti et al.), reactivity, and so on, have as deep engine the vibronic coupling. Besides, the vibronic paradigm makes the universe of molecular dynamic processes an elegant place, bringing qualitative symmetry reasons as background. If we adapt a bit of Aristotle’s philosophy about causal typologies to problems involving the dynamics of nuclei, we can assign the “causa materialis” to the considered particular system, interpret the trend to energy minimum as “causa finalis”, and take the symmetry reasoning as “causa formalis”, arriving at the vibronic interactions as “causa efficiens”, i.e. the driving force of the phenomena (bistabilities, distortion trends, dissociation of the reactive transition state). 344 4.6.2 4 Bond! Chemical Bond: Electronic Structure Methods at Work A Simple Approach of the H3 Prototypic System. Example for Reaction Potential Energy Surfaces and E ⊗ E-Type Jahn–Teller Effect The neutral trihidrogen molecule, H3, is one of the simplest systems to be considered for problems discussing profiles resembling a chemical reaction (whose transition state can be interpreted as a pseudo Jahn–Teller). Besides, the regular triangular H3 can be taken as Jahn–Teller instability of a degenerate state, as will be seen. It is nice that in this conjuncture one may also invoke a simple Valence Bond type model, known as London-Eyring-Polanyi-Sato (LEPS), (Eyring 1931; Polanyi and Wong 1969; Sato 1955) very conveniently placed in the continuation of the discussions from Sect. 4.4.1. Namely, we can take parameters describing the energy profiles of the H2 molecule as transferable for the account of interactions on edges of the H3 general triangle (or linear) configuration. The model is crude, but it retains enough quantum chemical truth. We will even rely on the quantities resulting from the very basic level of calculation, intentionally selected like this, for the sake of insight, namely the STO-3G based CASSCF curve for the ground singlet and triplet states of the H2 system. Then, for the H3 case we will not do new calculations, choosing instead a Heisenberg-Dirac-van Vleck (HDVV) spin Hamiltonian, whose principles were extensively discussed in Sect. 4.5.1, as a phenomenological version of a Valence Bond modeling (roughly satisfactorily, instead of the CASSCF approach). In a pragmatic approach we will take the data displayed in Fig. 4.16 for the ground potential energy curve (labeled 1Rg(0)) as a function of inter-atomic separation, R, fitting it with a Morse-type potential: ð4:45Þ R13 (Å) E13 (Hartree)  lðRÞ ¼ 2h0 þ D
ð1  Exp½a
ðR  R0 ÞÞ2 1 ; R12 (Å) R13 (Å) R12 (Å) Fig. 4.16 The contour plot (left side) and the 3D potential energy surface of the of the Ed ðR12 ; R13 ; R23 ¼ R12 þ R13 Þ ground state in the case of linear configuration of the H3 molecule, mimicking the landscape of a H + H2 ! H2 + H reaction 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 345 where h0 are the orbital energies of individual hydrogen atoms, D is the depth of the energy minimum with respect of the 2h0 plateau at large R, while R0 is the position of the minimum and a an adjustable parameter. For the energy of the triplet labeled 3 Ru in Fig. 4.16, we can take a simple exponential decay: mðRÞ ¼ 2h0 þ u
exp½b
R0 ; ð4:46Þ imposing that it tends also to the 2h0 limit, at large interatomic distances. We have h0 = −0.46 Hartree, D = 0.2186 Hartree, R0 = 0.7337 Å, a = 2.0642 Å−1, u = 4.1977 Hartree, and b = 3.2076 Å−1. We recall that we discard the issue of numeric accuracy. We should have obtained a h0 = −0.5 Hartree rigorous value, but this does not impinge upon the conceptual level. Otherwise, the R0 and D are relatively close to experiment or richer calculations. Now, we will consider that these curves are accounted by the Eq. (4.4), simplified to the case of no overlap, having then l(R) = 2h0 + Q(R) + J(R) for the ground curve and m(R) = 2h0 + Q(R) − J(R), for the excited one. The J(R) is negative in the entire domain, since effectively absorbed the terms due to the overlap, neglected in the explicit model, but acting in its background. Then, from Eqs. (4.45) and (4.46) one may find the distance dependence of the Coulomb and Exchange effective integrals: 1 QðRÞ ¼ ðlðRÞ þ mðRÞÞ  2h0 ; 2 ð4:47aÞ 1 JðRÞ ¼ ðlðRÞ  mðRÞÞ: 2 ð4:47bÞ These functions will be taken as ingredients for the HDVV modeling of the H3 molecule, fully equivalent with the above mentioned LEPS model (except the actual re-parameterization). With principles exposed in the above sections, one may easily form the Hamiltonian in the basis of the corresponding spin configurations, arriving at the following eigenvalues: Ed ðR12 ; R13 ; R23 Þ ¼ 3h0 þ QðR12 Þ þ QðR13 Þ þ QðR23 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi JðR12 Þ2 þ JðR13 Þ2 þ JðR23 Þ2  JðR12 ÞJðR13 Þ  JðR12 ÞJðR23 Þ  JðR13 ÞJðR23 Þ; ð4:48Þ for the two spin-doublet states, and Eq ðR12 ; R13 ; R23 Þ ¼ 3h0 þ QðR12 Þ þ QðR13 Þ þ QðR23 Þ  JðR12 Þ  JðR13 Þ  JðR23 Þ; ð4:49Þ for a spin quartet, discarded in the following discussion (since it is higher in energy and does not have interesting behavior). In the above formulas Rij are the distances 346 4 Bond! Chemical Bond: Electronic Structure Methods at Work between the H atoms of the H3 molecule, in any general circumstance. For the spin doublet states, the solution with minus sign in front of the square root is generally lower, the couple becoming degenerate at equilateral triangle R12 = R13 = R23 = R, when both energies are 3h0 + 3Q(R). Taking the convention of a linear placement of the atoms, e.g. imposing R23 = R12 + R13, (assuming the 3–1–2 ordering), one draws an energy surface containing the reaction coordinate of a H + H2 ! H2 + H substitution process. Reading this map qualitatively, one may figure the path: e.g. representing the atom #3 coming along one the valley of the decreasing R13 coordinate, progressively elongating the R12 bond of the initially stable H2 molecule, passing over a barrier representing the linearly symmetric activated complex, at about R12 = R13 * 1 Å, the process being completed by the mirror evolution, when atom #2 is departing along the growing R12 line, while the R13 parameter consolidates a new H2 molecule. Intuitively, for this process to go, the relative velocity of the initially separated H (#3) and H(#1)–H(#2) must incorporate enough kinetic energy of the nuclei (i.e. extra to those of the electrons in atomic and molecular orbitals). This must exceed the height of the encountered barrier. The compromise approach is to do molecular dynamics considering that the nuclei are moving by classical mechanics, combining the initial impetus of the moving parts with the energy gradients taken from a given potential energy surface as reactive forces, and drawing trajectories of the involved atoms. This can be done, in the frame of molecular dynamics, without charting a potential energy surface (or multi-dimensional hypersurface in larger systems, with many geometry variables), estimating gradients, e.g. computed by quantum mechanics procedures at an instantaneous position of nuclei, starting from an initial set of positions and velocities of the nuclei. To calculate gradients by quantum means can be rather costly, so that molecular dynamics calculations are done within approximations, as is the case of Car–Parrinello methods (Car and Parrinello 1985), assuming conventional mass to electrons, to be treated altogether with nuclei, in a DFT frame. Full and free molecular dynamics can drive the system in a situation problematic for the chosen electronic structure method, as would be a degenerate or near-degenerate state in a DFT approach. However, we will not debate here the domain of molecular dynamics. We just recall, without detailing, that, rigorously speaking, the nuclei must be in principle treated by wave functions and related probabilities of distribution. In quantum mechanics there is the so-called tunneling phenomenon, where a barrier can be passed even at lower kinetic energy than a faced potential barrier height. Such problems are not easy to tackle, and present-day computational chemistry is not able to approach them often, in ways sufficiently close to this conceptual level. Indulging ourselves in a joke, one may say that, whether quantum chemistry can offer rather good pictures of static molecules, the dynamics is yet in the very crude stage of the very first Lumière movies: not very accurate, with short and simple action, but yet nice and attractive. 347 R13 (Å) E (Hartree) E (Hartree) 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … R13 (Å) R12 (Å) R12 (Å) Fig. 4.17 Different views of the potential energy surface of the two doublet states (see Eq. 4.48) drawn as function of R12 and R13 edges, when the 3–1–2 angle is kept fixed at 60°. Note that on the main diagonal, with R12 = R13 = R23, the curves are crossing in an orbitally degenerate state having the case of Jahn–Teller instability Figure 4.17 aims to catch the Jahn–Teller effect occurring at equilateral triangle of the H3 molecule, fixing one angle at 60°, so that at the R12 = R13 the highest symmetry is achieved. One observes the crossing of the sheets describing the energies of spin doublet states, as a function of R12 and R13, while R23 = (R212 + R213 − R12R13)1/2, along the main diagonal of abscissas plane. In the right side panel, looking at the E versus R12 margin, while R13 is sufficiently large to assume a non-interacting third hydrogen atom, one observes the couple of the states originating from the ground singlet and excited triplet curve of the H2 molecule made of the #1 and #2 atoms. In order to see what is meant by conical intersection in the case of Jahn–Teller case, let us perform a transformation of coordinates. As sketched in Fig. 4.18, one may propose getting the relative positions of three atoms gliding them along three lines mutually rotated at 120°, with the R1, R2, and R3 distances from a fixed origin, allowing also negative values of the R parameters. There is no restriction in reproducing any mutual placement. For instance, keeping fixed the atoms #2 and #3, with a proper negative R1 shift, one may bring the #1 atom collinear with the QA QθE QεE R3 R1 R2 Fig. 4.18 Symmetrized coordinates for positions of the atoms in a general case of three atoms in a fixed plane, based on an equilateral triangle as reference 348 4 Bond! Chemical Bond: Electronic Structure Methods at Work other ones. Furthermore, we take now the collective movements expressed by the following coordinates: R1 R2 R3 QA ¼ pffiffiffi þ pffiffiffi þ pffiffiffi ; 3 3 3 rffiffiffi 2 R2 R3 QEh ¼ R1  pffiffiffi  pffiffiffi ; 3 3 6 R2 R3 E Qe ¼ pffiffiffi  pffiffiffi ; 2 2 ð4:50Þ graphically represented in Fig. 4.18. The combination coefficients are similar to those used in the design of sp2 hybrids, expressing the intention to have a reference based on the trigonal symmetry, since we aim to visualize the Jahn–Teller effect in the regular triangle. Indeed, one may see that, at regular triangle, namely equal radial Ri parameters, the QEh and QEe amplitudes are null. To make the scheme operational with the considered LEPS model, it is useful to have the Cartesian coordinates as a function of the introduced symmetrized modes: 0 x1 B @ x2 x3 1 0 R1 y1 C B  1 R2 y2 A ¼ @ 2 y3 0  12 R3 0 pffiffi 3 2 R2 pffiffi  23 R3 p1ffiffi QA 3 B B 1 A ¼B B  2pffiffi3 Q þ @  2p1 ffiffi3 QA þ þ 1 C A qffiffi 2 E 3Qh 1 ffiffi E p Q 2 6 h 1 ffiffi E p Q 2 6 h  þ 1 0 1 ffiffi E p Q 2 2 e 1 ffiffi E p Q 2 2 e 1 A 2Q  1 ffiffi E p Q 2 2 h  12 QA þ þ 1 ffiffi E p Q 2 2 h 1 2 þ qffiffi 3 E 2Q e 1 2 qffiffi 3 E 2Qe C C C: C A ð4:51Þ Thus, one may formulate the R12, R13, and R23, distances needed in the (4.48) equations as a function of the QA , QEh and QEe symmetrized amplitudes. In the D3h symmetry of a regular triangle, the couple constituting the spin-doublet orbital-quartet is labeled as 2E′, while the spin quartet resulting from three electron unpaired in s-type orbitals has the 4A2′ representation. The QA coordinate corresponds to the A1′ representation, its action not changing the point group. The two QE coordinates are belonging to the E′ representation. The individual action of the h coordinate, as defined previously, bring the triangle to isosceles with C2v point group, while gliding along e or a h − e combination leads to a Cs situation (in general, destroying completely the symmetry, down to C1, but it happens that the three-points case keeps inevitably the single molecular plane as symmetry element). In this case we met an E-type degenerate state, coupling with an E-set of distortion, the situation being denoted by E ⊗ E (or e ⊗ E), where the first labeling Fig. 4.19 The Jahn–Teller instability of the equilateral triangular configuration of the H3 molecule, represented as a function of symmetrized coordinates, within a LEPS derived model. The left side evidences the conical intersection at the QEh ¼ 0 and QEe ¼ 0 high symmetry point 349 E ( Hartree) 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … QθE QεE stands for the electronic state, while those after the circled product symbol (⊗) are representing the coordinates. Since the H3 molecule is not a bonded system (Mahapatra and Köppel 1998; Mistrík et al. 2001), it is maybe not the complete example of a Jahn–Teller manifestation, since we cannot find on the potential energy map of the states originating from the 2E′ terms, as function of QE coordinates, valleys with minima, where the system is stabilized after the distortion from the equilateral triangle. The system tends to distant H and H2 moieties. However, since the Jahn–Teller effect is about the instability itself, Fig. 4.19 illustrates clearly the case of the conical intersection, the acute point in the center of the lower energy sheet meeting the sharp bottom of the upper surface. The E ⊗ E type of effect is very popular in coordination chemistry, where it leads to distortions like tetragonally elongated octahedra, for Oh or O point group references, or compressed and elongated tetrahedra, in Td cases. Without entering the symmetry details, note two other prototypic cases, T1 ⊗ (E + T2) and T1 ⊗ (E + T2), saying that a triply degenerate state, T1 or T2 in octahedral or tetrahedral point groups, can couple with E and T2 distortions. Besides, if the symmetry carries parity labels, u and g, then, irrespective the index of the T representation for the electronic state, the reaction coordinates are always of even parity, namely Eg + T2g. Because the illustration of the E ⊗ E Jahn–Teller effect on the H3 case was a bit shadowed by the non-bonded status of the whole system, not catching minima on the potential energy surface, we expose now, without proof or deep discussion, the master Hamiltonian that led to the generic energy profiles shown in Fig. 4.20. Thus, in the basis of the two components of the E degenerate state, with parametric dependence on the QEh and QEe coordinates (whose concrete meaning is to be adapted to the problem at hand), the E ⊗ E Jahn–Teller Hamiltonian matrix is: HE E 0 ¼ @  2  2 QE þ QEe  h    2 2 QEh þ iQEe
V þ QEh þ QEe iQEh QEe
W 1 2k    1   2 2 QEh  iQEe
V þ QEh þ QEe þ iQEh QEe
W A;    2 2 1 QEh þ QEh 2k ð4:52Þ 350 4 Bond! Chemical Bond: Electronic Structure Methods at Work E E QεE QεE QθE QθE Fig. 4.20 Typical surfaces for the E ⊗ E Jahn–Teller effect. Left side: the so-called “Mexican hat” profile, obtained confining to first-order coupling parameter. Right side: the so-called “tricorn”, resulted including second-order non-diagonal vibronic coupling elements where k is the force constant associated with the average of state energies, V is the first-order vibronic coupling parameter, while W is the second-order companion. The formulation is the result of series expansion of the Hamiltonian, as a function of distortion coordinates, stopped to the second order, in both the diagonal and non-diagonal elements. One notes that at Q = 0 the two states are degenerate, with both energies conventionally imposed in the zero of the scale. The increasing distortion causes non-null non-diagonal elements, that are creating the two different energy solutions. An alternative expression is reached by converting the coordinate frame in a polar-alike formulation (QEh , QEe ) ! (q, u): HE E ¼  1 2 2 kq exp ðiuÞq
V þ exp ð2iuÞq2
W  exp ðiuÞq
V þ exp ð2iuÞq2
W : 1 2 2 kq ð4:53Þ The eigenvalues of this Hamiltonian are: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 e ¼ kq2  V 2 þ q2 W 2 þ q
V
W
Cosð3uÞ: 2 ð4:54Þ The linear pattern implied in the conical intersection profile is due to the predominance of the linear coupling parameter, V, in the vicinity of undistorted system, at q ! 0. If we take W = 0, the solutions are independent on the u coordinate (or on the QEh /QEe ratio), a minimum being reached at a circular profile with qmin = (2 V/k)1/2. This is the celebrated “Mexican hat” potential energy surface, shown in the left panel of Fig. 4.20. A vertical section yields the pattern from Fig. 4.15a. The higher order couplings are determining further warping of the potential energy surface, with the trigonal pattern due to the cos(3u) term in the 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 351 eigenvalues. Thus, one arrives at a series of three minima and three saddle points (minima along q coordinate, while maxima with respect of u), if we follow a circular profile nearby the qmin amplitude. 4.6.3 The Computational Approach of the Pseudo Jahn– Teller Effect (Second-Order Vibronic Coupling) In the following we will deal with the pseudo Jahn–Teller type of vibronic interactions. The master equation for this effect describes the curvature of a ground potential energy surface, K, with respect of a given nuclear coordinate, Q, in the spirit of second-order perturbation theory:  _ 2 _  X W0  @H  @2H  @2E @Q Wr   K¼ ; W 2 ¼ W0 @Q2 Er  E0 @Q2 0 r6¼0 ð4:55Þ the first term of the last member (expectation value of the second-order derivative from Hamiltonian) being called non-vibronic curvature, K0, while the summation over the excited states, labeled r, forms the vibronic contribution Kv. The vibronic instability appears when total negative curvature occurs, due to the predominance of Kv over K0 (positive, in the above terming). Computationally, this corresponds to detection of imaginary frequencies (since the frequency is proportional with the square root of K, having also the meaning of force constants for vibrations along the Q mode). At the same time, a non-vibronic and vibronic couple occurs in any force constant, even in the frequent situations where the K0 predominates, having locked a stability with respect of a Q-type distortion. Such analysis, as the calculation of molecular vibrations is tacitly conceived in situation of null gradient, ∂E/∂Q = 0. If we take the Hamiltonian literally, the dependence on the distortion coordinates  _  appears only in electron-nuclear terms and, then, the W @H =@QW matrix 0 r elements, called vibronic coupling parameters, are related with the electric field at nuclei. Correspondingly, K0 would be determined by the gradient of the electric field at the nuclei. However, with respect of computational practice, things are not so simple. The used Hamiltonians are self-consistent, incorporating the wave functions, which actually yield contributions to derivatives with respect of the moving nuclei. In other words, the above formula will be right only for the exact W0 ground state and excited states Wr of the same quality. In real practice, we must consider other formulations of second-order incidence of the vibronic coupling, to be consistent with a given sort of self-consistent calculation. The issue is caught in the following quotation from Pulay (1987): “Nuclear coordinates as perturbation parameters are different from other common perturbations, e.g. weak external fields”, who also noticed that “the basis functions must be coupled to the nuclei”. 352 4 Bond! Chemical Bond: Electronic Structure Methods at Work Working with atomic basis sets, we must consider that these are floating with nuclei, along a distortion coordinate, a fact ignored in a formulation like (4.55), or in many interpretations of potential energy calculations. The floating molecular orbitals are obtained keeping frozen the LCAO coefficients, while the AOs are allowed to move, tied to their nuclei. Taking Eq. (2.178) as the most general formulation of the ground energy, the non-vibronic term is what results applying the second-order derivative of the given coordinate to the molecular integrals: K0 ¼ M X M X l þ clm m @2 lj^ hjm @Q2 M X M X M X M 1X 2 l l0 m m0 @2 ðll0 jmm0 Þ; Cll0 jmm0 @Q2 ð4:56Þ while the density matrix coefficients are kept frozen, like in the Q = 0 point. This is adapted from a multi-configuration energy formula, but is in principle valid for any type of calculation. Note that the two-electron integrals undergo derivatives with respect of Q variable, even though the operator is entirely independent of nuclei, because of the orbital flotation (that affects the distance between the AOs giving the integrals). Without detailing a demonstration (retrievable from Chibotaru and Cimpoesu 1997), we give the vibronic component adapted to self-consistent calculation results: Kv ¼ 2 X r6¼0  2 _  1 @ @    W H Wr  E0 W W ; Er  E0 @Q 0 @Q 0 r ð4:57Þ where the excited states, labeled r, are taken as resulting from single-excitations with respect of the W0 wave function and its calculation method. Note that, to be distinguished from (4.55), where matrix elements from derivatives of the Hamiltonian are taken, the new formulation considers derivatives from the whole matrix elements. In this way, it is better suited for implementation in calculation methods, since, once the evaluation of matrix elements is known, their derivatives are technically approachable. The involved single-excited states are equivalent with the time-dependent (TD) (Kutzelnigg 1989) version of the given computation level. Thus, if the ground state is a Hartree–Fock result, we must use the complete set of TD-HF states. Similarly, vibronic analysis within the DFT method needs the TD-DFT states. For multi-configuration methods, e.g. CASSCF, there is a TD equivalent, implying single orbital promotion from core to the active space, from active space to virtuals, as well as from core to virtuals. The orbital promotions inside the active space will not change the solutions. The multi-configurational methods would be appropriate to adapt a discussion of Jahn–Teller case computations (first-order vibronic coupling), in the spirit of those sketched here for second-order effects. The fact that, for consistency, we must involve single-orbital promotions is in line with the situation 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 353 that the vibronic coupling follows a one-electron effective nature, paralleling the true one-electron nature of the electron–nuclear interaction component of the Hamiltonian. The very interesting fact is that the expansion of vibronic curvature over excited states from (4.57) is fully equivalent with a formulation in terms of orbital promotions, using quantities related with the so-called Coupled Perturbed theories (CP) (Yamaguchi et al. 1994) which, in short, define the algorithm to take derivatives with respect of different perturbation (molecular displacements, external fields) aside a self-consistent calculation method. Thus, in CP methods, in case of a perturbation by the Q distortion coordinate, the matrix U(Q) gives the transformation leading to the new self-consistent MOs, at a certain Q amplitude, starting from the floating orbitals (that kept the frozen LCAO coefficients, while the AOs were following the Q displacement, with the nuclei). In the given situation, we would need only the first-order expansion, U(Q) * I + Q
U(1), where I is the identity matrix. Taking a general level of electronic structure, expressible with the help of natural orbitals with ni occupation numbers, e.g. from a multi-configuration calculation, the vibronic curvature is done running over orbital couples: Kv ¼ 2 allMO X X i j i  nj  ni uð1Þ V ð1Þ : ji ij ð4:58aÞ Concretized to single-determinant methods, HF or DFT, the summation goes over excitations from occupied orbitals (formalized with ni = 1, in general unrestricted frame) to virtuals (nn = 0), and it becomes: KvHF=DFT ¼ 2 occ X virt X i n ð1Þ ð1Þ uni
Vin : ð4:58bÞ The u(1) coefficients represent the first-order U matrix from the CP techniques. The V(1) quantities, named orbital vibronic coupling elements, can be regarded as companions of the so-called B matrix from Coupled Perturbed methods. For the HF frame, where the B matrix elements are, in first order, defined as: ð1Þ ð1Þ ð0Þ ð1Þ Bi!n ¼ fni ei Sni þ occ occ X 1X ð1Þ S Ani ; 2 j l ni lj ð4:59Þ the newly introduced first-order orbital vibronic constants are: ð1Þ ð1Þ ð0Þ ð1Þ Vi!n ¼ fni  ei Sni  occ occ X 1X ð1Þ S Ani ; 2 j l ni lj ð4:60Þ where f(1) are the quantities resulting from first-order derivatives with respect of Q coordinate of all the integrals summing a Fock matrix element, S(1) is the 354 4 Bond! Chemical Bond: Electronic Structure Methods at Work first-order derivative of the overlap appearing along with floating orbitals (with identity matrix as zero-th order), and Ani lj ¼ 2ðiljnjÞ  ðiljjnÞ  ðijjlnÞ, a combination of two-electron integrals. For other methods, such as DFT or CASSCF, one may formulate vibronic coupling in a similar way, arriving at companions of the B matrix, by reverting the sign of the S(1) cofactors. In general, the vibronic part cumulates all the terms accounting for the relaxation of floating orbitals (with LCAOs frozen at Q = 0 undistorted reference) toward the new self-consistent canonical MOs of the given method. Since the CP methods are used in the analytic ab initio calculations of vibration modes and frequencies, their adaptation to vibronic analysis is naturally suited. The CP techniques can be reformulated as accomplishing the self-consistent status of a perturbed system, by the interaction of the ground level with excited states of TD type. We note again the remarkable fact that the formulation of vibronic coupling in terms of many-electron states from (4.57) is fully equivalent with one-electron orbital expansions from (4.58a) and (4.58b) equations. The following discussions will relate to the Hartree–Fock frame, where the vibronic analysis is done opening the black box of CPHF module of frequency calculations, in order to take the U matrix and produce the V companion of the B matrix. We did this entering the corresponding module of the GAMESS code. We remain at the HF level for the convenient handling of formulas and algorithmic adaptation, intending also to keep transparency at the conceptual level. We will take a simple problem, considering the pyramidal geometry of ammonia as a pseudo Jahn–Teller instability of the triangular planar configuration, having the D3h point group, i.e. the highest symmetry possible with the given NH3 composition. In general, the vibronic approach to stereochemistry consists in posing the question why a system having possible higher symmetry is found in lower ones, or presents certain floppiness. The ammonia shows also mobility, in “umbrella flipping” style, moving between conformations with nitrogen above and below the hydrogen made plane. Thus, in the D3h reference, the pyramidal NH3 can be interpreted as a (A1′ + A2″) ⊗ A2″ pseudo Jahn–Teller effect (Atanasov and Reinen 2001, 2002). This means that the A1′ ground state couples with one or more A2″ excited states, via the “umbrella” vibration mode, spanning the A2″ representation. The spectral terms are complementary to the orbital symmetries of the frontier couple: the HOMO has the a2″ representation, being in D3h an almost pure pz AO, perpendicular on the molecular plane, while LUMO is a1′ looking like a symmetric combination of NH anti-bonds. Formally, the excited state symmetries result from HOMO-LUMO orbital promotion A1′ = (a2″)2 ! A2″ = (a2″) (a1′). Figure 4.21 shows the profile of the “umbrella” mode, computed with the 6-311 ++G(2d,2p) basis set, by different methods: RHF, MP2, CCSD(T), and B3LYP. The relaxed energy profiles (done by tuning the out-of-planarity angle, a, and optimizing the bond lengths at each a point) look very similar for all the different methods, after translating the a = 0 point to the zero of the energy scale. This comparability is another reason for confining the vibronic analysis to the simplest 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … Fig. 4.21 The potential energy curve of the NH3 “umbrella” mode, computed tuning the deviation of NH line from the plane parallel with the H3 fragment, while the NH bond lengths are optimized. The results of different calculations (HF, MP2, CCSD(T), B3LYP/6-311++G(2d,2p)) are shifted with the planar configuration at zero energy HF MP2 355 CCSD(T) B3LYP 1.5 E(eV) 1 0.5 0 α(°) -0.5 -45 -30 -15 0 15 30 45 HF level, in the spirit of the “Ochkam razor” principle of avoiding unnecessary complexities: “Entia non sunt multiplicanda praeter necessitatem” (the thing should not be multiplied beyond strict necessity). The situation is due to the fact that the correlation energy is roughly parallel with the number and nature of electron pairs, being almost the same along a distortion, if no bond breaking occurs. Then, the HF simplest level contains already the basic mechanism deciding the molecular geometry. Such regularities may not hold in case of bond-breaking and formation processes. For instance, the reaction profile of the Cl(−)


H3CCl ! ClCH3


Cl(−) nucleophile substitution, can also be formulated as a (A1′ + A2″) ⊗ A2″ pseudo Jahn–Teller effect, with the symmetric [Cl


H3C


Cl](−) complex as D3h reference. However, in this case the energy profile of the system passing over the reaction barrier differs more with the computation method, because bond redistribution is implied. Returning now to the analysis of the pseudo Jahn–Teller effect in trigonal planar ammonia, Table 4.10 shows the predominant 2uni Vni terms from the vibronic curvature of A2″ pyramidalization mode (the first ten couples, ordered in decreasing absolute value of their contribution). Among the occupied orbitals, only the HOMO has the a2″ symmetry needed to realize the coupling with the A2″ deformation. However, in the range of virtual orbitals, in spite of the expectation to see a decisive role of the LUMO, one finds that the major contributions are collected from rather high virtuals spanning the a1′ symmetry, such as the #18 and # 38 components. These two contributions are giving a −0.452 mdyne/Å from the total vibronic curvature, Kv = −0.828 mdyne/Å. The rest is due to many small terms. The non-vibronic positive part is Kv = 0.252 mdyne/Å. One observes that the negative total curvature K, representing the unstable position on the top of the hill between the two minima of opposed orientations of the pyramidal geometries, is due to the predominance of the vibronic negative Kv over the K0. 356 4 Bond! Chemical Bond: Electronic Structure Methods at Work Table 4.10 The orbital couples with predominant contribution to the vibronic curvature Kv for the A2″ “umbrella” mode in planar NH3 Occupied MOs # MO 5 5 5 5 5 5 2 3 4 5 2 Virtual MOs Sym. ei(Hartree) a2″ a2″ a2″ a2″ a2″ a2″ a1′ e′ e′ a2′ a1′ −0.391 −0.391 −0.391 −0.391 −0.391 −0.391 −1.134 −0.660 −0.660 −0.391 −1.134 # MO 18 38 36 36 13 29 24 51 52 37 9 Vibronic couples Sym. en(Hartree) De(eV) a1′ a1′ a1′ a1′ a1′ a1′ a2″ e″ e″ a1′ a2″ 0.723 2.892 2.089 2.089 0.277 1.278 1.009 5.449 5.449 2.443 0.197 30.32 89.35 67.48 67.48 18.17 45.41 58.32 166.23 166.23 77.12 36.22 2 2 (b) (a) (c) 1 N H 0 2uni Vni (mdyne/Å) −0.284 −0.167 −0.082 −0.082 −0.052 −0.034 −0.031 −0.025 −0.025 −0.017 −0.012 1 N H 0 -1 -1 2 2 N H New bonding Fig. 4.22 Electronic density variations during the a2″ pyramidal distortion in the planar NH3, for a section through an N–H bond perpendicular to the molecular plane. a Nonvibronic part of density variation: the first derivative with respect of pyramidal deformation (floating AOs, frozen LCAOs). b Vibronic part of density variation: the first derivative giving the relaxation of the electronic cloud, toward the density of the distorted system. c The total first derivative of electronic density (vibronic plus non-vibronic). The dark areas denote density accumulation, the light ones show the depletion For a better insight, Fig. 4.22 offers a visual account of the non-vibronic, vibronic, and total density displacement, expressed as first derivative with respect of the a2″ pyramidalization coordinate. Note that, according to the so-called Wigner 2n + 1 theorem (Epstein 1974), the first derivative with respect of density or of MO functions is consistent with undertaking a second-order variation in the energy, as we considered for the pseudo Jahn–Teller theorem. In general, in the frame of self-consistent methods, the n-th derivative of density or LCAOs yields energy in the 2n + 1 order. Panel (a) from Fig. 4.22 shows the result of AO flotation, namely applying the atomic movements to the carried bases, seeing the density accumulation related to the upward movement of nitrogen and the downward shift of hydrogen atom of the 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 357 taken NH bond. The movement of both N and H sites is conventionally chosen to keep unmoved the center of gravity. The vibronic density derivative from panel (b) shows the opposed density flow, dominantly on the nitrogen atom (produced by the part due to derivatives of coefficients of the atomic orbitals in the molecular functions). This can be interpreted as analogue to Le Châtelier’s principle known in chemical equilibria, saying that a system undergoing a perturbation acts in the direction minimizing the induced stress, by counterpoise modification of its inner parameters. Thus, a part of the density moved upwards with the nitrogen atoms, seems to return back, in the vibronic relaxation movement. The part due to hydrogen is less visible, because this atom has lesser total density and also fewer atomic orbitals, as channels of reorganization. The total map of density flow shows in panel (c) a map of accumulation of density for both nitrogen and hydrogen atoms, looking like a build-up of an orbital overlap component. This is in line with the qualitative explanation by the so-called Walsh rules (1953), that the distortions are taking place in the sense of gaining more overlap encounters. Thus, in planar NH3, the pz orbital of the nitrogen, perpendicular to the molecular plane, cannot make bond with hydrogen, because this does not have a p-type capability in its valence shell. However, the pyramidalization entrains the pz in hybridization and turns it usable for the bonding interests. In heuristic manner, the vibronic approach can be presented as a combination of Walsh diagram interpretations (invoking orbital factors) and Nyholm–Gillespie electron pair repulsion models (as relaxation driven by as electron–electron effects). 4.6.4 The Vibronic Orbitals The exemplification of vibronic analysis for the “umbrella” mode in planar NH3 showed that the HOMO-LUMO couple is not found among the significant contributors to the vibronic curvature. It seems counter-intuitive to find significant vibronic coupling between orbitals separated by large energy gaps. In order to retrieve an effective simplicity, we introduce now the concept of vibronic orbitals. These functions will be designed by a transformation of the i occupied orbitals in a new set, a, bringing simultaneously the n virtuals to g modification: jai ¼ occ X jgi ¼ virt X i n r ia jii; ð4:61aÞ sng jni: ð4:61bÞ ð1Þ ð1Þ in a way concentrating the Kv term in as few uga
Vga terms as possible, in the (4.58b) type of summation. For this purpose we will use Lagrange multipliers 358 4 Bond! Chemical Bond: Electronic Structure Methods at Work imposing extrema conditions over Kv, within the constraint of orthogonality conservation: " #! occ occ X virt occ X X d X ð1Þ ð1Þ uga
Vga  k r ia dij r ja  1 ¼ 0; dr ia a g j i d dsng occ X virt X a ð1Þ uga ð1Þ
Vga g k " virt X virt X n m sng dmn smg  1 #! ¼ 0: ð4:62aÞ ð4:62bÞ These equations lead to the eigenvalue format: X r ja a X a smg " " virt  X ð1Þ uni n ð1Þ
Vnj þ ð1Þ unj ð1Þ
Vni #  kdij ¼ 0; # occ  X ð1Þ ð1Þ ð1Þ ð1Þ uni
Vmi þ umi
Vni  kdmn ¼ 0; i ð4:63aÞ ð4:63bÞ having the solutions equivalent with the diagonalization of the lij ¼ virt  X n ð1Þ ð1Þ ð1Þ ð1Þ uni
Vnj þ unj
Vni ð4:64aÞ matrix for the i, j belonging to the occupied orbital sets, concomitantly with the diagonalization of the mmn ¼ occ  X ð1Þ ð1Þ ð1Þ ð1Þ uni
Vmi þ umi
Vni i ð4:64bÞ block, in the basis of m, n virtuals. This procedure resembles the orbital localization, this time done in the view of vibronic analysis, arriving, by unitary transformations to the objects named vibronic orbitals (Mösch-Zanetti et al. 2000; Cimpoesu and Hirao 2003). Applying this transformation to the NH3 case, one finds most of vibronic curvature contained in a single u(1)
V(1) couple, amounting −0.775 mdyne/Å from the Kv = −0.828 mdyne/Å total. The occupied vibronic orbital is actually the same with the canonical HOMO, since in this set there are no other a2″ elements. Conversely, the virtual main vibronic orbital is sensibly different from LUMO, as the comparison drawn in Fig. 4.23 shows. In contrast with LUMO, this virtual vibronic orbital shows density contours covering the whole r skeleton, aside components resembling d-type AOs on the nitrogen atom. Then, the mechanism of pyramidal stereochemistry appears a bit more intricate than the simple hybridization idea, but it retains the same symmetry pattern. 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 359 LUMO , a2’’ Vibronic virtual, a2’’ HOMO, a1’ Vibronic occupied, a1’ Fig. 4.23 The canonical frontier MOs (left side) and vibronic orbitals (right side) of NH3 molecule considered as (a1′ + a2″) ⊗ A2″ effect in the D3h reference To deepen a bit the discussion dedicated to vibronic analysis, we will briefly consider the case of a floppy molecule, exemplified by the interesting silicon analogue of the allene. The allene is the linear H2C = C=CH2 rigid molecule, locked by successive double bonds. The silicon analogue is a floppy system (Ishida et al. 2003), due to the reluctance of this element to establish double bonds. Actually, the structure can be ascribed as {CH2((CH3)3C)2C}2Si = Si = Si{C(C(CH3) 3)2CH2}2, the silicon atoms being parts of C4Si rings, that form the terminal bonds at the ends of the Si3 row. The X-ray image of the molecule catches an averaged picture of the middle atom moving like a guitar string, taking opposite angular bends at different sides of the terminal Si


Si line, in two mutually perpendicular planes. The planes of movement can be described also as perpendicular to the CSiC planes at terminal atoms (equivalent to HCH moiety in allenes). The two CSiC planes are also mutually perpendicular, as formal consequences that the p bond planes of the concatenated double bonds are orthogonal. With the middle atom multiplied four times in the equatorial plane perpendicular to the marginal silicon atoms, the dynamic Si–Si–Si system takes the appearance of an elongated octahedron. Figure 4.24 shows the “brute force” approach to the problem, taking molecular dynamics simulation and the numeric experiment of relaxed potential energy surface at controlled bending of the Si–Si–Si line. These procedures are reproducing the floppiness of the molecule, but are not revealing explanations for this behavior. The dynamics is done with the electronic structure simulated at B3LYP/6-311G* level, and kinetic energy corresponding to room temperature, triggering the mobility of the atoms. The computed gradients at different geometry snapshots are deciding the trajectory steps, according to classical mechanics. The interval of 300 ps shown in the left upper corner of Fig. 4.24 caught about a quarter of the 360 4 Bond! Chemical Bond: Electronic Structure Methods at Work -1,178.80 E(Hartree) HF/6-311G* 0.01 B3LYP/6-311G* -1,178.83 -1,178.85 t (ps) -1,178.88 0 50 100 150 200 250 300 E(Hartree) 0 -0.01 -0.02 -0.03 -0.04 -120 -90 -60 -30 0 α(º) 30 60 90 120 Fig. 4.24 The computational illustration of the floppy Si–Si–Si moiety in the {CH2((CH3)3C)2C}2Si = Si = Si{C(C(CH3)3)2CH2}2 silicon allene, executing angular deformations around the linear configuration, in two perpendicular planes. The right side illustrates molecular (energy profile and snapshot of the molecular configurations at 0, 150, and 300 picoseconds, starting with linear form and ending with the maximal bending amplitude. The right side shows the energy profiles of relaxed (expressed as deviation from linearity in the Si–Si–Si chain) and B3LYP calculations (using the 6-311G* basis) shifted to a common origin of energy at the linear configuration periodic movement of the librational mode, from the maxima of the starting linear configuration, down nearby the minimum recorded at maximal bending amplitude. The variation of the energy is those represented with the appearance of the noise, the continuous line being the average of this trend. The rapid energy turns are due to the vibrations of the C–C and C–H skeleton, which, having frequencies higher than the Si–Si–Si bending, are repeated several times in the interval assigned to a quarter of the librational degree of freedom. The fact that shifting together the energy profiles of the HF and B3LYP calculations leads to comparable positions of the minima and heights of the hill (see right-side panel in Fig. 4.24), sustains the simplification of HF-based vibronic analysis. The results in this frame are briefly outlined in the following. In the S4 point group, the highest symmetry possible for the considered system, the instability is assignable to a (A + E) ⊗ E pseudo Jahn–Teller effect, implying in orbital part a couple of degenerate p-type MOs (in the occupied side, as will be shown immediately) degenerate and the E couple of equivalent bending in perpendicular planes. For both deformation components the vibronic analysis reveals the same quantities, since these are equivalent. Thus, the non-vibronic curvature is K0 = 0.315 mdyne/Å, encompassed in absolute value by the vibronic part, Kv = −0.483 mdyne/Å. In this case the HOMO-LUMO couple contributes rather significantly, with a u(1)
V(1) of −0.144 mdyne/Å. The orbital shapes from Fig. 4.25a.2, b.2 suggest that the active transition is from occupied p-type orbitals to a non-bonding r-type symmetry, represented for both conjugated modes by the LUMO function. This tendency for r–p mixing corresponds to the propensity of silicon to avoid the p bonding along a linear 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 361 (a) (a.1) (a.2) (b.1) (b.2) (a.3) (b) (b.3) Fig. 4.25 The angular deformation modes (libration) of the Si–Si–Si moiety in the {CH2((CH3)3C)2C}2Si = Si = Si{C(C(CH3)3)2CH2}2 silicon allene. The movements in perpendicular planes are represented in the (a) and (b) panels. The (a.1) and (b.1) represent the main characteristic of instability modes, as the off-center movement of the middle silicon atom. The (a.2) and (b.2) show the frontier orbitals coupling with the given modes (with HOMO as the components of a doubly degenerate couple and non-degenerate LUMO, common in both cases). The (a.3) and (b.3) show the couple of occupied and virtual conformation. The vibronic orbitals are concentrating in their maximal contributing couple the −0.352 mdyne/Å amount. Although not identical with the HOMO e-degenerate couple, the occupied virbronic orbitals are rather close, in shape, to this pair, while the virtual vibronic orbitals are quite different from the LUMO, illustrating a rich collection of excitations contributing to the instability. The vibronic orbitals are the functions designed as the most sensitive to a given molecular dynamic coordinate. Then, aside their utility in suggestively easing the vibronic analysis into a limited number of active couples, we can note other relevant practicalities of them. We suggest the vibronic orbitals associated to a CPHF level as useful preamble in setting the active orbitals in multi-configurational calculations in problems concerning potential energy curves or surfaces. Thus, taking the 362 4 Bond! Chemical Bond: Electronic Structure Methods at Work ammonia example, the canonical orbitals #18 and #38 are suggested as relevant if we want to consider the electron relaxation along the pyramidalization coordinate. However, in standard routines, for sure, these high orbitals will be ignored in setting the starting orbitals of a CASSCF procedure. Then, the vibronic analysis can offer important hints in setting active space in dedicated calculations. Another way to use vibronic orbitals resulting from the elementary HF-based analysis, just suggested here, without illustration attempts, concerns the detection of correlation effects. The setting will be as follows: perform geometry optimization at HF level, then by a desired correlated method (e.g. DFT), draw the atomic displacement relating the two geometries, considering it as a reaction coordinate subjected to a vibronic analysis. Then, the HF vibronic orbitals will reveal, in picturesque manner, which orbital excitations are serving as interaction channels for the correlation effects. Furthermore, vibronic orbitals can be conceived, together with the vibronic adaptation of modules performing Coupled Perturbed analytic frequency calculations for other methods, such as CASSCF or DFT frames. In this case we will have also two subsets, occupied core and virtuals, since the transformations inside active space are superfluous. The virtues of orbital formulation of vibronic analysis can be tentatively put in relation with the fact that a correct one-electron picture can in principle be achieved. First, one may invoke the Kohn–Sham orbitals, in the frame of DFT theorems. Then, a series of single excitations taken from a Slater determinant can be formulated as a new Slater determinant (Kryachko 1995), substantiating the use of orbital pictures in TD frames. Outside of DFT, one may invoke the so-called Dyson orbitals (Ortiz 1999). We did not discuss in this book the propagator theories that are backgrounding this concept. In short, the Dyson orbitals are functions expressed as the overlap between the ground states for a system in its N–1 and N electron count (their product, integrated over N–1 electrons) being therefore a one-electron function. In principle, a given state can be described as successive addition of electrons, starting from the bare nuclei, or from a collection of cores. 4.6.5 More on the Usage of Vibronic Modeling 4.6.5.1 Two State Models of Pseudo Jahn–Teller Effect Let us explore now the phenomenological side of the second-order vibronic coupling. The simplest model is a two-state Hamiltonian: 1  K þ Q2  E 2  V
Q   V
Q  ¼ 0; 1 2  2 K Q þ D  E ð4:65Þ 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 363 conceived on the basis of the states generically labeled W+ and W−, having different symmetries and not interacting in the Q = 0 reference point (where are separated by the D energy gap), but progressively coupled along the distortion coordinate, Q. For instance, the general notation (W+ + W−) ⊗ Q can stand for the (A + B) ⊗ b or (Ag+ Au) ⊗ au pseudo Jahn–Teller effects, if the C2 or Cs groups are respectively hosting a problem of this sort, or describe in effective manner the (a1′ + a2″) ⊗ A2″ situation in the D3h reference, discussed previously for the pyramidal stereochemistry of the ammonia molecule. It also can stand for the independent equivalent E displacements in a (A + E) ⊗ E mode describing the bending of a linear system, as in the previously discussed silicon allene, or of HOH molecule, if we enforce its linear configuration as the basis of discussion. Given the phenomenological intentional description, the vibronic coupling will be formalized as non-diagonal matrix element from the first-order derivative of the total Hamiltonian: _  dH  W ; V ¼ Wþ  dQ  ð4:66Þ (ignoring the technicalities exposed previously about orbital flotation, that are applied to a specific computational frame). On the diagonal of (4.65) matrix are acting the force constants associated to the harmonic vibrations along the Q mode: _ Kþ  d2 H  ¼ W þ  2 W þ ; dQ ð4:67aÞ _  d2 H  K ¼ W  2 W : dQ ð4:67bÞ Although the two quantities are different, one may suffice in practice to impose them equal: K+ = K− = K. With this assumption, the solutions of the above determinant equation are: E ðQÞ ¼  1 D þ KQ2 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D2 þ 4V 2 Q2 : ð4:68Þ For VQ << D, namely around the Q = 0 point, one may expand in series the square root term, i.e. by the (1 + x)1/2  1 + x/2 approximation, arriving at: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 D þ 4V Q ¼ D 1 þ 4V 2 Q2 =D2  D þ 2V 2 Q2 =D: The lowest energy is approximated as ð4:69Þ 364 4 Bond! Chemical Bond: Electronic Structure Methods at Work   1 2V 2 2 K E þ ðQÞ  Q þ


2 D ð4:70Þ the term in parenthesis being the new curvature of the potential energy function, made as sum of the positive former force constant with the negative vibronic part. Thus, the system is revealed as stable in symmetric form, at Q = 0, if K > 2 V2/D, while it acquires the negative curvature, signaling the distortion trend, for the K < 2 V2/D case. The instability is determined by a large coupling parameter or by a small D gap, that justifies naming the couple of involved states as quasi-degenerate. As function of the balance between the D, V, and K parameters, one may distinguish the three types of solutions depicted in Fig. 4.26. In the (a) case there is no pseudo Jahn–Teller effect, the vibronic coupling determining a smaller curvature of ground state, i.e. a parabola wider than those of the (1/2)KQ2 curve. The (b) case illustrates the pseudo Jahn–Teller effect and the realization of two minima, at ±Qmin, where the distortion amplitude is: Qmin sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V2 D2  ¼ ; K 2 4V 2 ð4:71Þ obtained from the condition dE+(Q)/dQ = 0. From the minima points, the vertical excitation energy (assigned to optical processes) is: wopt ¼ E ðQmin Þ  E þ ðQmin Þ ¼ 2 (a) (b) V2 : K ð4:72aÞ (c) E - (Q ) E - (Q ) E 0 {Ψ - (Q )} E 0 {Ψ + (Q )} E 0 {Ψ - (Q )} w opt E 0 {Ψ + (Q )} w opt ∆ ∆ E + (Q ) E + (Q ) Q Q Q -Q min w term +Q min E A (q A ) E B (q B ) Fig. 4.26 The situations of vibronic coupling in a model with two non-degenerate states. The dashed curves represent the (1/2)KQ2 harmonic oscillators of uncoupled ground and excited states. a The system is stable in symmetric state. b The active vibronic coupling determines negative curvature at the Q = 0 point and minima at distorted ±Qmin points. c Very strong vibronic coupling due to accidental degeneracy of ground and excited states (D ! 0) 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 365 Walking on the ground state, as presumable for thermally activated processes, the barrier formed at the Q = 0, between the minima, has the following height: wtherm ¼ jE þ ðQmin Þj ¼ ðKD  2V 2 Þ2 : 8KV 2 ð4:72bÞ For a practicable reaction coordinate, this barrier must be placed in the order of kBT quantity. The (c) panel in Fig. 4.26 shows the case of small D gap, when the minima get the appearance of two independent parabolas. From another perspective, it can be judged as similar to conical intersection from Jahn–Teller effect, but the symmetry conditions of this quasi-degenerate situation are, in general, essentially different from the true first-order coupling. 4.6.5.2 Vibronic Phenomenology of Mixed Valence Systems The vibronic modeling can be used to rationalize the issue of charge localization vs. delocalization in the so-called mixed valence systems. Mixed valence is a situation where the localization vs. delocalization balance of electrons (localized vs. itinerant electrons) is coupled with the movement of the nuclei (with deformations related with the change in the oxidation state of the involved coordination spheres). For exemplification, let us take a [L5M–X–ML5] binuclear, whose skeleton can be in principle symmetric, but the electron count corresponds to the formally asymmetric situation with the m oxidation state at one center and n, on the other. For the d-type transition metals, the oxidation states are usually differing by one unit (e.g. Fe(II) vs. Fe(III), Cu(I) vs. Cu(II) etc.) while for metals from the p block of the periodic table, the difference is of two electrons [e.g. Pb(II) vs. Pb(IV), Sb (III) vs. Sb(V)]. The symmetric geometry is realized in the case of the delocalized valence, when the oxidation state appears averaged over the two sites. [L5M(m+n)/2 – A L ]. This will be taken as reference geometry, irrespective whether, at X–M(m+n)/2 B 5 the end, it represents a stable or a metastable state. The geometry transformation accompanying the localization of the valence in two equivalent forms, [L5Mm A–X– L ] will be regarded as the ±Q opposed directions along MnBL5] and [L5MnA–X–Mm B 5 a normal coordinate (a vibration mode). More exactly, this is expected to be close to the anti-phase combination of the isotropic expansion-contraction movements of the two coordination sites (the breathing modes). Denoting by A and B the two sites and by q the mode approximated as the isotropic variation of the coordination bond lengths on each of them, the reaction coordinate is taken as follows: 1 Q  Q ¼ pffiffiffi ðqA  qB Þ: 2 ð4:73Þ Usually, the increase of the oxidation state results in shortening the metal-ligand bond lengths. A scheme of the localization–delocalization relationship, through the ±Q distortions, is seen in Fig. 4.27. The electronic states can be described as 366 4 Bond! Chemical Bond: Electronic Structure Methods at Work X X X +Q MAm X MBn MA(m+n)/2 MB(m+n)/2 Q=0 X -Q MAn MBm Fig. 4.27 Scheme of localized (central) versus delocalized (left and right extrema) valence in a homo-metallic [L5M-X-ML5] complex. symmetric and asymmetric superposition of the determinants concatenating the orbitals on the different sites, in the two versions of distributing the valence electrons: 1 W þ ¼ pffiffiffi ðjwA fmgwB fngj þ jwA fngwB fmgjÞ; 2 1 W ¼ pffiffiffi ðjwA fmgwB fngj  jwA fngwB fmgjÞ: 2 ð4:74aÞ ð4:74bÞ n In this definition, the pure localized form [L5Mm A–X–MBL5] results as the pffiffiffi n ðW þ þ W Þ= 2 combination, while the companion [L5MA–X–Mm B L5] comes from pffiffiffi the opposed modulation ðW þ  W Þ= 2. In general, each situation of the system can be described as a certain mixing of the defined states, the normal coordinates corresponding to the geometry distortion. The situations depicted in the (a), (b), and (c) panels of Fig. 4.26 correspond respectively to classes III, II, and I in the Robin and Day chemical classification of the mixed valence systems (Robin and Day 1967; Day et al. 2008). The Robin and Day classes are: Type I: systems in which oxidation states are clearly differentiated on the centers, so that the overall properties (optical, redox) act as a summation of the independent systems carrying coordination sites with the same collection oxidation states. (Examples: Pb3O4 with distinct PbII and PbIV sites, the composition being done as the PbO:PbO2 = 2:1 ratio; Fe3O4 with distinct FeII and FeIII sites, equivalent with the FeO:Fe2O3 = 1:1 composition; Sb2O4 with SbIII and SbV distinct oxidation states.) This corresponds to the (c) case from Fig. 4.26, where the strong coupling with the distortion trend leads to the firm localization in the individual minima. Type II: systems showing centers with differentiated crystallographic oxidation states, but their properties are different from those presented by chromophores separated in ideal distinct states. This corresponds to a mutual perturbation of the sites. (Example: the Prussian blue solid: Fe4[Fe(CN)6]3.4H2O, showing a lattice 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 367 with [FeIII(CN)6]3− and [FeII(NC)6]4− octahedral sites, the valence being partly delocalized along the FeIII–CN–FeII bridges.) This corresponds to the (b) scheme of vibronic coupling. The partial delocalization is due to a certain probability to go over the barrier, by thermal activation. For this reason, the magnetic and optical properties are tuned by external parameters. The system shows inter-valence transitions, corresponding to the processes marked wopt in the panel of Fig. 4.26b. Thus, the intense color of the Prussian blue is not due to local transitions on FeII or FeIII sites, but to a process that can be roughly described as the FeIII–CN– FeII ! FeII–CN–FeIII charge transfer. The associated vibronic paradigm can be taken also as representative for redox processes, according to Marcus-type models (Marcus 1956). Type III: the identity of oxidation states is lost in the overall delocalization, so that the centers can be described as having an average valence. [Examples: NaxWO3 metallic bronzes with x * 0.3 − 0.9 (the color turning from black to golden, through red, along with increasing x), with an itinerant electron delocalized, through and oxygen bridges, over the metal sites defined as WV and WVI averages; Fe4S4 or Fe4S3 cubane frames from ferredoxin enzymes, where the FeIII:FeII ratio can be tuned along the 1:3, 2:2, 3:1 sequence.] This corresponds to the type of Fig. 4.26a of vibronic schemes. Whether the case (a) from the panoply of Fig. 4.26 is apparently non-interesting in a general vibrational problem (i.e. having the case of a vibronic coupling that does not lead to pseudo Jahn–Teller activity), it is rather spectacular as the mixed valence type III, since it describes the apparently less intuitive case of averaged (possibly fractional) valence states. A very famous fully delocalized mixed valence system is the Creutz-Taube ion: [Ru2.5(NH3)5(Pyz)Ru2.5(NH3)5]5+ (Creutz and Taube 1969), where Pyz = pyrazine, the C4H4N2 aromatic ligand with nitrogen atoms placed in trans. The coordinate related with the differential variation of the oxidation state at one site can be perceived from the following information: (i) the distance Ru–N(Pyz) în [RuIII(NH3)5(Pyz)]3+ is by 0.07 Å smaller than in [RuII(NH3)5(Pyz)]2+; (ii) the Ru– N distance in [RuIII(NH3)6]3+ is by 0.04 Å smaller than in [RuII(NH3)6]2+. The reaction coordinate, corresponding to the shrinkage of RuIII sphere, concerted with the expansion of the RuII one has the B1u symmetry in the D2h point group of the symmetric molecule. Qualitatively, the coupling can be regarded as triggered by the HOMO!LUMO orbital excitation, whose direct product gets the symmetry of the Q coordinate: b3g ⊗ b3u = b1u. A representation of the frontier orbitals and distortion coordinate is given in Fig. 4.28. In fact, many other occupied–unoccupied pairs are participating in the coupling, other symmetry channels being: ag! b1u, b1u ! ag, au! b1g, b1g ! au, b2g ! b2u, b2u ! b2g, b3g ! b3u and b3u ! b3g. The Creutz-Taube complex presents an asymmetric optical band in near infrared at 1570 nm, assigned to the inter-valence charge transfer (IVCT), manifested at the same position in many solvents. This shows that the compound is certainly of type III in the Robin and Day scheme, the two metal centers appearing equivalent at the 10−13 s time scale of the vibrational spectroscopy. However, analogue complexes, 368 4 Bond! Chemical Bond: Electronic Structure Methods at Work x z y Ru HOMO (b3g) LUMO (b3u) N N Ru Q (B1u) Fig. 4.28 The frontier orbitals of the symmetric Creutz-Taube complex, [Ru2.5(NH3)5(Pyz) Ru2.5(NH3)5]5+, and qualitative scheme of the redox reaction coordinate (right side) as a function of the ligands, can behave as class II systems. For the type II, the inter-valence transition is dependent on temperature and solvent, since the dynamics over the barrier drives the vertical transition at different values, from the maximal wopt, taken at Qmin minima, to the minimal energy gap D, occurring from the hill of ground state at Q = 0 to the bottom of the excited state parabola. Conversely, in type III delocalized complexes, the optical transition is, in principle, always associated with the D gap, not influenced by external factors. 4.6.5.3 The Use of Vibronic Models to Fit Potential Energy Curves, Surfaces and Hyper-Surfaces In this subsection we will advocate the idea of using vibronic models as ancillary tools to construct accurate and meaningful fit functions to computed sets of potential energy data. This may have implications in various applications, such as getting handy effective potentials for the stereochemistry or reactivity of certain moieties, taken as parts of bigger systems, like proteins or other nano-scale edifices. For instance, having a good description of the energy landscape of the water molecule, as a function of its instantaneous geometry, can be relevant for further modeling of large aqua-assemblies in life sciences or in simulating the versatility of ice crystallization. Attempting to model the H2O molecule, we will consider its bent geometry as a consequence of pseudo Jahn–Teller instability of the linear configuration. We will take data from CASSCF(8,6) calculations with 6-311G* basis set. The active space aimed to formally contain the octet of electrons around the oxygen atom and six orbitals, counted as four hybrids from oxygen and two s-type hydrogen components. This setting should perform satisfactorily even in more demanding portions of general potential energy surfaces, such as bond dissociation, but the bending mode is rather non-problematic. The equal OH bond lengths of the symmetric water molecule are optimized at each considered HOH angle. It is convenient to represent the deviation from linearity as a = HOH–180°, taking positive and negative values, relative to the bending in opposite directions. The computed data are the marked points (circles) in Fig. 4.29, the same in all the panels. One of the most immediate options for fit is the polynomial expansions (or their multi-dimensional versions). However, a drawback of the polynomial approach is 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … (b) (a) 5000 -90 -60 0 0 -5000 -5000 -5000 -10000 -10000 -10000 -30 α (°) 0 30 60 90 120 α (°) -15000 -120 -90 -60 -30 0 30 60 90 120 (e) 5000 -60 0 0 -5000 -5000 -10000 -10000 -30 α (°) 0 30 60 90 120 -90 α (°) -15000 -120 -90 -12500 E (cm-1) -60 -30 0 30 60 90 120 -60 -30 E (cm-1) -13000 α (°) -15000 -120 E (cm-1) (f) 5000 E (cm-1) -15000 -90 5000 E (cm-1) 0 (d) -120 (c) 5000 E (cm-1) -15000 -120 369 0 30 60 90 120 -13500 65 75 α (°) 85 Fig. 4.29 Various fit models to the energy profile of bending the symmetric HOH molecule. The circle points are from a relaxed potential energy curve computed with CASSCF(8,6)/6-311G*. The abscissas denote the angle of deviation from linearity, a = HOH-180°, taken in opposed directions. a Fit by a {a2, a4} based polynomial, excluding the first three and last three points. b Fit by a {a2, a4} based polynomial, with all points. c Fit by a {a2, a4, a6} based polynomial, with all points. d Fit with simple vibronic model (based on K harmonic force constant, linear vibronic coupling V, and D gap), excluding the first three and last three points. e Fit with all the points, by simple vibronic model and a richer vibronic version (see the text), drawn, respectively, with green and blue continuous lines, mostly coincident. f A magnification around one minima for the same fit as in the previous panel the bad and uncontrolled behavior in the extrapolation regime. Thus, the simplest polynomial that can fit the double-well pattern is the u2a2 + u4a4 form, with u2 < 0, to yield the negative curvature at the symmetric maximum and u4 > 0, to acquire the two minima. To probe the behavior in extrapolation mode, we fit first (see panel (a) in the Fig. 4.29) excluding three points at the beginning and three others at the end of the computed sets. The considered data are reaching the minima at the symmetric sides, but do not feel much from the growing wall at extrema of the a scale. In this way, one finds that the extrapolated u2a2 + u4a4 polynomial does not recover well the trend along the missing data, rendering steeper walls. The fit with the same polynomial for the whole available dataset (namely, the computed energies for the a varying from −95° to +95° by a 5° step) yields rather good match, in spite of its simplicity (see panel (b) of the discussed figure). The average deviation on the given set is about 510 cm−1, i.e. 4.6% from the total barrier height, estimated at about 11,000 cm−1. Putting a higher order, in the attempt to improve the fit, i.e. the u2a2 + u4a4 + u6a6 polynomial, one gets the 140 cm−1 mean deviation, but with the cost of catastrophic behavior in the extrapolated extrema, as seen in panel (c) of Fig. 4.29. 370 4 Bond! Chemical Bond: Electronic Structure Methods at Work Thus, after probing the drawbacks of the polynomial procedures, we propose now the use of vibronic modeling as fitting tool. Taking the angle a as variable, we introduce a model Hamiltonian which is a richer parametric version of Eq. (4.65): ( eðaÞ ¼ Min Eig " ð2Þ ð4Þ k a a2 þ k a a4 V 1 a þ V 3 a3 V 1 a þ V 3 a3 ð2Þ 2 ð4Þ k b a þ k b a4 þ D !#) : ð4:75Þ Namely, the states are completed with quartic terms, taking different k(2) and k(4) couples for the two diagonal positions, while the linear vibronic coupling is enhanced with a cubic term. In this case we discard the excited state, taking only the ground level as minimum of the eigenvalues, as symbolized in the above equation. In first instance, we will consider the simplest version, taking only the second-order diagonal terms, equal in the two states, and the linear vibronic coupling. Namely it (2) is the (4.65) form, if we convert Q to a. Or, equivalently, take k(2) a = kb = (1/2) K and V1 = V, while cubic and quartic terms are fixed null. We fit with this model a series where the first and last three points from the computational experiment are taken out. Then, we see in panel (d) of Fig. 4.29 that the extrapolated curve copes well with this situation, intercepting rather closely the eliminated points at the left and right side extrema, definitely better than the polynomial test from panel (a). Panel (e) shows two fitted curves, superposed, basically indistinguishable at large scale. One curve (shown in green) corresponds to the vibronic fit with simplest parametric set, by K, V, and D. The achieved mean deviation is 35 cm−1, much better than those from fourth- and sixth-order polynomials, with the valuable advent of stable evolution at the extrapolated margins. The walls at large a absolute values are expected to grow steeply, since they correspond to enforce very small HOH angles, kicking the protons together. The blue curve (better distinguished in the magnification from panel (f) includes all the parameters in Eq. (4.75), reaching a deviation of about 6 cm−1. At the same time, this magnification shows that the simpler fit, with green line, goes also very close to the computed marked points. Now, we will illustrate the extension of this methodology in a bit more detail, taking into account the variation of the OH bonds, up to their breaking. In this view, we will add on the diagonal of the Eq. (4.75) a Morse-type dependence. Basically, there is no need to consider cross-terms formed as products of angular and bond length variables. It works for the whole hypersurface of the HOH molecule, namely if we take independent variation of OH bond lengths, but we will confine ourselves here to the illustration of the symmetric system, when both bonds vary in the same way, tuning the rOH value. Figure 4.30 shows the computed points and the fitted surface, remarking the success of this fit attempt, by the good match between star symbols (computed data) and circles (retrieved by model). In sections along the rOH variable, one notes the Morse-type asymmetric well, while the sections on a angle show the symmetric double-well discussed previously. At large distance, the double-well shaping is attenuated to a plateau, since distant hydrogen atoms are no longer tuning a molecular effect, the given model accounting well for this trend. 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 371 E (Hartree) rOH(Å) α(°) Fig. 4.30 The potential energy surface of symmetric HOH molecule, at the tuning of the HOH angle (representing the deviation from linearity, a = HOH-180°) and the rOH bond length. The calculations, done by CASSCF(8,6)/6-311G*, are represented by star symbols. The transparent surface and the points marked by empty circles correspond to a fit by a vibronic two-state model with a Morse-type curve incorporated on the diagonal, as a function of rOH. Note the good match of computed and fit data (the open circles retrieving well the positions of star symbols) The fact that the vibronic modeling behaves very well in the extrapolation regime and achieves a good precision, even in simple parameterization, shows that it hits the right phenomenology in the description of potential energy surfaces. Thus, in this sort of approach we do not focus on the pseudo Jahn–Teller effect in itself, but use it as a trick to account for the formation of a set of barriers and minima, considering that these should come, in last instance, from the interaction between states designed as building blocks carrying the initial minima (by the incorporated harmonic terms). With the pragmatic focus on the fitting of the ground state, we do not consider the companion excited state implied in the model, keeping it as a tacit side-support. The nominal parametric evolution of the excited state is merely a discarded by-product, but it keeps the conceptual scaffold, that the excited states are playing an effective role in stereochemical, dynamic, and reactivity problems. On the other hand, in other problems, needing the explicit account of many states, the vibronic modeling is the natural working frame. To illustrate a different case, but sharing the same modeling style, we took the interaction of a Cs atom with the ammonia molecule, showing the results in Fig. 4.31. It is not properly a chemical reaction, but it has a flavor of such process, the vibronic modeling being suitable to reproduce details of the potential energy 372 4 Bond! Chemical Bond: Electronic Structure Methods at Work E (Hartree) surface. The calculations, performed by the Coupled Cluster Single, Double, and partially Triple excitations method, CCSD(T) (with 6-311G* on ammonia and effective core potential on Cs), is a rather costly procedure. Then, the complete computational simulation may be cumbersome, having, for instance, only scattered data from different portions of the full scan, as suggested by marked points of Fig. 4.31. In such a situation, the virtues of vibronic modeling are revealed again. Complementing the calculations with a model Hamiltonian of vibronic offspring, one may recover the uncharted areas with rather good accuracy, having a very smooth behavior of the whole surface. A polynomial fit, with rkal terms will not be so well-tempered, showing undesired warping in areas far from available computed points. The surface of the considered Cs


NH3 reaction shows a Morse-type profile in the section along the rNCs variable. The section on the a out-of-planarity angle is asymmetric in the region of formed molecular complex, this being simulated with a linear term added on the diagonal of the two-state model. The asymmetry is due to the fact that, in the presence of the Cs atom at one side of the NH3 molecule, its “umbrella” type configurations are no longer equivalent, the minimum being in the case of nitrogen as the tip of the pyramid pointing toward the metal atom, establishing a sort of coordination and charge transfer bonding. The conformation with H atoms toward the metal, NH3


Cs forms the part with higher local minima, since the hydrogen moiety cannot act well for bonding contributions. At large rNCs separation, the a-driven profile trends to the symmetric double-well pattern, as can be seen on the right-side margin of Fig. 4.31. α(°) rNCs(Å) Fig. 4.31 Computed points (by CCSD(T) method, marked as star symbols) for the Cs–NH3 neutral system, with C3v symmetry, as function of N–Cs bond length and a angle measuring the deviation of NH3 moiety from planarity. The figured surface corresponds to a vibronic two-state model, noticing its excellent behavior in interpolating the scattered points and at extrapolating outskirts 4.6 Mobilis in Mobile: Electrons Moving Around Mobile … 373 Working in an imaginative manner, the vibronic models can be proposed as fitting functions of more complex potential energy surfaces, fulfilling in this way also certain heuristic needs of interpretation. The idea is that, if we have to account several minima on the energy landscape, it is technically convenient (and conceptually true) to see them as originating from the configuration interaction of several parabolas, assigned to a phenomenological set of states. As pointed out, the excited states are not of primary interest, but are ingredients improving sensibly the account of the ground surface. Being aware of this is like walking on the Earth, but being aware of the existence of the sky and the whole Cosmos, a plus of pure knowledge, paying off in the long run, starting from the very earthly pragmatic purposes. 4.7 Breaking Symmetry in Quantum Chemistry The wave–corpuscle duality in quantum mechanics is a well-known concept (see Chap. 1); however, in quantum chemistry it was the first side of the matter’s nature to prevail in bonding that is the wave function, i.e. the wave-nature of chemical bonding. However, in the wave function in chemistry it should be accompanied by a similar consolidated particle (appropriately called a bondon- see the Chap. 9 of the present book), at least for preserving the symmetry of the duality of the fundamental concept in chemistry; it is in such a context that quantum chemistry should predict the specific quantum particles carrying the chemical information in bonding, interaction, and reactivity. The present section makes such introduction, with further extension in Chap. 9. 4.7.1 The Symmetry Breaking of Chemical Field Generation The phenomenon of symmetry breaking may be naturally understood as the occurrence/emergence of stable state of the system lower than the symmetry of the interaction potential of the system (Ruckenstein and Berim 2010; Monin and Voloshin 2010). Such definition was successfully applied in the already consecrated disciplines of physics (Goldstone 1961; Nambu and Jona-Lasinio 1961; Goldstone et al. 1962; Higgs 1964a, b; Dirac 1978; Anderson 1963; Elitzur 1975; Englert and Brout 1964; Fradkin and Shenker 1979; Guralnik et al. 1964; Kibble 1967; Weinberg 1996; Frauendorf 2001), chemistry (Mikami and Yamanaka 2003; Terenziani et al. 2006; Sorkin et al. 2008; Putz 2016a; Putz and Ori 2015), and biology (Palmer 2004; Kuhn 2008). However, the last decade has witnessed the rise of graphene as the new cross-disciplinary material, worthy to be studied both for its intrinsic properties and for the emergence of new physical and chemical phenomena (Novoselov et al. 374 4 Bond! Chemical Bond: Electronic Structure Methods at Work 2012; Sheka 2012). This is because graphene brings the following (to date) unique properties for a nanomaterial (Sheka 2014): (i) the lightest material under the environmental conditions; (ii) the 2D perfection of the packed benzenoid units assuring the sp2 configuration; (iii) the mechanical strength due to the C–C high strength (1.41–1.47 Å). Nevertheless the so-called “Flagship graphene” program of realizations of “low performance” and “high performance” applications relies just on the degree to which the physical properties allow for chemical realizations (e.g. touch screens, rollable e-paper, foldable organic light emitting diodes—OLED, tunable sensors, solar cells, etc.) while just the chemistry of graphene resists physical applications (electronic devices like semimetal-semiconductors) on the other side, respectively. The explications are based on the peculiar chemical features of the above structural properties, namely (Putz et al. 2016, Putz 2016a,b): • The topochemical character of graphene mechanics (Bissett et al. 2013), rendering—for instance—the mechanical shocks by molecular anisotropic propagation (Long et al. 2015) according to the topological defects anisotropy induced by Stone-Wales waves (Ori et al. 2011); • The radical character of the graphene relates with the so-called “edge problem” for the carbon skeleton of the fixed structure at the level of odd electrons originating in the number of atoms in the structure added to those from its edge; moreover, towards obtaining graphenic nanoribbons (perfect 2D planar single-layer structures) by cutting graphenic sheets additional dangling bonds are generated and thus the unpaired electrons, all of these enhancing the graphenic radical properties, with the consequence of the chemical reaction proportion on the ribbon circumference with possible reorganization of the pristine graphene structure; this also relates with the atomic chemical susceptibility in relation to graphene magnetism (Sheka 2013) and conductivity (Sheka 2007); • The collective behavior character of the graphene is a direct consequence of the above “odd electrons”, as well as due to external actions of electric/magnetic fields along any chemical modification (e.g. hydrogenations (Elias et al. 2009), photoexcitation or mechanical loading; it is manifested by the sp2 ! sp3 transformation of the carbon valence electrons, destroying the 2D planar geometry. This effect is naturally accompanied by the C–C bond length redistribution, so preventing the knowledge about the atom/region (or quantum dot) caring the highest reactivity of graphene sheet and calling for the next deposition/interaction target; actually with C–C bonding growing up to 1.8–2 Å the densely packed sp2 honeycomb crystal changes due to the correlating electrons, due to the occurred radicalization so that the collective excited behavior is manifested (Staroverov and Davidson 2000); actually the benzenoid carcass of pristine graphene changes to cyclohexanoid units with high heteromorphic configuration—as observed with high resolution transmission electronic microscopy (HRTEM) providing more natural pictures than the recently proposed artificial “gigantic pseudomagnetic field” (Levy et al. 2010); nevertheless, such electronic p-correlation may be inhibited by deposition of monolayer graphene on substrates as such silicon carbide, boron nitride, and quartz 4.7 Breaking Symmetry in Quantum Chemistry 375 surfaces (Hwang et al. 2012), or by using the rarely yet regularly distributed nanoparticles grid (Wu et al. 2013), or even by designing the extreme case of graphene–substrate composites by adsorbed carbon monolayers of hexagon patterning (Hoffmann 2013), which assures the free standing of graphene by seeking the best surface science partners so inhibiting the carbon atoms radicalization and preserving the hexagonal-packed monolayers (Agnoli and Granozzi 2013). For the analytical developments, one starts with the working field potential under double-well form, adapted for the electronegativity-chemical hardness global potential nature of chemical interaction (Putz 2008, 2016a, b; Putz and Ori 2015; Putz et al. 2016): Vð/Þ ¼ l/2 þ 1 4 g/ ; 2 ð4:76Þ to generate bonding fields and its particles—the bondons—by changing from the upper (positive potential) branch to lower (negative potential) branch of the first-order particle ( /2 ) potential. Figure 4.32 represents graphically the parabolic dependence with v > 0, η > 0 (see the upper dashed curve) presenting the minimum zero potential for the vanishing field / = 0 (so for fermions when v1 identified with positive chemical potential of the system). Instead, a completely different picture is obtained if the same potential is considered with l < 0 (or −l = v > 0), l > 0 (so for the bosons for negative chemical potential) when two distinct non-zero minimum potential Fig. 4.32 The potential of Eq. (4.76) with +l (aka fermionic state for positive chemical potential) in dashed curve and with −v (aka bosonic states for negative potential) on continuous curve, illustrating the symmetry breaking of the parabolic (upper curve) to double-well potential (lower curve) that includes also the negative (vacuum, yet more stable) quantum states (Putz 2008, Putz 2016a, b; Putz and Ori 2015; Putz et al. 2016) 376 4 Bond! Chemical Bond: Electronic Structure Methods at Work values appear in its negative (vacuum) region for the chemical field acquiring the respective values (Putz 2008, Putz 2016a, b; Putz and Ori 2015; Putz et al. 2016): @V ¼ 0 ) /a;b ¼  @/ rffiffiffiffiffiffiffi l : g ð4:77Þ Solutions of Eq. (4.77) thus largely justify the bosonic appearance for the activation of the spontaneous symmetry breaking. It is worth going from the fermionic potential driven by the positive chemical potential (+l) to the bosonic potential driven by electronegativity (−l = v). Through such a “path” the chemical field naturally identifies with the chemical bonding field by shifting the minimum zero potential to its minimum negative range, in the quantum vacuum region from where the quantum particles are to be spontaneous created, namely the bondon and anti-bondon as the quantum particles of the chemical bonding fields. We arrive at the fundamental statements: the chemical bonding as described by the joint bonding–antibonding states has compulsory an associated bosonic character, the negative chemical potential, i.e. driven by positive electronegativity! Next, one likes employing this phenomenological analysis to analytically determine the bondonic mass through the quantum creation by the symmetry breaking mechanism (Putz 2008). To this aim, one considers the Lagrangian of the Schrödinger field ð/Þ (Putz 2016a): 2 h L ¼ i h/ /_  ðr/ Þðr/Þ  V/ /; 2m0 ð4:78Þ produced by the actual potential V Vð/Þ ¼ v/2 þ 1 4 g/ ; 2 ð4:79Þ by connecting the chemical field ð/Þ with the parabolic expansion of the chemical reactive/valence energy (Parr et al. 1978; Putz 2011). Accordingly, the stationary solutions of Eq. (4.77) become for the bondonic fields: /a;b rffiffiffi v : ¼ g ð4:80Þ It is worth remarking that the result (4.80) gives the basic indication that the chemical filed (and therefore also the related electronic density) directly depends on the chemical reactivity indices electronegativity and chemical hardness—so paving the route to analytically unfold the inverse density problem of chemical field theory and subsequent chemical bonding and reactivity modeling. 4.7 Breaking Symmetry in Quantum Chemistry 4.7.2 377 The Inverse Quantum Chemical Problem With the developments of the Density Functional Theory (DFT) of atoms and molecules, also the chemical reactivity in general and reactivity indices of electronegativity (Parr et al. 1978; Parr and Yang 1989; Kohn et al. 1996; Islam and Ghosh 2011; Putz 2016b, c):   dE½q v ¼ l ¼  ; dq q¼qðVÞ ð4:81Þ and chemical hardness (Pearson 1997): "  #    @v @ dE½q g¼ ¼ ; @N V @N dq q¼qðVÞ V ð4:82Þ were accordingly reformulated, at the conceptual level, as the density functionals relating the first- and the second-order derivatives of the total (or valence) energy functional respecting the electronic density itself, respectively. Moreover, in this framework, also the total (or valence) number of electrons also writes as the density functional through the integral: Z N ¼ qðrÞdr: ð4:83Þ Nevertheless, the main point is the apparent fundamental density as the basal concept and tool for the (observable or not) measures of chemical reactivity as electronegativity and chemical hardness of Eqs. (4.81) and (4.82). However, as the Hohenberg–Kohn Theorem stipulates (Hohenberg and Kohn 1964), the electronic density and applied (or effective) potential upon the system stay in a biunivoque relationship. Consequently, one may advance also the inverse density problem, according to which it is the effects of the chemical organization of electrons in atoms, molecules, or chemical bond; formally, such connection is apparent from the energetic functional derivative (Parr and Gázquez 1993): qðrÞ ¼  dE d VðrÞ  N : ð4:84Þ The fact that electronic density “follows/runs” for the best energetic in an electronic system, according with conceptual and computational DFT, is pregnant, for instance, in the constrained-search optimization formalism of Levy (1982). The Levy mechanism assures, via electronic density, the manifestation of the ground state energy (the same for valence case) out of the class of energies as density functionals mapped from the correspondent class of densities: 378 4 Bond! Chemical Bond: Electronic Structure Methods at Work EGS ¼ min min hWjðT þ Vee þ VÞjWi q W!q Z ¼ min min hWjðT þ Vee ÞjWi þ qðrÞVðrÞdr q W!q ¼ minðFHK ½q þ CA ½qÞ ð4:85Þ q ¼ minðE½qÞ: q Therefore one can explore also the potential measures determining the electronic density VðrÞ ! qðrÞ; ð4:86Þ which in terms of chemical reactivity (various orders) potentials, as electronegativity and chemical hardness are by their constructions, means looking for the analytical relationships of the general form qðrÞ ¼ f ðv; g; rÞN;V : ð4:87Þ For the practical electronic density–chemical reactivity indices relationship within the Feynman-Kleinert formalism (see Chap. 1, Sect. 1.6.2), one begins with the evaluation of the smeared out potential (1.227) for the general harmonic and anharmonic potential: Vjv;gi ð/Þ ¼ v/2 þ 1 4 g/ ; 2 ð4:88Þ leaving with the respective result: Va2 ð/0 Þ ð/0 Þ ¼ Zþ 1 1 yielding 1 Va2 jv;gi ð/0 Þ ¼ 2 d/ ð/  /0 Þ2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vð/Þ exp  2a2 ð/0 Þ 2pa2 ð/0 Þ " 3g a4 ð/0 Þ  2v/20 þ g/40 þ 2a2 ð/0 Þð6g/20  2vÞ # ð4:89Þ ð4:90Þ The potential (4.90) is to be further specialized within the Markovian limit (1.233) to become W1 ð/0 Þb!0 ffi Va2 ! b ð/0 Þ 3 2 12 1 2 2 g b  2v/ 0 7 1 6 48 7; ¼ 6 5 4 2 1 2 4 bð6g/0  2vÞ þ g/0 þ 12 ð4:91Þ 4.7 Breaking Symmetry in Quantum Chemistry 379 and, together with the whole approximated potential (1.239), will release the quantum statistical partition functions (1.231), respectively as: Zjv;gi Zþ 1 1 ¼ pffiffiffiffiffiffiffiffi d/0 exp½bW1 ð/0 Þ 2pb 1 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # 1 gb  4v bðgb  4vÞ2 K1 ¼ 4 p gb 4 64g "  # b 48v2  8vgb þ g2 b4  exp ; 192g ð4:92Þ being K[] the modified Bessel function of the second rank. With the help of the partition functions (4.92), in the Markovian limit, the associated electronic densities (1.230) work with the potential (4.79) as: qjv;gi ð/Þ ffi Zj1 v;gi /4 exp b v/ þ g 2  2  : ð4:93Þ Then, the expression for the total number of electrons, according to Eq. (4.83) is given: Njv;gi ð/Þ ¼ Zþ 1 1 ¼2 qjv;gi ð/Þd/ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2pb v  4v½Ha1  bg h 2i K14 b4gv K14 bð½Hagb4vÞ 64g  exp  ð4:94Þ 2 b2 ½Hað½Hagb  8vÞ ; 192 with the Hartree unit energy in electron-Volts, [Ha] = 27.11 eV so adjusting the correct units. Note that one can control the system’s temperature by the numerical casting of the inverse thermal energy parameter (NIST 2015): b¼ 1 11604:505 ½K/eV : ffi kB T T ½K ð4:95Þ 380 4 Bond! Chemical Bond: Electronic Structure Methods at Work Nevertheless, the fully chemical field representation of the electronic density (4.93) requires knowledge of the chemical field itself. The chemical field works as the first-order b-expansion, which interestingly (yet meaningfully) depends only on chemical hardness influence (Parr et al. 1978; Parr and Yang 1989; Kohn et al. 1996; Islam and Ghosh 2011; Putz 2016b, c): /¼e ar 0 pffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 2g½Ha 1 þ 3bge 4ar 0 pffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 2g½Ha : ð4:96Þ Note that the exponential factor in (4.96) was considered normalized to the nominal value of the Bohr radius, taken without dimension \a0 [ ¼ 0:529, however at the atomic scale (Ångstroms) so preserving the physical nano-bonding information inside. We have now solved analytically the inverse density problem in terms of electronegativity and chemical hardness. Next, we consider the number of electrons (4.94) and the chemical field density (4.93)-with-(4.96) to produce Figs. 4.33 and 4.34, for a working case of carbon–carbon simple bond in a saturated environment (Tudoran and Putz 2015), see the upper inset of Fig. 4.34, with respective chemical reactivity indices, v ¼ 8:9725 ½eV and g ¼ 7:31879 ½eV, yielding the following interpretation: • By inspection of Fig. 4.33 it is noted that while for asymptotic high temperature the number of electrons is restricted on valence states (driven by electronegativity and chemical hardness frontier chemical reactivity indices)—as in the top and bottom “N” representations—as the temperature decreases the electronic circulation increases (naturally) yet also the income/flux of electrons in the system increases as well, so turning the “aromaticity/fermionic” behavior into a collective/bosonic (bondonic) one (Putz 2010). • Figure 4.34 confirms on a temperature-radii representation (upper) and contour projection (bottom) the fermionic chemical bond manifested at the limited (still not asymptotic) high temperature limit yet confined (not aromatic) at certain radius limit from one reference carbon atom in carbon–carbon bonding, according with earlier analysis of bondonic representation of the chemical bond (Putz 2010). All in all, the maximum of two electrons limit in the bottom of Fig. 4.33 is still not acquired over a range of chemical reactivity electronegativity and chemical hardness indices for both simple and double carbon-carbon bondings in saturated adjacency valence electrons (Tudoran and Putz 2015). The present approach finely accords the quantum orbital spin-restriction, i.e. Pauli principle, with the fermionic branch of the general inharmonic potential for the chemical field, see Eq. (4.88) and 4.7 Breaking Symmetry in Quantum Chemistry 381 Fig. 4.33 The numbers of electrons (4.94), with numerical electronegativity and chemical hardness values v ¼ 8:9725 ½eV and g ¼ 7:31879 ½eV, associated to the open molecular electronic of single carbon– carbon bond in saturated valence adjacency (see the inset of Fig. 4.34), in the asymptotic temperature limit (1.233)—on top “N” representation, for decreasing the high temperature limit—in the middle “NN” representation, and on the carbon–carbon chemical bonding general range (Tudoran and Putz 2015) of electronegativity and chemical hardness dependency—in the bottom representation, respectively (Putz 2016b) discussion around Fig. 4.31. This situation can be accordingly fully turned into the bondonic-bosonic inverse problem of the chemical bond by reconsidering the sign of electronegativity in Eq. (4.88), according with the interpretation of Fig. 4.31 and with recent findings (Putz et al. 2016). 382 4 Bond! Chemical Bond: Electronic Structure Methods at Work Fig. 4.34 Identification of the chemical bonding apex and its close contour (bottom representation) along the radii of the single carbon–carbon bond in adjacency (top inset) and of the high temperature limit (4.92)–(4.93) for the electronic density in FKP (Feynman-Kleinert-Putz) chemical field approach driven by associated electronegativity (v ¼ 8:9725 ½eV) and chemical hardness (g ¼ 7:31879 ½eV) (Tudoran and Putz 2015; Putz 2016b) 4.8 Conclusions The chapter offered the following items of conceptual and technical know-how: • Explaining various computation methods. • Assessing the calculation results (e.g. on optimized molecular geometry). 4.8 Conclusions 383 • Comparison with experiment (e.g. molecular geometry from diffraction data). • Critical stand in comparative approach (experiment vs. simulation or inter-method outcome). • Writing input files for various calculation methods (HF, DFT, CASSCF, VB), with different codes: GAMESS, Gaussian, ADF, VB2000. • Explaining the action of several keywords, along with the exemplified methods and codes. • Hinting at non-usual controls: e.g. permutation of orbital list for ionization potentials, handling fractional population numbers in DFT with ADF, conducting a Broken Symmetry calculation, customizing the initial orbitals in CASSCF and VB procedures, driving the VB resonance structures in accordance with the chemical meaning. • Understanding the phenomenology of two-electron chemical bond as spin coupling. • Developing model Hamiltonians (e.g. Heisenberg–Dirac–van Vleck) spin coupling and use as ancillary tools in interpreting computational output. • Defining directed bonds, lone pairs, resonance structures in VB calculations. • Solving phenomenological VB calculations with graphic rules. • Interpreting calculation results with model Hamiltonians. • Solving the problems of the interaction between vibration and electronic degrees of freedom (vibronic coupling) and multiple facets of manifestations. • Introducing the terminology of vibronic coupling prototypic cases, the Jahn– Teller and pseudo Jahn–Teller effects as instability of a molecular configuration. • Enumerating the complex aspects of molecular dynamics and the incidence of vibronic coupling in rationalizing stereochemical and reactivity problems. • Exercising concrete vibronic modeling on simple cases. • Clarifying the computational implementation of vibronic analysis, in the frame of Coupled Perturbed (CP) techniques, using the time-dependent sets of excited states or equivalent formulation in terms of orbital promotion couples. • Introducing the concept of vibronic orbitals and showing illustrative examples. • Using vibronic language to rationalize different manifestations such as classification of mixed valence systems (Robin and Day scheme). • Hinting at the use of phenomenological vibronic modeling as an accurate and meaningful tool for fitting potential energy surfaces. • Conceptually inverting the hierarchy allowed in chemistry by the celebrated Hohenberg–Kohn density-applied potential bijection for ground/valence state. • Unfolding the previous idea in an analytical manner by extending the quantum statistical path integral formalism of Feynman and Kleinert to the first-order temperature parameter high expansion of the chemical field as appeared by breaking symmetry in chemical field generation. • Setting the chemical field by the quantum applied potential in an anharmonic manner, in accordance with the parabolic electronic system’s energy by electronegativity and chemical hardness. 384 4 Bond! Chemical Bond: Electronic Structure Methods at Work References ADF2012.01, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam [http://www.scm.com] Agnoli S, Granozzi G (2013) Second generation graphene: opportunities and challenges for surface science. Surf Sci 609:1–5 Anderson PW (1959) New approach to the theory of superexchange interactions. Phys Rev 115:2–13 Anderson PW (1963) Plasmons, gauge invariance, and mass. Phys Rev 130:439–442 Atanasov M, Reinen D (2001) Density functional studies on the lone pair effect of the trivalent group (V) elements: I. electronic structure, vibronic coupling, and chemical criteria for the occurrence of lone pair distortions in AX3 molecules (A = N to Bi; X = H, and F to I). J Phys Chem A 105:5450–5467 Atanasov M, Reinen D (2002) Predictive concept for lone-pair distortions: DFT and vibronic model studies of AXn-(n–3) molecules and complexes (A = NIII to BiIII; X = F–I to I–I; n = 3–6). J Am Chem Soc 124:6693–6705 Bencini A, Totti F, Daul CA, Doclo K, Fantucci P, Barone V (1997) Density functional calculations of magnetic exchange interactions in polynuclear transition metal complexes. Inorg Chem 36(22):5022–5030 Bersuker IB (1984) The Jahn–Teller effect and vibronic interactions in modern chemistry. Plenum Press, New York Bersuker IB (2001) Modern aspects of the Jahn–Teller effect theory and applications to molecular problems. Chem Rev 101:1067–1114 Bersuker IB (2013) Pseudo-Jahn–Teller effect: a two-state paradigm in formation, deformation, and transformation of molecular systems and solids. Chem Rev 113:1351–1390 Bissett MA, Konabe S, Okada S, Tsuji M, Ago H (2013) Enhanced chemical reactivity of graphene induced by mechanical strain. ACS Nano 7:10335–10343 Bode BM, Gordon MS (1998) Macmolplt: a graphical user interface for GAMESS. J Mol Graph Mod 16(3):133-138 Car R, Parrinello M (1985) Unified approach for molecular dynamics and density-functional theory. Phys Rev Lett 55:2471–2474 Ceulemans A, Chibotaru LF, Cimpoesu F (1997) Intramolecular charge disproportionation and the band structure of the A3C60 superconductors. Phys Rev Lett 78:3275–3728 Chibotaru LF, Cimpoesu F (1997) Vibronic instability of molecular configurations in the Hartree– Fock-Roothan method. Int J Quant Chem 65:37–48 Cimpoesu F, Hirao K (2003) The ab initio analytical approach of vibronic quantities: application to inorganic stereochemistry. Adv Quant Chem 44:370–397 Cooper DL, Thorsteinsson T, Gerratt J (1999) Modern VB representations of CASSCF wave functions and the fully-variational optimization of modern VB wave functions using the CASVB strategy. Adv Quantum Chem 32:51–67 Creutz C, Taube H (1969) Direct approach to measuring the Franck-Condon barrier to electron transfer between metal ions. J Amer Chem Soc 91:3988–3989 Day P, Hush NS, Clark JR (2008) Mixed valence: origins and developments. Phil Trans R Soc A 366:5–14 Dirac, PAM (1978) Mathematical foundations of quantum theory. In: Marlow A (ed). Academic Press, New York Dutuit O, Tabche-Fouhaile A, Nenner I, Frohlich H, Guyon PM (1985) Photodissociation processes of water vapor below and above the ionization potential. J Chem Phys 83:584–596 Elias DC, Nair Mohiuddin TMG, Morozov SV, Blake P, Halsall MP, Ferrari AC, Boukhvalov DW, Katsnelson MI, Geim AK, Novoselov KS (2009) Control of graphene’s properties by reversible hydrogenation: evidence for graphene. Science 323:610–613 Elitzur S (1975) Impossibility of spontaneously breaking local symmetries. Phys Rev D 12:3978–3982 Englert F, Brout R (1964) Broken symmetry and the mass of gauge vector mesons. Phys Rev Lett 13:321–323 Epstein ST (1974) The variation method in quantum chemistry. Academic Press, New York References 385 Eyring H (1931) The energy of activation for bimolecular reactions involving hydrogen and the halogens, according to the quantum mechanics. J Am Chem Soc 53:2537–2549 Flükiger P, Lüthi HP, Portmann S, Weber J (2000-2002) MOLEKEL version 4.3. Swiss Center for Scientific Computing, Manno (Switzerland) Fonseca Guerra C, Snijders JG, te Velde G, Baerends EJ (1998) Towards an order-N DFT method. Theor Chem Acc 99:391–403 Fradkin E, Shenker SH (1979) Phase diagrams of lattice gauge theories with Higgs fields. Phys Rev D 19:3682–3697 Frauendorf S (2001) Spontaneous symmetry breaking in rotating nuclei. Rev Mod Phys 73:463–514 Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Scalmani G, Barone V, Mennucci B, Petersson GA, Nakatsuji H, Caricato M, Li X, Hratchian HP, Izmaylov AF, Bloino J, Zheng G, Sonnenberg JL, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Vreven T, Montgomery JA Jr, Peralta JE, Ogliaro F, Bearpark M, Heyd JJ, Brothers E, Kudin KN, Staroverov VN, Kobayashi R, Normand J, Raghavachari K, Rendell A, Burant JC, Iyengar SS, Tomasi J, Cossi M, Rega N, Millam JM, Klene M, Knox JE, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Martin RL, Morokuma K, Zakrzewski VG, Voth GA, Salvador P, Dannenberg JJ, Dapprich S, Daniels AD, Farkas Ö, Foresman JB, Ortiz JV, Cioslowski J, Fox DJ (2009) Gaussian 09. Gaussian Inc, Wallingford CT Goldstone J (1961) Field theories with superconductor solutions. Nuovo Cim 19:154–164 Goldstone J, Salam A, Weinberg S (1962) Symmetry groups in nuclear and particle physics. Phys Rev 127:965–970 Guralnik GS, Hagen CR, Kibble TWB (1964) Global conservation laws and massless particles. Phys Rev Lett 13:585–587 Heisenberg WZ (1928) Zur theorie des ferromagnetismus. Z f Phys 49:619 Higgs PW (1964a) Broken symmetries, massless particles and gauge fields. Phys Lett 12:132–133 Higgs PW (1964b) Broken symmetries and the masses of gauge bosons. Phys Rev Lett 13:508– 509 Hirao K, Nakano H, Nakayama K, Dupuis M (1996) A complete active space valence bond (CASVB) method. J Chem Phys 105:9227–9239 Hoffmann R (2013) Small but strong lessons from chemistry for nanoscience. Angew Chem Int Ed 52:93–103 Hohenberg P, Kohn W (1964) Inhomogeneous electronic gas. Phys Rev 136:864–871 Hwang C, Siegel DA, Mo SK, Regan W, Ismach A, Zhang Y, Zettl A, Lanzara A (2012) Fermi velocity engineering in graphene by substrate modification. Sci Rep 2:590/1–4 Ishida S, Iwamoto T, Kabuto C, Kira M (2003) A stable silicon-based allene analogue with a formally sp-hybridized silicon atom. Nature 421 (6924):725-727 Islam N, Ghosh DC (2011) The electronegativity and the global hardness are periodic properties of atoms. J Quantum Info Sci 1:135–141 Kibble TWB (1967) Symmetry breaking in non-abelian gauge theories. Phys Rev 155:1554–1561 Kirchner MT, Das D, Boese R (2008) Cocrystallization with acetylene: molecular complex with methanol. Cryst Growth Des 8:763–765 Kohn W, Becke AD, Parr RG (1996) Density functional theory of electronic structure. J Phys Chem 100:12974–12980 Kryachko ES (1995) On wave-function correction to Hellmann-Feynman force: Hartree–Fock method. Int J Quant Chem 56:3–7 Kuhn H (2008) Origin of life—symmetry breaking in the universe: emergence of homochirality. Curr Opin Colloid Interface Sci 13:3–11 Kutzelnigg W (1989) Ab initio calculation of molecular properties. J Mol Struct (Theochem) 202:11–61 Levy M (1982) Electron densities in search of Hamiltonians. Phys Rev A 26:1200–1208 Levy N, Burke SA, Meaker KL, Panlasigui M, Zettl A, Guinea F, de Castro Neto AH, Crommie MF (2010) Strain-induced pseudo-magnetic fields greater than. Science 329:544–547 386 4 Bond! Chemical Bond: Electronic Structure Methods at Work Li J, Duke B, McWeeny R (2007) VB2000 Version 2.0, SciNet Technologies, San Diego, CA Li J, McWeeny R (2002) VB2000: pushing valence bond theory to new limits. Int J Quantum Chem 89:208–216 Long X, Zhao F, Liu H, Huang J, Lin Y, Zhu J, Luo S-N (2015) Anisotropic shock response of stone wales defects in graphene. J Phys Chem C 119:7453–7460 Mahapatra S, Köppel G (1998) Quantum mechanical study of optical emission spectra of Rydberg-excited H3 and its isotopomers. Phys Rev Lett 81:3116–3120 Marcus RA (1956) On the oxidation-reduction reactions involving electron transfer. J Chem Phys 24:966–978 McWeeny R (2001) Methods of molecular quantum mechanics. Academic Press, London Mikami K, Yamanaka M (2003) Strain-induced pseudo-magnetic fields greater than 300 Tesla in graphene nanobubbles. Chem Rev 103:3369–3400 Mistrík I, Reichle R, Helm H, Müller U (2001) Predissociation of H3 Rydberg states. Phys Rev A 63:042711 Mitani M, Mori H, Takano Y, Yamaki D, Yoshioka Y, Yamaguchi K (2000) Density functional study of intramolecular ferromagnetic interaction through m-phenylene coupling unit (I): UBLYP, UB3LYP and UHF calculation. J Chem Phys 113:4035–4050 Monin A, Voloshin MB (2010) Spontaneous and induced decay of metastable strings and domain walls. Ann Phys 325:16–48 Mösch-Zanetti NC, Roesky HW, Zheng W, Stasch A, Hewitt M, Cimpoesu F, Schneider TR, Prust J (2000) The first structurally characterized aluminum compounds with terminal acetylide groups. Angew Chem Int Ed 39:3099–3101 Nakano H, Sorakubo K, Nakayama K, Hirao K (2002) In: Cooper DL (ed) Valence Bond Theory. Elsevier Science, Amsterdam, pp 55–77 Nambu J, Jona-Lasinio G (1961) Dynamical model of elementary particles based on an analogy with superconductivity. Phys Rev 122:345–358 NIST (2015) National Institute for Science and Technology (NIST) on Constants, Units, and Uncertainty. http://physics.nist.gov/cuu/Constants/index.html Noodleman L (1981) Valence bond description of antiferromagnetic coupling in transition metal dimmers. J Chem Phys 74:5737–5743 Noodleman L, Davidson ER (1986) Ligand spin polarization and antiferromagnetic coupling in transition metal dimmers. Chem Phys 109:131–143 Novoselov KS, Fal’ko VI, Colombo L, Gellert PR, Schwab MG, Kim K (2012) A roadmap for graphene. Nature 490:192–200 Onishi T, Takano Y, Kitagawa Y, Kawakami T, Yoshioka Y, Yamaguchi K (2001) Theoretical study of the magnetic interaction for M–O–M type metal oxides: comparison of broken-symmetry approaches. Polyhedron 20:1177–1184 Ori O, Cataldo F, Putz MV (2011) Topological anisotropy of stone-wales waves in graphenic fragments. Int J Mol Sci 12:7934–7949 Ortiz JV (1999) Toward an exact one-electron picture of chemical bonding. Adv Quantum Chem 35:33–52 Palmer AR (2004) Symmetry breaking and the evolution of development. Science 306:828–833 Parr RG, Donnelly RA, Levy M, Palke WE (1978) Electronegativity: the density functional viewpoint. J Chem Phys 68:3801–3807 Parr RG, Gázquez JL (1993) Hardness functional. J Phys Chem 97:3939–3940 Parr RG, Yang W (1989) Density functional theory of atoms and molecules. Oxford University Press, New York Pauling L, Wheland GW (1933) The nature of the chemical bond. V. The quantum-mechanical calculation of the resonance energy of benzene and naphthalene and the hydrocarbon free radicals. J Chem Phys 1362 Pearson RG (1997) Chemical hardness. Wiley-VCH, Weinheim Polanyi JC, Wong WH (1969) Location of energy barriers. I. Effect on the dynamics of reactions A + BC. J Chem Phys 51:1439–1450 References 387 Pulay P (1987) Analytical derivative methods in quantum chemistry. In: Lawley KP (ed) Ab Initio Methods in Quantum Chemistry. John Wiley, New York Putz MV (2008) The chemical bond: spontaneous symmetry-breaking approach. Symmetry: Cult Sci 19:249–262 Putz MV (2010) The bondons: the quantum particles of the chemical bond. Int J Mol Sci 11:4227– 4256 Putz MV (2011) Electronegativity and chemical hardness: different patterns in quantum chemistry. Curr Phys Chem 1:111–139 Putz MV (2016a) Quantum nanochemistry: a fully integrated approach. Vol 3: quantum molecules and reactivity. Apple Academic Press & CRC Press, Toronto Putz MV (2016b) Chemical field theory: the inverse density problem of electronegativity and chemical hardness for chemical bond. Curr Phys Chem 7(2):133-146. June 2017 doi: 10.2174/ 1877946806666160627101209 Putz MV (2016c) Quantum nanochemistry: a fully integrated approach. Vol II: quantum atoms and periodicity. Apple Academic Press & CRC Press, Toronto Putz MV, Ori O (2015) Predicting bondons by Goldstone mechanism with chemical topological indices. Int J Quantum Chem 115:137–143 Putz MV, Ori O, Diudea M, Szefler B, Pop R (2016) Bondonic chemistry: spontaneous symmetry breaking of the topo-reactivity on graphene. In: Ali Reza Ashrafi, Mircea V Diudea (eds) Chemistry and physics: distances, symmetry and topology in carbon nanomaterials. Springer, Dordrecht, pp 345–389 Raimondi M, Cooper DL (1999) Ab initio modern valence bond theory. In: Surján PR, Bartlett RJ, Bogár F, Cooper DL, Kirtman B, Klopper W, Kutzelnigg W, March NH, Mezey PG, Müller H, Noga J, Paldus J, Pipek J, Raimondi M, Røeggen I, Sun JQ, Surján PR, Valdemoro C, Vogtner S (eds) Topics in current chemistry: localization and delocalization, vol 203. Reidel, Dordrecht, pp 105–120 Robin MB, Day P (1967) Mixed-valence chemistry: a survey and classification. Adv Inorg Chem Radiochem 10:247–422 Ruckenstein E, Berim GO (2010) Contact angle of a nanodrop on a nanorough solid surface. Adv Colloid Interface Sci 154:56–76 Ruiz E, Cano J, Alvarez S, Alemany P (1999) Broken symmetry approach to calculation of exchange coupling constants for homobinuclear and heterobinuclear transition metal complexes. J Comp Chem 20:1391–1400 Sato S (1955) On a new method of drawing the potential energy surface. J Chem Phys 23:592 Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, Koseki S, Matsunaga N, Nguyen KA, Su SJ, Windus TL, Dupuis M, Montgomery JA (1993) General atomic and molecular electronic structure system. J Comput Chem 14:1347–1363 Sheka EF (2007) Odd electrons in molecular chemistry, surface science, and solid state magnetism. Int J Quantum Chem 107:2935–2955 Sheka EF (2012) Computational strategy for graphene: insight from odd electrons correlation. Int J Quantum Chem 112:3076–3090 Sheka EF (2013) In: Hetokka M, Brandas E, Maruani J, Delgado-Barrio G (eds) Progress in theoretical chemistry and physics, vol 27, pp 249–284 Sheka EF (2014) The uniqueness of physical and chemical natures of graphene: their coherence and conflicts. Int J Quantum Chem 114:1079–1095 Song L, Chen Z, Ying F, Song J, Chen X, Su P, Mo Y, Zhang Q, Wu W (2012) XMVB 2.0: an ab initio non-orthogonal valence bond program. Xiamen University, Xiamen 361005, China Song L, Mo Y, Zhang Q, Wu W (2005) XMVB: a program for ab initio nonorthogonal valence bond computations. J Comput Chem 26:514–521 Sorkin A, Iron MA, Truhlar DG (2008) Density functional theory in transition-metal chemistry: relative energies of low-lying states of iron compounds and the effect of spatial symmetry breaking. J Chem Theory Comput 4:307–315 Staroverov VN, Davidson ER (2000) Distribution of effectively unpaired electrons. Chem Phys Lett 330:161–168 388 4 Bond! Chemical Bond: Electronic Structure Methods at Work te Velde G, Bickelhaupt FM, van Gisbergen SJA, Fonseca Guerra C, Baerends EJ, Snijders JG, Ziegler T (2001) Chemistry with ADF. J Comput Chem 22:931–967 Terenziani F, Painelli A, Katan C, Charlot M, Blanchard-Desce M (2006) Chromophores: symmetry breaking and solvatochromism. J Am Chem Soc 128:15742–15755 Thorsteinsson T, Cooper DL (1998) Modern valence bond descriptions of molecular excited states: an application of CASVB. Int J Quant Chem 70:637–650 Tudoran MA, Putz MV (2015) Molecular graph theory: from adjacency information to colored topology by chemical reactivity. Curr Org Chem 19:358–385 van Vleck JH, Sherman A (1935) The quantum theory of valence. Rev Mod Phys 7:167–228 Walsh, AD (1953). The electronic orbitals, shapes, and spectra of polyatomic molecules. Part IV. Tetratomic hydride molecules, AH3. J Chem Soc 2296–2301 Weinberg S (1996) The quantum theory of fields, vol 2: modern applications. Cambridge University Press, Cambridge Wolfram S (2003) The mathematica book, 5th edn. Wolfram-Media, Champaign, Illinois Wu Q, Wu Y, Hao Y, Geng J, Charlton M, Chen S, Ren Y, Ji H, Li H, Boukhvalov DW, Piner RD, Bielawski CW, Ruoff RS (2013) Chem Commun 49:677–679 Yamaguchi Y, Osamura Y, Goddard JD, Schaeffer H III (1994) A new dimension to quantum chemistry: analytic derivative methods in ab-initio molecular electronic structure theory. Oxford University Press, Oxford