Charge localization in a diamine cation provides a test of energy functionals and self-interaction correction

Density functional theory (DFT) is widely applied in calculations of molecules and materials. Yet, it suffers from a well-known over-emphasis on charge delocalization arising from self-interaction error that destabilizes localized states. Here, using the symmetric diamine N,N′-dimethylpiperazine as a model, we have experimentally determined the relative energy of a state with positive charge localized on one of the two nitrogen atoms, and a state with positive charge delocalized over both nitrogen atoms. The charge-localized state was found to be 0.33 (0.04) eV higher in energy than the charge-delocalized state. This provides an important test of theoretical approaches to electronic structure calculations. Calculations with all DFT functionals commonly used today, including hybrid functionals with exact exchange, fail to predict a stable charge-localized state. However, the application of an explicit self-interaction correction to a semi-local functional identifies both states and gives relative energy in excellent agreement with both experiment and CCSD(T) calculations.

C harge transfer (CT) between two or more charge centres and charge delocalization across separate parts of a molecular system are of great importance in many areas of chemistry as they often define molecular structure and reactivity [1][2][3] . Their relevance further extends to biological polymers such as DNA and to light-harvesting processes using artificial photosynthesis [3][4][5][6][7] . Localized electronic states and separation of charges are, furthermore, an important aspect of semiconductor devices and solar cells [8][9][10][11][12] . Yet, the calculation of charge localization and charge delocalization is challenging, especially for computational methods that are applicable to extended systems.
The commonly used Kohn-Sham density functional theory (DFT) functionals are known to suffer from bias towards delocalized states because of self-interaction error [13][14][15] . Hartree-Fock (HF), on the other hand, is overly biased towards localized states. Hybrid functionals that mix the two approaches are widely used in all branches of chemistry and are increasingly used in condensed matter calculations, but the balance between localized and delocalized states depends strongly on the proportions used in the mixing. This mixing ratio is sometimes treated as an empirical, system-dependent parameter. A more fundamental, parameter-free approach is, however, needed to accurately predict the delicate balance between localized and delocalized electronic states, especially in extended systems. Over 30 years ago, an explicit self-interaction correction approach was proposed by Perdew and Zunger (PZ-SIC) 16 . Although it was applied early on in studies of the electronic structure of molecules (see, for example, ref. 17), it has not become a commonly used approach. A variational, self-consistent implementation of this PZ-SIC using complex valued orbitals has recently been applied successfully in studies of molecules and solids, including Rydberg excited states [18][19][20] . A discussion of the self-interaction error and our implementation of PZ-SIC can be found in Supplementary Note 1. We show here that although all the commonly used DFT functionals, including hybrid functionals, fail to produce the localized charge state of the system studied, calculations using PZ-SIC give results in excellent agreement with the experimentally determined relative energy of a localized and delocalized electronic state.
To achieve this comparison, experimental benchmarks for chargelocalized and charge-delocalized states are required. This is challenging because in most cases electronic systems exhibit just one particular charge distribution, making comparison between states with different charge distributions difficult. In a recent study, however, we have found that in N,N'-dimethylpiperazine (DMP), both a charge-localized state and a charge-delocalized state can be observed 21 . DMP has previously served as a prototype for exploring electron lone pair interactions and CT between nitrogen atoms [22][23][24] . In its cationic state, the positive charge can be localized on one of the nitrogen atoms (DMP-L þ ), or delocalized over the two equivalent nitrogen atoms (DMP-D þ ). The charge-localized ion DMP-L þ , however, has not been observed until the recent ultrafast timeresolved experiment 21 . It is possible to identify the different charge states because they have distinct spectral and temporal signatures. Upon optical excitation, the charge-localized state is initially generated. A subsequent CT then leads to the charge-delocalized state. This discovery of two states with distinct charge distributions has now laid the foundation for an experimental approach for measuring their relative energy. The approach is based on photoionization from Rydberg states, whose binding energies have been found to be remarkably dependent on the nuclear arrangements and the charge distribution of the molecular ion core. Further, the Rydberg electron-binding energy is independent of a molecule's internal energy, making it ideally suited for the exploration of high-energy charge states that are populated when the molecules are highly energized [25][26][27][28] . Because the Rydberg electron has a small effect on the bonding and molecular structure, the configuration of the molecular ion cores of Rydberg states closely resemble those of the cationic states. The binding energies of the Rydberg states, therefore, yield information about the chargelocalized and charge-delocalized cationic states.
The relative energy of a charge-localized state and a chargedelocalized state of a given molecule has not been determined previously, as far as we know, even though it is a critically important parameter for calibrating theoretical approaches. In the current study, the relative energy was determined using a newly devised experimental approach that takes advantage of an equilibrium established on a picosecond time scale upon excitation to the Rydberg states. Calculations using conventional DFT, self-interaction-corrected DFT and ab initio methods were tested against the experimental measurements. We aim to evaluate the capability of each method to properly describe the charge localization and delocalization in the ground electronic state of the molecular cation.

