Diagnosis of two evaluation paths to density-based descriptors of molecular electronic transitions

TL;DR

This paper evaluates the reliability of two computational methods for analyzing electronic transitions in molecules, comparing their accuracy and convergence across various molecules, basis sets, and functionals.

Contribution

It introduces a diagnostic framework to assess the accuracy of population analysis relative to numerical integration for density-based descriptors.

Findings

Population analysis can deviate from numerical integration results.

Eight diagnostic tests effectively identify when population analysis is unreliable.

The reliability varies with basis sets and exchange-correlation functionals.

Abstract

In this paper we discuss the reliability of two computational methods (numerical integration on Cartesian grids, and population analysis) used for evaluating scalar quantities related to the nature of electronic transitions. These descriptors are integrals of charge density functions built from the detachment and attachment density matrices projected in the Euclidean space using a finite basis of orbitals. While the numerical integration on Cartesian grids is easily considered to be converged for medium-sized density grids, the population analysis approximation to the numerical integration values is diagnosed using eight diagnostic tests performed on fifty-nine molecules with a combination of fifteen Gaussian basis sets and six exchange-correlation functionals.

Figures

Tables

Table 1. Table 1: Summary of the eight detachment/attachment PA diagnosis strategies. Molecules for d 1 subscript 𝑑 1 d_{1} and d 2 subscript 𝑑 2 d_{2} are taken from Ref. etienne_toward_2014 ; molecules for d 3 subscript 𝑑 3 d_{3} to d 5 subscript 𝑑 5 d_{5} are taken from Ref. peach_excitation_2008 ; and PP is taken from Ref. pastore_unveiling_2017 . The sketches of the molecules are reported in figure S1 of the Supplementary Information.

Diagnosis	Molecule	System size	CT nature	BS	xc-functional
$d_1$	nIII, nV-X, nVI-c	$\surd$	$\surd$	$\surd$	$\times$
$d_2$	nV-X	$\times$	$\surd$	$\surd$	$\times$
$d_3$	N2, CO, HCl, H2CO	$\times$	$\times$	$\times$	$\surd$
$d_4$	acenes	$\surd$	$\times$	$\surd$	$\surd$
$d_5$		$\times$	$\times$	$\surd$	$\times$
$d_6$	PP	$\times$	$\surd$	$\times$	$\times$
$d_7$		$\times$	$\times$	$\surd$	$\times$
$d_8$		$\times$	($\surd$)	$\times$	$\surd$

Table 2. Table 2: Descriptors values and their PA-approximated deviations, Δ Δ \Delta ( q 𝑞 q ) = q absent 𝑞 =q (NI ) − q )-q (PA), obtained for the three first excited states of the oligo-acenes series with the B3LYP/cc-pVTZ level of theory (part of d 4 subscript 𝑑 4 d_{4} ).

$n$	E-S	${\varphi }_S$	$\Delta$(${\varphi }_S$)	$\phi$	$\Delta$($\phi$)	$\psi$	$\Delta$($\psi$)
1	1	0.92	-0.06	0.26	0.13	0.82	-0.09
	2	0.98	0.00	0.15	0.00	0.90	0.00
	3	0.88	-0.07	0.32	0.11	0.78	-0.08
2	1	0.90	0.08	0.29	0.17	0.80	-0.12
	2	0.98	-0.01	0.15	0.02	0.90	-0.01
	3	0.91	-0.06	0.27	0.09	0.81	-0.07
3	1	0.90	-0.09	0.27	0.10	0.81	-0.12
	2	0.97	-0.01	0.15	0.03	0.90	-0.02
	3	0.90	-0.06	0.30	0.07	0.79	-0.05
4	1	0.90	0.10	0.27	0.18	0.81	-0.13
	2	0.89	-0.05	0.30	0.05	0.79	-0.04
	3	0.97	-0.01	0.15	0.03	0.90	-0.02
5	1	0.90	-0.09	0.26	0.17	0.82	-0.12
	2	0.90	-0.06	0.33	0.06	0.78	-0.05
	3	0.89	-0.05	0.35	0.09	0.76	-0.07

Table 3. Table 3: Detachment/attachment spatial overlap integral values and their PA-approximated deviations, Δ Δ \Delta (X) = ϕ S absent subscript italic-ϕ 𝑆 =\phi_{S} (X) − ϕ S ℓ superscript subscript italic-ϕ 𝑆 ℓ -\phi_{S}^{\ell} (X), calculated for the different molecules from the nV-X molecular test-set (part of d 1 subscript 𝑑 1 d_{1} ). Those values were obtained using PBE0/6-311++G(2d,p) excited-state calculations, based on PBE0/6-311G(d,p) geometries.

$n$	E-S	${\varphi }_S$(NH)	$\Delta$(NH)	${\varphi }_S$(O)	$\Delta$(O)	${\varphi }_S$(S)	$\Delta$(S)	${\varphi }_S$(Se)	$\Delta$(Se)
1	1	0.69	0.02	0.68	0.03	0.70	0.03	0.72	0.04
	2	0.48	0.36	0.48	0.36	0.48	0.39	0.48	0.38
	3	0.34	0.14	0.52	0.40	0.51	0.41	0.54	0.03
2	1	0.57	0.01	0.60	0.03	0.58	0.02	0.68	0.02
	2	0.49	0.42	0.48	0.35	0.76	0.25	0.52	0.25
	3	0.40	0.03	0.71	0.04	0.72	0.00	0.60	0.12
3	1	0.48	0.01	0.53	0.02	0.51	0.02	0.62	0.02
	2	0.73	-0.01	0.78	-0.02	0.75	-0.01	0.79	-0.01
	3	0.63	0.01	0.69	0.00	0.70	-0.01	0.74	-0.01
4	1	0.39	0.00	0.47	0.02	0.52	0.01	0.57	0.01
	2	0.60	0.00	0.77	-0.02	0.75	-0.02	0.78	-0.01
	3	0.71	-0.03	0.66	0.00	0.70	-0.01	0.74	-0.01
5	1	0.34	0.01	0.42	0.01	0.47	0.00	0.53	0.00
	2	0.52	0.00	0.75	-0.02	0.73	-0.02	0.76	-0.02
	3	0.71	-0.04	0.65	-0.01	0.68	-0.02	0.72	-0.02

Table 4. Table 4: Descriptor values for the first three excited states of PP obtained with B3LYP and CAM-B3LYP xc-functionnals and six BSs out of the fifteen investigated in this contribution. Two dihedral angles have been choosen (0 ∘ and 80 ∘ ). K 𝐾 K is the number of basis functions.

angle 0∘
BS	K	B3LYP		CAM-B3LYP		B3LYP		CAM-B3LYP		B3LYP		CAM-B3LYP
		${\varphi }_S$	${\varphi }_S^{\ell }$	${\varphi }_S$	${\varphi }_S^{\ell }$	$\phi$	${\phi }^{\ell }$	$\phi$	${\phi }^{\ell }$	$\psi$	${\psi }^{\ell }$	$\psi$	${\psi }^{\ell }$
3-21G	135	0.33	0.43	0.74	0.71	0.86	0.82	0.60	0.63	0.23	0.31	0.57	0.54
6-31G	135	0.32	0.42	0.73	0.69	0.86	0.83	0.62	0.65	0.23	0.30	0.55	0.52
6-31+G	187	0.31	0.41	0.70	0.66	0.87	0.83	0.63	0.66	0.22	0.30	0.53	0.50
6-31G(d)	213	0.31	0.42	0.73	0.69	0.86	0.83	0.61	0.65	0.22	0.30	0.56	0.52
6-311+G(d)	313	0.31	0.40	0.71	0.66	0.87	0.84	0.63	0.66	0.22	0.29	0.54	0.50
6-311++G(2d,p)	414	0.31	0.41	0.70	0.66	0.87	0.83	0.63	0.66	0.22	0.29	0.53	0.50
angle 80∘
BS	K	B3LYP		CAM-B3LYP		B3LYP		CAM-B3LYP		B3LYP		CAM-B3LYP
		${\varphi }_S$	${\varphi }_S^{\ell }$	${\varphi }_S$	${\varphi }_S^{\ell }$	$\phi$	${\phi }^{\ell }$	$\phi$	${\phi }^{\ell }$	$\psi$	${\psi }^{\ell }$	$\psi$	${\psi }^{\ell }$
