Benzene and its radical cation

We studied the $^1 A_{1g}$ ground state of the benzene molecule by using a very simple one particle basis set: for the AGP, a 2s1p DZ set centered on the carbon atoms and a 1s SZ on the hydrogen, instead for the 3-body Jastrow, a 1s1p DZ-GTO set centered only on the carbon sites. $C_6H_6$ is a peculiar molecule, since its highly symmetric ground state, which belongs to the $D_{6h}$ point group, is a resonance among different many-body states, each of them characterized by three double bonds between carbon atoms. This resonance is responsible for the stability of the structure and therefore for its aromatic properties. We started from a non resonating 2-body Jastrow wave function, which dimerizes the ring and breaks the full rotational symmetry, leading to the Kekulé configuration. As we expected, the inclusion of the resonance between the two possible Kekulé states lowers the VMC energy by more than 2 eV. The wave function is further improved by adding another type of resonance, that includes also the Dewar contributions connecting third nearest neighbor carbons.

Table 3.5: Bond lengths ($r$ ) for the two lowest $^2 B_{2g}$ and $^2 B_{3g}$ states of the benzene radical cation. The angles $\alpha$ are expressed in degrees, the lengths in $a_0$ . The carbon sites are numerated from 1 to 6.

	$^2 B_{2g}$	$^2 B_{3g}$	Computational method
	acute	obtuse
$r(C_1-C_2)$	2.616	2.694	B3LYP/cc-pVTZ (4)
	2.649	2.725	BLYP/6-31G* (3)
	2.659(1)	2.733(4)	SR-VMC 3.1
$r(C_2-C_3)$	2.735	2.579	B3LYP/cc-pVTZ (4)
	2.766	2.615	BLYP/6-31G* (3)
	2.764(2)	2.628(4)	SR-VMC 3.2
$\alpha(C_6 C_1 C_2)$	118.4	121.6	B3LYP/cc-pVTZ (4)
	118.5	121.5	BLYP/6-31G* (3)
	118.95(6)	121.29(17)	SR-VMC 3.1

As reported in Tab. 3.4, the gain with respect to the simplest Kekulé wave function amounts to 4.2 eV, but the main improvement arises from the further inclusion of the three-body Jastrow factor, which allows to recover the $89 \%$ of the total atomization energy at the VMC level. The main effect of the three body term is to keep the total charge around the carbon sites to approximately six electrons, thus penalizing the double occupation of the $p_z$ orbitals. The same important correlation ingredient is present in the well known Gutzwiller wave function already used for polyacetylene (71,72). Within this scheme we have systematically included in the 3-body Jastrow part the same type of terms present in the AGP one, namely both $g^{a,b}$ and $\lambda^{a,b}$ are non zero for the same pairs of atoms. As expected, the terms connecting next nearest neighbour carbon sites are much less important than the remaining ones because the VMC energy does not significantly improve (see the full resonating + 3-body wave function in Tab. 3.4). A more clear behaviour is found by carrying out DMC simulations: the interplay between the resonance among different structures and the Gutzwiller-like correlation refines more and more the nodal surface topology, thus lowering the DMC energy by significant amounts.

Figure 3.1: Electron density (atomic units) projected on the plane of $C_6H_6$ . The surface plot shows the difference between the resonating valence bond wave function, with the correct $A1g$ symmetry of the molecule, and a non-resonating one, which has the symmetry of the Hartree Fock wave function.

Therefore it is crucial to insert into the variational wave function all these ingredients in order to have an adequate description of the molecule. For instance, in Fig. 3.2 we report the density surface difference between the non-resonating 3-body Jastrow wave function, which breaks the $C_6$ rotational invariance, and the resonating Kekulé structure, which preserves the correct $A_{1g}$ symmetry: the change in the electronic structure is significant. The best result for the binding energy is obtained with the Kekulé Dewar resonating 3 body wave function, which recovers the $98,6\%$ of the total atomization energy with an absolute error of 0.84(8) eV. As Pauling (73) first pointed out, benzene is a genuine RVB system, indeed it is well described by the JAGP wave function. Moreover Pauling gave an estimate for the resonance energy of 1.605 eV from thermochemical experiments in qualitative agreement with our results. A final remark about the error in the total atomization energy: the latest frozen core CCSD(T) calculations (74,62) are able to reach a precision of 0.1 eV, but only after the complete basis set extrapolation and the inclusion of the core valence effects to go beyond the psudopotential approximation. Without the latter corrections, the error is quite large and even in the CCSD approach it is 0.65 eV (74). In our case, such an error arises from the fixed node approximation, whose nodal error is not compensated by the difference between the atomic and the molecular energies, as already noticed in the previous subsection.