Results
Experimental determination of the relative energy. To measure the energy of the cationic state, we first experimentally determined the relative energy of 3sD and 3sL Rydberg states with the charge-delocalized ion core DMP-D þ and the charge-localized ion core DMP-L þ , respectively, by measuring the equilibrium composition of the Rydberg states as a function of the excitation energy. As shown in Supplementary Fig. 1, because the binding energy, that is, the energy difference between the cationic state and the Rydberg state, is measured in the experiment, the relative energy of the cationic states is known once the relative energy of the 3s Rydberg states is known. DMP was excited from its ground state to the 3p or the 3s Rydberg state using pump photons with wavelengths tuned in the range of 193.0-240.8 nm. The probe photon monitored the time-dependent dynamics by ionizing the Rydberg-excited molecules. The kinetic energy of the ejected photoelectrons was measured to determine the binding energy of the Rydberg states. The 3p state, which was reached with the shorter wavelengths, decayed by internal conversion into 3s within several hundred femtoseconds. Details of the experimental setup and parameters are given in Supplementary Note 2. As we have shown previously, upon optical excitation to the localized charge state, an equilibrium between the localized and delocalized states is reached after several picoseconds 21 .
The complete set of the time-resolved photoelectron spectra at several pump wavelengths are given in Supplementary Fig. 2. Because the energy is conserved in internal conversion, pump photons of different wavelengths deposit different amounts of energy into the vibrational manifolds. For pump wavelengths between 193.0 and 240.8 nm, the molecule has effective vibrational temperature between 565 and 980 K after relaxation, assuming the energy is distributed across all vibrational modes. The detailed calculation of the effective temperature is discussed in the Supplementary Note 6.
Two 3s peaks (Fig. 1), located at 2.70 (0.03) eV and at 2.81 (0.04) eV, are assigned to 3sD and 3sL, respectively 21 . The 3sL dominates at short delay times (Fig. 1a,b), whereas 3sD dominates at longer delay times (Fig. 1a,c,d), indicating an early population in 3sL and a lower energy for 3sD. As the temperature decreases from 980 to 565 K with the pump wavelength increasing from 193.0 nm to 240.8 nm, the intensity of 3sL decreases almost to the baseline in the spectra with the molecules at equilibrium (Fig. 1d) Because the change of the Gibbs free energy (DG) scales linearly with the natural logarithm of the equilibrium constant (K), the enthalpy (DH) and entropy (DS) change in the transition can be determined from the logarithm of the equilibrium constant as a function of inverse temperature: where R is the gas constant and T is the temperature. The spectra at each point in time were fitted using two Lorentzians with variable peak centres to derive the equilibrium constants in the 3s states, as described in Supplementary Fig. 2. The logarithm of the equilibrium constant is indeed found to depend approximately linearly on the estimated reciprocal temperature. A fit using equation (1) gives À 20.9 (3.7) kJ mol À 1 or À 0.22 (0.04) eV and À 17.7 (4.6) J K À 1 mol À 1 for DH and DS of the transition from 3sL to 3sD. Because the binding energy difference between 3sL and 3sD is 0.11 (0.01) eV, the energy of the DMP-L þ ion is determined to be 0.33 (0.04) eV higher than that of the DMP-D þ ion. A schematic cut through the energy surface, deduced from these measurements, is shown in Fig. 1e.
Test of theoretical methods. The ability of various theoretical approaches to describe the charge localized and delocalized states can now be assessed by comparison with these experimental results. First, calculations were carried out to determine the optimal molecular geometries of the two states. The DMP-L þ and DMP-D þ ion structures were optimized with the Gaussian 09 (refs 29,30) and NWChem software 31    cation calculations 32,33 . The cation structures were also optimized using PZ-SIC with the GPAW software [34][35][36] , where a real space grid over a cubic simulation cell of 20 Å side length and 0.13 Å mesh was used. The PZ-SIC was applied to the PBE semi-local functional. Although in the neutral molecule each nitrogen atom assumes a pyramidal structure with the methyl group in equatorial position, the nitrogen becomes pseudo-axial and pseudo-planar in the cation. The Cartesian coordinates of selected optimized structures are listed in Supplementary  Tables 4 and 5. DFT calculations with any one of the available hybrid functionals implemented in the Gaussian 09 software failed to provide a stable localized state, except for the BHandHLYP functional that has previously been found to describe CT interactions well 37 while giving generally poor results for other molecular properties such as total energy (and therefore not commonly used) 38 . When DFT calculations were started from MP2 or HF-optimized structures for the localized state, the minimization of the energy resulted in a conversion to the charge delocalized state. Most surprisingly, the M06-HF functional, which contains 100% HF exchange and is therefore widely deemed to be particularly appropriate for CT [39][40][41] , also fails to localize the charge in this case. However, the PZ-SIC calculation gives a stable localized state with similar structure as that obtained from MP2 and CCSD. The bond lengths are typically 0.02-0.04 Å shorter than those obtained from MP2 and CCSD. The minimum energy path between the two states calculated using the nudged elastic band method 42 and the PZ-SIC as well as the M06-HF functional is shown in Fig. 3. An energy barrier of 0.2 eV for the transition from the localized to the delocalized state is obtained with PZ-SIC, whereas no barrier is obtained in the M06-HF calculations. Table 1 lists the calculated energy difference between the optimized DMP-L þ and DMP-D þ structures obtained using HF, MP2, CCSD, DFT with selected functionals and PZ-SIC. CCSD(T) calculations were carried out to obtain the singlepoint energy for each one of these structures. Although calculations using HF, MP2 and CCSD produce both the DMP-L þ and the DMP-D þ states, the relative energy of these states is poorly estimated, especially by HF, which gives lower energy for the localized state than the delocalized state. Single-point CCSD(T) calculations using the MP2 or CCSD geometries give a relative energy in good agreement with the experimental measurements. The PZ-SIC calculation also yields a relative energy that is close to the experimental results, see Table 1 and  Supplementary Table 3. Supplementary Note 7 also reports a satisfactory agreement of the experimental and computed entropy differences.
The close agreement between the relative energy obtained from the high-level CCSD(T) calculations and the experimental results confirms the validity of the interpretation of the experimental data. To further cement the correspondence of experimental and computational results, we have calculated the Rydberg electronbinding energy with the optimized DMP-L þ and DMP-D þ structures using the equation of motion CCSD. For comparison, the binding energy was also calculated using PZ-SIC, which has previously been shown to give good estimates of Rydberg binding energy of both molecules and molecular clusters 18,27,28,43 . The total energy of the Rydberg excited states using PZ-SIC was obtained using the delta self-consistent field method 44 and the binding energy obtained by subtracting the total energy of the excited state from that of the ion. As listed in Table 2, the calculated binding energy is in good agreement with the experimentally measured values for both the DMP-L þ and DMP-D þ structures, supporting the assignment of the observed spectroscopic features.
The Rydberg orbitals and the associated spin densities are shown in Fig. 4. The 3sL Rydberg orbital (Fig. 4a) anchors on the planar nitrogen, whereas the 3sD Rydberg orbital (Fig. 4b) is centred symmetrically between the two nitrogen atoms. Both orbitals are extended and comprise the whole molecule, as expected. The spin densities (shown in Fig. 3c,d), which were generated by subtracting the spin-down density from the spin-up density, illustrate the charge distributions of the localized and  delocalized cations. The charge is localized on the planar nitrogen in the DMP-L þ (Fig. 3c) but delocalized between the two nitrogen atoms in DMP-D þ (Fig. 3d). Intriguingly, there is also net spin density between the two intermediate C-atoms, indicating a through-bond-interaction 22 in the chargedelocalized ion.