3-21G	135	0.53	0.54	0.37	0.34	0.81	0.82	0.90	0.93	0.37	0.37	0.25	0.22
6-31G	135	0.55	0.55	0.33	0.30	0.81	0.82	0.92	0.94	0.38	0.38	0.22	0.20
6-31+G	187	0.50	0.51	0.30	0.26	0.84	0.85	0.93	0.95	0.35	0.34	0.20	0.17
6-31G(d)	213	0.54	0.54	0.33	0.30	0.81	0.83	0.92	0.94	0.37	0.37	0.22	0.20
6-311+G(d)	313	0.49	0.50	0.30	0.26	0.84	0.86	0.93	0.95	0.34	0.33	0.20	0.17
6-311++G(2d,p)	414	0.49	0.50	0.30	0.27	0.84	0.85	0.93	0.94	0.34	0.33	0.20	0.18

Equations

$$\text{r}⟼n(\text{r})=\sum _{r=1}^L\sum _{s=1}^L(\mathbf{\gamma }{)}_{rs}{\chi }_r(\text{r}){\chi }_s^{\ast }(\text{r}).$$

$${\chi }_r(\text{r})=\sum _{\mu =1}^K(\text{C}{)}_{\mu r}{\varphi }_{\mu }(\text{r})(1\le r\le L)$$

$$n(\text{r})=\sum _{\mu =1}^K\sum _{\nu =1}^K(\text{P}{)}_{\mu \nu }{\varphi }_{\mu }(\text{r}){\varphi }_{\nu }^{\ast }(\text{r})$$

$$\text{P}=\text{C}\mathbf{\gamma }{\text{C}}^{†}.$$

$$\forall (\mu ,\nu )\in [[1,K{]]}^2,(\text{S}{)}_{\mu \nu }:={\int }_{{\mathbb{R}}^3}d\text{r}{\varphi }_{\mu }^{\ast }(\text{r}){\varphi }_{\nu }^{\ast }(\text{r})$$

$${\int }_{{\mathbb{R}}^3}d\text{r}n(\text{r})=\mathrm{tr}\left(\text{PS}\right).$$

$${\int }_{{\mathbb{R}}^3}d\text{r}n(\text{r})=\mathrm{tr}\left({\text{S}}^x\text{P}{\text{S}}^{1-x}\right).$$

$$\text{r}⟼f_{\pm }(\text{r})=\left|f(\text{r})\right|\pm f(\text{r}).$$

$$f_{+}(\text{r})-f_{-}(\text{r})=f(\text{r});f_{+}(\text{r})+f_{-}(\text{r})=\left|f(\text{r})\right|.$$

$$G[f]$$

$$:=\frac{1}{4}{\int }_{{\mathbb{R}}^3}d\text{r}\sum _{\tau \in \left\{-1,1\right\}}\underset{{2f}_{\tau }(\text{r})}{\underbrace{\left|f(\text{r})\right|+\tau f(\text{r})}}.$$

$$\sum _{\tau \in \left\{-1,1\right\}}\tau {\int }_{{\mathbb{R}}^3}d\text{r}f(\text{r})=0⟺G[f]=\frac{1}{2}{\int }_{{\mathbb{R}}^3}d\text{r}\left|f(\text{r})\right|.$$

$${\int }_{{\mathbb{R}}^3}d\text{r}f(\text{r})=0⟺{\int }_{{\mathbb{R}}^3}d\text{r}{f}_{+}(\text{r})={\int }_{{\mathbb{R}}^3}d\text{r}{f}_{-}(\text{r}),$$

$$H[\text{Q}]:=\frac{1}{4}\sum _{\tau \in \left\{-1,1\right\}}\mathrm{tr}\left[\sqrt{{\left(\text{Q}\circ I_m\right)}^2}+\tau \left(\text{Q}\circ I_m\right)\right]$$

$${\left(\text{Q}\circ \text{E}\right)}_{ij}=(\text{Q}{)}_{ij}(\text{E}{)}_{ij},$$

$${\left(\text{Q}\circ I_m\right)}_{ij}={\left(\text{Q}\right)}_{ij}{\delta }_{ij},{\left(\sqrt{{\left(\text{Q}\circ I_m\right)}^2}\right)}_{ij}=\left|{\left(\text{Q}\right)}_{ij}\right|{\delta }_{ij},$$

$$\sum _{\tau \in \left\{-1,1\right\}}\tau \mathrm{tr}\left(\text{Q}\circ I_m\right)=0⟺2H[\text{Q}]=\mathrm{tr}\sqrt{{\left(\text{Q}\circ I_m\right)}^2}.$$

$${\mathbf{\gamma }}^{\Delta }={\mathbf{\gamma }}^q-{\mathbf{\gamma }}^0\in {\mathbb{R}}^{L\times L},$$

$$\exists (\text{M},\text{K})\in {\left({\mathbb{R}}^{L\times L}\right)}^2,{\text{M}}^{†}{\mathbf{\gamma }}^{\Delta }\text{M}=\mathrm{diag}(K_i{)}_{i\in [[1,L]]}.$$

$${\text{K}}_{\pm }=\frac{1}{2}\left(\sqrt{{\text{K}}^2}\pm \text{K}\right).$$

$$\text{M}{\text{K}}_{-}{\text{M}}^{†}={\mathbf{\gamma }}^d;\text{M}{\text{K}}_{+}{\text{M}}^{†}={\mathbf{\gamma }}^a.$$

$$n_d(\text{r})$$

$$=\sum _{\mu =1}^K\sum _{\nu =1}^K(\underset{\text{D}}{\underbrace{\text{C}{\mathbf{\gamma }}^d{\text{C}}^{†}}}{)}_{\mu \nu }{\varphi }_{\mu }(\text{r}){\varphi }_{\nu }^{\ast }(\text{r})$$

$$n_a(\text{r})$$

$$=\sum _{\mu =1}^K\sum _{\nu =1}^K(\underset{\text{A}}{\underbrace{\text{C}{\mathbf{\gamma }}^a{\text{C}}^{†}}}{)}_{\mu \nu }{\varphi }_{\mu }(\text{r}){\varphi }_{\nu }^{\ast }(\text{r}).$$

$$\vartheta$$

$$={\int }_{{\mathbb{R}}^3}d\text{r}{n}_a(\text{r})=\mathrm{tr}\left({\text{S}}^x\text{A}{\text{S}}^{1-x}\right)$$

$$n_{\Delta }(\text{r})=n_a(\text{r})-n_d(\text{r})\in \mathbb{R};\eta (\text{r})=\sqrt{n_d(\text{r})n_a(\text{r})}\in {\mathbb{R}}^{+},$$

$$n_{\Delta }(\text{r})=n_n(\text{r})-n_0(\text{r}),$$

$${\varphi }_S:={\vartheta }^{-1}{\int }_{{\mathbb{R}}^3}d\text{r}\sqrt{n_d(\text{r})n_a(\text{r})}=2{\vartheta }^{-1}G[\eta ]\in [0,1]$$

$$\chi :=G[n_{\Delta }];\vartheta \ge \chi ⟹\exists !\phi \in [0,1],\phi ={\vartheta }^{-1}\chi .$$