Figure 3.2: Surface plot of the charge density projected onto the molecular plane. The difference between the non-resonating (indicated as HF) and resonating Kekulé 3-body Jastrow wave function densities is shown. Notice the corresponding change from a dimerized structure to a $C_6$ rotational invariant density profile.

The radical cation $C_6 H_6^+$ of the benzene molecule has been the subject of intense theoretical studies(4,3), aimed to focus on the Jahn-Teller distorted ground state structure. Indeed the ionized $^2 E_{1g}$ state, which is degenerate, breaks the symmetry and experiences a relaxation from the $D_{6h}$ point group to two different states, $^2 B_{2g}$ and $^2 B_{3g}$ , that belong to the lower $D_{2h}$ point group. In practice, the former is the elongated acute deformation of the benzene hexagon, the latter is its compressed obtuse distortion. We applied the SR structural optimization, starting from the $^2 E_{1g}$ state, and the minimization correctly yielded a deformation toward the acute structure for the $^2 B_{2g}$ state and the obtuse for the $^2 B_{3g}$ one; the first part of the evolution of the distances and the angles during those simulations is shown in Fig.3.3. After this equilibration, average over 200 further iterations yields bond distances and angles with the same accuracy as the all-electron BLYP/6-31G* calculations reported in Ref. (3) (see Tab.  3.5).

Figure 3.3: Plot of the convergence toward the equilibrium geometry for the $^2 B_{2g}$ acute and the $^2 B_{3g}$ obtuse benzene cation. Notice that both the simulations start form the ground state neutral benzene geometry and relax with a change both in the $C-C$ bond lengths and in the angles. The symbols are the same of Tab. 3.5.

As it appears from Tab. 3.6 not only the structure but also the DMC total energy is in perfect agreement with the BLYP/6-31G*, and much better than SVWN/6-31G* that does not contain semi empirical functionals, for which the comparison with our calculation is more appropriate, being fully ab-initio.

The difference of the VMC and DMC energies between the two distorted cations are the same within the error bars; indeed, the determination of which structure is the real cation ground state is a challenging problem, since the experimental results give a difference of only few meV in favor of the obtuse state and also the most refined quantum chemistry methods are not in agreement among themselves (3). A more affordable problem is the determination of the adiabatic ionization potential (AIP), calculated for the $^2 B_{3g}$ state, following the experimental hint. Recently, very precise CCSD(T) calculations have been performed in order to establish a benchmark theoretical study for the ionization threshold of benzene (4); the results are reported in Tab. 3.7. After the correction of the zero point energy due to the different structure of the cation with respect to the neutral molecule and taken from a B3LYP/cc-pVTZ calculation reported in Ref. (4), the agreement among our DMC result, the benchmark calculation and the experimental value is impressive. Notice that in this case there should be a perfect cancellation of nodal errors in order to obtain such an accurate value; however, we believe that it is not a fortuitous result, because in this case the underlying nodal structure does not change much by adding or removing a single electron.

Table: Total energies for the $^2 B_{2g}$ and $^2 B_{3g}$ states of the benzene radical cation after the geometry relaxation. A comparison with a BLYP/6-31G* and SVWN/6-31G* all-electron calculation (Ref. (3)) is reported.

	VMC	DMC	BLYP/6-31G*	SVWN/6-31G*
$^2 B_{2g}$	-231.4834(15)	-231.816(3)	-231.815495	-230.547931
$^2 B_{3g}$	-231.4826(14)	-231.812(3)	-231.815538	-230.547751

Therefore we expect that this property holds for all the affinity and ionization energy calculations with a particularly accurate variational wave function as the one we have proposed here. Nevertheless DMC is needed to reach the chemical accuracy, since the VMC result is slightly off from the experimental one just by few tenths of eV. The AIP and the geometry determination for the $C_6 H_6^+$ are encouraging to pursue this approach, with the aim to describe even much more interesting and challenging chemical systems.

Table: Adiabatic ionization potential of the benzene molecule; our estimate is done for the $^2 B_{3g}$ relaxed geometries of the benzene radical cation, with an inclusion of the zero point motion correction between the $^2 B_{3g}$ state and the $^1 A_{1g}$ neutral molecule ground state, calculated in Ref. (4) at the B3LYP/6-31G* level.

	VMC	DMC	CCSD(T)/cc-pV$\infty$ Z (4)	experiment (75)
AIP	8.86(6)	9.36(8)	9.29(4)
$\Delta ZPE_{ad}$	-0.074	-0.074	-0.074
best estimate	8.79(6)	9.29(8)	9.22(4)	9.2437(8)