Discussion
The present study advances experimentally the state-of-the-art in exploring charge localization and delocalization in systems with multiple charge centres. The wavelength-dependent Rydberg electron-binding energy spectroscopy enables us to experimentally determine the energy and entropy change of the CT reaction. It is the experimental tool of choice to probe the CT process in the presence of large vibrational energy, as is needed to observe the higher-energy state before CT. Indeed, previous time-resolved spectroscopic studies on DMP radical cations prepared by a photo-induced electron transfer technique only observed the DMP-D þ ion 23,24 possibly because of the low temperature in the system. Together with the binding energy information obtained from the spectra, the present experiment creates an important benchmark to evaluate various theoretical approaches.
The MP2 and CCSD methods are found to give good estimates of the molecular structure as the single-point CCSD(T) calculation gives good agreement with the experimental binding energies of the 3sL and 3sD states as well as the relative energy of the two cation states. The computational effort of these methods, however, scales unfavourably with system size and they are thus limited to small systems. The computational effort of DFT calculations increases slower with size as N 3 , where N is the number of electrons, and is the only viable approach for many problems involving large molecules and condensed phase systems. Conventional DFT functionals are, however, found here not to predict the metastable localized state of DMP. The explicit inclusion of the self-interaction correction, as proposed by Perdew and Zunger, implemented in a variational and selfconsistent way with complex-valued orbitals, can remedy these shortcomings of the DFT approach. The PZ-SIC calculations give similar results for the binding energies of the 3sL and 3sD states as the equation of motion CCSD calculations and for the relative energies of the DMP-L þ and DMP-D þ states as the CCSD(T) calculations. The computational effort in this approach is larger than conventional DFT but still scales with size as N 3 . These results are expected to guide future improvements to energy functionals describing electronic systems.