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.

Charge 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 reactivity123. Their relevance further extends to biological polymers such as DNA and to light-harvesting processes using artificial photosynthesis34567. Localized electronic states and separation of charges are, furthermore, an important aspect of semiconductor devices and solar cells89101112. 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 error131415. 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 states181920. 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 charge-localized 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 observed21. DMP has previously served as a prototype for exploring electron lone pair interactions and CT between nitrogen atoms222324. 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 time-resolved experiment21. 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 energized25262728. 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 charge-localized and charge-delocalized cationic states.

The relative energy of a charge-localized state and a charge-delocalized 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 picoseconds21.

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, respectively21. 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)

Figure 1

Photoelectron spectra of DMP.

(a) The time-resolved spectrum of DMP with 231.5 nm pump photon. The colour bar represents the logarithmic intensity scale. (b) The 0–3 ps time-integrated spectrum of a. (c) The 20–200 ps time-integrated spectrum of a. (d) The 20–200 ps time-integrated spectra of DMP at five selected pump wavelengths. The relative populations of the charge-localized (3sL) and the charge-delocalized (3sD) states can be determined as a function of temperature from this data, thus providing an estimate of the relative energy of the two states, which turns out to be 0.33 eV. (e) A schematic cut of the potential energy surface for DMP+. The red and blue lines illustrate the vibrational states of DMP-L+ and DMP-D+, respectively.

Because the change of the Gibbs free energy (ΔG) scales linearly with the natural logarithm of the equilibrium constant (K), the enthalpy (ΔH) and entropy (ΔS) 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 Notes 4 and 5. Details of the fits are shown in Supplementary Notes 4 and 5, Supplementary Figs 3 and 4 and Supplementary Tables 1 and 2. The results are shown in 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 ΔH and ΔS 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.

Figure 2

Temperature dependence of the equilibrium constants.

Measured values of the equilibrium constant for the 3sL to 3sD states of the DMP molecule are shown as a function of reciprocal temperature estimated from the photon energy. The red line shows a linear best fit providing an estimate of the energy and entropy difference between the two states.

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 software31 at various levels of theory including HF, MP2 (Møller–Plesset perturbation theory, truncated at the second order), DFT with all the commonly used functionals (complete list available in Supplementary Note 2), and CCSD (coupled cluster method with single and double excitations). Unless specified otherwise, the aug-cc-pVDZ basis set was used in the Rydberg state calculations and the cc-pVTZ basis set in the cation calculations3233. The cation structures were also optimized using PZ-SIC with the GPAW software343536, 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 well37 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 CT394041, 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 method42 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.

Figure 3

Calculated minimum energy path between the localized and delocalized state of the DMP cation.

The energy of images, Eimg, in the nudged elastic band calculations is given with respect to the energy of the localized state, E1, as a function of the accumulated displacement of the atoms, ΔR. The red dots show results of a PZ-SIC calculation where a barrier of 0.2 eV separates the metastable, localized state from the delocalized state. The green dots show results of calculations using the M06-HF functional where MP2 optimized structures are used for the end points. In the M06-HF calculations, the energy barrier is not present and a structure optimization starting from the localized state converges on the delocalized state. Similar results were obtained for all other commonly used DFT functionals.

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 single-point 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.

Table 1

RE of the DMP-L+ and DMP-D+ states obtained using various computational methods.

Method	RE (eV)
HF	−0.53
MP2	0.81
DFT
B3LYP	—a
M06	— a
M06-2X	— a
M06-HF	— a
PBE0	— a
BHandHLYP	0.19
PZ-SIC	0.34
CCSDb	0.23
MP2_CCSD(T)-SP	0.39
CCSD_CCSD(T)-SP	0.38
Experiment	0.33 (0.04)

CCSD, coupled cluster method with single and double excitations; DFT, density functional theory; DMP, N,N'-dimethylpiperazine; HF, Hartree–Fock; PZ-SIC, Perdew and Zunger self-interaction correction; RE, relative energy.

Single-point energy calculations were carried out with the CCSD(T) method, in one case using a structure obtained with MP2 and in the other case using a structure obtained with CCSD.

Zero point energy correction has not been applied but an estimate based on ground vibrational states would reduce the calculated relative energy by 0.07 eV, see Supplementary Note 3. As the molecule is at high temperature, the full, ground vibrational state zero point correction is an overestimate.

aNo value shown because the DMP-L+ state was not stable at this level of theory.

bThe CCSD optimizations were carried out with the aug-cc-pVDZ basis set because calculations with the cc-pVTZ were found to be too demanding for our computational resources.

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 electron-binding 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 clusters18272843. The total energy of the Rydberg excited states using PZ-SIC was obtained using the delta self-consistent field method44 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.

Table 2

Calculated Rydberg binding energy (in eV) of 3sD and 3sL states for MP2 optimized DMP-L+ and DMP-D+ structures using the EOM-CCSD and PZ-SIC methods.

Method	3sD	3sL
PZ-SIC	2.71	2.87
EOM-CCSD	2.65	2.72
Experiment	2.70 (0.03)	2.81 (0.04)

EOM-CCSD, equation of motion coupled cluster method with single and double excitations; PZ-SIC, Perdew and Zunger self-interaction correction.

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-interaction22 in the charge-delocalized ion.

Figure 4

The Rydberg orbitals and the associated spin densities.

(a,b) Calculated 3sL and 3sD Rydberg orbitals, respectively, rendered at 0.001 Å−3/2 isovalues. (c,d) Calculated spin density of the DMP-L+ and DMP-D+ ion, respectively, at isovalue of 0.2 electron per Å−3.

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+ ion2324 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 N3, 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 self-consistent 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 N3. These results are expected to guide future improvements to energy functionals describing electronic systems.

Additional information

How to cite this article: Cheng, X. et al. Charge localization in a diamine cation provides a test of energy functionals and self-interaction correction. Nat. Commun. 7:11013 doi: 10.1038/ncomms11013 (2016).