$$\psi :={\left(\pi /2\right)}^{-1}\arctan \left(\frac{{\varphi }_S}{\phi }\right)\in [0,1[$$

Full text

 

Keywords: Excited states ; electronic-structure reorganization ; one-body reduced density matrices ; orbital population analysis.

I Introduction

The possibility of analyzing the charge transfer character of electronic transitions of complex molecular systems has been of interest for decades, and still attracts considerable attention in the excited-state electronic structure community luzanov_interpretation_1980 ; head-gordon_analysis_1995 ; furche_density_2001 ; tretiak_density_2002 ; martin_natural_2003 ; batista_natural_2004 ; tretiak_exciton_2005 ; dreuw_single-reference_2005 ; mayer_using_2007 ; surjan_natural_2007 ; wu_exciton_2008 ; luzanov_electron_2009 ; li_time-dependent_2011 ; plasser_analysis_2012 ; plasser_new_2014 ; plasser_new_2014-1 ; bappler_exciton_2014 ; ronca_charge-displacement_2014 ; ronca_density_2014 ; etienne_toward_2014 ; etienne_new_2014 ; etienne_probing_2015 ; pastore_unveiling_2017 ; pluhar_visualizing_2018 ; mewes_communication_2015 ; li_particlehole_2015 ; etienne_transition_2015 ; plasser_statistical_2015 ; poidevin_truncated_2016 ; wenzel_physical_2016 ; plasser_entanglement_2016 ; savarese_metrics_2017 ; plasser_detailed_2017 ; etienne_theoretical_2017 ; mai_quantitative_2018 ; mewes_benchmarking_2018 ; skomorowski_real_2018 ; park_low_2018 ; barca_excitation_2018 ; etienne_comprehensive_2018 . In this context, a variety of quantities have been designed in order to unveil the nature of electronic transitions with scalars that provide – through a simple number – particular information related to the light-induced electronic structure reorganization luzanov_interpretation_1980 ; peach_excitation_2008 ; luzanov_electron_2009 ; le_bahers_qualitative_2011 ; garcia_evaluating_2013 ; guido_metric_2013 ; etienne_toward_2014 ; etienne_new_2014 ; guido_effective_2014 ; etienne_probing_2015 ; plasser_statistical_2015 ; wenzel_physical_2016 ; poidevin_truncated_2016 ; plasser_entanglement_2016 ; savarese_metrics_2017 ; pastore_unveiling_2017 ; etienne_theoretical_2017 ; plasser_detailed_2017 ; mai_quantitative_2018 ; barca_excitation_2018 ; mewes_benchmarking_2018 ; etienne_charge_2018 ; campetella_quantifying_2019 . An example of information that one might wish to seek when considering a molecular system undergoing a photoinduced intramolecular charge transfer is the locality of this charge transfer, i.e., the extent to which the electronic cloud has been perturbed and polarized etienne_toward_2014 ; etienne_new_2014 ; etienne_probing_2015 ; etienne_theoretical_2017 ; etienne_charge_2018 . To such probing of the charge transfer locality one can for example join the extent to which the detached and attached charges contribute to the net charge displacement occurring during an electronic transition etienne_probing_2015 ; etienne_theoretical_2017 ; etienne_charge_2018 . These two pieces of information are contained in two separate quantities, named $\phi_{S}$ and $\varphi$, respectively.

Whilst the possibility of using computational resources for reaching such theoretical insights has already been used in theoretical and applied research, three different implementations were mentioned for evaluating these two density-based descriptors etienne_probing_2015 . Two of them were based on the population analysis mulliken_electronic_1955 ; mulliken_electronic_1955-1 ; mulliken_criteria_1962 ; davidson_electronic_1967 ; reed_natural_1985 ; huzinaga_extended_1990 ; parrondo_natural_1994 ; clark_density_2004 ; mayer_lowdin_2004 ; bruhn_lowdin_2006 ; bachrach_population_2007 ; saha_are_2009 ; jacquemin_what_2012 of the so-called detachment and attachment one-body reduced density matrices head-gordon_analysis_1995 , and the third one was the direct numerical integration of the related one-body charge density functions using Cartesian grids. When first published etienne_probing_2015 , the population analysis approximations to our density-based descriptors were compared to the numerical integration method as a reference solely for organic push-pull molecules using a single level of theory, for the sake of demonstrating its hypothetical usefulness, given that its implementation allows an extremely fast evaluation of the $\phi_{S}$ and $\varphi$ descriptors.

In this contribution we wish to extend this diagnosis to a much wider set made of various types of molecules, and bring a rigorous and extensive diagnosis of the population analysis as an evaluation tool for our density-based descriptors. For this sake, we will use and combine in this paper multiple diagnosis criteria, and provide a precise motivation for the choice of the population analysis to perform, as well as a discussion based on the understanding of the possible preliminary transformations that are applied to the basis used before the theoretical evaluation of the descriptors is performed using the population analysis we selected. These principles will be involved in the critical discussion of the reliability of this evaluation path, which might be discarded in certain situations. With this we hope to provide the reader and user with the possibility to rationalize one’s choice among the different possible derivations of the descriptors.

This report is organized as follows: after bringing a short reminder on density matrices and density functions, as well as the possibility to rewrite a density matrix in Löwdin’s symmetrically-orthogonalized basis, we will bring few mathematical tools that will be useful for writing the remaining expressions in the text in the most compact way possible. This technical part will also be the opportunity to justify certain choices for writing the descriptors in previous contributions. Afterward, the $\phi_{S}$ and $\varphi$ descriptors (among others) will be introduced, together with their population-analysis approximate value. Then, a general comment about the choice of the rules to diagnose in this contribution will be given, before the detailed diagnosis strategy is fully exposed in a dedicated section. Subsequently, the results of this diagnosis are introduced and discussed, before we conclude with the recommendations one could make to the future user of the related descriptors, with regards to the choice of the relevant evaluation method to use according to the system of interest. With this in hand, and taking into account certain restrictions related to the systems themselves, one should be able to consider the population-analysis approximation as being self-diagnosed — the population-analysis value of the descriptors is a partial measure of the reliability of the population-analysis approximation itself.

II Theoretical introduction

In this contribution, we will bring few modifications to the nomenclature used until now in Refs. etienne_toward_2014 ; etienne_new_2014 ; etienne_probing_2015 ; pastore_unveiling_2017 ; etienne_theoretical_2017 ; etienne_charge_2018 . These amendments are fully detailed in Appendix VIII.1 with explanations about the reason why they were used then, and why they will no longer be employed.

II.1 Brief reminder on density matrices and density functions

In a finite, orthonormal basis $\left(\chi_{r}\right)_{r\in\llbracket 1,L\rrbracket}$ of one-particle spinorbital spatial parts, from a one-particle reduced density matrix (1–RDM cioslowski_many-electron_2000 ; coleman_structure_1963 ; davidson_reduced_1976 ; lowdin_quantum_1955 ; mcweeny_recent_1960 ; mcweeny_methods_1992 ; mcweeny_density_1959 ; mcweeny_density_1961 ) $\bm{\gamma}$, one can unequivocally give the expression of the corresponding one-body charge density (1–CD) function $n$ defined on $\mathbb{R}^{3}$

[TABLE]

If the spatial part of every spinorbital function is expanded in the normalized but non-orthogonal basis of atomic functions $\left(\phi_{r}\right)_{r\in\llbracket 1,K\rrbracket}$ as

[TABLE]

where the $r^{\mathrm{th}}$ column of C stores the so-called LCAO (standing for Linear Combination of Atomic Orbitals) coefficients of the $r^{\mathrm{th}}$ spinorbital, the 1–CD function can be recast as

[TABLE]

with

[TABLE]

The so-called overlap matrix S stores the spatial overlap between pairs of atomic functions, and its knowledge allows the straightforward evaluation of the integral of the 1–CD function over the whole space: Writing

[TABLE]

implies

[TABLE]

For any real number $x$, we have $\textbf{S}=\textbf{S}^{1-x}\textbf{S}^{x}$ and, according to the invariance of the trace of a matrix product upon cyclic permutations,

[TABLE]

In particular for ($x=1/2$), we have that $\textbf{S}^{1/2}\textbf{P}\textbf{S}^{1/2}$ corresponds to the 1–RDM in the so-called Löwdin symmetrically-orthogonalized basis of atomic orbitals lowdin_nonorthogonality_1950 ; kashiwagi_generalization_1973 ; scofield_note_1973 ; scharfenberg_new_1977 ; aiken_lowdin_1980 ; attila_szabo_modern_1996 ; mayer_lowdins_2002 ; mayer_simple_2003 ; mayer_lowdin_2004 ; bruhn_lowdin_2006 ; szczepanik_several_2013 .

II.2 Tools for writing the density-based descriptors

In order to simplify the notations in this contribution, we will define two functionals, that will be used throughout this paper.

The first functional is constructed by considering a real-valued function $f$ defined on $\mathbb{R}^{3}$. Its value in every point of space is used for constructing two others, $f_{-}$ and $f_{+}$, according to the sign of the $f$ entries:

[TABLE]

Note that, for every r in $\mathbb{R}^{3}$,

[TABLE]

Averaging the integrals of $f_{+}$ and $f_{-}$ over $\mathbb{R}^{3}$ gives the $G[f]$ functional

[TABLE]

This is basically how four of our density-based descriptors were defined in Refs. etienne_new_2014 ; etienne_probing_2015 ; pastore_unveiling_2017 (their relation to those reported here is discussed in Appendix VIII.2). We can see however that $G[f]$ can be reduced to a simple integral, according to

[TABLE]

The reason why those descriptors were computed as in Eq. (9) is that, although the functions considered were satisfying

[TABLE]

in practice the integrals of $f_{+}$ and $f_{-}$ were averaged together in order to decrease the numerical imprecision inherent to their practical evaluation.

For the second functional used in this report, we will consider any square, real matrix $\textbf{Q}\,\in\,\mathbb{R}^{m\times m}$, and define

[TABLE]

where this time the integration has simply been replaced by a trace. The “$\circ$” symbol stands for the Hadamard product (sometimes called the entrywise product), defined as

[TABLE]

and $I_{m}$ is the $m\times m$ identity matrix, so that $\textbf{Q}\circ I_{m}$ is a diagonal matrix whose nonzero elements are the diagonal elements of Q only, and

[TABLE]

so that in practice, $H[\textbf{Q}]$ simply isolates the diagonal of Q in order to average its positive and negative diagonal entries sums. For the same reason as for $G$ this procedure was used in Ref. etienne_probing_2015 to compute most of the matrix-algebraic approximations to the descriptors (their link to those derived here is highlighted in Appendix VIII.2) though, again, we see that we can simplify its expression:

[TABLE]

In this contribution, $G$ and $H$ will be used to revisit and generalize various descriptors of electronic transitions and to express their derivation using two types of procedures: numerical integration on a Cartesian grid ($G$) and matrix-algebraic approximation ($H$).

Finally, note that the $G$ and $H$ functionals are similar mappings in a continuous and discrete basis, respectively.

II.3 The target density-based descriptors

We start by considering the so-called detachment and attachment 1–RDMs. Those can be seen as the hole and particle density matrices, in a one-hole/one-particle representation of an electronic transition. They are constructed by diagonalizing the one-body difference density matrix (1–DDM, $\bm{\gamma}^{\Delta}$), i.e., the matrix obtained, in the spinorbital space, by subtracting the ground electronic state 1–RDM ($\bm{\gamma}^{0}$) to the one corresponding to the $q^{\mathrm{th}}$ excited state ($\bm{\gamma}^{q}$)

[TABLE]

[TABLE]

where the columns of M are called natural difference orbitals. This operation also produces eigenvalues, which are often referred to as transition occupation numbers, and can be sorted with respect to their sign:

[TABLE]

This operation is similar to those involved in the construction of $G$ and $H$. Backtransformation provides the $L\times L$ real-valued detachment (\textcolorblack$\bm{\gamma}^{d}$) and attachment (\textcolorblack$\bm{\gamma}^{a}$) 1–RDMs

[TABLE]

Once projected in the Euclidean space, they give the detachment and attachment 1–CD functions:

[TABLE]

and

[TABLE]

which are both positive or zero everywhere in $\mathbb{R}^{3}$. We have, for every real number $x$,

[TABLE]

This $\vartheta$ represents the quantity of charge that has been involved in the electronic transition. It is this integral of the detachment and attachment densities that will be used as a normalization factor for the quantities defined below. Note that in previous publications etienne_toward_2014 ; etienne_new_2014 ; etienne_probing_2015 ; etienne_theoretical_2017 ; pastore_unveiling_2017 ; etienne_charge_2018 the integral of the detachment and attachment densities was averaged to give $\vartheta$, for the same reasons as those described in the penultimate paragraph of Appendix VIII.1.

Note that here “A” and “D” should not be confused with electron “Acceptor” and “Donor”.

We can now write two functions involving the detachment/attachment densities

[TABLE]

where

[TABLE]

and write the detachment/attachment overlap integral etienne_toward_2014

[TABLE]

as well as the non-normalized etienne_new_2014 ($\chi$) and normalized etienne_probing_2015 ($\varphi$) detachment/attachment-based charge-displacement descriptor

[TABLE]

$\varphi$ measures the fraction of detachment/attachment density contributing to the net transferred charge. $\phi_{S}$ and $\varphi$ were then combined in Ref. etienne_probing_2015 into the normalized, general $\psi$ descriptor

[TABLE]

jointly describing the range of a charge separation, and the amount of charge transferred during an electronic transition.

II.4 Alternative approach: the population analysis

Since it is not possible to solve analytically the $\phi_{S}$ and $\varphi$ integrals reported above, only approximate values can be obtained. In this paper we extend the matrix-algebraic procedure proposed in Ref. etienne_probing_2015 , where numerical integrations on Cartesian grids usually used for computing the density-based descriptors $\phi_{S}$ and $\varphi$ were approximated by performing a detachment/attachment population analysis (these two methods were loosely called direct- and Hilbert-space derivations in Ref. etienne_probing_2015 , the latter being employed because one-body orbital wave function integrals were part of the derivation procedure): If the atomic functions are centered on atoms, a known procedure is to approximate the atomic population, i.e., the partial charge that an atom bears for a given electronic distribution, through a Mulliken population analysis by multiplying a density matrix with the basis set overlap matrix S, and to consider the $k^{\mathrm{th}}$ diagonal entry of this product of matrices as the contribution of the $k^{\mathrm{th}}$ atomic orbital to the electronic population of the system. The actual approximation then consisted in considering the fact that the basis of atomic functions was local, with the Gaussian basis functions grossly covering a given region of space, so that manipulating (in particular, subtracting or multiplying) the diagonal entries of the matrix product of the atomic-space detachment and attachment density matrices with S was an approximation to the manipulation of elementary fractions of a charge density in a given region of space. Another possible approximation reported in Ref. etienne_probing_2015 is the Löwdin detachment/attachment population analysis, where the density matrices are contracted to the left and to the right by the square root of S to orthogonalize the basis set while keeping as much as possible its local character. In fact, any combination of real-valued exponents for the left and right multiplication can arbitrarily be used, as long as the sum of the exponents is equal to unity, in order to satisfy Eq. (6).

Though the choice of the $x$ exponent seems quite arbitrary, it was demonstrated in Appendix VIII.3 that only a Löwdin-like population analysis produces atomic populations that are not susceptible of being negative or greater than the maximum allowed occupancy, which would be unphysical. Therefore, though the matrix-algebraic definition of the approximate descriptors are given below for any value of $x$, in the next section of this paper we will solely use the Löwdin-like scheme (denoted by the “$\ell$” symbol), unless the contrary is explicitly stated.

The limitations of the approximation made by using symmetrically-orthogonalized orbitals populations as if we were mapping these charges to point charges in space when computing our descriptors is diagnosed and discussed in this report.

According to Eq. (6), since the integral of a 1–CD function is unequivocally equal to the trace of the product of S with the density matrix corresponding to the charge distribution when this matrix is known, the $\vartheta$ index is computed without approximation as the trace of the unrelaxed difference density matrix multiplied by S, using $x=0$ for sparing the effort of diagonalizing S.

The population analysis approximations to the detachment/attachment spatial overlap $\phi_{S}$ are

[TABLE]

with $x\in\mathbb{R}$, and $\bm{\Omega}\in\mathbb{R}^{K\times K}$, with, for every $(\mu,\nu)$ in $\llbracket 1,K\rrbracket^{2}$,

[TABLE]

where the last factor is called “Iverson’s brackets”. This factor takes a value of zero when the condition between the brackets is not true, and a value of one when the condition is true. Note that the $\phi_{S}^{x}$ functional should not be confused with $H$ in Eq. (15). Note also that the Iverson brackets are imposed in Eq. (29) for every matrix element, though this condition is necessary only for the diagonal elements, because the non-diagonal entries are removed by the Hadamard product in the definition of the $H$ functional. This justification is transferrable to the derivation of each other quantity below.

Due to the fact that solely the Löwdin analysis produces populations that cannot be negative, the Iverson brackets above are always equal to one for the diagonal elements of $\displaystyle{\textbf{S}^{1/2}\textbf{D}\textbf{S}^{1/2}}$ and $\displaystyle{\textbf{S}^{1/2}\textbf{A}\textbf{S}^{1/2}}$, so the Löwdin-derived approximation to $\phi_{S}$, written here $\phi_{S}^{\ell}$, reads

[TABLE]

If we are now interested in the computation of $\varphi$ with matrix algebra, we find that the $H$ functional defined in Eq. (15) becomes very useful, since we can write the generalised (i.e., for any real $x$) approximation to $\varphi$ as

[TABLE]

with, for every $(\mu,\nu)$ in $\llbracket 1,K\rrbracket^{2}$,

[TABLE]

where $\bm{\Delta}$ is the unrelaxed difference density matrix in the atomic space ($\textbf{C}\bm{\gamma}^{\Delta}\textbf{C}^{\dagger}$). When the Löwdin analysis is concerned for approximating $\varphi$, we see that the two Iverson brackets are equal to unity, leading to

[TABLE]

II.5 Beyond the unrelaxed picture of electronic transitions

If the description of the excited state electronic structure involves a pseudo-orbital relaxation through the addition of a Z-vector ($\bm{\gamma}^{\mathrm{Z}}$ in the spinorbital space) to the 1–DDM, one obtains the “relaxed” difference, detachment, and attachment 1–RDMs. Those are obtained using the same procedure as in Eqs. (16) to (19). The corresponding 1–CDs $n^{\prime}_{\Delta}$, $n^{\prime}_{d}$ and $n^{\prime}_{a}$ are also constructed using the same rules as in Eqs. (20), (21), and (23). Finally, the integral of the relaxed detachment/attachment 1–CDs, computed in a way similar to the one reported in Eq. (22), will be written $\vartheta^{\prime}$.

In a recent contribution pastore_unveiling_2017 we were interested in describing the nature of this pseudo-orbital relaxation. For this purpose, we designed three quantities: the first one is the difference between the rearranged (detached/attached) charge in the relaxed and unrelaxed schemes,

[TABLE]

We also have looked at the relaxation-induced detachment/attachment rearranged charge ($\bm{\nu}\in\left\{\textbf{D},\textbf{A}\right\}$)

[TABLE]

($\alpha^{\ddagger}_{\bm{\Gamma}}$ and $\alpha^{\ddagger}_{\bm{\Lambda}}$ in the original publication, see Appendix VIII.1), which corresponds to the separate impact of the Z-vector on the detachment and attachment densities. Finally, the third quantity used for unveiling the nature of the post-linear response pseudo-orbital relaxation, is the charge effectively displaced by the relaxation process:

[TABLE]

with

[TABLE]

This $n_{\mathrm{Z}}$ density function actually corresponds to the projection of the Z-vector in the Euclidean space.

If we now turn to the population analysis computation of these unnormalized descriptors we see that the $\zeta^{\mathrm{Z}}$ quantity reported in Eq. (37) can also be approximated using $H$. First, we rewrite the $\bm{\Delta}^{\mathrm{Z}}$ we met in Eq. (37)

[TABLE]

where $\bm{\Delta}^{\prime}$ is the atomic-space relaxed 1–DDM, $q$ is the excited state number (as in Eq. (16)), and each P is an atomic-space 1–RDM, related to the spinorbital-space 1–RDM by Eq. (4). With this, we find

[TABLE]

with, for every $(\mu,\nu)$ in $\llbracket 1,K\rrbracket^{2}$,

[TABLE]

Our considerations about Löwdin’ scheme for population analysis lead us to write

[TABLE]

Finally, as far as the $\alpha^{\ddagger}_{\textbf{D}}$ and $\alpha^{\ddagger}_{\textbf{A}}$ descriptors are concerned, we give here the following approximation to their value:

[TABLE]

with, for every $(\mu,\nu)$ in $\llbracket 1,K\rrbracket^{2}$,

[TABLE]

which is given for $\bm{\nu}\in\left\{\textbf{D},\textbf{A}\right\}$. The Löwdin scheme returns

[TABLE]

II.6 Diagnosis

Each descriptor that has been given an approximation with population analysis above is either derivable exactly ($\vartheta$, $\vartheta^{\prime}$, and $\lambda^{\mathrm{Z}}$), or by using the $H$ rule ($\varphi^{\ell}$, $\zeta^{\mathrm{Z}}_{\ell}$, $\alpha^{\ddagger}_{\bm{\nu},\ell}$) or the $F$ one ($\phi_{S}^{\ell}$). Therefore, the diagnosis performed for this report actually concerns the $H$ and $F$ rules. In particular, for the $H$ rule we selected to investigate the $\varphi^{\ell}$ approximation, since the $\varphi^{\ell}$ number can be combined to $\phi_{S}^{\ell}$ to produce $\psi^{\ell}$ without performing any additional matrix manipulation. It is assumed here that since the evaluation of $\zeta^{\mathrm{Z}}_{\ell}$ and the two $\alpha^{\ddagger}_{\bm{\nu},\ell}$ quantities follows computational schemes identical to the one of $\varphi^{\ell}$, our diagnosis on the latter also holds for the former.

Finally, note that due to Iverson’s brackets in the expression of the matrix elements of $\bm{\Omega}_{x}$, $\bm{\mu}_{x}$, $\bm{\Xi}_{x}$, and $\bm{\varpi}_{x}$, negative populations are ignored in the computation of the related descriptors, so that there might be occurrences in which the trace of the matrices on the elements of which Iverson’s brackets build conditions is not preserved since these negative populations can be overbalanced in other (positive) diagonal entries that are not ignored by Iverson’s brackets. This also means that these negative entries, which actually contribute to the value of the relevant 1–CD in every point of space, are conveniently ignored as they would appear below a square root in $H$. Fortunately, such case never occurs with the Löwdin scheme we use in practice and in this paper (see Appendix VIII.3).

III Diagnosis strategy

Our diagnosis strategy is summarized in Table 1. All the calculations reported in this paper were performed in vacuum, using the Gaussian16 (revision A03) package g16 . A set of fifty-nine metal-free molecules of varying size (from diatomics to complex molecular systems) was used in order to perform our diagnosis. The Cartesian grids were generated using the Cubegen utility from Gaussian16, and the electronic-structure analyses (numerical integration - NI - of the Cartesian grids, and population analyses - PA) were performed using a homemade code, to be released soon.

In order to determine to which extent both methods are consistent, different variables have been chosen for comparing NI and PA values: the system size, the charge transfer (CT) nature, the basis set (BS) size and nature, and the type of exchange-correlation (xc) functional used in time-dependent density functional theory.

Through this study, the influence of the Cartesian grid density of points has been evaluated by comparing medium (6 points/bohr for each direction) grids to fine (12 points/bohr for each direction) ones.

Our first molecular test-set (S1) is taken directly from Ref. etienne_toward_2014 and is constituted by a series of conjugated molecules with nitro (-NO2) acceptor and dimethylamino (-NMe2) donor terminal groups, separated by a varying number (ranging from one to five) of bridge subunits. The level of theory used here for our computations is exactly the one from the reference cited above.

One of the spacers is an oligo-ethylene for the nIII group of molecules. The oligo-ethylene might separate two donor groups for what we will call the nIII-NMe2-NMe2 series here, or two nitro groups (nIII-NO2-NO2), or one of each (nIII-NMe2-NO2).

If the donor and the acceptor are separated by a varying number of phenyl moieties, we have the nVI-c series.

Finally, in the nV-X series, the dimethylamino and nitro groups are connected through an heterocyclic oligomer (with five-member rings containing one oxygen, sulfur, or selenium) or through oligo-pyrrols.

The equilibrium geometries of the molecules belonging to the S1 set were obtained using the PBE0/6-311G(d,p) level of theory. Vibrational frequencies have been calculated to verify that the optimized geometries are actual minima. The excited-state calculations were done at the PBE0/6-311++G(2d,p) level of theory.

The particular interest of this set lies in the possibility to investigate electronic transitions of different nature through a series, with molecules of different size, hence calculations with basis sets of different size, so that the only stable variable is the xc-functional through the diagnosis. When switching the heteroatom in the nV-X series from oxygen to sulfur or selenium, one can also see how an increase of the number of electrons and basis functions for molecules of comparable size influences the precision of the PA approximation. This will be our first diagnosis, denoted $d_{1}$ in the results section, as well as in Table 1 below.

In addition to this, excited-state calculations were performed on the nV-X series using a small- (6-31G) and intermediate-sized (6-31+G(d,p)) BS in order to check the BS convergence for the PA values at a given geometry for this type of push-pull molecule ($d_{2}$).

The second test-set (S2) originates from Ref. peach_excitation_2008 and will be investigated using the same level theory as in the original paper, unless the contrary is explicitly stated. In particular, we will investigate the diatomic (N2, CO, and HCl) and formaldehyde molecules from S2 in order to diagnose the reliability of $F$ and $H$ rules for very small molecules ($d_{3}$). From S2 we will also focus on the acenes ($n+1$ fused benzene rings, with $n$ ranging from one to five) series, as we know that for such molecules the detachment and attachment densities will exhibit a substantial spatial overlap, no matter the number of subunits. This diagnosis will allow us to know whether, for such an extreme case of transition nature, there might be an influence from the size of the system on the accuracy of the PA approximation ($d_{4}$). Note that we also performed excited-state calculations on the acenes series using the B3LYP/6-31G and B3LYP/6-311++G(2d,p) levels of theory in order to assess the basis set convergence on the reference and approached descriptor values ($d_{5}$).

The last system that was used for this report is the pyridinium phenolate (PP) that was already investigated in Ref. pastore_unveiling_2017 . We used a rigid scan of the central dihedral angle, scanning the angle from zero to ninety degrees with a one-degree step, in order to see the evolution of the nature of the donor-acceptor intramolecular charge transfer. Such an approach allowed us to investigate a system of fixed size, at a given level of theory, with tunable CT character, for a one-variable diagnosis (d6) on a medium-sized molecular system. This study was also repeated including a variation of the level of theory by performing this ninety-one-point scan for a combination of fifteen split-valence Pople basis sets ditchfield_selfconsistent_1971 ; frisch_selfconsistent_1984 and six xc-functionals (the hybrid B3LYP becke_densityfunctional_1993 , PBE0 adamo_toward_1999 , and M06-2X zhao_m06_2008 , the long-range corrected CAM-B3LYP yanai_new_2004 and LC-$\omega$HPBE vydrov_importance_2006 , and the range-separated, dispersion-corrected $\omega$B97X-D chai_long-range_2008 ), so that we could select some geometries and evaluate the separate influence of the basis set size (and nature), as well as the type (hybrid/range-separated hybrid) of xc-functional (d7 and $d_{8}$). Configuration Interaction with single excitations (CIS hirata_configuration_1999 ) and Time-Dependent Hartree-Fock (TDHF hirata_configuration_1999 ) calculations were also performed on PP, combined with the fifteen BSs used in TDDFT hirata_configuration_1999 above.

The sketch of the molecules investigated in this report is given in \textcolorblackfigure S1.

IV Results

Due to our practical use of the descriptors, we report herafter their values rounded to two decimal places and, when comparing two evaluation paths, we will consider the consistency between the two methods as achieved when the results are consistent within two decimal places. Moreover, consistently with the discussion initiated in Ref. etienne_probing_2015 , we will consider a tolerance for the absolute deviation of $\pm$ 0.05 arbitrary units for the descriptors as the maximum allowed deviation between two methods.

IV.1 Influence of the NI grid density of points

When studying the whole set of molecular systems for this report, we quasi-systematically used the two density of points (6 and 12 points/bohr in each direction) for the sake of comparison and, according to our convergence/consistency criterion enonciated above, we noticed that the results on the descriptors are strictly converged for medium-density grids (6 points/bohr in each direction) and can be considered exact for the level of theory used, which means that the user who wishes to compute our descriptors can reliably spare the cost of producing fine grids eight times denser than the medium-density one. We drag the attention of the reader on the fact that when the numerical integration is performed, one has always to carefully check whether the integral of the detachment and attachment densities are strictly equal. We report in SI (see \textcolorblackTables S29 to S34) our data for the nIII group of molecules, and show that for the nIII-NMe2-NMe2 series and one excited state of the nIII-NMe2-NO2 series the standard grid sizes lead to a significant fraction of charge lost in the attachment. The integrals obtained with finer grids are also provided in order to prove that the fraction of charge lost comes from the grid size instead of the grid density of points. In order to overcome such issue the user should increase the grid size with the same density of points to get the right detachment/attachment integral. In this regard one can use the exactness of the trace of DS and AS to diagnose the quality of the integration grid, as it was done in \textcolorblackTables S29, S31, and S33.

IV.2 The $F$ and $H$ rules

Again, when considering the total set of data that we report in Supplementary Information for the fifty-nine molecules studied in this paper, we observed that, according to our accuracy/consistency criterion mentioned at the beginning of this section, when the PA-derived $\phi_{S}$ values were found to be reliable, so were the $\varphi$ values, though they were not obtained from the exact same rule. Therefore, from now on in this section we solely will discuss the $\phi_{S}$ and $\phi_{S}^{\ell}$ values for the sake of brevity, though the data in the tables and in the Supplementary Information will also contain the $\varphi$ and $\psi$ diagnosis.

IV.3 The oligo-acenes series

As expected, when studying the oligo-acenes series with our descriptors, we noticed that for the three xc-functionals used, i.e., PBE, B3LYP and CAM-B3LYP, in the $d_{4}$ diagnosis (see Table 1), the electronic cloud remains almost unperturbed by the electronic transitions, leading to $\phi_{S}$ values systematically greater than 0.9, which is quite an extreme case of detachment/attachment spatial overlap. This statement holds irrespective of the size of the molecule and of the BS (diagnosis $d_{5}$), or of the xc-functional used. Results in Table 2 and in \textcolorblackTables S1 to S5 show that in such an extreme case the reliability of the PA to approximate the value of $\phi_{S}$ cannot be firmly established. Indeed, we saw among the results both exact NI/PA matching and deviations beyond our tolerance criterion, which means that for such extreme situations the use of the PA-approximation $\phi_{S}^{\ell}$ is strongly discouraged. These conclusions also hold for the poly-acetylene series from the same S2 set of molecules (\textcolorblacksee Tables S38 to S40).

IV.4 Influence of the size of the system

When the nV-X and nVI-c sets of molecular systems are concerned, we globally observe that the molecular composition influences the nature of the CT, and that when performing the $d_{1}$ diagnosis the $\phi_{S}^{\ell}$ approximation can become (but is not systematically) problematic both for small systems and for electronic transitions characterized by a small $\phi_{S}^{\ell}$.

For the nV-X and nVI-c series (see Table 3, Tables S27 and S28, S35 to S37, and S44 to S47) we observe that for medium- to large-sized systems ($n>3$) the $\phi_{S}^{\ell}$ descriptor reproduces very accurately the $\phi_{S}$ value. Conversely, for small systems ($n=1$ or $n=2$) the PA-approximation is no longer reliable for $\phi_{S}$ values lower or equal to 0.52. This was also observed for the oligo-peptides of the S2 set of molecules (\textcolorblacksee Tables S6 to S8, S12 to S14, and S48 to S50). Note that in the S2 set of molecules the geometrical origin of the NI-PA deviation also applies for higher $\phi_{S}$ values, for example when the CO, N2, HCl, and H2CO molecules (see Tables S9 to S11, and S18 to S26), which are too small for considering the $\phi_{S}^{\ell}$ approximation as reliable, are concerned (diagnosis $d_{3}$).

The $d_{2}$ diagnosis (\textcolorblacksee Tables S27 and S28, S36 and S37, and S44 to S47) allowed us to conclude that the BS convergence was rapidly achieved for the $\phi_{S}$ value with medium-sized BS, and that when comparing the converged $\phi_{S}$ values to their PA-approximated values, the accuracy of the PA approximation was independent from the BS size. In other words, if the level of theory contains a BS of intermediate or large size, increasing the size of the BS did not improve the accuracy of the PA, \textcolorblackas it was already observed in the $d_{5}$ diagnosis. Following the same idea, we notice in Table 3 that replacing oxygens by sulfur or selenium might alter the $\phi_{S}$ value but not the accuracy of its PA-approximation, although in such substitutions we considerably modify the number of basis functions centered on the heteroatom positions before the symmetric orthogonalization.

IV.5 Study of medium-sized molecules

The accurate $\phi_{S}^{\ell}$ results for the S2-belonging DMABN and N-phenylpyrrole molecules (\textcolorblacksee Tables S15 to S17 and S41 to S43) should be considered with care. Indeed, we have seen when investigating the PP molecule that tuning the CT character of electronic transitions for a given molecule has an impact on the reliability of the PA-approximation ($d_{6}$ diagnosis).

The PP molecule excited states, when computed at the B3LYP/6-311++G(2d,p) level of theory, exhibit avoided crossing between the first and second singlet excited states and a conical intersection between the second and third excited states. This behaviour is very well reproduced by both the $\phi_{S}$ and $\phi_{S}^{\ell}$ computation (see the left part of figure 1), which lead us to state that one should not systematically discard $\phi_{S}^{\ell}$ results for the only reason that it is characterized by a sizeable gradient. However, one sees in figure 1 that when the $\phi_{S}^{\ell}$ value drops below 0.4 the values might or might not be reliable. For instance the first electronic transition at this level of theory is quite well described by $\phi_{S}^{\ell}$ when tuning down the twist angle from 90 to ca. 55 degrees while the trends completely differ afterward, i.e., when $\phi_{S}^{\ell}$ is dropping below 0.4. Significant divergences also occur for the second and third excited state as long as $\phi_{S}^{\ell}$ is maintained below 0.4, while non-negligible but tolerable deviations are observed when $\phi_{S}^{\ell}$ is lying between 0.4 and 0.5 for the second excited state.

Whilst there is neither avoided crossings nor conical intersections in the PP excited-state energy landscape for the first three excited states when computed at the CAM-B3LYP/6-311++G(2d,p) level of theory, the observations reported for the B3LYP xc-functional strictly apply to the CAM-B3LYP computations (see the right part of figure 1). Finally, the results in Table 4 show that at given geometries, since the xc-functional influences the CT character, it indirectly influences the accuracy of the PA-approximation. Indeed, while the xc-functional does not explicitly play any role in the definition of the $F$ rule, since we established that the $\phi_{S}^{\ell}$ approximation reliability was CT nature-dependent, overtuning the detachment/attachment overlap down by changing the xc-functional at certain geometries might result in $\phi_{S}$ values that cannot be approximated by Löwdin’ symmetrically-orthogonalized orbitals population analysis. Note however that when turning to the $d_{8}$ diagnosis the conclusions related to the CAM-B3LYP computation of PP excited states are transferrable to the other two range-separated xc-functionals, LC-$\omega$HPBE and $\omega$B97X-D (as well as the CIS and TDHF results). On the other hand the B3LYP results also hold for the other hybrid xc-functional PBE0 while M06-2X reproduces the electronic-structure analysis of range-separated hybrid xc-functionals, which is a known behaviour already reported in Ref. etienne_fluorene-imidazole_2016 . All the data for the six xc-functionals and the data related to the CIS and TDHF excited-state computations can be found in \textcolorblackTables S51 to S170 and figures S3 to S378. B3LYP and CAM-B3LYP data are also reported in Table 4 where two PP twist angles, 80∘ and 0∘, have been selected, as they allow us to highlight the particular behaviour of B3LYP for this molecule with a low dihedral angle.

Finally, as it is reported in Table 4 for six basis sets, and in SI for the nine others that were used in this investigation ($d_{7}$), again, once the convergence on the basis set is achieved for the $\phi_{S}$ value, the accuracy of the $\phi_{S}^{\ell}$ approximation is stable. This allows us to conclude that if the BS used is reliable for the excited-state computation, its size will not influence the accuracy of the $\phi_{S}^{\ell}$ approximation.

V Discussion

Our approximation here is to consider the Löwdin symmetrically-orthogonalized orbitals as being centered in one point and multiplying or subtracting diagonal entries of $\textbf{S}^{1/2}\textbf{D}\textbf{S}^{1/2}$ and $\textbf{S}^{1/2}\textbf{A}\textbf{S}^{1/2}$ as if we were manipulating point charges with precise locations (centered on the atomic positions), so that Löwdin’ symmetric orthogonalization tails are ignored in such a model. However, as we saw, there are few cases in which our approximation is limited.

In a diatomic molecule for example: If the atomic orbitals are indeed centered on atomic positions, the orthogonalization tail for the transformation of the orbital centered on the first atom, will be centered on the position of the second one (see figure in chapter 3 of Ref. mayer_simple_2003 ). More generally, Löwdin’ symmetrically-orthogonalized orbitals have orthogonalization tails centered on the position of other atoms, and can no longer be considered as strictly centered on one position, unlike the original atomic orbitals.

Therefore, for small molecular systems, the presence of these orthogonalization tails introduces a bias in the picture of the electronic distribution given by the population analysis we wish to use. Since the number of atoms is reduced, amongst the atomic orbitals there is a large percentage of them that is susceptible of having a significant overlap, and the orthogonalization tails are distributed on few atoms in a close region of space, so the deviation from our “local population” model (i.e., ignoring the orthogonalization tails), introduces a bias of significant relative importance in the population analysis.

The bias can also be significant for very high values of $\phi_{S}^{\ell}$ (greater than 0.90), since in that case there are numerous contributions from many pairs of atomic orbitals that are to be considered, which increases the number and importance of the orthogonalization tails that are ignored in the model, or more generally there is a collective effect observed from the sum of the individual small contributions to the deviation from the exact integral value.

For low values of $\phi_{S}^{\ell}$ (less than 0.5) the problem is quite different: Such low value for the detachment/attachment Löwdin population overlap means that only few diagonal product amplitudes are non-negligible. Since these contributions are summed, introducing a bias on the value of a few number of entries of small amplitude can introduce a strong deviation on the total outcome of the evaluation of $F$.

VI Conclusion

We reported the study of two computational methods for evaluating two types of integrals related to the nature of molecular electronic transitions. The first type of integral is the spatial overlap between two density functions, while the second type separately sums the positive and negative entries of a difference-density function. The first integral evaluation method diagnosed is the projection of the one-body reduced density matrices corresponding to the density functions in the Euclidean space using Cartesian grids, and discrete summations. The second evaluation path is the population analysis of the related density matrices in the basis of symmetrically-orthogonalized orbitals. The diagnosis is related to two rules ($F$ and $H$) mapping the orbital space to a scalar quantity.

Our first conclusion concerns the first integral evaluation method: We observed that the precision of the results of our numerical integrations on Cartesian grids is converged for medium-density grids with six points per bohr for the three directions of space, and that the exactness of the trace of some relevant matrices allows the quality of the grid to be diagnosed.

When the diagnosis of the population analysis is concerned, we first noticed that the behaviour of the $F$ and $H$ rules was very similar, and decided to focus the analysis on the $F$ rule leading to the Löwdin PA approximation ($\phi_{S}^{\ell}$) to the detachment/attachment spatial overlap integral $\phi_{S}$.

We then deduced from an extensive diagnosis of the $F$ rule that its accuracy depends on the size of the target molecular systems, and on the nature of the electronic transitions considered. Basically, the PA-approximation to $\phi_{S}$ is limited to medium- and large-sized molecular systems. In addition to that, $\phi_{S}^{\ell}$ values below 0.5 (for small and medium-sized molecules) and above 0.9 for any type of molecule should not be trusted, so that numerical integrations should be performed in such occurrences.

We also noticed that large $\phi_{S}^{\ell}$ gradients are not necessarily corresponding to large errors in the population analysis, and that globally the accuracy of the PA-approximation is independent from the size of the basis set used, apart for very small basis sets that in any case would not be used for the accurate computation of excited-state properties in general.

Supplementary Information

Tables S1 to S170, and Figures S1 to S378, are available at the \textcolorbluearxiv.org/abs/1902.05840 address using the “Ancillary files” section.

VII Tables and Figures

VIII Appendices

VIII.1 Nomenclature and comments

The modifications related in this Appendix concern the nomenclature of objects, functions, and quantities, taken from Refs. etienne_toward_2014 ; etienne_new_2014 ; etienne_probing_2015 ; pastore_unveiling_2017 ; etienne_theoretical_2017 ; etienne_charge_2018 .

Since we are not dealing anymore with centroids of charge, the use of “$\bm{\xi}=(\xi_{1},\xi_{2},\xi_{3})$” for indicating the position vector $\textbf{r}=(x,y,z)$ is no longer necessary for writing compact formulas, so any three-dimensional position vector has been denoted as “r” in this paper. Density functions were denoted with an “$n$” here rather than the “$\varrho$” used in Refs. etienne_toward_2014 ; etienne_new_2014 ; etienne_probing_2015 .

When detachment and attachment 1–RDMs are used in the atomic space, those are no longer written $\bm{\Gamma}$ and $\bm{\Lambda}$, but simply D and A respectively. The previous nomenclature was used in order to avoid any confusion with the terms “Donor” and “Acceptor”, commonly used in molecular science for talking about fragments that can donate or accept electron(s).

Note also that in the previously cited contributions, some quantities were given a definition corresponding to their practical evaluation, like the $\tilde{\chi}$ and $\tilde{\varphi}$ descriptors, which used to take a $\sim$ when derived from the detachment/attachment density difference, in contrast to their evaluation from the excited/ground state density difference. Though in practice the values might differ due to some numerical imprecision in their computation, which explains the use of this symbol to highlight how the quantity was actually computed, this symbol has been dropped according to the fact that the quantities derived from the two derivation schemes are, by construction, identical in principle.

Finally, when relaxed 1–DDMs were used in Ref. pastore_unveiling_2017 , the related quantities were labeled with an “R” while the unrelaxed ones were labeled with a “U”. In this contribution, the “relaxed” quantities were simply labelled with a prime while the unrelaxed ones were bearing no distinctive mark.

VIII.2 Correspondence between some descriptors related to G and H

Around G

The four density-based descriptors previously derived using a $G$–like functional were $\zeta^{\mathrm{Z}}$ (unchanged here), as well as $\tilde{\chi}$, $\tilde{\varphi}$, and $\alpha^{\ddagger}_{\bm{\nu}}$ with $\bm{\nu}\in\left\{\bm{\Gamma},\bm{\Lambda}\right\}$, renamed respectively ${\chi}$, ${\varphi}$, and $\alpha^{\ddagger}_{\bm{\nu}}$ with $\bm{\nu}\in\left\{\textbf{D},\textbf{A}\right\}$ in this contribution.

Around H

The $\tilde{\chi}_{1}$, $\tilde{\varphi}_{1}$, and $\phi_{S}^{\Omega}(1)$ in Ref. etienne_probing_2015 have been renamed respectively $\chi^{\ell}$, $\varphi^{\ell}$, and $\phi_{S}^{\ell}$ here, while $\tilde{\chi}_{2}$, $\tilde{\varphi}_{2}$ and $\phi_{S}^{\Omega}(2)$ are not explicitely mentioned in this paper, but correspond to $\chi^{(x=0})$, $\varphi^{(x=0})$, and $\phi_{S}^{(x=0)}$ according to the current nomenclature.

VIII.3 The particular case of Löwdin population analysis

In this appendix we demonstrate that when generalizing the orbital population analysis, only the Löwdin model prevents us from deriving unphysical atomic orbital populations (negative or greater than the maximal allowed occupancy).

For this purpose we will show that for a two-electron system with a two-orbital model, any other population analysis can, in principle, lead to an unphysical atomic population. In this appendix, any matrix element $(\textbf{L})_{ij}$ will be written $L_{ij}$ for the sake of readability. This applies to the matrix elements of C, P, and S.

Let us consider a basis of two real-valued normalised but not orthogonal atomic orbitals, $\left(\phi_{1},\phi_{2}\right)$. Their overlap matrix reads

[TABLE]

If $S$ was allowed to take the $-1$ or the $1$ value, the atomic orbitals would no longer be linearly independent.

Let us now consider two real-valued molecular orbitals, $\varphi_{1}$ and $\varphi_{2}$, that form a complete and orthonormal basis with respect to the space spanned by the atomic orbitals basis set. Their spatial parts, $\chi_{1}$ and $\chi_{2}$, expand as

[TABLE]

The expansion coefficients are real (note that the present derivation can be generalised easily to complex-valued atomic and molecular orbitals). The matrix of coefficients,

[TABLE]

satisfies a generalised orthonormality relationship with respect to the descriptor defined by S, expressed as

[TABLE]

The molecular orbitals are assumed to be linearly independent, hence C is invertible and

[TABLE]

S is a Gram matrix. Therefore, it is positive definite. Its eigenvalues are $1\pm S>0$. A possible eigendecomposition of it is

[TABLE]

From this, we can define, for any real $x$,

[TABLE]

where $(1\pm S)^{x}>0$. In particular,

[TABLE]

[TABLE]

and

[TABLE]

Among the possible square roots of S, the only positive definite one, which is unique and termed “principal” square root of S, is chosen here, and written $\textbf{S}^{1/2}$.

The generalised orthonormality relationship fulfilled by C can be recast as follows (Löwdin symmetric orthonormalisation),

[TABLE]

Thus, C being square and invertible,

[TABLE]

Note that this would not hold if the molecular orbitals set were not complete (i.e., if C were rectangular). $\textbf{S}^{1/2}\textbf{C}$ is a real unitary matrix (i.e., orthogonal) and can be chosen as a rotation matrix (this parametrisation is not unique),

[TABLE]

where $\theta$ is real. As a consequence,

[TABLE]

so that

[TABLE]

Let us now consider a two-electron closed-shell wave function in the form a single Slater determinant,

[TABLE]

where “s” collects the four spin-spatial variables. The corresponding 1–CD reads

[TABLE]

where, again, $\chi_{1}$ is the spatial part of $\varphi_{1}$. The $n$ function satisfies

[TABLE]

The corresponding density matrices, P and the duodempotent $\bm{\gamma}$, satisfy

[TABLE]

[TABLE]

The Mulliken gross orbital populations are defined as

[TABLE]

Thus,

[TABLE]

From the expression of $\bm{\gamma}$ and C given above, and according to Eq. (4) we get

[TABLE]

[TABLE]

This can be formulated using expansions in terms of Pauli matrices

[TABLE]

as

[TABLE]

where we identify

[TABLE]

We can also rewrite

[TABLE]

Generalised orbital-population analyses can also be considered, based on $\textbf{S}^{x}\textbf{PS}^{1-x}$ with $0\leq x\leq 1$. We get

[TABLE]

[TABLE]

with

[TABLE]

The diagonal entries (orbital populations), thus read

[TABLE]

and

[TABLE]

For $x=0$ or $x=1$ (Mulliken), the populations satisfy

[TABLE]

They can be negative, depending on the values of $S$ and $\theta$. Large values of $|S|$ make this situation likely to occur (see the top of figure S2). It might also happen that the orbital population exceeds the maximal allowed orbital population (see the bottom of figure S2). These two situations can occur also for any other value of $0\leq x\leq 1$ except for $x=1/2$. In the latter case, $\textbf{S}^{1/2}\textbf{PS}^{1/2}$ is positive semidefinite and its diagonal entries are all nonnegative. In particular here,

[TABLE]

This justifies that only Löwdin-like population analysis was used in our calculations.