Rivo Mamy

Dianzinga

N -representable density matrix perturbation theory

Keywords: Matrice densité, réponses électroniques, champ auto-cohérent, croissance linéaire N-representable density matrix perturbation theory Density matrix, electronic response functions, self-consistent field, linear scaling

Titre : Théorie des perturbations en matrice densité Nreprésentable

Résumé :

Alors que les approches standards de résolution de la structure électronique présentent un coût de calcul à la puissance 3 par rapport à la complexité du problème, des solutions permettant d'atteindre un régime asymptotique linéaire, O(N), sont maintenant bien connues pour le calcul de l'état fondamental. Ces solutions sont basées sur la "myopie" de la matrice densité et le développement d'un cadre théorique permettant de contourner le problème aux valeurs propres. La théorie des purifications de la matrice densité constitue une branche de ce cadre théorique. Comme pour les approches de type O(N) appliquées à l'état fondamental, la théorie des perturbations nécessaire aux calculs des fonctions de réponse électronique doit être révisée pour contourner l'utilisation des routines coûteuses. L'objectif est de développer une méthode robuste basée uniquement sur la recherche de la matrice densité perturbée, pour laquelle seulement des multiplications de matrices creuses sont nécessaires. Dans une première partie, nous dérivons une méthode de purification canonique qui respecte les conditions de N-representabilité de la matrice densité à une particule. Nous montrons que le polynôme de purification obtenu est auto-cohérent et converge systématiquement vers la bonne solution. Dans une seconde partie, en utilisant une approche de type Hartree-Fock, nous appliquons cette méthode aux calculs des tenseurs de réponses statiques non-linéaires pouvant être déterminés par spectroscopie optique. Au delà des calculs à croissance linéaire réalisés, nous démontrons que les conditions Nrepresentabilité constituent un prérequis pour garantir la fiabilité des résultats.

Résumé

Le Chapitre 1 résume les postulats et notions de la mécanique quantique nécessaires à l'introduction de la matrice densité à une particule comme variable fondamentale dans la résolution des équations de champ moyen auto-cohérent (self-consistent field -SCF), héritées de l'approximation mono-déterminantale de la fonction d'onde électronique. Nous montrons qu'il est possible de contourner la résolution de l'équation de Schrödinger -qui revient à résoudre un problème aux valeurs propres du type :

HC -CE = 0, (H, C) ∈ R M ×M , H = H t , C t C = CC t = I, E = diag{ϵ 1 ϵ 2 • • • ϵ M }-
par la résolution directe de l'équation de Liouville-von Neumann du type : HD -DH = 0, D ∈ R M ×M , D = D t , Tr{D} = N , D 2 = D, 0 < N < M , dont la seule variable est la matrice densité D. La détermination de la matrice densité via les vecteurs et valeurs propres, respectivement C et E, de la matrice hamiltonienne, H, est numériquement coûteuse puisqu'elle nécessite une étape de diagonalisation. Les ressources de calcul nécessaires à la réalisation de cette étape présentent une croissance à la puissance 3 par rapport à la complexité du problème, cette dernière étant généralement définie par le nombre d'états occupés, N , ou le nombre d'atomes, ou encore la taille des matrices M . La détermination directe de la matrice densité -sans passer par le calcul des états propres-peut être explicitée sur la base d'un principe de minimisation -sous contrainte-de l'énergie du système, qui dans le cadre de cette thèse, correspond à l'énergie Hartree-Fock. En se référant aux procédures standards de minimisation lagrangienne, il est montré que les propriétés d'idempotence et de conservation de la trace de la matrice densité sont nécessaires et suffisantes pour garantir l'unicité de la solution. Ces propriétés sont regroupées sous le terme générique de conditions de N -representabilité de la matrice densité à une particule. Dans ce travail, tous les calculs sont réalisés sur la base de l'approche Hartree-Fock semiempirique de Pariser, Parr et Pople (PPP). Les techniques d'accélération de la procédure SCF, notamment l'interpolation à paramètre constant (damping) et l'extrapolation par inversion directe du sous-espace des itérations (direct inversion of the iterative subspace -DIIS) sont également discutées. Un exemple d'application montre clairement les avantages de la méthode DIIS. xii Le Chapitre 2 présente les méthodes les plus couramment utilisées pour la résolution directe de la matrice densité. Elles peuvent être regroupées en deux familles : (i) les méthodes de minimisation, et (ii) les méthodes de purification, chacune s'appuyant partiellement sur les conditions de N -representabilité évoquées précédemment. Dans le premier cas, le polynôme de purification de McWeeny dérive du principe de minimisation des moindres carrés de la contrainte d'idempotence, couplé à un algorithme de descente de gradient. Une alternative proposée par Li, Nunes et Vanderbilt (LNV) est basée sur la minimisation de la fonctionnelle de l'énergie sous une contrainte faible d'idempotence, couplée à un algorithme de gradient conjugué. Dans leur formulation grand canonique, ces deux approches sont en mesure de garantir les conditions de N -representabilité si et seulement si le potentiel chimique est connu à l'avance. En d'autres termes, les énergies correspondant au dernier état occupé et premier état inoccupé doivent être déterminées au préalable. Notons que pour un taux d'occupation, θ = N/M = 50%, le potentiel chimique peut être évalué avec une certaine précision à partir des limites supérieure et inférieure du spectre des valeurs propres. Une autre catégorie de purification, qualifiée de canonique, est également présentée. Dans cet ensemble, il n'est plus nécessaire d'évaluer les valeurs propres internes du spectre de l'hamiltonien. Néanmoins, leur application implique de prendre en compte plusieurs facteurs de stabilité, ce qui peut limiter leur efficacité et par conséquent, aussi complexifier les algorithmes de calcul. C'est dans ce cadre qu'est développée la première originalité de notre travail : en introduisant une méthode de purification canonique simple et robuste qui s'affranchit des considérations heuristiques de ces prédécesseurs. Cette nouvelle variante est basée sur la reformulation lagrangienne du principe de minimisation de l'idempotence de McWeeny en introduisant une contrainte explicite sur la trace de la matrice densité de faible idempotence. De cette façon, la méthode de purification est auto-cohérente -l'ajustement a posteriori du pôlynome n'est plus nécessaire-et vérifie les conditions de N -représentabilité à chaque itération. Dans le cadre de l'approximation des liaisons fortes, une étude détaillée des différentes méthodes de purification canoniques est réalisée. Il est prouvé que les performances de cette nouvelle approche sont comparables aux méthodes heuristiques tout en montrant des propriétés intéressantes de convergence monotonique et variationnelle. Toujours dans ce même chapitre, après avoir résumé les approximations permettant l'application de l'algèbre linéaire creuse aux méthodes de minimisation et purification, des calculs SCF-HF-PPP à croissance linéaire sont réalisés sur une série de nanotubes de carbone. Les difficultés liées à l'application des méthodes de troncatures numérique et radiale sont étudiées. xiii Dans le Chapitre 3, une généralisation des méthodes de variation/perturbation de la matrice densité à une particule est proposée, puis étendue à la résolution des équations couplées-perturbées du champ auto-cohérent (coupled-perturbed self-consistent field -CPSCF). Il en résulte trois approches : (i) la méthode standard basée sur la résolution spectrale de la matrice de Fock qui sera considérée par la suite comme la référence, (ii) l'approche proposée par Kussmann et Ochsenfeld basée sur les relations de commutations généralisées et l'utilisation du gradient conjugué, et (iii) la méthode de Niklasson basée sur le développement par récursion de l'opérateur de Fermi-Dirac perturbé. Dans le cadre du troisième formalisme, deux méthodes de purification sont utilisées, dont notre pôlynome canonique auto-cohérent. Il est à noter que pour les deux premières approches, le calcul des fonctions de réponse d'ordre supérieur nécessite la connaissance des fonctions de réponse associées aux ordres inférieurs, alors que dans les cas des purifications toutes les matrices densités perturbées sont calculées simultanément. En d'autres termes, l'ordre zéro (non-perturbé) et l'ordre supérieur (cible), ainsi que tous les ordres intermédiaires (si nécessaires), sont déterminés au cours du même processus de purification. Pour toutes les approches mentionnées la résolution des équations CPSCF est accélérée par l'algorithme de la dérivée du DIIS (D-DIIS) introduit par Weber et Daul. En fin de chapitre, nous généralisons cet algorithme pour n'importe quel ordre de perturbation.

Le Chapitre 4 est consacré à l'application des différentes méthodes de perturbation de la matrice densité pour le calcul des fonctions de réponse électronique induites par l'application d'un champ électrique externe et statique. Ces fonctions de réponse permettent, entre autres, de déterminer des grandeurs accessibles via des mesures d'optique non-linéaire, comme la polarisabilité, et/ou la première et seconde hyperpolarisabilité. Dans ce même chapitre, le principe des méthodes de différence de champ fini (finite field difference -FFD) est résumé. Dans un premier temps, toutes les méthodes sont appliquées aux calculs des propriétés optiques d'une série de petites molécules π-conjuguées en utilisant l'approche CPSCF-HF-PPP. L'accord remarquable observé démontre la fiabilité des différentes implémentations. Dans un deuxième temps, les calculs sont étendus à des systèmes facilement réplicables tels que des hydrocarbures insaturés. L'erreur de chaque méthode de perturbation basée sur la matrice densité (ii-iii) par rapport à la méthode de référence (i) est évaluée en fonction de la taille du polymère. Les résultats démontrent que la purification canonique est la seule à conserver une précision remarquable, et ce quelque soit l'ordre de perturbation. Ce résultat, qui constitue la deuxième originalité de ce travail, est directement relié au respect des précieuses conditions de N -representabilité discutées dans les chapitres précédents. Toutefois il est à noter que pour les approches ne vérifiant pas explicitement ces contraintes, l'erreur observée reste acceptable puisqu'elle [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF]5) using a simple approach, the CDA and the DIIS optimization. In (a) is displayed the energy minimization during the iterative process. In (b) is represented the energy error during the iterative process using a logarithmic scale. . . TC2-CPSCF and CG-CPSCF, at zero (D), first (D (1) = D (x) ), second (D (2) = D (xx) ) and third (D (3) = D (xxx) ) orders. Results obtained for PPV+ of Set 1 with n = 28. . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Idempotency of the density matrices during the SCF iterations for the HPCP-CPSCF, TC2-CPSCF and CG-CPSCF, at zero (D), first (D (1) = D (x) ), second (D (2) = D (xx) ) and third (D (3) = D (xxx) ) orders. Results obtained for PPV+ of Set 1 with n = 28. . . . . . . . . . . . . . . . . . . 4.8 Frobenius norm of the error vector [cf. Eq. (3.69)] of the density matrices during the SCF iterations for the HPCP-CPSCF, TC2-CPSCF and CG-CPSCF, at zero (D), first (D (1) = D (x) ), second (D (2) = D (xx) ) and third (D (3) = D (xxx) ) orders. This example is for PPV+ of Set 1 with n = 28. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 First form of illustration for the converged density matrices using the HPCP, from zero to third order. nnz is the number of non zero elements at 10 -3 . This example is for PPV+ of Set 1 with n = 28. D (1) = D (x) , D (2) = D (xx) , D (3) = D (xxx) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Second form of illustration for the converged density matrices using the HPCP, from zero to third order. nnz is the number of non zero elements at 10 -3 . This example is for PPV+ of Set 1 with n = 28. D (1) = D (x) , D (2) = D (xx) , D (3) = D (xxx) . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Histogram of number of SCF iterations with respect to the size of the systems, for the four density matrix perturbation methods at first and second orders. The results are obtained for TPA+ and PPV+ of Set 1. . 4.12 Calculation time for the polarizability α xx and first hyperpolarizability β xxx as a function of number of atoms. Results are obtained for TPA+ (Set 1). For each density matrix perturbation method, the calculation time is pictured along with its fit. . . . . . . . . . . . . . . . . . . . . . .

Introduction

Sustained by the fast increase of investments in high performance computing technologies, quantum mechanics accuracy, as found in standard electronic structure methods, is about to reach the mesoscale within the next century. This implies the development of advanced parallelized programs including adapted theoretical frameworks and efficient numerical algorithms. At the same time, spectroscopies are improving in resolution and increasing in complexity with respect to the size of the probed systems, causing new difficulties for interpreting spectra, and new challenges for the theoreticians. There is no doubt in the importance of probing structure of matter at the atomic scale to establish clear relationships between the macroscopic properties and the atoms' arrangement in a sample. For that purpose, electromagnetic spectroscopies are banal experiments to obtain molecular fingerprint in physico-chemical analyses. Beyond their standard use, when we have no a priori knowledge on the system -but let's say the chemical composition for a material sample-analysis of the spectra may become tedious, or very difficult when dealing with disordered, amorphous or soft matter. Even if spectrocospies are continuously increasing their possibilities in resolution, support of theoretical prediction remains of primary interest for spectrum assignment and structure elucidation. On this way, the accuracy of the theoretical methods and the size of the investigated system consitute interrelated bottlenecks that need to be addressed. Accurate quantum methods based on the explicit resolution of the many-electron Schrödinger equation, where dynamic and/or static electron correlation can selectively be accounted for remain limited to a few dozens of atoms or less depending on the level of theory. If now, we are willing to sacrifice accuracy in order to resolve electronic structure for larger systems where the number of electrons is above a few hundreds, single-determinant theories, such as Hartree-Fock [1] (HF) or Kohn-Sham (KS) density functional theory [2] are, until now, the only relevant methods. In that case, the manyelectron Schrödinger equation is reduced to a mean-field one-electron equation, whose variational solutions are obtained by solving an eigenvalue problem. As a result, solutions are obtained by minimizing the electronic energy by means of the self-consistent field (SCF) procedure. Whatever the basis set used to expand the wave functions, SCF methods have a common limitation on the size of the problem in such way that they require calculation of the full or partial set of eigenstates of the Hamiltonian matrix at each iteration of the SCF. The computational effort related to the direct/iterative diagonalization techniques [3,4] or state-by-state conjugate gradient (CG) algorithm, [5,6] increases with the cube of the number of electrons. Over the last two decades, alternative methods which scale linearly with the size of the problem were proposed as solution to these standard energy minimizations. [7][8][9][START_REF] Kussmann | Linear-scaling selfconsistent field methods for large molecules[END_REF] These methods are based on the Kohn's principle of electronic structure nearsightedness, [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF][START_REF] Prodan | Nearsightedness of electronic matter[END_REF] which under certain conditions, eg. non-vanishing electronic gap, shows an exponential decay of the density matrix (DM) elements with respect to the distance. [START_REF] Des | Energy bands and projection operators in a crystal: analytic and asymptotic properties[END_REF][START_REF] Ismail | Locality of the density matrix in metals, semiconductors, and insulators[END_REF] On exploiting this natural property, that is by enforcing sparsity of the matrices using a truncation scheme, O(N ) can be achieved by replacing the diagonalization step with DM solvers along with sparse-matrix multiply (SpMM) algorithms. [4,[START_REF] Challacombe | A general parallel sparse-blocked matrix multiply for linear scaling scf theory[END_REF][START_REF] David R Bowler | Parallel sparse matrix multiplication for linear scaling electronic structure calculations[END_REF] Predictions of spectroscopic observables for molecules and solids rely on the solid approximation that the strengths of the electromagnetic radiations are negligible with respect to magnitude of the electron bonding allowing the safely use of the Rayleigh-Schrödinger wave function perturbation theory to compute the electronic response at any order. The first applications of this theory to SCF methods based on molecularorbital (MO) wave functions were introduced during the 60s for the computation of molecular properties such as magnetic susceptibility, [START_REF] Stevens | Perturbed hartree-fock calculations. i. magnetic susceptibility and shielding in the lih molecule[END_REF] static polarizabilities and force constants, [START_REF] Gerratt | Force constants and dipole-moment derivatives of molecules from perturbed hartree-fock calculations[END_REF][START_REF] Gerratt | Force constants and dipole-moment derivatives of molecules from perturbed hartree-fock calculations. ii. applications to limited basis-set scf-mo wavefunctions[END_REF] which are all related to second-order energy derivatives through the calculation of the first-order change of the wavefunctions with respect to the small perturbation. Similarly to the unperturbed case, variational solutions of the perturbed MOs are obtained by solving the so-called coupled-perturbed self-consistent field (CPSCF) equations. [START_REF] Thomsen | Calculation of molecular one-electron properties using coupled hartree-fock methods: I. computational scheme † part of this work was performed at the max-planck-institut für physik und astrophysik, münchen[END_REF][START_REF] Ditchfield | Self-consistent perturbation theory of diamagnetism: I. a gaugeinvariant lcao method for nmr chemical shifts[END_REF][START_REF] J_A Pople | Derivative studies in hartree-fock and møller-plesset theories[END_REF] These early developments based either on the perturbed MOs or mixed perturbed AOs-MOs are well-known to involve cumbersome matrix transformations. [START_REF] Frisch | Direct analytic scf second derivatives and electric field properties[END_REF]24] In 1962, McWeeny had already introduced the elegant formalism of the density matrix perturbation theory (DMPT), [25] which was extended to the CPSCF equations resolution by Diercksen and McWeeny[26] for the evaluation of π-electron polarizabilities using the Pariser-Parr-Pople model. This work has first inspired Moccia to generalize the McW-CPSCF equations resolution to non-orthogonal basis. [27,28] Perturbationdependent non-orthogonal basis implementation was then proposed by Dodds, McWeeny, Sadlej and Wolinski [START_REF] Dodds | Self-consistent perturbation theory: Generalization for perturbation-dependent non-orthogonal basis set † this research was partly supported by the institute of low temperatures and structure research of the polish academy of sciences[END_REF][30][31] for the calculation of atomic (hyper)-polarizabilities using HF method and gaussian-type orbitals. The advantages of the McWeeny's approach over MO/AO-CPSCF have been clearly outlined, for instance, in the seminal article of Wolinski, Hinton, and Pulay [32] dealing with the calculation of magnetic shieldings as measured in nuclear magnetic resonance (NMR) spectroscopy.

All the methods mentioned above have in common two limitations which narrow their applicability to few hundred atoms system at most, which are: (i) the unperturbed eigenstates are required prior the evaluation of the perturbed quantities, (ii) CPSCF equations resolution involves dense matrix multiplications which scale as M 3 . Considering the specific case of AO/GTO basis, the construction of the effective Hamiltonian matrix -and the corresponding derivatives-should also be considered as an additional ratelimiting step, although robust linear scaling methods are nowadays well recognized. [33][START_REF] Challacombe | Linear scaling computation of the fock matrix[END_REF][START_REF] Schwegler | Linear scaling computation of the fock matrix. ii. rigorous bounds on exchange integrals and incremental fock build[END_REF][START_REF] Ochsenfeld | Linear and sublinear scaling formation of hartree-fock-type exchange matrices[END_REF] Disregarding by now this specific feature, linear scaling can be achieved only if the two following conditions are fullfilled: (i) the theoretical framework involves the density matrix as the unique variable, that is no wavefunctions enter anymore in the formalism, (ii) perturbed density matrices must preserve some locality pattern allowing for SpMM algebra. Whereas, to the authors knowledge, analytical demonstration of the former point has not yet been proposed, raw numerical analysis have already shown that first-order perturbed density matrices for insulating systems present an approximate exponential decay of the elements [START_REF] Ochsenfeld | A reformulation of the coupled perturbed self-consistent field equations entirely within a local atomic orbital density matrix-based scheme[END_REF][START_REF] Weber | Ab initio linear scaling response theory: Electric polarizability by perturbed projection[END_REF] which apply also, to a lesser extent, to higher orders. [START_REF] Weber | Higher-order response in o (n) by perturbed projection[END_REF] Current methods dealing with condition (i) are intrinsically related to the DM solver used for the unperturbed case. Concerning the schemes, Ochsenfeld and Head-Gordon first reformulated the CPSCF equations in terms of the density matrix only [START_REF] Ochsenfeld | A reformulation of the coupled perturbed self-consistent field equations entirely within a local atomic orbital density matrix-based scheme[END_REF] (referred as CG-CPSCF by the authors) starting from the Li-Nunes-Vanderbilt (LNV) unconstrained energy functional [START_REF] Li | Density-matrix electronic-structure method with linear system-size scaling[END_REF] where the McWeeny purification polynomial [START_REF] Mcweeny | The density matrix in self-consistent field theory. i. iterative construction of the density matrix[END_REF][START_REF] Mcweeny | Some recent advances in density matrix theory[END_REF] is used as input DM. Later, Kussmann and Ochsenfeld recognized important deficiencies in this initial version which were corrected in the alternative derivation of Ref. [START_REF] Kussmann | Linear-scaling method for calculating nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals within hartree-fock and density-functional theory[END_REF][START_REF] Kussmann | A density matrix-based method for the linear-scaling calculation of dynamic second-and third-order properties at the hartree-fock and kohn-sham density functional theory levels[END_REF].

In this manuscript we present the necessary and sufficient materials for developing a one-particle density matrix solver for electronic structure perturbation theory, especially within the framework of single-determinant theory. Chapter 1 introduces the main quantum mechanical fondations necessary to approach the density as the main variable. Chapter 2 presents the current methods applied to solve for the density matrix. In Chapter 3, a comprehensive description of the density matrix perturbation theory is presented. In Chapter 4, a large set of numerical experiments are performed to compare the various schemes.

Chapter 1

Density matrices and electronic structure 1.1 Density operator and stationary condition

We shall start from an ensemble of particles within some external potential. From quantum mechanics first postulate, information on this ensemble -eg. positions and momenta-is completely specified by a mathematical object: the wavefunction Ψ, which relates the probability amplitude of finding the system in a state -usually symbolized by the ket |Ψ⟩-to its physical observation. The physical observation can only be realized through the scope of an operator Ô, for which, when applied to |Ψ⟩ and integrated over the space of the possibilities, results in the most probable value, that is, the expectation value of the observable (operator). Obviously, depending on what we want to observe, the operator is chosen accordingly. Whatever is this observable, at the end, the process is always the same, that is, bring |Ψ⟩ to the space specified by the operator, ie. | ÔΨ⟩, and integrate over that space, ⟨Ψ| ÔΨ⟩. The time evolution of the ket is governed by the time-dependent Schrödinger equation [START_REF] Schrödinger | An undulatory theory of the mechanics of atoms and molecules[END_REF]

iℏ ∂ |Ψ(t)⟩ ∂t = Ĥ(t) |Ψ(t)⟩ (1.1)
where Ĥ is the Hamiltonian operator which describes the energy of the system. For an unperturbed and closed system, where the Hamiltonian operator does not depend explicitely on time, the time dependence of the wave function can be separated assuming

|Ψ(t)⟩ = e -i Ĥt/ℏ |Ψ⟩ (1.2)
The expectation value of the time-independent Hamiltonian, ie. the energy E, is then given by

E = ⟨Ψ(t)| Ĥ|Ψ(t)⟩ ⟨Ψ(t)|Ψ(t)⟩ = ⟨Ψ| Ĥ|Ψ⟩ ⟨Ψ|Ψ⟩ (1.3)
where the second equality emphasizes that E is independent on the time, that is, the eigenstate {E, Ψ} is stationary. Given a properly normalized state,

⟨Ψ|Ψ⟩ = 1 (1.4)
evaluating E from the definition of Eq. (1.3), requires to solve an eigenvalue problem: the time-independent Schrödinger equation

Ĥ |Ψ⟩ = E |Ψ⟩ (1.5)
1.1 Density operator and stationary condition

7

Let us now introduce another object, the density operator, which from a mathematical point of view represents an (orthogonal) projection from a vector space to the same vector space. The density operator is defined according to

D := |Ψ⟩ ⟨Ψ| (1.6)
such that, from the above definition and the normalization condition (1.4), D verifies the following properties:

hermicity: D = D † idempotency: D2 = D normalization: Tr{ D} = 1 (1.7a) (1.7b) (1.7c)
For the last equality, we made use of the property of the projection operators: the trace of the projection matrix is equal to the inner product of its constitutive eigenvectors.

From here, we can search for the equation of motion of D retaining the Schrödinger picture of Eq. (1.1). This gives rise to the Liouville-von Neumann equation [START_REF] Von | Mathematical foundations of quantum mechanics[END_REF]

iℏ ∂ D(t) ∂t = Ĥ(t), D (t) (1.8)
where [•, •] denotes a commutator. Again, if we consider a conservative system, the solution of Eq. (1.8) is found to be D(t) = e -i Ĥt/ℏ D e i Ĥt/ℏ (1.9)

In that case, the stationary condition (1.5) can be recast in an operator form, according to Ĥ D = D Ĥ (1.10)

We may call it the time-independent Liouville-von Neumann equation. On multipying on the left (or on the right) by D and assuming that conditions (1.7) are respected, we obtain

D Ĥ D = D2 Ĥ ⇔ Ĥ D2 = D Ĥ D |Ψ⟩ E ⟨Ψ| = D Ĥ ⇔ Ĥ D = |Ψ⟩ E ⟨Ψ| E = Tr{ D Ĥ} ⇔ Tr{ Ĥ D} = E (1.11a) (1.11b) (1.11c)
where me made explicit the evaluation of the energy of the stationary state. Therefore, we found that the expectation value of Ĥ does not necessarily require to resolve the Schrödinger equation (1.5). An alternative route is offered by the calculation of the Density matrices and electronic structure matrix representation of the density operator. It is worth noticing that there is a one-to-one correspondence between (E, Ψ) and ( Ĥ, D) such that, for non-degenerate cases, the unique solution of the Schrödinger equation leads to a unique definition of the density operator. This remark can, in principle, be extended to the calculation of (time-dependent/independent) properties based on Rayleigh-Schrödinger perturbation theory.

If associated with |Ψ⟩ we defined an (abstract) vector space as being a separable Hilbert space of elements {|u i ⟩} ∞ i=1 , such that, ⟨u i |u j ⟩ = δ ij , the ket is expanded into this basis according to

|Ψ⟩ = i ⟨u i |Ψ⟩ |u i ⟩ (1.12)
On inserting the above definition into Eq. (1.6), the density operator transforms to

D = i,j ⟨u i |Ψ⟩ |u i ⟩ ⟨u j | ⟨Ψ|u j ⟩ = i,j |u i ⟩ ⟨u i |Ψ⟩ ⟨Ψ|u j ⟩ ⟨u j | = i,j |u i ⟩ D ij ⟨u j | with: D ij := ⟨u i | D|u j ⟩ (1.13a) (1.13b) (1.13c) (1.13d)
As a result, D can be expressed as a superposition of basis projectors. In that case, it is easy to show that the definition (1.13c) also verifies the properties (1.7), the trace of the density operator being defined by

Tr{ D} = i D ii (1.14)
If now, instead of enforcing the system to be described by a single pure state, we allow for a statistical description. On introducing a probability distribution over all the possible pure states, the ensemble (S ) density operator turns to be

DS := i p i |Ψ S i ⟩ ⟨Ψ S i |
subject to:

p i ≥ 0 i p i = 1 (1.15a) (1.15b) (1.15c)
1.2 Density matrix for fermion systems

9

where p i is the probability of the system being found in the microstate |Ψ S i ⟩. Depending on the statistical ensemble, eg. canonical (S = NVT ) or grand canonical (S = µVT ), the particle number of each |Ψ S i ⟩ may vary. Within the NVT ensemble, subject to the constraints of Eq. (1.15) we can show that DS verifies the hermiticity and the normalization conditions of Eq. (1.7), whereas idempotency is lost. This observation allows to distinguish a mixed state from a pure state. By induction, we can state that a pure state is well-defined within the microcanonical ensemble (S = NVE) at T = 0.

Density matrix for fermion systems 1.2.1 Generalities

In this work, we are mainly interested in computing the energy of an ensemble of N electrons within the external potential created by K nuclei, and latter in the manuscript, its variation(s) with respect to some external perturbation(s). Since electrons are fermions, we must insure that the ket |Ψ⟩ -as defined in some vector space-respects the antisymmetry principle. As a consequence, the space of representation of |Ψ⟩ is reduced to the anti-symmetric Hilbert space, such that, given the state vectors {|u i ⟩} ∞ i=1 allowing to define the N -particle symmetric state,

|u 1 u 2 • • • u N ) := |u 1 ⟩ |u 2 ⟩ • • • |u N ⟩ (1.16)
any of the N -particle anti-symmetric state is obtained from the following definition:

|u 1 u 2 • • • u N ⟩ := A |u 1 u 2 • • • u N ) with: A := 1 √ N ! p∈S N (-1) p P p (1.17)
where A is the anti-symmetrization operator, P is the permutation operator of two particles, and (-1) p relates the parity of the permutation. 1 It is customary in chemistry to solve the time-independent Schrödinger equation [Eq. (1.5)] -or the time-independent Liouville-von Neumann [Eq. (1.10)]-within the position-spin space. In that context, the N -particle wavefunction is expressed in an abstract basis of continuous position 10 Density matrices and electronic structure vectors, according to

Ψ(x 1 x 2 • • • x N ) = ⟨x 1 x 2 • • • x N |A |Ψ⟩ (1.18)
where {|x i ⟩ := |r i , σ i ⟩} stands for the space and spin coordinates. In this new vector space, the matrix representation of the density operator as given in Eq. (1.13c) reads

D(x 1 x 2 • • • x N , x ′ 1 x ′ 2 • • • x ′ N ) = Ψ(x 1 x 2 • • • x N )Ψ * (x ′ 1 x ′ 2 • • • x ′ N ) (1.19)
Note that by using Eq. (1.18) and the resolution of identity in a continuous basis,

dx |x⟩ ⟨x| = I (1.20)
where dx := drdσ, we can easily demonstrate that the hermiticity and imdempotency properties still hold for the density matrix defined in Eq. (1.19). Looking especially to the trace, we obtain

Tr{D(x 1 • • • x N , x ′ 1 • • • x ′ N )} = dx 1 • • • dx N D(x 1 • • • x N , x 1 • • • x N ) = dx 1 • • • dx N Ψ(x 1 • • • x N )Ψ * (x 1 • • • x N ) (1.21a) (1.21b)
On multiplying Eq. (1.19) by infinitesimal space-spin elements centered on each particle coordinate, keeping only the diagonal elements,

D(x 1 x 2 • • • x N , x 1 x 2 • • • x N ) dx 1 dx 2 • • • dx N (1.22)
we obtain the probability of an electron is in the space-spin volume element dx 1 located at x 1 with spin state s 1 , while simultaneously another electron is in dx 2 at x 2 with spin state s 2 and so on.

Reduced density matrices

Since one-and two-electron operators are necessary to fully describe the electronic Hamiltonian, the generalized expression of Eq. (1.19) describing the so-called N th-order density matrix can be reduced following the reduced density matrix theory. For the second-order reduced density matrix D 2 this gives

D 2 (x 1 x 2 , x ′ 1 x ′ 2 ) := N (N -1) 2 dx 3 • • • dx N D(x 1 x 2 • • • x N , x ′ 1 x ′ 2 • • • x N ) (1.23)
1.2 Density matrix for fermion systems 11 and for the first-order density matrix D 1 ,

D 1 (x 1 , x ′ 1 ) := N dx 2 • • • dx N D(x 1 x 2 • • • x N , x ′ 1 x 2 • • • x N ) (1.24)
It is worth to emphasize that D 2 fully determines D 1 . This is apparent by noting that

D 1 (x 1 , x ′ 1 ) = 2 N -1 dx 2 D 2 (x 1 x 2 , x ′ 1 x 2 ) (1.25)
From the above definitions, we found that the trace of the second-order density matrix leads to the number of particle pairs,

Tr{D 2 (x 1 x 2 , x ′ 2 x ′ 2 )} = dx 1 dx 2 D 2 (x 1 x 2 , x 1 x 2 ) = N (N -1) 2 (1.26)
whereas, the first-order density matrix is normalized to the number of particles

Tr{D 1 (x 1 , x ′ 1 )} = dx 1 D 1 (x 1 , x 1 ) = N (1.27)
Indeed by integrating out the spin variable, the diagonal of D 1 is recognized as the one-electron density function, -usually identified as ρ(r)-, that is, the probability of finding one electron in dr at position r assuming that the others are anywhere else, the indistinguishability of the fermions being properly accounted for in Eq. (1.27). In this work, we have reduced the scope of our investigations to closed-shell systems where the total (electron) spin momentum is zero. Thereafter, we shall integrate out the spin variables. Within the context of continuous reduced density matrices, the expectation values of any one-particle ( Ô1 ) and two-particle ( Ô2 ) operators are given by

Tr{O 1 D} = dr 1 dr ′ 1 O 1 (r 1 , r ′ 1 )D 1 (r ′ 1 , r 1 ) Tr{O 2 D} = dr 1 dr 2 dr ′ 1 dr ′ 2 O 2 (r 1 r 2 , r ′ 1 r ′ 2 )D 2 (r ′ 1 r ′ 2 , r 1 r 2 ) (1.28) (1.29)
If we assume that those operators are local, ie.

O 1 (r 1 , r ′ 1 ) = O 1 (r 1 )δ(r 1 -r ′ 1 ) O 2 (r 1 r 2 , r ′ 1 r ′ 2 ) = O 2 (r 1 r 2 )δ(r 1 -r ′ 1 )δ(r 2 -r ′ 2 ) (1.30) (1.31)
Eqs. (1.28) and (1.29) simplify to

Tr{O 1 D} r 1 =r ′ 1 = dr 1 O 1 (r 1 )D 1 (r 1 , r 1 ) Tr{O 2 D} r 1 =r ′ 1 ,r 2 =r ′ 2 = dr 1 dr 2 O 2 (r 1 r 2 )D 2 (r 1 r 2 , r 1 r 2 ) (1.32) (1.33)
The non-relativistic time-independent Hamiltonian operator for an ensemble of N fermions within the external potential (V ext ) created by K nuclei -considered fixed in positions {R A } K A=1 with charges {Z A } K A=1 -can be expressed as the sum over the electron kinetic energy operator ( T ), the electron-electron Coulomb interaction ( Ĝ), and the aforementioned external potential

Ĥ := T + Vext + Ĝ (1.34)
with the following explicit definitions2

T = - 1 2 i ∇ 2 i , Ĝ = i<j 1 |r i -r j | and Vext = i v(r i ) with: v(r i ) = - A Z A |r i -R A | (1.35)
where i runs over the electrons, (i < j) the electron pairs, and A the nuclei. By recognizing that T and Vext are one-electron operators, and Ĝ is a two-electron operator, in virtue of Eqs. (1.32) and (1.32), the electronic energy is given by

E [D 2 ] = T [D 1 ] + V ext [D 1 ] + G [D 2 ] (1.36)
according to the following definitions

T [D 1 ] := Tr{ T D} = dr 1 δ(r 1 -r ′ 1 ) - 1 2 ∇ 2 r D 1 (r 1 , r ′ 1 ) V ext [D 1 ] := Tr{ Vext D} = dr 1 v(r 1 )D 1 (r 1 , r 1 ) G [D 2 ] := Tr{ ĜD} = dr 1 dr 2 D 2 (r 1 r 2 , r 1 r 2 ) |r 1 -r 2 | (1.37) (1.38) (1.39)
where, for the kinetic energy component, we made explicit the fact that, first, the Laplacien is applied to D 1 (r 1 , r ′ 1 ), and then, integration over space coordinate is performed for r ′ = r. The square braket [•] in Eqs. (1.37)- (1.39) indicates the functional dependence of the energy contribution over the density matrix. As previously noted, since D 1 is determined by D 2 , it is necessary and sufficient to know the second-order reduced density matrix. Starting from an initial guess for D 2 without any prior knowledge of the electronic wavefunction, and solving a variational principle related to Eq. (1.36) in order to evalute the energy of an ensemble of interacting electrons is an extraordinarily difficult task that we leave to the specialists of the reduced density matrix theory (RDMT). [START_REF] Löwdin | Quantum theory of many-particle systems. i. physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction[END_REF][START_REF] Löwdin | Quantum theory of many-particle systems. ii. study of the ordinary hartree-fock approximation[END_REF][START_REF] Löwdin | Quantum theory of many-particle systems. iii. extension of the hartree-fock scheme to include degenerate systems and correlation effects[END_REF][START_REF] Mcweeny | The density matrix in many-electron quantum mechanics. i. generalized product functions. factorization and physical interpretation of the density matrices[END_REF][START_REF] Mcweeny | The density matrix in many-electron quantum mechanics. ii. separation of space and spin variables; spin coupling problems[END_REF][START_REF] John | Reconstructive approaches to one-and two-electron density matrix theory[END_REF] At this stage, it is important to note that the RDMT embrasses, at some point, the formalism of the (orbital-free) density functional theory [2,[START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Valentin | Issues and challenges in orbital-free density functional calculations[END_REF][START_REF] Xia | Can orbital-free density functional theory simulate molecules[END_REF][START_REF] Chen | Petascale orbital-free density functional theory enabled by small-box algorithms[END_REF] (DFT) and Kohn-Sham DFT (KS-DFT), [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] in the sense where, for a given V ext [D 1 ], both of them are trying to approximate T [D 1 ] and G [D 2 ] -in terms of D 1 (r 1 , r 1 ) for DFT and D 1 (r 1 , r ′ 1 ) for KS-DFT-, without requiring the support of Ψ, which is also a tedious challenge.

Interestingly for our work, addressing the difficulties mentioned aboved involved introducing a set of constraints which must be fullfiled during the search of the solution to guarantee that, at convergence, the density matrix corresponds to an acceptable anti-symmetrized wavefunction. These constraints are called the N -representability conditions [START_REF] John | Structure of fermion density matrices[END_REF][START_REF] Kummer | n-representability problem for reduced density matrices[END_REF][START_REF] Aj Coleman | Necessary conditions for n-representability of reduced density matrices[END_REF] for RDMT, in addition to the v-representability conditions, [2,[START_REF] Levy | Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem[END_REF][START_REF] Elliott | Density functionals for coulomb systems[END_REF] a specificity of the DFT.

Density matrix for a single determinant

Let us now consider a more standard Hilbert space built from a discrete set of square integrable functions {|ψ i ⟩} ∞ i=1 , such that, ψ i (r) = ⟨r|ψ i ⟩ and ⟨ψ i |ψ j ⟩ = δ ij . If we impose that the N -wavefuntion is approximated by a single anti-symmetrized product [Eq. (1.17)] -also called a Slater determinant-of a subset of these functions, that is {|ψ i ⟩} N i=1 , the first-order density matrix in the coordinate representation, reads

D 1 (r 1 , r ′ 1 ) = N i=1 ψ i (r 1 )ψ * i (r ′ 1 ) (1.40)
We can show that there exists a one-to-one mapping between D 1 and a single antisymmetrized product of the form (1.17). As a consequence of the density operator properties: (i) the trace of D 1 is equal to the number of electrons [Eq. (1.27)] and (ii) D 1 is idempotent, ie. dr ′′ D 1 (r, r ′′ )D 1 (r ′′ , r ′ ) = D 1 (r, r ′ ), constitute necessary and sufficient conditions for the first-order density matrix to correspond to a pure state approximated by a single Slater determinant. [2] In that context, the second-order density matrix can be defined in terms of D 1 according to

D 2 (r 1 r 2 , r ′ 2 r ′ 2 ) = 1 2 (D 1 (r 1 , r ′ 1 )D 1 (r 2 , r ′ 2 ) -D 1 (r 1 , r ′ 2 )D 1 (r 2 , r ′ 1 )) (1.41)
On introducing the Eqs. (1.26) and (1.41) into Eq. (1.36), we obtain the definition of the Hartree-Fock (HF) energy expressed in density matrix form:

E HF [D 1 ] := T [D 1 ] + V ext [D 1 ] + G [D 1 ] with: G [D 1 ] := J [D 1 ] + K [D 1 ] (1.42)
such that,

J [D 1 ] := + 1 2 dr 1 dr 2 D 1 (r 1 , r 1 )D 1 (r 2 , r 2 ) |r 1 -r 2 | K [D 1 ] := - 1 2 dr 1 dr 2 D 1 (r 1 , r 2 )D 1 (r 2 , r 1 ) |r 1 -r 2 | (1.43) (1.44)
where the electron-electron energy is a sum over the classical Coulomb repulsion J, and the quantum exchange energy K arising from the anti-symmetry principle of Eq. (1.17).

In textbooks it is also common to find the following condensed expression for the Hartree-Fock energy, 3

E HF [ρ(r, r ′ )] = r=r ′ dr - 1 2 ∇ 2 + v(r) ρ(r, r ′ ) + 1 2 drdr ′ ρ(r)ρ(r ′ ) |r -r ′ | - ρ(r, r ′ )ρ(r ′ , r) |r -r ′ | (1.45)

Density matrix representation in finite non-orthogonal basis

Practical calculation of the HF energy and other properties requires a closed and separable Hilbert subspace, although the accuracy of result is closely linked with the dimension of such space through the variational principle. In chemistry, it is rather natural to use local atomic orbitals (AO) to describe chemical bonds in molecules. When properly parametrized, these orbitals permit to reach a high level of accuracy with a limited set of variational parameters. Unfortunatly, these orbitals form a non-orthogonal basis for the representation of the operators, increasing the computational complexity.

3 Using the following substitutions: ρ(r) := D 1 (r, r) and ρ(r, r ′ ) := D 1 (r, r ′ ).

Density matrix for fermion systems 15

Let us introduce a set of M non-orthogonal atomic-like basis functions {|ϕ µ ⟩} M µ=1 to expand the one-electron states4 {|ψ a ⟩} M a=1 -the sp-called molecular orbitals (MO). According to Eq. (1.12), we have

|ψ a ⟩ = M µ=1 ⟨ϕ µ |ψ a ⟩ |ψ a ⟩ (1.46)
which, in coordinate representation, transforms to

⟨r|ψ a ⟩ = µ ⟨ϕ µ |ψ a ⟩ ⟨r|ϕ µ ⟩ ⇔ ψ a (r) = µ c µa ϕ µ (r) (1.47)
where the inner products, {⟨ϕ µ |ψ a ⟩} µ , are generally identified as the linear combination of atomic orbital (LCAO) coefficients c µa µ . These coefficients constitute the set of variational parameters to be optimized. By making use of Eqs. where i runs over the N occupied states, that is, in restricted Hartree-Fock theory:

N = N/2.
The set of elements {D µν } constitutes the HF one-particule density matrix in the atomic orbitals basis. In coordinate representation Eq. (1.48) transforms to

D(r, r ′ ) := ⟨r| D|r ′ ⟩ = µ,ν ϕ µ (r)D µν ϕ ν (r ′ ) (1.49)
Note in passing, for such kind of basis set we need to introduce the overlap matrix S to respect the idempotency relation, ie.

D(r, r ′ ) = dr ′′ dr ′′′ D(r, r ′′ )S(r ′′ , r ′′′ )D(r ′′′ , r ′ ) (1.50)
For a properly normalized basis, the elements of the S matrix are determined by

S µν = ⟨µ|ν⟩ = dr ϕ * µ (r -R µ )ϕ ν (r -R ν )    S µµ = 1 0 < |S µν | < 1 Tr{S} = M    (1.51) 16
Density matrices and electronic structure

In order to make the link between the density matrix formalism described above and programming, we shall recast Eq. (1.48) using a more suitable form. Assuming that we already know the LCAO coefficients for each |ψ a ⟩, we generally defined a coefficient matrix, in our case C ∈ R M ×M , using the format described below:

C ≡                   • • • • • • ψ 1 ψ N ψ N +1 ψ M occupied unoccupied (1.52)
where C collects the M coefficients (sorted in column) of the M eigenvectors (sorted in row). Consequently, the one-particle density matrix is expressed as

D = COC † (1.53)
where O is the matrix of the particle occupation numbers, where Ō is the matrix of the hole occupation numbers,

O := diag{I N , 0 N } ≡                               occupied 1 0 0 . . .
Ō := diag{0 N , I N } ≡                               occupied 0 0 0 . . . 0 0 0 1 0 unocuppied . . . 0 1 (1.56)
It is worth to emphasize that O and Ō are the matrix representations of the one-particle and one-hole density operator, respectively, in the molecular orbitals basis. Using this representation, the idempotency property of Eq. (1.50) writes:

D = DSD = COC † SCOC † subject to: C † SC = I (1.57)

Restricted Hartree-Fock energy

In restricted Hartree-Fock (RHF) theory applied to closed-shell systems, the electronic energy of Eq. (1.42) can be recast in the following form

E HF = Tr{D(2h + G)} = Tr{D(h + F )} (1.58) (1.59)
with h and G the one-electron and two-electron contributions. In Eq. (1.59) we have expressed the HF energy in terms of the Fock matrix, F := h + G (cf. Section 1.5). The matrix elements of the one-particle Hamiltonian (also called the core hamiltonian), within the AO basis, are defined according to

h µν = ⟨µ| ĥ|ν⟩ with: ĥ = - 1 2 ∇ 2 - K A=1 Z A |r -R A | (1.60a) (1.60b)
where we recognize the kinetic energy operator and the external potential created by the K nuclei, respectively. The matrix elements of G are given by:

G µν = ⟨µ| VHF |ν⟩ with: VHF = η,κ D ηκ 2 dr 2 ϕ * η (r 2 )ϕ κ (r 2 ) |r 1 -r 2 | -dr 2 ϕ * η (r 2 )P 12 ϕ κ (r 2 ) |r 1 -r 2 | (1.61a) (1.61b)
where P 12 is the 12-permutation operator5 already introduced in Eq. (1.17). The first and second term in Eq. (1.61b) are easily recognized as the Coulomb and exchange operators:

Ĵ := η,κ dr 2 ϕ * η (r 2 )D ηκ ϕ κ (r 2 ) |r 1 -r 2 | K := η,κ dr 2 ϕ * η (r 2 )D ηκ P 12 ϕ κ (r 2 ) |r 1 -r 2 | (1.62a) (1.62b) such that: VHF = 2 Ĵ - K
We emphasize that the Coulomb operator of Eq. (1.62a) can be re-written in terms of the one-particle electron density ρ(r), according to

2 Ĵ = dr 2 ρ(r 2 ) |r 1 -r 2 | with: ρ(r) = 2 η,κ D ηκ ϕ * η (r)ϕ κ (r) (1.63a) (1.63b)
If we look for expressions of G matrix elements, by inserting Eqs. (1.62a) and (1.62b) into Eq. (1.61b), we obtain: G µν = 2J µν -K µν with:

J µν = ⟨µ| Ĵ|ν⟩ = η,κ dr 1 dr 2 ϕ * µ (r 1 )ϕ * η (r 2 )D ηκ ϕ κ (r 2 )ϕ ν (r 1 ) |r 1 -r 2 | K µν = ⟨µ| K|ν⟩ = η,κ dr 1 dr 2 ϕ * µ (r 1 )ϕ * η (r 2 )D ηκ ϕ κ (r 1 )ϕ ν (r 2 ) |r 1 -r 2 | (1.64a) (1.64b)
We shall introduce some of the commonly used notations in Chemistry for the matrix elements of the Coulomb and exchange operators. They are expressed in terms of the electron repulsion integrals (ERI). Following the Chemists' notation, the formal expression of an ERI is given according to

(µν|ηκ) := dr 1 dr 2 ϕ * µ (r 1 )ϕ ν (r 1 )ϕ * η (r 2 )ϕ κ (r 2 ) |r 1 -r 2 | (1.65)
where the coordinate of the electron 1 and 2 appear side by side with respect to the vertical bar in (µν|ηκ). If we identify Eq. (1.65) to be the Coulomb integral as found in Eq. (1.64a), the exchange integral appearing in Eq. (1.64b) writes

(µκ|ην) := dr 1 dr 2 ϕ * µ (r 1 )ϕ κ (r 1 )ϕ * η (r 2 )ϕ ν (r 2 ) |r 1 -r 2 | (1.66)
As a result, the G matrix elements are defined according to

G µν = η,κ D ηκ [2(µν|ηκ) -(µκ|ην)] (1.67)
On including the one-electron Hamiltonian in Eq. (1.67), we obtain the following expression for the Fock matrix elements:

F µν = h µν + η,κ D ηκ [2(µν|ηκ) -(µκ|ην)] (1.68)
We shall briefly review other commonly found definitions of the electronic energy since it can be confusing in literature. From Eq. (1.58), the RHF energy can be alternatively expressed as

E = Tr{D(2H + G)} = 2Tr{D(H + 1 2 G)} (1.69)
By introducing the bond-order matrix, P := 2D, we may write

E = Tr{P (H + 1 2 G)} = 1 2 Tr{P (H + F )} with: G µν = η,κ P ηκ (µν|ηκ) - 1 2 (µκ|ην) (1.70a) (1.70b)
If we briefly review this Section, we observed that the construction for the Fock matrix can be a serious bottleneck for large scale calculations. Since the AOs extend over the whole molecule, the amount of information in each AO is proportional to N . [START_REF] Bowler | Recent progress in linear scaling ab initio electronic structure techniques[END_REF] Based on a naive analysis of Eqs. (1.60a), (1.64a) and (1.64b), we found that building the one-core Hamiltonian involved integrating the product of two AOs scales as O(N 2 ).

Density matrices and electronic structure

While with the product of four AOs, the Coulomb and the exchange matrices scale as O(N 4 ). Therefore, the cost of the Fock matrix construction is basically O(N 4 ). For large molecular systems, relying on some distance criteria, the scaling can be reduced to reach asymptotically N 2 . In order to calculate the Coulomb and exchange integrals in full linear scaling, specific numerical techniques for GTOs have been developped, [START_REF] Ochsenfeld | Linear-scaling methods in quantum chemistry[END_REF]8] such as the continuous fast multipole method [33] for the Coulomb integrals, and the LinK method [START_REF] Ochsenfeld | Linear and sublinear scaling formation of hartree-fock-type exchange matrices[END_REF][START_REF] Ochsenfeld | Linear scaling exchange gradients for hartree-fock and hybrid density functional theory[END_REF] for the exchange integrals.

Pariser-Parr-Pople method

The entire implementation and all the applications performed during this thesis were based on a semi-empirical method derived for the calculation of energetics and properties of πconjugated systems. In regard to the all-electron HF approaches based on non-orthogonal extended local basis sets, this choice permits to focus our efforts mainly on algorithmic developments related to the density matrix solvers, the computational ressources used for evaluating the matrices being negligeable in that case. It has also the merit to overcome the intricacies encountered when trying to modify some of the routines found in standard quantum chemistry packages. Using such kind of simplified HF model, we were able to performed a fast and fair comparison between various methods, while making sure that coding was optimized for all of them. For these reasons, we have considered the most simple semi-empirical HF model, the Pariser-Parr-Pople method [START_REF] Pariser | A semi-empirical theory of the electronic spectra and electronic structure of complex unsaturated molecules. i[END_REF][START_REF] Pariser | A semi-empirical theory of the electronic spectra and electronic structure of complex unsaturated molecules. ii[END_REF][START_REF] Ja Pople | Electron interaction in unsaturated hydrocarbons[END_REF] (PPP), which was originally developed for treating conjugated hydrocarbons.

Zero-differential-overlap approximation

For atomic-like basis functions centered on atoms, ie. which are explicilty dependent on the nucleus coordinates, such as Slater-type orbital (STO) or Gaussian-type orbitals (GTO) -commonly used in quantum chemistry-, the ERI of Eq. (1.65) explicitly writes as:

dr 1 dr 2 ϕ * µ (r 1 -R µ )ϕ ν (r 1 -R ν )ϕ * η (r 2 -R η )ϕ κ (r 2 -R κ ) |r 1 -r 2 | (1.71)
where

{R µ , R ν , R η , R κ }
|ϕ µ (r 1 -R µ )| 2 |ϕ η (r 2 -R η )| 2 |r 1 -r 2 | (1.73)
Note that the Kronecker functions in Eq. (1.72), implies that the basis set is orthonormal, such that:

⟨µ|ν⟩    1 if µ = ν 0 if µ ̸ = ν (1.74)
On introducing Eq. (1.72) into (1.68), we obtain:

F µν = h µν + η,κ D ηκ [2(µµ|ηη)δ µν δ ηκ -(µµ|ηη)δ µκ δ ην ] (1.75)
which, without loss of generality, simplifies to

F µν = h µν + 2 η D ηη (µµ|ηη)δ µν - η,κ D νν (µµ|ηη)δ µκ δ ην for: µ = ν F µµ = h µµ + 2 ν D νν (µν|µν) -D µµ (µµ|µµ) for: µ ̸ = ν F µν = h µν -D νµ (µν|µν) (1.76a) (1.76b) (1.76c)
In literature, we identify: Γ µµ := (µµ|µµ)

Γ µν := (µν|µν)
which are the one-center integral corresponding to an energy constant, and the twoelectron repulsion integral which is parametrized with respect to a set of internal geometric parameters. As a result, the ZDO approach greatly simplifies the problem at the cost of a parametrization, which reduces the ab initio character of the Hartree-Fock method.

Pariser-Parr-Pople model parameterization

Several parameterizations of the PPP model can be found in literature. [START_REF] Mataga | Electronic structure and spectra of nitrogen heterocycles[END_REF][START_REF] Ohno | Some remarks on the pariser-parr-pople method[END_REF][START_REF] Schulten | Correlation effects in the spectra of polyenes[END_REF][START_REF] Zhang | Excitation energy calculation of conjugated hydrocarbons: A new pariser-parr-pople model parameterization approaching caspt2 accuracy[END_REF] In this work we haved used the Ohno's parameterization [START_REF] Ohno | Some remarks on the pariser-parr-pople method[END_REF] using standard parameters [START_REF] Sony | A general purpose fortran 90 electronic structure program for conjugated systems using pariser-parr-pople model[END_REF][START_REF] Kondayya | A fortran 90 hartree-fock program for one-dimensional periodic π-conjugated systems using pariser-parr-pople model[END_REF] collected in Table {1 distance of the first nearest neighbour(s). Then, the electron-nuclei interaction is added to the TB part in order to give the final one-electron Hamiltonian matrix elements h µν .

In the Table {1.1}, the electron-nuclei interaction is the second term in h µµ . Finally, the two-electron contribution is constructed with the density matrix and the elements Γ µµ and Γ µν , which added to the core Hamiltonian, leads to the Fock matrix elements of Eq. (1.76).

Minimization of the Hartree-Fock energy

L HF [D] := 2Tr{D(h + 1 2 G(D))} -2 Tr{Λ(D 2 -D)} + µ(Tr{D} -N ) (1.83)
where the matrix of Lagrange multipliers Λ and the scalar µ have been introduced to constrained idempotency and trace conservation, respectively. From here, we shall search for minimizing such functional with respect to D. This yields to solve:

∇L HF [D] = 2(h + G(D) -ΛD -DΛ + Λ -µI) = 0 (1.84)
where we have used the following functional derivative properties,

8 ∇Tr{DA} = A † ∇Tr{D 2 A} = (DA + AD) † ∇Tr{DG(D)} = 2G(D) †
In the last statement, we used Tr{XG(Y )} = Tr{Y G(X)} [START_REF] Liu | The extremum method in quantum chemistry: Part i. extremum of tr (btpc) and the trace algebra[END_REF][START_REF] Liu | The extremum method in quantum chemistry: Part ii. introduction to some applications[END_REF] for X = Y = D. On recalling that all the operators are Hermitian -more specifically in this work all the matrices are symmetric-, the following working equation is found

F (D) -µI = ΛD + DΛ -Λ (1.85)
where, F = h + G, is the Fock matrix already introduced in Section 1.3. It should be outlined that despite the appealing form of this equation, to our knowledge, there is only one paper dealing with its resolution. [START_REF] Adhikari | Augmented lagrangian method for order-n electronic structure[END_REF] To demonstrate that solving Eq. (1.85) leads to an unique solution decribing a pure state within the NVE ensemble, we may try to recover the famous Roothaan-Hall equation [START_REF] Carel | New developments in molecular orbital theory[END_REF][START_REF] George | The molecular orbital theory of chemical valency. viii. a method of calculating ionization potentials[END_REF] widely used in quantum chemistry. By taking the commutator of Eq. (1.85) with respect to D, we obtain:

[D, F -µI] = DΛD + D 2 Λ -DΛ -ΛD 2 -DΛD + DΛ (1.86)
If the density matrix is exactly idempotent, the above equation reduced to:

F D = DF (1.87)
which is the single-determinant time-independent Liouville-von Neumann equation in matrix form. For instance, multiplying on the right by the coefficient matrix C and using the definition of Eq. (1.53), we have

F DC = DF C F COC † C = COC † F C F CO = COE (1.88) (1.89) (1.90)
where, for an orthogonal basis set, E is a diagonal matrix containing the M eigenvalues of F . 9 The presence of occupation number matrix O indicates that Eq. (1.90) gives access only to the eigenvalues of occupied states. Based on symmetry considerations, it is easily proved that Eq. (1.87) holds as well for the one-hole density matrix [cf. Eq. (1. [START_REF] Xia | Can orbital-free density functional theory simulate molecules[END_REF]]. As a result, the eigenvalues of the unocuppied states can be obtained by solving

F C Ō = C ŌE (1.91)
On assembling Eqs. (1.90) and (1.91), we obtain the condensed matrix form

F (D)C -CE = 0 subject to: C † C = I (1.92a) (1.92b)
which will be referred as the Roothaan-Hall equation. 10 On multiplying Eq. (1.92a) from the right by C † and using Eqs. (1.90) and (1.91), the Fock matrix reads

F = COEC † + C ŌEC † (1.93)
Hence, the spectrum of the Fock matrix can be resolved according to Eq. (1.13c), that is with i and j running over the energy-weighted projectors for the occupied and unoccupied subspace, respectively.

F = i ϵ i D i +

The self-consistent field procedure

Since the eigenvalue problem of Eqs. (1.92) is non-linear with respect to the density matrix, the Roothaan-Hall equations are resolved iteratively using a two-step approach: (i) solve a linear eigenvalue problem for a fixed Fock matrix, (ii) update the new Fock matrix from the previous solutions. Iterations are repeated until some convergence criteria is met. This procedure, called self-consistent field (SCF), is illustrated in Figure {1.1}. The update of the Fock matrix at iteration n + 1, from the Fock matrix at iteration n,

D 0 F [D n ] F [D n ]
Diag.

---→ C n+1 n = n + 1 The intial guess is one of the key steps of the SCF procedure which may have a strong impact on the convergence rate. A poor or wrong initial guess can slow down the convergence, or even worse, to a divergence. The former is generally related to oscillations when approaching the final state so that it can not be reached with a reasonable number of iterations. The latter indicates that the initial guess has no physical significance or is too far away from the expected solution. There are different ways to define the starting guess. In this work we merely start from the solutions of a tight-binding calculation, as decribed in Section 1.4.2. The fact still remains that, even a good initial guess does not prevent convergence instabilities. For that reason, numerous suggestions [START_REF] Konstantin N Kudin | A black-box selfconsistent field convergence algorithm: One step closer[END_REF][START_REF] Hu | Accelerating self-consistent field convergence with the augmented roothaan-hall energy function[END_REF][START_REF] Garza | Comparison of self-consistent field convergence acceleration techniques[END_REF][START_REF] Cancès | Can we outperform the diis approach for electronic structure calculations[END_REF][START_REF] Konstantin | Converging self-consistent field equations in quantum chemistry-recent achievements and remaining challenges[END_REF] have been made to solve these issues. Many of those were combined into hybrid methods 1.6 The self-consistent field procedure 27 attempting to overcome the weakness of the standalone models [START_REF] Hu | Accelerating self-consistent field convergence with the augmented roothaan-hall energy function[END_REF][START_REF] Garza | Comparison of self-consistent field convergence acceleration techniques[END_REF]. For our purpose here, we recall two schemes which consist in optimizing the Fock matrix construction in order to stabilize and accelerate the SCF convergence. The constant damping algorithm [START_REF] Dr Hartree | The calculation of atomic structures[END_REF] (CDA) is a simple method based on a linear interpolation between the past and futur events, that is, the Fock matrix at iteration n is defined according to Fn :

D n+1 = C n+1 OC † n+1 D n+1 vs. D n D ∞ No Yes SCF cycle

Constant damping algorithm

D 0 F [D n ] F [D n ] = F [D n-1 ] + λ (F [D n ] -F [D n-1 ]) F [D n ] Diag. ---→ C n+1 n = n + 1 D n+1 = C n+1 OC † n+1 D n+1 vs. D n D ∞ No Yes SCF cycle
= λF n + (1 -λ)F n-1 , (1.95)
where λ, referring to a damping factor, is a constant chosen freely in the interval [0, 1]. The damping step is outlined by the green chart in Figure {1.2} which can be compared to algorithm of Figure {1.1}. The main drawback with the CDA, is that the convergence rate is now fixed by the value of λ. In other words, having chosen heuristically λ to initiate the SCF processus, if for some reasons convergence problem persists, one has to stop and restart the processus with another damping parameter. It is true that some approaches were proposed to dynamically optimize the damping factor during the SCF [START_REF] Michael | A dynamical damping scheme for converging molecular scf calculations[END_REF][START_REF] Cancès | Can we outperform the diis approach for electronic structure calculations[END_REF]. It remains that using a constant or dynamic approach, the CDA is not always successful for solving convergence issues.

Direct inversion of the iterative subspace extrapolation

Another more general technique to speed up and stabilize the SCF convergence is the direct inversion of the iterative subspace (DIIS) extrapolation proposed by Pulay [START_REF] Pulay | Convergence acceleration of iterative sequences. the case of scf iteration[END_REF][START_REF] Pulay | Improved scf convergence acceleration[END_REF][START_REF] Tracy | Direct inversion in the iterative subspace (diis) optimization of open-shell, excited-state, and small multiconfiguration scf wave functions[END_REF]. This method and several improvements [START_REF] Hu | Accelerating self-consistent field convergence with the augmented roothaan-hall energy function[END_REF][START_REF] Garza | Comparison of self-consistent field convergence acceleration techniques[END_REF] have shown to be very efficient [START_REF] Hu | Accelerating self-consistent field convergence with the augmented roothaan-hall energy function[END_REF][START_REF] Garza | Comparison of self-consistent field convergence acceleration techniques[END_REF] and it constitutes one of our ingredient for efficient calculation of the response properties presented in Chapter 3. The idea of the DIIS is to extrapolate the Fock matrix at iteration n, from a linear combination of Fock matrices taken in the history of the SCF procedure according to

Fn := n i=n-m c i F i (1.96)
where m is the size of the set of historical Fock matrices and n is the iteration from which the DIIS is switched on. The coefficients of the linear combination are determined using a set of Pulay's error vectors,

{e i }, corresponding to {F i , D i }, according to e i = [F i , D i ].
The DIIS approach assumes that this linear combination is a good approximation of the

D 0 F [D n ] e n min n i,j=n-m c i c j (e i • e j ) , n i=n-m c i = 1 F [D n ] = n i=n-m c i F [D i ] F [D n ] Diag. ---→ C n+1 n = n + 1 D n+1 = C n+1 OC † n+1 D n+1 vs. D n D ∞ No Yes SCF cycle
c i (F ∞ + e i ) = F ∞ n i=n-m c i + n i=n-m c i e i (1.97)
From this equation, in order to have Fn = F ∞ , we need to minimize the norm of the second term while requiring the sum in the first term to be normalized, which gives

11 min    ∥ n i=n-m c i e i ∥ 2 n i=n-m c i = 1    (1.98)
where the coefficients {c i } are assumed real, ie. c * i = c i . As for the minimization of the Hartree-Fock energy developped in Section 1.5, we can use the Lagrange multiplier technique to minimize the expression of Eq. (1.98). This yields to

L DIIS := n i,j=n-m c i c j B ij -λ   n i=n-m c i -1   (1.99)
where λ is the Lagrange multiplier, and B ij = (e i • e j ) . The first derivative of L DIIS with respect to the coefficient c l , gives

∂L DIIS ∂c l = i,j ∂c i ∂c l c j B ij + i,j c i ∂c j ∂c l B ij -λ i ∂c i ∂c l = i,j δ il c j B ij + i,j c i δ jl B ij -λ i δ il = j c j B lj + i c i B il -λ (1.100)
where, for the last step, we used the following properties

∂c i ∂c j = δ ij , i δ ij = 1, i c i δ ij = c j (1.101)
Density matrices and electronic structure Substituting j by i in the first sum of Eq. (1.100), and using the fact that the matrix B of elements {B ij } is symmetric, Eq. (1.100) simplifies to

∂L DIIS ∂c l = 2 i c i B il -λ (1.102)
As a result, solutions of eq. (1.98) translates to the minimization of the Lagragian L DIIS such that

i c i B il -λ = 0 (1.103)
The coefficients {c i } are finally given by the Eq. (1.103) which corresponds to a system of (m + 1) linear equations [START_REF] Berrut | Barycentric lagrange interpolation[END_REF][START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF][START_REF] Meijering | A chronology of interpolation: from ancient astronomy to modern signal and image processing[END_REF][START_REF] Quarteroni | Scientific computing with matlab[END_REF]. Once the coefficients are found, the Fock matrix is updated following the linear expansion of Eq. (1.96).

The DIIS steps are defined by the orange chart in the diagram of Figure {1.3}. In the particular case of m = 2, the DIIS corresponds to the damping scheme of Eq. (1.95) but without any restriction on λ, which means that the DIIS is a dynamical extension of the CDA where the coefficients are optimized on-the-fly. As a concrete example, in Figure {1.4} are presented the SCF convergence profiles obtained for a non-optimized procedure (the simple SCF), the CDA interpolation using two different values of the damping parameter, and the DIIS extrapolation. It is quite clear that the convergence with the DIIS is the most efficient.

Chapter 2

Density matrix purifications and minimizations

The orthodox resolution of the SCF equations [cf. Section 1.6] as obtained from a single Slater determinant is based on the diagonalization of the Fock matrix at each iteration. The diagonalization step is well known for being an expensive task, which becomes rapidly the limiting step for large scale calculations. For this reason, it has been suggested that one can solve the SCF equations relying only on the one-particle density matrix. [START_REF] Mcweeny | The density matrix in self-consistent field theory. i. iterative construction of the density matrix[END_REF] Even if, in their native forms, density matrix (DM) methods also present an asymptotical cubic scaling, they constitute the first ingredient towards linear scaling regime. Density matrix solvers can be classified following the physical motivations they originated from. These solvers are: (i) the iterative density matrix functional minimizations [START_REF] Li | Density-matrix electronic-structure method with linear system-size scaling[END_REF][START_REF] Andrew D Daniels | Semiempirical methods with conjugate-gradient density-matrix search to replace diagonalization for molecular-systems containing thousands of atoms[END_REF][START_REF] John | Linear scaling conjugate gradient density matrix search as an alternative to diagonalization for first principles electronic structure calculations[END_REF][START_REF] Andrew | What is the best alternative to diagonalization of the hamiltonian in large scale semiempirical calculations[END_REF][START_REF] Bowler | Density matrices in o (n) electronic structure calculations: theory and applications[END_REF][START_REF] Challacombe | A simplified density matrix minimization for linear scaling self-consistent field theory[END_REF][START_REF] Helgaker | Direct optimization of the ao density matrix in hartree-fock and kohn-sham theories[END_REF][START_REF] Larsen | Direct optimization of the atomic-orbital density matrix using the conjugate-gradient method with a multilevel preconditioner[END_REF] where for one-determinant SCF theories, the HF or KS energy functional is minimized with respect to an auxiliary density matrix used in place of the conventional fixed DM built from the eigenvectors, (ii) the recursive density matrix polynomial expansion where the Fermi-Dirac ground state DM at the zero electronic temperature limit is obtained by a recursive application of projection polynomials -also referred to as purifications. [START_REF] Goedecker | Efficient linear scaling algorithm for tight-binding molecular dynamics[END_REF][START_REF] Baer | Chebyshev expansion methods for electronic structure calculations on large molecular systems[END_REF][START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF][START_REF] Bates | Comparison of conjugate gradient density matrix search and chebyshev expansion methods for avoiding diagonalization in large-scale electronic structure calculations[END_REF][START_REF] Liang | Improved fermi operator expansion methods for fast electronic structure calculations[END_REF][START_REF] Pino | Purification of the first-order density matrix using steepest descent and newton-raphson methods[END_REF][START_REF] David | Towards idempotent reduced density matrices via particle-hole duality: Mcweeny's purification and beyond[END_REF][START_REF] Németh | Linear scaling density matrix search based on sign matrices[END_REF][START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF][START_REF] Anders Mn Niklasson | Implicit purification for temperature-dependent density matrices[END_REF][START_REF] Anders Mn Niklasson | Trace resetting density matrix purification in o (n) self-consistent-field theory[END_REF][START_REF] Kohalmi | Idempotency-conserving iteration scheme for the one-electron density matrix[END_REF][START_REF] Sałek | Linearscaling implementation of molecular electronic self-consistent field theory[END_REF]

Density matrix minimization principle

In continuation of Section 1.5, where we have shown that the minimization of the Hartree-Fock (HF) energy functional can be expressed in terms of the one-particle density matrix only, we can either try to minimize the HF Lagrangian of Eq. (1.83) by releasing some of the constraints, or by introducing another objective functional to minimize.

Idempotency error functional minimization

In order to obtain an exactly idempotent density matrix from a roughly idempotent initial guess, McWeeny has proposed to minimize the sum of the squares of the idempotency errors, that is, ∥D 2 -D∥ 2 , using a steepest gradient descent method. [START_REF] Mcweeny | The density matrix in self-consistent field theory. i. iterative construction of the density matrix[END_REF][START_REF] Mcweeny | Some recent advances in density matrix theory[END_REF] This is fully equivalent to minimize the following functional:

Ω McW := Tr{(D 2 -D) 2 } (2.1)
Using the trace algebra summarized in Section 1.5, the gradient for this functional is given according to

∇Ω McW = 2(2D 3 -3D 2 + D) (2.2)

Density matrix minimization principle 35

The optimal step length γ for the line search of the steepest descent can be derived from the Cauchy relation:

γ := min γ Tr{(D 2 γ -D γ ) 2 } (2.3)
where

D γ := D -γ∇Ω McW (2.4)
Working on Eq. ( 2.3), to the second order in γ, the optimum value is found to be

γ = Tr{(D 2 -D) 2 } Tr{(D 2 -D) 2 (2D -I) 2 } (2.5)
On substituting D := D ′ + δ, where D ′ is trully idempotent, [START_REF] Mcweeny | The density matrix in self-consistent field theory. i. iterative construction of the density matrix[END_REF] and expanding Eq. ( 2.5), it can be easily shown that γ ≃ 1. As a result, for a fixed step length gradient descent, Eq. (2.4) reduces to

D = 3D 2 -2D 3 (2.6)
It is drawn from Eq. (2.6) that the fixed step gradient descent gives rise to an alternative approach to obtain an idempotent DM relying on the following recursive formula:

D n+1 = 3D 2 n -2D 3 n (2.7)
where n is the iteration index. This relation is the so-called McWeeny purification. We note that, in line with Section 1.1, the term purification clearly indicates that repeated application of the polynomial (2.7) to a mixed state -the initial guess (vide infra) for the density matrix-transforms it into a pure state: the idempotent one-particle density matrix. It is worth emphasizing that solution to the minimization problem of Eq. (2.1) is not restricted to the (fixed step) steepest descent. One can also consider other gradient descent based algorithms such as the conjugate gradient (CG) or Newton-Raphson method, each of them coming with their own pros and cons.[107]

Energy functional minimization

Another way to find the density matrix is to minimize an energy functional using a conjugate gradient routine. [START_REF] Li | Density-matrix electronic-structure method with linear system-size scaling[END_REF][START_REF] Andrew D Daniels | Semiempirical methods with conjugate-gradient density-matrix search to replace diagonalization for molecular-systems containing thousands of atoms[END_REF][START_REF] John | Linear scaling conjugate gradient density matrix search as an alternative to diagonalization for first principles electronic structure calculations[END_REF][START_REF] Andrew | What is the best alternative to diagonalization of the hamiltonian in large scale semiempirical calculations[END_REF][START_REF] Bowler | Density matrices in o (n) electronic structure calculations: theory and applications[END_REF][START_REF] Challacombe | A simplified density matrix minimization for linear scaling self-consistent field theory[END_REF] Within the tight-binding (TB) framework, Li, Nunes and Vanderbilt (LNV) have proposed to minimize the grand potential functional [START_REF] Li | Density-matrix electronic-structure method with linear system-size scaling[END_REF] at the zero temperature limit as defined below,

Ω µ := E [D] -µN (2.8)
where µ is the chemical potential and E[D] stands for one-eletron energy functional of the one-particle density matrix. Instead of working on D directly, which would have led to consider a set of Lagrange multipliers in Eq. (2.8) -cf. Eq. (1.82) of Section 1.5 and discussion therein-, LNV have considered the auxiliary DM of Eq. (2.7) allowing to introduce variational degrees of freedom within the grand potential functional of Eq. (2.8).

Considering the Fock matrix as input, the LNV energy functional reads:

Ω LNV := Tr{F (3D 2 -2D 3 )} -µTr{3D 2 -2D 3 }
with:

∇Ω LNV = 3(DF ′ + F ′ D) -2(D 2 F ′ + DF ′ D + F ′ D 2 ) (2.9) (2.10)
and F ′ := F -µI. Latter, Xu and Scuseria [START_REF] Hui | An o (n) tight-binding study of carbon clusters up to c 8640: the geometrical shape of the giant icosahedral fullerenes[END_REF] (XS) have proposed a slight modification of the LNV functional minimization by a damping method based on updating the chemical potential value between consecutive conjugate gradient iterations. They reported an improvement in the convergence of the CG minimization.

Unfortunately, since Ω µ (or Ω LNV ) is only well-defined within the µVT ensemble, unconstrained minimization of Eq. (2.9) is not expected to yield the correct number of particles unless the chemical potential is known exactly. This poses sever problems for unsymmetric cases, ie. when µ is not in (or close to) the middle of the eigenvalue spectrum, in other words, when the one-electron one-orbital picture is abandoned (vide infra). To correct this drawback, a more general approach was introduced by Millam and Scuseria [START_REF] John | Linear scaling conjugate gradient density matrix search as an alternative to diagonalization for first principles electronic structure calculations[END_REF] (MS), where the update of the chemical potential is constrained via the trace of the gradient. The advantage for this functional is to explicitly calculate the chemical potential during the conjugate gradient iterations so that the electron number is preserved. Nevertheless, this scheme implies that the chemical potential must be zero at convergence, [START_REF] John | Linear scaling conjugate gradient density matrix search as an alternative to diagonalization for first principles electronic structure calculations[END_REF] which clearly restricts its domain of applicability unless modifications within the working equations are derived.

The other major issue of LNV functional minimization is that the idempotency of the density matrix is not guaranteed. For this reason, the density matrix has to be purified by the McWeeny polynomial of Eq. (2.7) outside the conjugate gradient. [START_REF] Bowler | Density matrices in o (n) electronic structure calculations: theory and applications[END_REF] Extensive analysis of the LNV shortcomings and solutions are given in Refs. [START_REF] Challacombe | A simplified density matrix minimization for linear scaling self-consistent field theory[END_REF][START_REF] Bowler | Density matrices in o (n) electronic structure calculations: theory and applications[END_REF]116]. Algorithms related to the XS, MS and the original LNV are given in Appendix B, along with the CG routine used in this work. [START_REF] Charles | Global convergence properties of conjugate gradient methods for optimization[END_REF]

Density matrix polynomial expansion

The other class of density matrix solver is based on the Fermi-Dirac (FD) operator expansion. [START_REF] Anders Mn Niklasson | Density matrix methods in linear scaling electronic structure theory[END_REF] Within the µVT (or NVT ) ensemble [cf. Section 1.1 and Eq. (1.15)] at non-zero electronic temperature, the one-particle density operator for a single antisymmetrized product of one-electron function [cf. Section 1.2.3 and Eq. (1.40)] is given by: D

= i η i |ψ i ⟩ ⟨ψ i | (2.11)
where {ψ i } are the set of molecular orbitals (MO), such that, For the sake of demonstration, we shall assume that our molecular system can be equated with an ensemble of weakly interacting particles at the thermodynamic equilibrium and obeying the Fermi-Dirac (FD) statistic. [START_REF] Paul | On the theory of quantum mechanics[END_REF] As a result, for a given fermion temperature T and chemical potential µ, we may associate the occupation numbers {η i } of Eq. (2.11) with the occupation probabilities of the single-particle energy states {ϵ i } following:

⟨ψ i |ψ j ⟩ = δ ij [cf. Sec- tion 1.2.
η µ,T (ϵ i ) = 1 1 + e β(ϵ i -µ) , with β = 1 k B T (2.

Density matrix purifications and minimizations

Note that the one-particle density matrix of Eq. (2.11) with the occupation probabilities as defined in Eq. (2.12) must respect the N -representability conditions:

Tr{ D} = i η i = N, and η i ∈ [0, 1] (2.13)
Within the MO basis, the matrix representation of Eq. (2.11) is readily recognized as the matrix of the particle occupation numbers [cf. Eq. (1.54)] which, at non-zero temperature, reads:

O µ,T = diag{η 1 η 2 • • • η M },
where M is the size of the basis set. Therefore, in the atomic orbitals (AO) representation, the DM is expressed as:

D µ,T (F ) = 1 + e β(F -µI) -1 (2.14)
The equation above indicates that there exists a one-to-one non-linear correspondence between the Fock and the density matrix. In other terms, at a given temperature, the statistical density matrix is determined by µ and F . By expanding the right-hand-side of Eq. (2.14) using appropriate polynomials, one can expect to obtain the N -representable one-particle density matrix, provided that conditions (2.13) are fulfilled. This constitutes the framework of the density matrix polynomial expansion (DMPE) theory. [START_REF] Anders Mn Niklasson | Density matrix methods in linear scaling electronic structure theory[END_REF] Since the pioneering works of Goedecker and Colombo, [START_REF] Goedecker | Efficient linear scaling algorithm for tight-binding molecular dynamics[END_REF] several variations of the DMPE have been proposed. These variations can be differentiated by: (i) the polynomials used for the expansion, and (ii) the statistical ensemble chosen for describing the system. In any case, solving Eq. (2.14), implies to proceed by iteration. Within the µVT or NVT canonical ensemble, the statistical mixture of one-electron states described in Eq. (2.14) can be purified following three different ways, depending on the variable we choose to operate on.

• In µVT , for fixed chemical potential and temperature: the FD distribution of Eq. (2.14) is expanded in terms of Chebyshev polynomials. [START_REF] Goedecker | Tight-binding electronic-structure calculations and tight-binding molecular dynamics with localized orbitals[END_REF][START_REF] Baer | Chebyshev expansion methods for electronic structure calculations on large molecular systems[END_REF]. Despite the good performance of the approach, [START_REF] Liang | Improved fermi operator expansion methods for fast electronic structure calculations[END_REF][START_REF] Liang | Fast methods for resumming matrix polynomials and chebyshev matrix polynomials[END_REF] its application requires a precise knowledge of the chemical potential, that is for isolated molecular system, knowledge of interior eigenvalues. Linear scaling algorithms addressing such task along with some improvements were proposed, [START_REF] Emanuel | Determination of the chemical potential and homo/lumo orbitals in density purification methods[END_REF][START_REF] Emanuel | Computation of interior eigenvalues in electronic structure calculations facilitated by density matrix purification[END_REF][START_REF] Emanuel | Interior eigenvalues from density matrix expansions in quantum mechanical molecular dynamics[END_REF] but one may wonder if this step is really necessary.

• In NVT , for a fixed number of particles N : we might try to cool down the system towards the zero temperature limit, as depicted on Figure {2.1}. This was first proposed by Daw in his seminal paper [START_REF] Murray S Daw | Model for energetics of solids based on the density matrix[END_REF] where he demonstrated that the McWeeny purification of Eq. (2.7) is related to a "temperature-driven" density matrix equation of motion. The notion of statistical ensembles and the possibility of solving Eq. (2.14) in the NVT were rationalized by Palser and Manolopoulos [START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF] (PM) -without establishing direct relationships with Daw's proposal. In the same work, PM introduced a way of enforcing the McWeeny polynomial to preserve the N -representability conditions throughout the purification process. Later in the manuscript, the PM approach will be referred to as canonical purification (CP).

• In µ(N )V T , at the zero temperature: Niklasson [START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF] proposed to approach the FD distribution at zero temperature, ie. the Heaviside step function, [START_REF] Ernst | Heaviside's Operational Calculus as Applied to Engineering and Physics: as applied to engineering and phsyics[END_REF] by varying the number of occupied states around the exact N , that is, adjusting the polynomial dynamically during the recursion without enforcing requirements of Eq. (2.13), such that the N -representable ground state DM is obtained only at convergence. As a result, the Niklasson's method implicitly assumes that the system is coupled to a bath of particles, which are added or withdrawn with respect to the target value. This family of polynomials constitutes the basis for trace-correcting (TC) purification.

It should be emphasized that the purification methods mentioned in the last two points can be easily adapted to the grand canonical ensemble. [START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF][START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF] Calculations of the ground state density matrix using the CP and TC methods are based on the recursive application of projection polynomials {P n } to evaluate the step function, Θ(µI -F ), centered at the (unknown) chemical potential. This can be formally written as follows:

Θ(µI -F ) = lim n→∞ P n ( P n-1 (...P 2 ( P 1 ( D 0 (F ; µ) ) )...) ) with D 0 = α 1 I -α 2 (µI -F ) (2.15) (2.16)
The density matrix purification is initialized by performing the linear transformation of Eq. (2.16) where {α 1 , α 2 } are parameters judiciously chosen to: (i) map the eigenvalues of the Fock matrix into the [0, 1] interval, and (ii) depending on the purification method, to verify: Tr{D 0 } = N . An example of purification process for a (2 × 2) mixed state is presented in Figure {2.2}. At this stage, the problem of the chemical potential remains a serious bottleneck since, for a given chemical potential or number of occupied states, fixed inflection point polynomials, eg. the McWeeny polynomial of Eq. (2.7), yield the N -representable ground state density matrix, if and only if, the inflection point is located at the middle of the eigenvalue spectrum as depicted in Figure {2.1}. As a result, there exist two possibilities: (i) determine the value of µ to shift the eigenvalues of the initial

(x) = 3x 2 -2x 3 . The initial (2 × 2) mixed state, with occupation numbers {η 1 , η 2 } n=0 = {0.3, 0.7}, is purified to eventually reach {η 1 , η 2 } n=4 = {0, 1}.
mixed state of Eq. (2.16) towards the left and the right of µ, respectively, and purify according to Eq. (2.15), or (ii) define an approximate guess without the support of µ and purify using flexible inflection point polynomials. This work is dealing with the second approach.

As discussed by Niklasson in Ref. [START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF] performances of the purification methods, that is the number of iterations, depends upon the location of µ in the [ϵ min , ϵ max ] interval, where ϵ min and ϵ max are the lower and upper bounds of the eigenvalue spectrum, or equivalently, on the value of the filling factor θ = N/M , where M is the number of available states. Typical values of θ are about 1/2 when dealing with one-electron one-orbital many-electron systems, and around 1/20 for calculations based on extended basis set where, for instance, there are 10 basis functions per electron. The influence of the filling factor over the performances of the purification methods can be qualitatively understood from the fact that the preconditioning of Eq. (2.16) leads to a clustered set of eigenvalues around θ, the range of this cluster being inversely proportional to the gap of the system. As a result, for extreme values of θ -let us say θ < 0.1 (or equivalently for θ > 0.9)-where the initial DM eigenvalues are located around 0.1, the polynomials must be flexible enough to send a few of the eigenvalues towards the upper bound (η = 1) of the DM eigenspectrum, whereas all the others must be kept around the lower bound (η = 0), and purified accordingly. In the next section, we shall present the most popular purification polynomials.

Canonical purification

The Palser and Manolopoulos canonical purification [START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF] (PMCP) is based on the introduction of a flexible inflection point within the McWeeny polynomial of Eq. (2.7) that allows to address the issues mentioned above, and moreover, to preserve the Nrepresentability properties of the initial guess throughout the recursive process. The density matrix is purified according to the following algorithm:

D n+1 =        - 1 1 -c n D 3 n + 1 + c n 1 -c n D 2 n + 1 -2c n 1 -c n D n if c n ≤ 1 2 - 1 c n D 3 n + 1 + c n c n D 2 n if c n > 1 2
(2.17) with:

c n = Tr{D 2 n -D 3 n } Tr{D n -D 2 n } (2.18)
The polynomials of Eq. (2.17 in the [0, 1] interval. For c = 1/2, where the inflexion point is located at x = 1/2, we found that both functions behave like the McWeeny polynomial, with stationary points at x = 0 and x = 1. For the extreme value c = 1 (c = 0), the inflexion point is located at x = 1 (x = 0) with the stationary end point outside the maximum (minimum) bound of the spectrum, at x > 1 (x < 0), whereas the other stationary point remains fixed at x = 0 (x = 1). The N -representable initial guess as introduced by PM in Ref. [START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF] is generated from the following normalization relation:

D 0 = α(μI -F ) + θI α = min N εmax - μ, M -N μ -εmin μ = Tr{F }/M (2.19a) (2.19b) (2.19c)
where εmax and εmin are estimates of the highest (ϵ max ) and lowest (ϵ min ) eigenvalues of the Fock matrix, respectively. These values are usually accessed, at low cost, using the Gershgorin's formulas [START_REF] Geršgorin | [END_REF][START_REF] Hardy Wilkinson | The algebraic eigenvalue problem[END_REF]:

εmax = max i    F ii + M j̸ =i |F ij |    εmin = min i    F ii - M j̸ =i |F ij |    (2.

Trace-correcting and trace-resetting purifications

To circumvent this issue, Niklasson [START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF] has proposed an alternative method where the N -representability constraints are alleviated and the density matrix purification is performed using the following trace-correcting polynomials:

D n+1 =    P (a) m (D n ) = I -(I -D n ) m (I + mD n ) if Tr{D n } ≤ N P (b) m (D n ) = D m n (I + m(I -D n )) if Tr{D n } > N (2.21)
where (m + 1) gives the order of the polynomial, and n is the index of the recursion. Note that for the special case of m = 1, the conditions on the trace in Eq. (2.21) have to be swapped. The set of polynomials {P (a) m } and {P (b) m } are plotted in Figure {2.4} for m = {1, 2, 3}. For m = 1 (2nd order polynomials), the stationary points are fixed for x = 0 and x = 1, for P Consequently, polynomials with higher flexibililty (m > 2) might be costly in ressources (for large scale systems) compared to lower order polynomials if the total number of purifications needed to reach convergence is not reduced. The initial guess for the TC family of purifications is given by: where β m ∈ [0, 1] is the stable fixed point such that: P (a) m (β m ) = β m , and

D 0,m = 1 -2β m εmax -εmin (ε max I -F ) + β m I (2.
P (b) m (1 -β m ) = 1 -β m .
It is worth emphasizing that normalization of Eq. (2.22) does not enforce -indeed it must not-the trace of D 0 to be equal to the correct value of N . Owing to its simplicity and efficiency, the second-order trace-correcting polynomials (TC2) is the most popular. [START_REF] Mj Cawkwell | Computing the density matrix in electronic structure theory on graphics processing units[END_REF][START_REF] Emanuel | Interior eigenvalues from density matrix expansions in quantum mechanical molecular dynamics[END_REF][START_REF] Mj Cawkwell | Computation of the density matrix in electronic structure theory in parallel on multiple graphics processing units[END_REF] By setting m = 1 in Eq. (2.21), it writes:

D n+1 =    D 2 n if Tr{D n } ≥ N 2D n -D 2 n if Tr{D n } < N with D 0 = (ε max I -F )/(ε max -εmin ) (2.23) (2.24)
An alternative solution to correct the CP deficiences at low or high filling factor, was brought by Niklasson, Tymczak and Challacombe through the trace-resetting mechanism. [START_REF] Anders Mn Niklasson | Trace resetting density matrix purification in o (n) self-consistent-field theory[END_REF] This is a hybrid method involving both trace-correcting and tracepreserving polynomials. The authors proposed to subsitute the robust PMCP by a more flexible, ie. efficient, trace-preserving polynomials for which, if appearing, instabilities are controlled via a resetting option based on the TC2 projections. The working equations of the trace-resetting (TRS) density matrix purification are the following,

D n+1 =        F(D n ) + γ n G(D n ) for γ n ∈ [γ min , γ max ] D 2 n if γ n < γ min 2D n -D 2 n if γ n > γ max with: γ n = N -Tr{F(D n )} Tr{G(D n )}
, γ min = 0, and γ max = 6

(2.25)

(2.26)
where the recipe for initializing the density matrix is identical to the TC2 purification [cf. Eq. ( 2. 24)]. The parameter γ n is analogous to the flexibe inflection point c n [cf.

Eq. (2.17)] but for the composite polynomial of Eq. (2.25) given by: The combination of the these two quartic polynomials (TRS4) results to a trace conserving purifiation for γ bounded in the [0, 6] interval. As mentioned by Niklasson et al., adding the function γG to F and increasing γ continuously change the TRS4 polynomial from F(x) at γ = 0 to the mirror function 1 -F(x) at γ = 6. [START_REF] Anders Mn Niklasson | Trace resetting density matrix purification in o (n) self-consistent-field theory[END_REF] The variations of the TRS4 function with respect to the γ value are plotted in Figure {2.5}, along with the TC2 polynomials. Convergence is achieved when γ → 3, that is, when the TRS4 polynomial transforms to the McWeeny purification: D n+1 = 3D 2 n -2D 3 n . Since the TRS4 function presents (at most), two inflection points, runaway solutions may appear when γ / ∈ [0, 6]. In that case the trace-resetting mechanism supplied by the TC2 projection takes over from the TRS4 and remaps the density matrix within the trace-conserving domain. Minute details including performances of the method can be found in Ref. [START_REF] Anders Mn Niklasson | Trace resetting density matrix purification in o (n) self-consistent-field theory[END_REF].

F(D n ) = D 2 n (4D n -3D 2 n ) G(D n ) = D 2 n (I -D n ) 2
The Table {2.1} summarizes the key points of the density matrix Fermi-Dirac polynomial expansion including the number of MMs performed for each recursive call. We emphasize that the number of MMs is governed by the order of the polynomials. For comparison, the characteristics of LNV density matrix minimization are also reported. In this case, it should be mentioned that there is an additional cost related to the use of the conjugate-gradient routine, which is indicated by the number in parenthesis. This number corresponds to the number of MMs performed during the CG line search. Algorithms used in this work are provided in Appendix C.

Hole-particle canonical purification

It is worth to mention that beyond their mathematical characteristics, some of the density matrix purifications presented above were based (more or less) on physical motivations. Since the early work of McWeeny in 1956, the increasing complexity in attempting to derive more robust and efficient polynomials has reached a stationary point, with for example the TRS4 or other approaches introduced for instance by Mazziotti, [START_REF] David | Towards idempotent reduced density matrices via particle-hole duality: Mcweeny's purification and beyond[END_REF] Kryachko [START_REF] Kryachko | Generalized idempotency purification transform in linear scaling self-consistent field theory[END_REF] or Holas. [START_REF] Holas | Transforms for idempotency purification of density matrices in linearscaling electronic-structure calculations[END_REF] Intuitively, the NVT ensemble appears as the natural framework to derive density matrix purification for isolated system, although hybrid methods (TC2 or TRS4) relying on a ill-defined µ(N )V T ensemble are already very efficient. The PMCP was, so far, the only strictly canonical purification, in the sense that it conserves N -represensability conditions throughout the iterative process and converges systematically as the order of the recursion increases. Nevertheless, Palser and Manalopoulos did not provide any physical interpretation nor insight, for explaining their formulation. We also outlined that, all the DMPE methods presented in this chapter invoked a conditional statement with respect to the trace of the density matrix in order to adjust the polynomial accordingly.

In the work presented below, by re-considering the original proposition of McWeeny described in Section 2.1.1, we introduce a constrained minimization principle where the N -representability conditions are fulfilled from the early steps to the end of the recursion. The very simple purification polynomial emerging from it, called hole-particle canonical purification (HPCP), is given by:

D n+1 = (1 -2c n )D n + 2(1 + c n )D 2 n -2D 3 n , with: c n = Tr{D 2 n -D 3 n } Tr{D n -D 2 n } (2.29)
where c n is the flexible inflexion point already introduced by Palser and Manolopoulos [cf.

Eq. (2.18)]. In terms of both, the one-particle and one-hole density matrix, Eq. (2.29) can be recast as

D n+1 = D n + 2(D 2 n Dn -c n D n Dn ), with: c n = Tr{D 2 n Dn } Tr{D n Dn } (2.30)
From a conditioned D 0 , the recursion relation Eq. (2.29) or Eq. (2.30) is able to deliver the exact ground-state density matrix without the need of correcting the trace, nor adjusting the polynomial during the purification process. A a consequence, the polynomial of Eq. (2.30) is self-consistent. The initial guess suitably conditioned for the HPCP approach is defined according to

D 0 = αD min + (1 -α) D max (2.31)
where α ∈ [0, 1] is the mixing coefficient between D min and D max , which are evaluated from the following recipe:

D min = λ o (µI -F ) + θI D max = λ q (µI -F ) + θI µ = Tr{F } M λ o = min {λ 1 , λ 2 } λ q = max {λ 1 , λ 2 } λ 1 = N M (ϵ max -µ) λ 2 = M -N M (µ -ϵ min ) (2.32a) (2.32b) (2.32c) (2.32d) (2.32e) (2.32f) (2.32g)
For the extended comparison presented in the following section, HPCP denotes for the HPCP associated with the initial guess described above. A Lagrangian formulation for the constrained search for the N-representable one-particle density matrix based on the McWeeny idempotency error minimization is proposed, which converges systematically to the ground state. A closed form of the canonical purification is derived for which no a posteriori adjustment on the trace of the density matrix is needed. The relationship with comparable methods is discussed, showing their possible generalization through the hole-particle duality. The appealing simplicity of this self-consistent recursion relation along with its low computational complexity could prove useful as an alternative to diagonalization in solving dense and sparse matrix eigenvalue problems. C 2016 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4943213] As suggested 60 years ago, 1 the idempotency property of the density matrix (DM) along with a minimization algorithm would be sufficient to solve for the electronic structure without relying on the time consuming step of calculating the eigenstates of the Hamiltonian matrix. The celebrated McWeeny purification formula 2 has inspired major advances in electronic structure theory based on (conjugategradient) DM minimization 3-8 (DMM) or DM polynomial expansion 9,10 (DMPE), where the DM is evaluated by the recursive application of projection polynomials (commonly referred to as purification). DMPE resolution includes the Chebyshev polynomial recursion, 9-15 the Newton-Schultz sign matrix iteration, 16-18 the trace-correcting 19 and the trace-resetting 20 purification (TCP and TRS, respectively), and the Palser and Manolopoulos canonical purification (PMCP). 21 They constitute, with sparse matrix algebra, the principal ingredient for efficient linear-scaling tight-binding (TB) and self-consistent field (SCF) theories. 22,23 Since all these methods were originally derived within the grand canonical ensemble, 24 for a given total number of states (M), none of them are expected to yield the correct number of occupied states (N) unless the chemical potential (µ) is known exactly. As a result, their implementation to the canonical ensemble involves heuristic considerations, where the value of µ 12 or the polynomial expansion 19 is adapted a posteriori to reach the correct value for N, which adds irremediably to the computational complexity. Despite the remarkable performances of the DMPE approaches for solving for sparse 6,25 and dense 26-28 DMs, it remains desirable to develop an approach that overcomes the use of the chemical potential while respecting the canonical requirement of constant-N.

Communication: Generalized canonical purification for density matrix minimization

In this letter, we derive a rigorous and variational constrained search for the one-particle density matrix which a) Electronic mail: lionel.truflandier@u-bordeaux.fr does not rely on ad hoc adjustments and respects the N-representability constraint throughout the minimization process. We shall start from the McWeeny unconstrained minimization of the error in the idempotency of the density matrix, 1 given by minimize

D→ D µ Ω McW {D; (H , µ)}, (1a)
with

Ω McW = Tr{(D 2 -D) 2 }, (1b)
where for a given fixed Hamiltonian 29 H and chemical potential µ, the density matrix D µ is the ground-state for that Hamiltonian and chemical potential. The initial guess (D 0 ) is generally constructed as a function H , suitably scaled,

D 0 = β 1 I + β 2 (µI -H ), (2)
where β 1 and β 2 stand for preconditioning constants such that the eigenvalues of D 0 lie within a predefined range. The double-well shape of the McWeeny function with 3 stationary points: 2 minima at x p = 1 and x p = 0 and 1 local maximum at x m = 1 2 (see Fig. 1(a), red curve), are important features in developing robust DMM algorithms. Finding the minimum of Ω McW would be easily performed by stepwise gradient descent, 1 where the DM is updated at each iteration n,

D n+1 = D n -σ n ∇Ω McW , (3a)
with

∇Ω McW = 2 2D 3 n -3D 2 n + D n , (3b)
and σ n ≥ 0 represents the step length in the negative direction of the gradient. Considering an optimal fixed step length descent (σ = 1/2), on inserting Eq. (3b) into Eq. (3a), the McWeeny purification formula appears,

D n+1 = 3D 2 n -2D 3 n , (4)
where the right-hand side of the equation above can be view as an auxiliary DM. For a well-conditioned D 0 , i.e., λ(D 0 ) ∈ [-1 2 , 3 2 ], repeated application of the recursion identity [Eq. ( 4)] naturally drives the eigenvalues of D n+1 towards 0 or 1. For basic TB Hamiltonians where the occupation factor (θ = N/M) is close to 1/2 and µ can be determined by symmetry 21 or when the input DM is already strongly idempotent, the minimization principle (1a) is able, on its own, to deliver the correct N-representable ground-state DM (D). Beyond these very specific cases, we have to enforce the objective function (1b) to keep N constant during the minimization. From Eq. ( 4), a sufficient condition would be to impose the trace of the auxiliary DM to give the correct number of occupied states. This leads us to solve a constrained optimization problem which can be formulated in terms of the McWeeny Lagrangian (L McW ) by minimize

{D→ D|Tr{D}=N } γ L McW {D, γ; (H ), N }, (5a)
with

L McW = Ω McW -γ Tr{3D 2 -2D 3 } -N , ( 5b
)
where γ is the constant-N Lagrange multiplier. The McWeeny Lagrangian can be minimized using

∇L McW = ∇Ω McW -6γ D -D 2 , (6a)
∂ γ L McW = Tr{3D 2 -2D 3 } -N. (6b)
Taking trace Eq. (6a) we obtain the expression for γ,

γ = 1 3 - 2 3 c - 1 6 d, (7a)
with c = Tr{D 2 -D 3 } Tr{D -D 2 } , (7b)
d = Tr{∇L McW } Tr{D -D 2 } . (7c)
Then, Eqs. (6a) and (7a) are updated at each iteration by requiring Tr{∇L McW } = 0, that is d = 0, for all D. As a result, given D 0 such that Tr{D 0 } = N and [H , D 0 ] = 0, from the fixed-step gradient descent minimization described above, we obtain a recursion formula, Taking advantage of the closure relation,

D n+1 = D n - 1 2 ∇L McW {D n ; γ n } , (8)
D + D = I, (9)
where D stands for the hole density matrix, 30 a more appealing form for the McWeeny canonical purification [Eq. (8)] can be derived by reformulating Eqs. (6a) and (7b) in terms of D and D,

D n+1 = D n + 2 ( D 2 n Dn - Tr{D 2 n Dn } Tr{D n Dn } D n Dn ) . (10)
Notice that since at convergence D D = 0, Tr{D D} must be chosen as the termination criterion in the recursion of Eq. ( 10) to avoid numerical instabilities when approaching the minima. The closed-form of this recurrence relation is remarkable: providing N and H used to build D 0 [Eq. ( 2)], we have a self-consistent purification transformation which should converge to D without any support of heuristic adjustments. Indeed, Eq. ( 10) can also be derived from the PMCP relations by working on both D and D and enforcing relation ( 9) at each iteration (see the Appendix). Consequently, we can also demonstrate 31 that the hole-particle canonical purification (HPCP) of Eq. ( 10) converges quadratically on D as shown in Fig. 2(b).

To assess the efficiency and limitations of the HPCP, we have investigated the dependence of the number of purifications (p) on the occupation factor (θ) and the energy gap (∆ϵ gap = ϵ N +1ϵ N ), defined by the higher-occupied (ϵ N ) and lower-unoccupied (ϵ N +1 ) states. Similarly to the protocol of Niklasson, 15,19 sequences of M × M dense Hamiltonian matrices (M = 100) with vanishing off-diagonal elements were generated, having eigenvalues randomly distributed in the range [-2.5, ϵ N ] ∪ [ϵ N +1 , 2.5] for various ∆ϵ gap ∈ [10 -7 , 1.0]. As a first test, results are compared to the PMCP, 21 along with the original initial guess [Eq. ( 2)], where β 1 = θ and β 2 = min β, β , with

β = θ  H max -µ , β = θ µ - H min , µ ≃  µ = Tr{H } M , (11)
and θ = 1θ = N/M, N being the number of unoccupied states. The lower and upper bounds of the Hamiltonian eigenspectrum (  H min and  H max , respectively) were estimated from to the Geršgorin's disc theorem. 32 The preconditioning of D 0 given in Eq. [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF] guarantees that the DM eigenvalues lie in the interval [0, 1] and gives rise to the following additional constraints:

Tr{D 0 } = N, (12a)
Tr{D 0 } > Tr{D 2 0 } > Tr{D 3 0 }, (12b)
Tr{D 3 0 } > 2Tr{D 2 0 } -Tr{D 0 }, (12c)
which are also necessary and sufficient conditions for c ∈ [0, 1] at the first iteration. Convergence was achieved with respect to the idempotency property, such that Tr{D n Dn } ≤ 10 -6 for all the calculations. Additional tests on the Frobenius norm 33 and the eigenvalues of the converged density matrix (D ∞ ) were performed, using

∥D ∞ ∥ F -  Tr{D ∞ } < 10 -6 , ( 13a
)
∥D ∞ ∥ F -N < 10 -6 , ( 13b
)
∥diag{D ∞ } -diag{I N , 0 N }∥ F < 10 -6 , (13c)
which ensures that, at convergence, the representation of D ∞ is orthogonal, and D ∞ corresponds to D. The variation of the average number of purifications ( p) with respect to θ and ∆ϵ gap is displayed in Fig. 2(a) using a color map for p ∈ [START_REF] Kussmann | Linear-scaling selfconsistent field methods for large molecules[END_REF][START_REF] Mcweeny | The density matrix in many-electron quantum mechanics. i. generalized product functions. factorization and physical interpretation of the density matrices[END_REF]. For a given energy gap, the HPCP shows a net improvement over the PMCP approach regarding moderate low and high occupation factors. Nevertheless, as previously noted by Niklasson and Mazziotti, 19,30 the extreme values of θ remain pathological for the original canonical purification and to a lesser extent for the HPCP. One solution would be to break the symmetry of the McWeeny function by moving x m towards x p or x p depending on the θ value. Basically, this requires a higher polynomial degree for Ω McW , i.e., Tr{(D n -D) 2 } n>2 , resulting in a higher computational complexity. Assuming optimal programming, we emphasize that the PMCP and HPCP involved only two matrix multiplications per iteration. As already proved in Ref. [START_REF] Ditchfield | Self-consistent perturbation theory of diamagnetism: I. a gaugeinvariant lcao method for nmr chemical shifts[END_REF] and highlighted by the energy convergence profiles in Fig. 2(b), the PMCP and HPCP approach E monotonically.

The dependence of p on the band gap plotted in Fig. 2(c) confirms the early numerical experiments, 19,25 where p increases linearly with respect to ln(1/∆ϵ gap ). The influence of θ is clearly apparent if we compare the minimum number of purifications as required for the wider band gap ( y-axis intercept), where for example, with θ = 0.5, both canonical purifications reach the ideal value of about 10 purifications, whereas for θ = 0.05, pHPCP = 23 and pPMCP = 37.

Let us consider how to improve the performance of the canonical purifications by working on the initial guess, regarding the hole-particle equivalence (or duality 30 ). Instead of searching for D, we may choose to purify D, which simply requires replacing D with D in relation [START_REF] Kussmann | Linear-scaling selfconsistent field methods for large molecules[END_REF]. In that case, the initial hole density matrix, satisfying λ( D0 ) ∈ [0, 1], would be given by Eqs. ( 2) and [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF], with β 1 = θ and β 2 = -max β, β . Then, intuitively, the guess for the particle density matrix should be improved by using this additional information. Therefore, a more general preconditioning is proposed,

D + 0 = αD 0 + (1 -α)(I -D0 ), (14)
where α can be viewed as a mixing coefficient. 34 Results obtained with this new preconditioning are plotted in Fig. 2 (notated PMCP+ and HPCP+). As evident from Fig. 2(a), the naive value of α = 0.5 leads to a net improvement of the PMCP and HPCP performances over the range 0.3 < θ < 0.7, inside of which the number of purifications becomes independent of θ. Outside this interval, runaway solutions were encountered due to the ill-conditioning of c, where either of the constraints in Eq. ( 12b) or (12c) is violated. The solution to this problem is to perform a constrained search of α in Eq. ( 14), such that the first inequality of Eq. ( 12b) is respected, that is,

search 0≤α≤1 δ >0      Tr{D 2 0 } =      N -δN, if θ < (1 -δ) N -δ N, if θ > (1 -δ)      , (15)
which leads to solve a second-order polynomial equation in α, at the extra cost of only one matrix multiplication. Obviously, the parameter δ has to be carefully chosen such that the second equality of Eq. (12b) and condition (12c) are also respected. We found δ ≃ 2/3 as the optimal value. 31 From Fig. 2, the benefits of this optimized preconditioning are clear when focussing within the range [0.0, 0.3] ∪ [0.7, 1.0], albeit with one or two extra purifications around the poles θ = {0.3, 0.7}. These benefits are even clearer in Fig. 2(c), where we also show the plots of p as a function of ln(1/∆ϵ gap ) for the test case θ = 0.01. At the intercept, we find pPMCP ≃ 38 compared to pHPCP ≃ 21, showing the improvement bring by the hole-particle equivalence. We have also compared our method against the most efficient of the trace updating methods, TRS4, 20 and find that for non-pathological fillings, the two are comparable in efficiency. For the pathological cases, where TRS4 adjusts the polynomial, we found it more efficient, but at the expense of non-variational behaviour in the early iterations.

To conclude, we have shown how, by considering both electron and hole occupancies, the density matrix for a given system can be found efficiently while preserving Nrepresentability. This opens the door to a more robust, stable ground state minimisation algorithm, with application to standard and linear scaling DFT approaches. L.A.T. would like to acknowledge D. Hache for his unwavering support and midnight talks about how to move beads along a double-well potential.

APPENDIX: ALTERNATIVE DERIVATION OF THE HOLE-PARTICLE CANONICAL PURIFICATION

We demonstrate that by symmetrizing the Palser and Manolopoulos equations with respect to D, the closed-form of Eq. ( 10) appears naturally. Let us start from Eq. ( 16) of Ref. [START_REF] Ditchfield | Self-consistent perturbation theory of diamagnetism: I. a gaugeinvariant lcao method for nmr chemical shifts[END_REF],

for c n ≤ 1 2 , (A1a)
D n+1 = - 1 1 -c n D 3 n + 1 + c n 1 -c n D 2 n + 1 -2c n 1 -c n D n , for c n > 1 2 , D n+1 = - 1 c n D 3 n + 1 + c n c n D 2 n , (A1b)
with c n given in Eq. (7b). We may search for purification relations dual to Eq. (A1), i.e., function of D. We obtain

for cn ≥ 1 2 , ( A2a
) Dn+1 = - 1 1 -cn D3 n + 1 + cn 1 -cn D2 n + 1 -2 cn 1 -cn Dn , for cn < 1 2 , Dn+1 = - 1 cn D3 n + 1 + cn cn D2 n , (A2b)
with cn = 1 -c n . Instead of purifying either D or D, we shall try to take advantage of the closure relation [Eq. ( 9)] in such a way that, if we choose to work within the subspace of occupied states, the purification of D [Eq. (A1)] is constrained to verify D = I -D. By inserting this constraint in Eq. (A2), we obtain

for c n ≤ 1 2 , D n+1 = I - ( - 1 c n (I -D n ) 3 + 2 -c n c n (I -D n ) 2 - 1 -2c n c n (I -D n ) ) , ( A3a
) for c n > 1 2 , D n+1 = I - ( - 1 1 -c n (I -D n ) 3 + 2 -c n 1 -c n (I -D n ) 2 ) . ( A3b
D n+1 = D n + 2 D 2 n Dn -c n D n Dn . ( A4a
)

54

Density matrix purifications and minimizations

Extended comparison of density matrix purifications

For this extended comparison of purification methods we shall consider the same set of Hamiltonian matrices used for the numerical experiment of Section 2.2.3 [cf. Figure 2 of the article and text therein]. The set of 32 Hamiltonians, ie. H ∈ R M ×M with M = 100, was constructed using the following recipe:

1.

H := random{N × N }, such that: H ij ∈ [-2, +2]
2. H := (H + H t )/2, such that:

H ij = H ji 3. H ij := H ij /|i -j| 2 ∀i ̸ = j, such that: lim |i-j|→0 H ij → 0 4. (E, C) := diagonalize{H}, such that: E = diag{ϵ 1 ϵ 2 • • • ϵ M }

Density matrix polynomial expansion 55

As represented on Figure {2.7}a, this protocol was repeated by varying the filling factor (the x-dimension) in the range ]0, 1[, and the energy gap in the range [10 -7 , 1.0] (the y-dimension). This figure clearly highlights the performances of the HPCP with respect to the TC2, TRS4 and PMCP, where each pixel represents the average number of purifications over 32 randomized Hamiltonians. The number of purifications required to achieve convergence for the TC2 is much more higher than for all the other methods (about a factor 2). The TRS4 and HPCP have similar color maps, which are constant along the θ direction, although at the very low or high θ, HPCP looses in efficiency. As discussed in the original article, the impact of the energy gap over the number of purifications is clearly observed. In order to emphasize on the computational performance of the purification methods, the number of matrix multiplications (MMs) is a parameter to consider. From an optimal programming perspective, the numbers of MMs per iteration summarized in Table {2.1} allow to compute the total number of MMs realized to reach convergence. Results obtained for a set of 256 random test Hamiltonians -for two different values of the band gap-as a function of θ are plotted in Figure {2.8}. We found that, while the HPCP and TRS4 are comparable for θ ∈ [0.1, 0.9], TRS4 is more efficient when dealing with pathological cases. Note that the TC2 also presents the inversed bell-shape distribution. Regarding the influence of the gap, we should note that for the HPCP, the standard deviation is almost negligeable, indicating that the HPCP performances do not depend on the eigenvalue distributions. For wide and low gaps, the TC2 remains the most performant purification method when θ ∈ [0.2, 0.8].

Linear scaling strategies

2.3 Linear scaling strategies

The density matrix for insulator contains naturally small elements which can be considered as zero with respect to some threshold. In order to accelerate the calculation, one can get rid of zero elements and work only with significant (also referred to as non-zero elements, nnz). Removing the zero elements involves to truncate the density matrix. Working only with non-zero elements requires a structure representation for the truncated density matrix. Performing density matrix minimizations and purifications with these two ingredients enables to achieve a linear scaling calculation.

Density matrix truncations

The most popular density matrix truncations that we use in this thesis are the following:

Numerical truncation

Since there are small elements in the density matrix, a simple and direct scheme to truncate the density matrix is to drop these small elements by setting them to zero if the absolute value is below a predefined threshold τ [START_REF] John | Linear scaling conjugate gradient density matrix search as an alternative to diagonalization for first principles electronic structure calculations[END_REF][START_REF] Paul E Maslen | Locality and sparsity of ab initio one-particle density matrices and localized orbitals[END_REF]. With respect to the chosen numerical threshold, the density matrix elements are said to be filtered such as

D = FILTER( D , τ ) =    D ij if |D ij | > τ 0 otherwise (2.33)
D is the truncated density matrix. In a minimization or purification algorithm, we apply the truncation (2.33) after each MM.

Radial truncation

The magnitude of density matrix elements depends on the distance between the basis function centers of atom centers. Therefore, one way to truncate the density matrix is to neglect the matrix elements that correspond to distances between basis function centers larger than a predefined cutoff radius R c [START_REF] Li | Density-matrix electronic-structure method with linear system-size scaling[END_REF][START_REF] Larsen | Direct optimization of the atomic-orbital density matrix using the conjugate-gradient method with a multilevel preconditioner[END_REF][START_REF] Shao | Curvy steps for density matrix-based energy minimization: Application to large-scale self-consistent-field calculations[END_REF][START_REF] Bates | Comparison of conjugate gradient density matrix search and chebyshev expansion methods for avoiding diagonalization in large-scale electronic structure calculations[END_REF][START_REF] Liang | Improved fermi operator expansion methods for fast electronic structure calculations[END_REF][START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF]. Hence, the filtered matrix elements are such as

D = FILTER( D , R c ) =    D ij if |r i -r j | < R c 0 otherwise (2.34)
In a minimization or purification algorithm, we apply the truncation (2.34) only on the density matrix D as shown in Algorithm {10} (step ( 9)). We apply the truncation above all else (any polynomial, function or gradient) depending on the density matrix.

Sparse matrix representations

Using the truncation reinforces the matrix sparsity. The access to only non-zero elements is a great advantage for the minimizations and purifications algorithms. Instead of using the standard algebra for matrices in dense format, it is better to employ sparse matrices algebra. In the latter, the sparse matrix is compressed into some matrix data structures. A more detailed discussion of data structures for sparse matrices can be found in Ref. [START_REF] Rubesson | Matrix Algebra for Quantum Chemistry[END_REF]. There are several matrix data structures for sparse matrices [START_REF] Emanuel H Rubensson | A hierarchic sparse matrix data structure for large-scale hartree-fock/kohn-sham calculations[END_REF][START_REF] Emanuel H Rubensson | Sparse matrix algebra for quantum modeling of large systems[END_REF][START_REF] Langr | Adaptive-blocking hierarchical storage format for sparse matrices[END_REF][START_REF] Hackbusch | A sparse matrix arithmetic based on\ cal h-matrices. part i: Introduction to {\ Cal H}-matrices[END_REF][START_REF] David R Bowler | Parallel sparse matrix multiplication for linear scaling electronic structure calculations[END_REF].

For instance, the tool package SPARSKIT [START_REF] Saad | Sparskit: A basic tool kit for sparse matrix computations[END_REF] presents a multitude of sparse matrices formats such as the compressed sparse column (CSC) format. The CSC representation is the simplest representation for the sparse matrices, that we use in this thesis.

Applications to carbon nanotubes

In this section we shall compare the performances of O(N ) density matrix energy functional minimizations and polynomial expansions presented in Sections 2.1.2 and 2.2, with respect to the truncation schemes. For that study, the Pariser-Parr-Pople (PPP) semi-empiral method described in Section 1.4 will be applied to a set of π-conjugated systems: the carbon nanotubes. We recall that, within the SCF-PPP framework, the ideal case of the one-electron one-orbital picture is imposed. For neutral π-conjugated network, this implies that the filling factor is automatically set up to 1/2.

Carbon nanotubes

The π-conjugated systems considered in this part are the carbon nanotubes (CNTs). A CNT is a compound which has the shape of a cylindrical tube made of a rolled single layer of carbon atoms. This single layer of carbon atoms is known as graphene. [START_REF] Ah Castro Neto | The electronic properties of graphene[END_REF][START_REF] Konstantin | Graphene: status and prospects[END_REF] A detailed description on CNTs and their properties can be found, for instance, in the Refs. [START_REF] Saito | Physical properties of carbon nanotubes[END_REF][START_REF] Reich | Carbon nanotubes: basic concepts and physical properties[END_REF][START_REF] Smalley | Carbon nanotubes: synthesis, structure, properties, and applications[END_REF]. The rolling of the graphene to form CNTs is modulated by the chirality vector, ⃗ C h , reproduced on Figure {2.9a}. This vector has two components, such that: where {⃗ a 1 , ⃗ a 2 } are the unit vectors of the graphene sheet. The values of the pair of indices (p, q) define the way the graphene sheet is wrapped. From the definition of ⃗ C h and the symmetry of translation, three types of CNTs can be generated:

⃗ C h = p ⃗ a 1 + q ⃗ a 2 , (p, q) ∈ Z (2.35)
• the armchair type, where p = q,

• the zigzag type, where p = 0 or q = 0,

• the chiral type, where p ̸ = q.

They are represented in Figure {2.9}. Note that p and q are generally given as positive integers with p > q. The indices (p, q) do not only determine the carbon atoms' arrangement in the tube, but they can also provide information on the properties of the CNT. [START_REF] Saito | Physical properties of carbon nanotubes[END_REF] For instance:

• the armchair CNTs are all metallic,

• the zigzag CNTs are metallic when p is a multiple of 3,

• the chiral CNTs are metallic when (p -q) is a multiple of 3. semi-empirical and the first-principles approaches, with differences which do not exceed 5% (for the Ohno parametrization only). This demonstrates that the major part of the two-electron π interactions is well reproduced by the PPP model.

CNTs # atoms

Numerical truncation for SCF calculations

SCF calculation based on density matrix solvers is controlled via two threshold parameters: (i) the first one is related to the density matrix (DM) convergence inside each SCF cycle, ie. τ D = ∥D n+1 -D n ∥,1 (ii) the second one is related to the convergence of the SCF procedure itself, and based on an energy criterion, ie. τ SCF = |E n+1 -E n |. For the results presented in the next section, we have used τ D = 10 -3 and τ SCF = 10 -6 . In order to achieve the linear scaling regime, the density matrix purifications and minimizations must be supported by techniques which reinforce the density matrix sparsity (cf. Section 2.3.1). Using the numerical truncation scheme of Eq. (2.33), the sparsity of the DM is controlled by the numerical threshold τ . In Table {2.4} are reported the SCF energies obtained for the set of CNTs described in Table {2.3} using the hole-particle canonical purification (HPCP) as density matrix solver. Different values of τ were considered within the range [10 -5 , 10 -7 ]. ∆E is the error between the exact SCF energy obtained (without truncation) and the truncated density matrix purification. N SCF is the number of SCF iterations achieved to reach convergence, and nnz is the number of non-zero DM elements. As expected, by reading Table {2.4} (along the rows), for all the CNTs, we observe that the energy approaches the exact value as τ decreases, ie. the density of the non-zero elements increases. Interestingly, by reading the same table (along the columns), for a given threshold, and given the π-π * gaps reported in Table {2.3}, we note that the sparsity of the DM decreases for increasing energy gap. For reasonable values of the threshold (τ ≥ 10 -6 ), the convergence is roughly achieved after a number of 6 cycles, whatever the CNT. This is to be compared with the trend obtained for highly sparse DM where, depending on the CNT, twice as many SCF cycles are required. To illustrate the influence of the numerical truncation on the density matrix during the SCF, sparsity patterns were represented in Figure {2.10} for the CNT [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF]5). These sparsity patterns correspond to an average of the density matrices over the six SCF cycles summarized into four sequences. In Figure {2 Whereas the numerical truncation operates directly on the DM elements without consideration on the density matrix topology, the radial truncation of Eq. (2.34) assumes that the relevant elements are localized within spheres of a certain radius centered on each atom. As a result, controlling the radius of the spheres, controls the number of non-zero elements which, compared to numerical truncation, has fixed positions within 66 Density matrix purifications and minimizations the DM network. As a result, the cutoff radius R c fixes the superior limit beyond which the density matrix elements are enforced to be zero. In this work, the upper bound of R c is defined by the lattice parameter a in the longitudinal-like direction, such that: R c ≤ a/2 [cf. Figure {2.12}]. It is clear that larger is the radial cutoff, higher is the number of significant nnz of the density matrix. We have investigated the convergence of the SCF energy with respect to R c . Results are displayed in Figure {2.13}. For each type of CNT -zigzag or chiral-we observe a variational-like convergence in agreement with the earlier works of Refs. [START_REF] Hernández | Linear-scaling density-functional-theory technique: the density-matrix approach[END_REF][START_REF] Otsuka | Accuracy of order-n density-functional theory calculations on dna systems using conquest[END_REF]. However, the purifications used with the radial truncation can be problematic for some cutoff radii. For instance, in Figure {2.14} is represented the convergence of the energy within a single SCF cycle (ie. for a fixed Fock matrix). For the variational trace-conserving purification method, we observe an increase of the energy after the 5th iteration [blue curve on panel (b)], which indicates that the (local) minimum can not be reached. Note that the relation E n+1 > E n was already proposed by Palser and Manolopoulos as the termination criterion when radial Density matrix purifications and minimizations Density matrix purifications and minimizations truncation is applied. [START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF] As a result, unlike the numerical truncation which drops small elements disregarding their locations in the density matrix, and can be used with any of the density matrix solver, the radial truncation seems to be, for some of the cases treated here, incompatible with the purifications.

Radial truncation for SCF calculations

As an attempt of solution, we have performed the following test: when E n+1 > E n , the radial truncation is switched off until E n ′ +1 < E n ′ , where then, the truncation is switched back on. Result (denoted by Radial * ) is dispalyed in red on the Figure {2.14}. In that case, we observed that the monotonic convergence is not interrupted at the crossing point corresponding to n = 5. Note that compared to the LNV minimization, the output energy (as obtained at the end of the purification) is lower. Nevertheless, for all the cases, we found that the LNV is more robust than the HPCP with regard to the radial truncation. In the following, only results obtained for the purifications and radial truncation at non-problematic cutoff radii are discussed. We now illustrate the radial truncation using the CNT [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF]5) at R c = 10 Å, as already used for the numerical truncation in Section 2.4.2. The evolution of the density matrix sparsity during the SCF iterations is summarized in Figures {2.15} and {2.16}. Since the number of non-zero elements is predifined by the radial truncation, compared to Figures {2.10} and {2.11}, we do not observe the "let-it-grow" evolution of the numerical scheme. The progression of the red color region for the non-zero elements is limited.

Linear scaling SCF calculations and conclusion

Applications of the radial truncation to achieve linear scaling on large systems are presented in Figure {2.17}. For that purpose, we have replicated the CNT [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF]5) in the longitudinal direction up to about 10,000 atoms. The figure displays the variation of the CPU time as a function of the system size. Calculations were performed using the diagonalization (Diag) and the following density matrix solvers: the standard LNV minimization, McWeeny purification (McW), trace correcting (TC2), canonical purification (PMCP), hole particle canonical purification (HPCP), and the trace resetting (TRS4). Figure {2.17(a)} shows that purifications are more efficient than the LNV minimization. In this ideal case of half-filling, the "fixed point" McWeeny polynomials presents the best performance. These results are in agreement with other comparative studies. [START_REF] Rudberg | Assessment of density matrix methods for linear scaling electronic structure calculations[END_REF][START_REF] Daniel | Comparison of two genres for linear scaling in density functional theory: Purification and density matrix minimization methods[END_REF] The influence of the density matrix sparsity on the linear scaling behavior can be appreciated in Figure {2.17(b)}, where the calculation time is reduced by more than a factor 2 when the truncation radius is decreased from 50 to 10 Å. Note in passing that, all the density matrix methods are cleary proved more efficient than the diagonalization. Using the two previous forms of matrix illustration [cf. Section 2.4.2 or Section 2.4.3], we display in Figure {2.18} the profile of the C and D matrices obtained at the end of SCF convergence, for the CNT [START_REF] Kohn | Density functional and density matrix method scaling linearly with the number of atoms[END_REF]5). C is completely full, while D is sparse with decaying profile. That justifies the net difference between the diagonalization scaling and density matrix methods scaling.

In this chapter we have shown standard diagonalization can be circumvented using density matrix based solvers, which, when combined with sparse matrix algebra can lead to linear scaling regime. We found that, enforcing sparsity by the numerical approach is more robust compared to the radial truncation. Nevertheless, in that case, the number of non-zero elements is controlled a posteriori. For all the density matrix purifications, we found that radial truncation may cause convergence problems. Obviously, definitive conclusions on the advantages and drawbacks of the various schemes are beyond the scope of this work. A more systematic comparison using non-orthogonal basis sets of various sizes, along with a browder set of systems should be beneficial to identify and analyse limitations in accuracy and stability of O(N ) methods. The main aim of this dissertation is to derive and test linear scaling algorithms for density matrix perturbation theory, which will be introduced in the next chapter.

Density matrix perturbation theory

In Chapters 1 and 2 of this manuscript, we have presented the resolution of timeindependent Liouville-von Neumann equation in order to determine the energy of a molecular system. However, to relate the results from quantum chemical calculations to experiment, it is essential to compute quantities that are directly comparable to measurements. For this purpose, the density matrix of the ground state obtained from the resolution of the Liouville-von Neumann equation is not sufficient. It is therefore necessary to compute further quantities that characterize the molecular system of interest. Theses quantities can be classified as follow:

1. Energy differences, such as reaction energies, dissociation energies, that involve energy information at different points on the Born-Oppenheimer potential energy surface.

2. Molecular properties, like dipole moment, polarizabilities, vibrational frequencies, nuclear magnetic resonance parameters, that require information of perturbed electronic states at a single point on the potential surface.

3. Transitions energy between different electronic states, as for instance, electronic excitation energies, radiative life times, that involve information for electronic states coupling.

The concept of spectroscopy refers to the observation of a physical phenomenon onto an energy scale or any quantity related to an energy, like frequency or wave length. Until nowadays, the spectroscopy principle is greatly expanding in many research fields such as astronomy, biophysics, chemistry, acoustics. The basic idea of the spectroscopy consists to subject to radiation the matter and measure the response. Comparing the original radiation with the response, one can extract information related to some of its intrinsic properties, for instance: structural, electronic or magnetic. A specific radiation, ie. wavelength range, allows to probe specific properties of the system. For the second class of properties enumerated above, the perturbation theory is necessary (and sufficient) to access various spectroscopic observables. This Chapter takes a look beyond the timeindependent Liouville-von Neumann equation, to derive density matrix time-independent perturbation theory. The density matrix methods such as the minimizations and the purifications, highlighted in Chapter 2, are used for solving the following unperturbed SCF equations

F (D n ) D n+1 -D n+1 F (D n ) = 0 subject to: Tr{D n } = N D 2 n = D n , ∀ n (3.4a) (3.4b) (3.4c)
Let us now consider that this system is disturbed by an external time-independent perturbation, which takes the form of a matrix W . The perturbed Fock and one-particle density matrices formally read

F λ := F + λW F D λ := D + λW D (3.5a) (3.5b)
where we have introduced a scaling parameter λ ∈ [0, 1] to modulate the strength of the perturbation. Consequently, the analogue of the non-perturbed SCF equations (3.4), for

78

Density matrix perturbation theory the perturbed case, are given by

F λ (D λ,n ) D λ,n+1 -D λ,n+1 F λ (D λ,n ) = 0 subject to: Tr{D λ,n } = N, D 2 λ,n = D λ,n , ∀ n (3.6a) (3.6b) (3.6c)
Since we are concerned about the response of the system to the perturbation, we need to evaluate the variation of the energy with respect to an infinitesimal variation of the perturbation strength. It is worth to note that the perturbation theory is valid if the response is small compared to the eigenspectrum of F . There exist two ways for evaluating the energy response: (i) the numerical solution based on finite differences [cf.

Appendix E], or (ii) the analytical solution based on the perturbation expansion of the relevant quantities. Density matrix perturbation theory presented in this Chapter is obviously dealing with the second possiblity. In that context, one can write the matrices F λ and D λ down as power series in λ following the standard Taylor expansion:

D λ = D + λD (1) + λ 2 D (2) + ... + λ k D (k) F λ = F + λF (1) + λ 2 F (2) + ... + λ k F (k) (3.7a) (3.7b)
with D (k) and F (k) the shorthand notations for the kth order variation of the density and Fock matrix. More explicitly, they should read

D (k) := 1 k! d k D dλ k F (k) := 1 k! d k F dλ k (3.8a) (3.8b)
From here, the idea is to apply the perturbation theory on the stationary conditions of Eq. (3.6). On inserting the expansion of Eq. (3.7a) into (3.6c), and equating the perturbation orders of the left and right sides, we obtain the following perturbed idempotency relations:

DD (1) + D (1) D = D (1)
DD (2)

+ (D (1) ) 2 + D (2) D = D (2)
DD (3) + D (1) D (2) + D (2) D (1) + D (3) D = D (3) . . .

DD (k) + D (1) D (k-1) + D (2) D (k-2) + ... + D (k) D = D (k)
or more generally: [F, D (1) ] + [F (1) , D] = 0 [F, D (2) ] + [F (1) , D (1) ] + [F (2) , D] = 0

D (k) = k l=0 D (l) D (k-l) (3.9a) (3.9b) (3.9c) (3.9d) (3.
[F, D (3) ] + [F (1) , D (2) ] + [F (2) , D (1) ] + [F (3) , D] = 0 . . .

[F, D (k) ] + [F (1) , D (k-1) ] + ... + [F (k) , D] = 0 or more generally:

k l=0 [F (l) , D (k-l) ] = 0 (3.10a) (3.10b) (3.10c) (3.10d) (3.10e)
The perturbed idempotency constraints (3.9) and SCF conditions (3.10) are both the background of the density matrix perturbation theories which are presented in the next Section. From these relations, we may also emphasize that the evaluation of the perturbed density matrix at order (k) is based on the knowledge of the order (k -1). In other words, evaluation of the perturbed quantities requires -in principle-to proceed order by order.

Wavefunction coupled perturbed self-consistent field formulation

The standard density matrix perturbation method is the atomic orbitals coupled perturbed self-consistent-field [25,[START_REF] Stevens | Perturbed hartree-fock calculations. i. magnetic susceptibility and shielding in the lih molecule[END_REF]26,[START_REF] T_ | Self-consistent perturbation theory for conjugated molecules[END_REF] (AO-CPSCF). In order to derive it, we simply need to decompose the non-perturbed Fock matrix in terms of its eigenvalues and (one-state) projection operators built from the eigenvectors.

The density matrix perturbation theory as developed by Diercksen and McWeeny [26] is based on the partitioning of D (k) into four distinct contributions and their resolutions. Using the closure relation (3.1c), any operator X can be recast into the following projected components:

X = (D + D)X(D + D) = DXD + DX D + DXD + DX D = X oo + X ov + X vo + X vv (3.11a) (3.11b) (3.11c)
where the subscripts (oo) and (vv) designate the occupied-occupied, and virtual-virtual contributions related to the original orthogonal subspaces H occ and H virt , for the occupied states and for the unoccupied states, respectively. Likewise, (ov) and (vo) stand for the perturbation induced by coupling terms associated with subspaces H occ-virt and H virt-occ , respectively. The spectral resolution of Eq. (1.94b) already introduced in Section 1.5 allows any projected matrices of Eq. (3.11c) to be decomposed into a sum of single-projected components, following:

X oo = DX D = i,j D i XD j X ov = DX D = i,j D i X Dj X vo = DX D = i,j
Di XD j

X vv = DX D = i,j Di X Dj (3.12a) (3.12b) (3.12c) (3.12d)
Note that the spectral resolution (3.12) implies the eigenstates to be known, which irremediably involves the diagonalization of the unperturbed Fock matrix.

First-order response

Let us start with the first-order of perturbation. On applying the projection decomposition of Eq. (3.11) to both sides of Eq. (3.9a), we obtain:

2D (1) oo + D (1) ov + D (1) vo = D (1) oo + D (1) ov + D (1) vo + D (1) vv (3.13) Proceeding by identification, it can be deduced that D (1) oo = 0, D (1) vv = 0 (3.14)

3.2 Wavefunction coupled perturbed self-consistent field formulation 81 and D (1) = D (1) ov + D (1) vo = D (1) ov + D (1) † ov (3.15a)

(3.15b)
For the last statement (3.15b), we relied on the symmetry property of the perturbed density matrix, that is, D (1) † = D (1) . The relation (3.15b) demonstrates that the determination of D (1) involves only the evaluation of the occupied-virtual transition matrix. For that purpose, multiplying Eq. (3.10a) by D from the left, and by D from the right, we obtain

F (1) ov (D (1) ov,n ) = [F, D (1)
ov,n+1 ] subject to: Tr{D (1) ov,n } = 0, ∀n

(3.16a) (3.16b)
where we have re-introduced the iteration indice n of the self-consistent resolution, and the dependence of the perturbation matrix F (1) ov over D (1) ov . 1 This relation constitutes the first-order -coupled perturbed self-consistent field (CPSCF)-density-matrix equations, analogous to the non-perturbed Eq. (3.4). Given an initial guess for the first-order perturbed density matrix D (1) ov [cf. Eq. (3.14)], we project the first-order perturbed Fock matrix F (1) , to build F (1) ov , and finally solve Eq. (3.16). The iterative process is repeated until convergence is achieved. By substituting Eq. (1.94a) into (3.16), one obtains ov,ij = D i D (1) Dj

F (1) ov,ij = D i F (1) Dj (3.17a) (3.17b) (3.17c)
which yields to the well-known sum-over-states (SOS) first-order equation:

i,j D (1) ov,ij = i,j F (1) ov,ij (ϵ i -εj ) ⇔ D (1) ov = i,j F (1) ov,ij (ϵ i -εj ) (3.18)

Second-order response

The second-order equation can be derived by applying the resolution of identity to both side of Eq. (3.9b). Keeping the notations of Eq. (3.11), we have 2D (2) oo + D (2) ov + D (2) vo + (D (1) D (1) ) oo + (D (1) D (1) ) vv + (D (1) D (1) ) ov + (D (1) D (1) ) vo = D (2) oo + D (2) ov + D (2) vo + D (2) vv (3.19) By resolving the product of first-order perturbed density matrices according to

D (1) D (1) = D (1)
ID (1) = D (1) (D + D)D (1) = D (1) (D 2 + D2 )D (1) we obtain (D (1) D (1) ) oo = D (1) oo D (1) oo + D (1) ov D (1) vo (D (1) D (1) ) vv = D (1) vv D (1) vv + D (1) vo D (1) ov (D (1) D (1) ) ov = D (1) ov D (1) vv + D (1) oo D (1) ov (D (1) D (1) ) vo = D (1) vo D (1) oo + D (1) vv D (1) vo

(3.21a) (3.21b) (3.21c) (3.21d)
On inserting the right-hand side (rhs) of Eqs. (3.21) into (3.19), and using the properties (3.14), we have 2D (2) oo + D (2) ov + D (2) vo + D (1) ov D (1) vo + D (1) vo D (1) ov = D (2) oo + D (2) ov + D (2) vo + D (2) vv (3.22) Therefore, it comes D (2) oo = -D (1) ov D (1) vo , D (2) vv = +D (1) vo D (1) ov (3.23) Unlike the 1st-order perturbation, the diagonal components of the 2nd-order perturbed density matrix are likely to be non zero and can be computed from the 1st-order perturbed density matrix. Relying furthermore on the symmetry of the perturbed density, it leaves only the occupied-virtual coupling matrix to evaluate. On resolving the 2nd-order perturbed Fock matrix using Eq. (3.10b), we obtain 1) , D (1) ] ov (3.24) 3.2 Wavefunction coupled perturbed self-consistent field formulation

F (2) ov = [F, D (2) ov ] + [F (

83

Using the spectral resolution of the non-perturbed Fock matrix and the perturbed density matrix, Eq. (3.24) transforms as

F (2) ov = i,j D (2)
ov,ij (ϵ i -εj ) + [F (1) , D (1) ] ov,ij (3.25) which leads to the 2nd-order SOS equation

D (2) ov = i,j (ϵ i -εj ) -1 F (2)
ov -[F (1) , D (1) ] ov ij (3.26) The final 2nd-order perturbed density matrix is obtained summing over D (2) ov , its conjugatetransposed, D (2) vo , and the diagonal contributions of Eq. (3.23).

Third-order response

Using the same route than for the first-and second-order, the third-order response equations are derived from Eq. (3.9c). This yields to 2D (3) oo + D (3) ov + D (3) vo + (D (1) D (2) ) oo + (D (1) D (2) ) vv + (D (1) D (2) ) ov + (D (1) D (2) ) vo + (D (2) D (1) ) oo + (D (2) D (1) ) vv + (D (2) D (1) ) ov + (D (2) D (1) ) vo = D (3) oo + D (3) ov + D (3) vo + D (3) vv (3.27) where (D (1) D (2) ) oo = D (1) oo D (2) oo + D (1) ov D (2) vo (D (1) D (2) ) vv = D (1) vv D (2) vv + D (1) vo D (2) ov (D (1) D (2) ) ov = D (1) ov D (2) vv + D (1) oo D (2) ov (D (1) D (2) ) vo = D (1) vo D (2) oo + D (1) vv D (2) vo (D (2) D (1) ) oo = D (2) oo D (1) oo + D (2) ov D (1) vo (D (2) D (1) ) vv = D (2) vv D (1) vv + D (2) vo D (1) ov (D (2) D (1) ) ov = D (2) ov D (1) vv + D (2) oo D (1) ov (D (2) D (1) ) vo = D (2) vo D (1) oo + D (2) vv D (1) vo 2D (3) oo + D (3) ov + D (3) vo + D (1) ov D (2) vo + D (2) ov D (1) vo + D (1) vo D (2) ov + D (2) vo D (1) ov = D (3) oo + D (3) ov + D (3) vo + D (3) vv (3.29) Density matrix perturbation theory Consequently, it follows

(3.28a) (3.28b) (3.28c) (3.28d) (3.28e) (3.28f) (3.
D (3) oo = -D (1)
ov D (2) vo + D (2) ov D (1) vo D (3) vv = + D (1) vo D (2) ov + D (2) vo D (1) ov

(3.30a) (3.30b)
Once again, these last equations show that the 2nd-order density matrix is necessary and sufficient to compute the 3rd order diagonal components. Again, at this point, we emphasize that only the occupied-virtual transition matrix needs to be evaluated selfconsistently since the perturbed density matrix is (at least) Hermitian. Using Eq. (3.10c), the 3rd-order perturbed Fock matrix reads ) , D (2) ] ov + [F (2) D (1) ] ov (3.31) Relying on the spectral resolution, this relation transforms to ) , D (2) ] ov + [F (2) D (1) ] ov ij (3.32) which leads to the 3rd-order SOS equation

F (3) ov = [F, D (3) ov ] + [F (1
F (3) ov = i,j D (3) ov,ij (ϵ i -εj ) + [F (1
D (3) ov = i,j (ϵ i -εj ) -1 F (3)
ov -[F (1) , D (2) ] ov -[F (2) , D (1) ] ov ij (3.33)

kth-order response

As a matter of fact, from the expansions (3.9) and (3.10), we can generalize the response equations to any kth-order, by proceeding in the same way and using the components from the lower orders. Then, the (diagonal) fixed components are defined according to

D (k) oo = - k-1 i=1 D (i) oo D (k-i) oo + D (i) ov D (k-i) vo (3.34) D (k) vv = k-1 i=1 D (i) vv D (k-i) vv + D (i) vo D (k-i) ov (3.35)
whereas the (off-diagonal) transition matrices, which involve a self-consistent resolution, are defined according to

D (k) ov = i,j (ϵ i -εj ) -1 F (k) ov -M (k)
ov ij (3.36) 3.3 Density matrix coupled perturbed self-consistent field formulation 85 (3.37) using the shorthand notations

M (k) ov,ij = k-1 l=1 (D (k) F (k-l) ) ov -(F (k) D (k-l) ) ov ij
F (k) ov,ij = D i F (k) Dj (D (k) F (k-l) ) ov,ij = D i D (k) F (k-l) Dj , (F (k) D (k-l) ) ov,ij = D i F (k) D (k-l) Dj
The diagonal components are used as appropriate guess to initiate the perturbed density matrix whereas the off-diagonal components are taken as the iterative part using the eigenvectors of the non perturbed Fock matrix.

Density matrix coupled perturbed self-consistent field formulation

In the AO-CPSCF formalism, calculation of the perturbed density matrix D (k) requires the eigenstates of the unperturbed Fock matrix, which implies a high computational effort for systems of increasing size. In this work, we also consider the density matrix-based perturbation formalism originally proposed by Oschenfeld and Head-Gordon, [START_REF] Ochsenfeld | A reformulation of the coupled perturbed self-consistent field equations entirely within a local atomic orbital density matrix-based scheme[END_REF] and reformulated by Oschenfeld and Kussmann. [START_REF] Kussmann | Linear-scaling method for calculating nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals within hartree-fock and density-functional theory[END_REF][START_REF] Kussmann | A density matrix-based method for the linear-scaling calculation of dynamic second-and third-order properties at the hartree-fock and kohn-sham density functional theory levels[END_REF][START_REF] Ochsenfeld | Ab initio nmr spectra for molecular systems with a thousand and more atoms: A linear-scaling method[END_REF] The approach relies on solving -self-consistently-the commutation relations (3.10) using a conjugate-gradient-based minimization. It will be referred in this manuscript to as CG-CPSCF. The simple derivation [START_REF] Kussmann | Linear-scaling method for calculating nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals within hartree-fock and density-functional theory[END_REF] of the CG-CPSCF formalism from the background equations (3.9) and (3.10) is presented below.

First-order response

The idea behind the Oschenfeld and Kussmann method is to constrain Eq. (3.10a) to commute with the unperturbed density matrix. That is, multiplying Eq. (3.10a) from the left and from the right separately by D, and substracting, leads to F, D, D (1) + 2DF (1) D -D, F (1) = 0 (3.39) This correponds to the first-order CG-CPSCF equation. Relying on the symmetry properties for the perturbed density and Fock matrices, this equation is also Hermitian. It is worthwhile to note that on multiplying the CG-CPSCF equation (3.39) from the left by D, and from the right by D, the AO-CPSCF equation Eq. (3.16) is recovered. This clearly demonstrates the relationships between the two formalisms. Nevertheless, unlike the AO-CPSCF, the resolution of the first-order CG-CPSCF equation yields directly to the first-order perturbed density matrix D (1) .

Second-and third-order response

If we want to write down the second-order CG-CPSCF equations using only the perturbed density matrices, that is, without relying on the resolution of the identity given in Eq. (3.11)

D (2) oo = -(D (1)
D (1) ) oo = -DD (1) D (1) D D (2) vv = +(D (1) D (1) ) vv = + DD (1) D (1) D

If we want to circumvent the use of D to define D (2) vv and use instead D, first, we should perform the following transformations: D (2) vv = DD (1) D (1) D = (D (1) D (1) ) vv = D (1) vv D (1) vv + D (1) vo D (1) ov [using Eq. (3.21b)] = D (1) vo D (1) ov [using Eq. (3. Also, by noting that:

D (2) vv = DD (2) D = D( DD (2) D) D = D(D (2) vv ) D (3.41)
we can deduce the expression for 2nd order CG-CPSCF initial guess:

D (2) oo = -DD (1) D (1) D, D (2) vv = D (1) DD (1) (3.42)
As for the first-order response, the commutator between the non-perturbed density matrix and Eq. (3.9b) yields to the second-order CG-CPSCF equation,2 following: F, D, D (2) + 2DF (2) D -D, F (2) = D, D (1) , F (1) (3.43)

Using Eqs. (3.14) and (3.28) Eq. (3.28f), the diagonal contributions for the 3rd-order response can be recast as: (1) ) oo + (D (1) D (2) ) oo = -(DD (2) D (1) D + DD (1) D (2) D) D (3) vv = + (D (2) D (1) ) vv + (D (1) D (2) ) vv = D(D (2) D (1) ) D + D(D (1) D (2) ) D

D (3) oo = -(D (2) D
Applying the transformations (3.40) and (3.41), we obtain:

D (3) vv = DD (3) vv D = D(D (2) DD (1) ) D + D(D (1) DD (2) ) D (3.45)
The 3rd order CG-CPSCF initial guess is therefore given by

D (3) oo = -D(D (2) D (1) + D (1) D (2)
)D D (3) vv = D (2) DD (1) + D (1) DD (2) Finally, 3rd-order CG-CPSCF equation reads:

F, D, D (3) + 2DF (3) D -D, F (3) = D, D (1) , F (2) + D, D (2) , F (1) (3.47)

Compared to AO-CPSCF equations (3.31), the solution of this equation should lead to D (3) ov , if we perform a spectral resolution of D (3) .

kth-order response

Consequently, we can easily generalize the CG-CPSCF equations at any kth-order, the diagonal components D (k) oo and D (k) vv defining the initial guess, according to

D (k) oo = -D k-1 i=1 D (k-i) D (i) D (3.48) D (k) vv = k-1 i=1 D (k-i) DD (i) (3.49) 88
Density matrix perturbation theory with the kth-order transition matrices being the solutions of the following equation

F, D, D (k) + 2DF (k) D -D, F (k) = D, k i=1 D (k-i) , F (i) (3.50)
In Eq. (3.50), the right-hand side changes with respect to the density and the Fock matrices at lower orders, whereas the left-hand side (lhs) depends only on the kth order density matrix to be determined. As a result, solving the CG-CPSCF equation is analogous to solve a linear system of equations AX = B. The description of the CG-CPSCF equation resolution using a conjugate-gradient algorithm is given in Appendix D.

Perturbed projection by trace-correcting purification

Derived from the generalized equations for the SCF conditions (3.10) and for the idempotency constraints (3.9), the AO-CPSCF method uses the eigenstates of the unperturbed Fock matrix, while the CG-CPSCF method employs uniquely the density matrix. Both approaches require to proceed order by order. In other words, given a perturbed Fock matrix we are solving a linear problem at each order. On the other hand, it is possible to compute, the perturbed density matrix at a desired order without prior calculation of the lower terms, ie. all the perturbed density matrices of lower orders, all of them, being computed on-the-fly during the CPSCF processus. In other words, we are solving a non-linear set of equations. This approach proposed by Weber and co-workers [START_REF] Weber | Ab initio linear scaling response theory: Electric polarizability by perturbed projection[END_REF][START_REF] Weber | Higher-order response in o (n) by perturbed projection[END_REF], is based on inserting the perturbative expansions (3.7) obtained for the density and Fock matrices within the trace correcting purification (TCP) formalism. [START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF] By doing so, they show that the approach is capable to purify the kth and lower orders perturbed density matrices within the SCF loop. This perturbed purification method is called the perturbed projection by the trace-correcting purification (TC2-CPSCF). The unperturbed TCP aims to purify simultaneously each perturbed density matrix following the relation:

D λ,n+1 =    D 2 λ,n if Tr{D λ,n } ≥ N 2D λ,n -D 2 λ,n if Tr{D λ,n } < N (3.51)
using the following initial guess:

D λ,0 = (ϵ max I -F λ )/(ϵ max -ϵ min ) (3.
D n+1 =    ( 3 k=0 D (k) n ) 2 if Tr{ 3 k=0 D (k) n } ≥ N 2( 3 k=0 D (k) n ) -( 3 k=0 D (k) n ) 2 if Tr{ 3 k=0 D (k) n } < N D 0 = ϵ max I - 3 k=0 F (k) / (ϵ max -ϵ min ) (3.54)
Considering the trace, Tr{(D n + D (1) n + D (2) n + D (3) n )} = Tr{D} + Tr{D (1) n } + Tr{D (2) n } + Tr{D (3) n } = N + δN (1) + δN (2) + δN (3) If we start from a well-conditioned initial guess, such as: δN (k) < τ, ∀k ≥ 1, where τ is some threshold parameter (typically about 10 -3 ), we may expect the perturbed density matrix to preserve the trace, that is, Tr{(D n + D (1) n + D (2)

n + D (3) n )} ≃ Tr{D n } (3.56)
This constraint means that the perturbation does not create nor annihilate particles. By developing and assembling terms by perturbation order, Eq. (3.53) can be recast in the following form

Tr{D n } ≥ N                D n+1 = D 2 n D (1) n+1 = D n , D (1) n D (2) n+1 = D (1) n 2 + D n , D (2) n D (3) n+1 = D n , D (3) n + D (1) n , D (2) n (3.57)
Tr{D n } < N                  D n+1 = 2D n -D 2 n D (1) n+1 = 2D (1) n -D n , D (1) n D (2) n+1 = 2D (2) n -D (1) n 2 + D n , D (2) n D (3) n+1 = 2D (3) n -D n , D (3) n + D (1) n , D (2) n (3.58)

90

Density matrix perturbation theory

The initial guesses are defined accordingly,

D = (ϵ max I -F )/(ϵ max -ϵ min )
D (1) = -F (1) /(ϵ max -ϵ min )

D (2) = -F (2) /(ϵ max -ϵ min ) D (3) = -F (3) /(ϵ max -ϵ min )
We can even generalize the perturbed TCP at any order k (≥ 0)

D (k) n+1 =        k l=0 D (l) n D (k-l) n if Tr{D n } ≥ N 2D (k) n -k l=0 D (l) n D (k-l) n if Tr{D n } < N (3.60)
with the initial guess as

D (k) = -F (k) /(ϵ max -ϵ min ) (3.61)
In principle, the perturbed projection can be derived from any purification approach. In this work, we have combined this method with the hole-particle canonical purification (HPCP).

Perturbed projection by hole-particle canonical purification

Regarding its polynomials, the TCP is the simplest method, compared to other purifications. However, let us recall that, using the TCP the density matrix trace reaches the correct value only at convergence [cf. Section 2.2.4], whereas the HPCP maintains the N -representability conditions throughout the purification process. As a result, our objective is to compare the performances of the two polynomials in terms of efficiency and reliability. From Eqs. (2.29), (2.31) and (2.32)], the HPCP perturbed projection equations are given below up to the 3rd-order.

0th order

D n+1 = (1 -2c n )D n + 2(1 + c n )D 2 n -2D 3 n D = αD min + (1 -α)D max D min = λ o (µI -F ) + θI, D max = λ q (µI -F ) + θI
n+1 = (1 -2c n )D (1) n + 2(1 + c n ) D n , D (1) n -2 D n D (1) n D n + D 2 n , D (1)
n D (1) = αD (1) min + (1 -α)D (1) max D

(1) min = -λ o F (1) , D (1) max = -λ q F (1) (3.63) 2nd order

D (2) n+1 = (1 -2c n )D (2) n + 2(1 + c n ) D (1) n 2 + D n , D (2) n -2 D n D (2) n D n + D (1) n D n D (1) n + D 2 n , D (2) n + D n , D (1) n 2 D (2) = αD (2) min + (1 -α)D (2) max D (2) min = -λ o F (2) , D (2) max = -λ q F (2) (3.64) 3rd order D (3) n+1 = (1 -2c n )D (3) n + 2(1 + c n ) D n , D (3) n + D (1) n , D (2) n -2 D (1) n 3 + D 2 n , D (3) n + D n , D (1) n , D (2) n -2 D n D (3) n D n + D (1) n D n D (2) n + D (2) n D n D (1) n D (3) = αD (3) min + (1 -α)D (3) max D (3) min = -λ o F (3) , D (3) max = -λ q F (3) (3.65)
Finally, the generalization of the perturbed HPCP at any order k (k ≥ 1) reads

D (k) n+1 = (1 -2c n )D (k) n + 2(1 + c n ) k l=0 D (l) n D (k-l) n -2 k l,j=0 D (l) n D (j) n D (k-l-j) n D (k) = αD (k) min + (1 -α)D (k) max D (k) min = -λ o F (k) , D (k) max = -λ q F (k) (3.66a) (3.66b) (3.66c)
We emphasize that, the constants c n in Eq. (2.18) and µ in Eq. (2.32c) are not expected to change since the constraint is applied on the trace of the unperturbed density matrix. Compared to the CG-based minimization [cf. Section 3.3], the generalized idempotency constraints (3.9e) constitute the kernel of the polynomials for the perturbed projections.

The SCF conditions of Eq. (3.10e) are used in this case to accelerate the perturbed projection.

Derivative of direct inversion of the iterative subspace

The simultaneous SCF computation for all the perturbed density matrices (from the 1st up to 3rd order) can be accelerated by means of the D-DIIS [START_REF] Weber | Improved coupled perturbed hartree-fock and kohn-sham convergence acceleration[END_REF] (derivative direct inversion in the Iterative subspace), which is an extension of the DIIS method presented in Chapter 1 and used for the calculation of the non-perturbed density matrix [cf. Section 1.6.2]. For perturbation orders k ≥ 1, we have implemented the D-DIIS inside the perturbed projection algorithms. The analogue of the DIIS update [cf. Eq. (1.96)] for the derivatives of the Fock matrix is given by

F (k) n := n i=n-m c (k) i F (k) i (3.67)
The coefficients {c

(k)
i } of the linear combination are obtained by minimizing the norm of the error vectors, {e

i }, defined as the commutators between perturbed density matrices and the corresponding Fock matrices [cf. Section 1.6.2]. For the D-DIIS algorithm [START_REF] Weber | Improved coupled perturbed hartree-fock and kohn-sham convergence acceleration[END_REF], these error vectors read e (1)

i := F i , D (1) i + F (1) i , D i e (2) i := F i , D (2) i + F (1) i , D (1) i + F (2) i , D i (3.68a) (3.68b)
which can be recognized as the SCF conditions of Eq. (3.10a) and (3.10b), respectively. As a result, from Eq. (3.10e), the error vector can be generalized at any kth-order, following:

e (k) i := k l=0 F (l) i , D (k-l) i (3.69)
The D-DIIS procedure is identical to the DIIS, the only difference being the definition of the error vectors. The minimization of the norm of the error is performed following the constraint of normalization of the coefficients [cf. Section 1.6.2, Eq. (1.98)], such as min

   f (D-)DIIS , n i=n-m c (k) i = 1    (3.70) with f D-DIIS (c (k) n-m , ..., c (k) n ) := n l,p=n-m c (k) l c (k) p (e (k) l • e (k) p ) (3.71)

Discussions

93

The solution to the problem (3.70) is given by the Euler-Lagrange equation (1.103), which corresponds to a system of (m + 1) linear equations, corresponding to the m coefficients {c (k) i } to be determined, including the Lagrange multiplier λ.

  B (k) 1 t 1 0     c (k) λ   =   0 1  
with: 0 = (0, ..., 0) and 1 = (1, ..., 1)

c (k) = (c (k) n-m , ..., c (k) n ) B (k) lp = (e (k) l • e (k) p ) (3.72) (3.73) (3.74) (3.75)
The pseudo-algorithm using the DIIS/D-DIIS is presented in the Appendix A.

Discussions

In this Chapter we have derived three density matrix-based perturbation methods: (i) the standard AO-CPSCF based on the diagonalization, (ii) the conjugate-gradient based method CG-CPSCF, and (iii) the perturbed projections TC2-CPSCF and HPCP-CPSCF. For all the methods, the density matrix sparsity can be controlled by means of sparse matrix algebra [Section 2.3.2] and a truncation scheme [Section 2.3.1]. This allows, in principle, to perform linear scaling CPSCF calculation. However, we noted that the linear scaling regime already achieved was based only on the numerical truncation. [START_REF] Weber | Higher-order response in o (n) by perturbed projection[END_REF][START_REF] Weber | Ab initio linear scaling response theory: Electric polarizability by perturbed projection[END_REF][START_REF] Kussmann | Linear-scaling method for calculating nuclear magnetic resonance chemical shifts using gauge-including atomic orbitals within hartree-fock and density-functional theory[END_REF][START_REF] Kussmann | A density matrix-based method for the linear-scaling calculation of dynamic second-and third-order properties at the hartree-fock and kohn-sham density functional theory levels[END_REF]. In Table {3.1} is reported the number of matrix multiplication (MM) involved at each CG-Method 1st order 2nd order 3rd order kth-order

AO-CPSCF 5 8 14 5 + 1.5 × k!, k ≥ 2 CG-CPSCF 4 7 11 5 + k! + (6), k ≥ 2 TC2-CPSCF 1 2 4 2 k-1 , k ≥ 1 HPCP-CPSCF 4 8 16 2 k+1 , k ≥ 1
Table 3.1 Number of matrix multiplication for the density matrix methods at different perturbation orders.

step or perturbed projection for all the methods, at any kth-order. For the AO-CPSCF and CG-CPSCF, this number includes the iterative and non-iterative parts. The number given in parenthesis for the CG-CPSCF at the kth-order correponds to the number of MM during the line search of the conjugate gradient. Using these methods, the k order perturbed density matrix is computed and will be used for the evaluation of molecular properties, such as the static non-linear optical properties.

Chapter 4

Applications to non-linear optical properties of π-conjugated systems

Non linear optical properties

The Chapter 3 has discussed the density matrix perturbation methods. However, this discussion was rather general in regards of perturbation to be considered. The molecular properties, as measured by spectroscopy, are the observable responses of the molecular system to an external perturbation. In other terms, the considered perturbartion allows to probe specific properties of the molecular system. We focus here on the response(s) induced by a pertubative external static electric field(s), which define the static optical properties of the molecular system.

Perturbed energy expression for the PPP model

On considering the molecular system in the static electric field ⃗ E , the interaction between the system and the field is pictured by the electric molecular dipole. [START_REF] Schwartz | Principles of electrodynamics[END_REF]1] The classical molecular dipole is defined from the sum of the nuclear contributions, added to the sum over the ponctual electronic charge times the electron position operator. Since in the present formalism nuclei position are not quantized, and we are dealing with a minimal HF-PPP approach [cf. Section 1.4], only the contributions from the π-electrons are relevant. The perturbation operator describing the coupling of the electric field (λ := ⃗ E ) with the electronic dipole ⃗ p, is defined by

∆λ := -⃗ p • ⃗ E (4.1)
with

⃗ p := - Ne i r i (4.2)
Since ∆λ is a one-electron operator, when added to the electronic Hamiltonian [cf.

Eq. (1.34)], its contributions appear only within the one-electron contribution of Eq. (1.60b), according to

h λ µν := dr 1 ϕ * µ (r 1 ) - 1 2 ∇ 2 1 -r 1 • ⃗ E - M A=1 Z A r 1A ϕ ν (r 1 ) (4.3)
The latter can be decomposed into two terms:

h λ µν = h µν + ∆ λ µν (4.4)

Non linear optical properties

97

where h is the original non-perturbed one-electron matrix and ∆ λ is the dipole-electric field coupling matrix, whose elements are given by:

∆ λ µν = ⟨µ|r • ⃗ E |ν⟩ = a∈{x,y,z} E a ⟨µ|a|ν⟩ (4.5)
In the equation above, E a is the component of the electric field along the cartesian direction a. The additional term in the one-electron core hamiltonian due to the presence of the electric field is incorporated during the whole SCF procedure. Consequently, the quantities such as the density and the Fock matrices along with the energies are necessarily modified. As a result, we may express the perturbed Fock matrix and electronic energy by

F λ = h λ + G (D λ ) E λ = Tr{D λ (h λ + F (D λ ))} (4.6a) (4.6b)
where we made explicit the dependence of the two-electron contribution to the perturbed density matrix D λ . Referring to Chapter 1, the expectation value of an operator is the trace of the product between this operator and the density matrix. As a result, for each cartesian component a ∈ {x, y, z}, we have:

⟨p (a) ⟩ := -Tr{Dh (a) } (4.7)
where h (a) is the dipole moment matrix of elements ⟨µ|a|ν⟩. Within the HF-PPP model, it nicely simplifies to

h (a) = ⟨a⟩δ µν (4.8)
This corresponds to a diagonal matrix, with for elements, the position vector component along the direction a.

Energy and response expansions

The conventional expansion [START_REF] Sp | Frequency dependent nonlinear optical properties of molecules: formulation and implementation in the hondo program[END_REF] of the energy for a system perturbed by an external electric field ⃗ E is given by

E( ⃗ E ) = E - a µ a E a - 1 2! a,b α ab E a E b + - 1 3! a,b,c β abc E a E b E c - 1 4! a,b,c,d γ abcd E a E b E c E d -... ( 4.9)
where (a, b, c, d) ∈ {x, y, z}, and E is the non perturbed total electronic energy. The response tensors µ (1st rank ≡ vector), α (2nd rank), β (3rd rank) and γ (4th rank) refer to the dipole moment, the polarizability, the first hyperpolarizability and the second hyperpolarizability, respectively. The density matrix and the Fock matrix are also expanded in terms of the electric field, [START_REF] Sp | Frequency dependent nonlinear optical properties of molecules: formulation and implementation in the hondo program[END_REF] according to

D( ⃗ E ) = D + a D (a) E a + 1 2! a,b D (ab) E a E b + + 1 3! a,b,c D (abc) E a E b E c + 1 4! a,b,c,d D (abcd) E a E b E c E d + ... (4.10) F ( ⃗ E ) = F + a F (a) E a + 1 2! a,b F (ab) E a E b + + 1 3! a,b,c F (abc) E a E b E c + 1 4! a,b,c,d F (abcd) E a E b E c E d + ... ( 4.11)
where (F (a) , F (ab) , ..., F (k) ) are the perturbed Fock matrices and correspond to kth derivative of Eq. (4.6a), which are defined, at any order, by

F (k) =    h (k) + 2J[D (k) ] -K(D (k) ), k = 1 (= a) 2J[D (k) ] -K(D (k) )
k > 1 (= ab, abc, abcd, ...) (4.12)

where h (1) is the dipole moment matrix as defined above. J and K are respectively the Coulomb matrix and the exchange matrix, defining the bi-electronic term G and depending on the perturbed density matrix D (k) . On taking the first derivative of Eq. (4.9) with respect to the electric field components, leads to the dipole moment.

p (a) ( ⃗ E ) = -µ a - 1 2! b α ab E b + - 1 3! b,c β abc E b E c - 1 4! b,c,d γ abcd E b E c E d -... ( 4.13)
Inserting firstly Eq. (4.10) into (4.7), then comparing the resulting equations with Eq. (4.13), leads to the following definitions for the response tensors,

µ a = Tr{h (a) D} α ab = Tr{h (a) D (b) } β abc = Tr{h (a) D (bc) } γ abcd = Tr{h (a) D (bcd) } (4.14a) (4.14b) (4.14c) (4.14d)
That system gives the basic definition of the response tensors (up to 4th order) for a molecular system within a static electric field. Each tensor is defined by a perturbed density matrix, as already outlined for the density matrix perturbation methods discussed in Chapter 3. The alternative to compute these quantities is based on differentiating the electronic energy E with respect to the applied field components at the zero field limit. Formally this writes:

µ a = - ∂E( ⃗ E ) ∂E a ⃗ E =0 , α ab = - ∂ 2 E( ⃗ E ) ∂E a ∂E b ⃗ E =0 , β abc = - ∂ 3 E( ⃗ E ) ∂E a ∂E b ∂E c ⃗ E =0 , γ abcd = - ∂ 4 E( ⃗ E ) ∂E a ∂E b ∂E c ∂E d ⃗ E =0 (4.15)
The Eq. (4.1.2) has actually presented two definitions of the response tensors for a molecular system under a static electric field: (i) the analytic definition of Eq. (4.14) using the density matrix perturbation methods and, (ii) the numerical definition of Eq. (4.15) where the response tensors are evaluated by a numerical differentation using the finite field difference method (FFD). Detailed description of the FFD applied to the calculation of static NLO properties is presented in Appendix E. It is worth to mention that the idempotency relation for the density matrix in a FFD is not verified. As a result, the perturbed electronic energy can not be calculated using a purification or a minimization method. In the algorithm that we have implemented for the finite difference method [Appendix E], the perturbed density matrix is calculated by means of the diagonalization.

Outline of the implementation

Before presenting the results of our calculations of optical properties, this section oultines our implementation. We have modified a pre-existing code, [START_REF] Sony | A general purpose fortran 90 electronic structure program for conjugated systems using pariser-parr-pople model[END_REF] written in Fortran90. The Intel Math Kernel library [START_REF] Wang | Intel math kernel library[END_REF] (MKL) routines were used for the diagonalization and CSC algebra. We note that any other package/routine can be easily adapted to the code.

The various steps for the computation of the unperturbed density matrix is outlined in Figure {4.1}. Given the atomic coordinates of the system, the distance matrix is calculated up to a predefined radial cutoff. From the chosen PPP parameterization, the hopping terms t µν and the PPP parameters Γ µν are evaluated accordingly. In Figure {4.1}, D 0 is the converged density matrix obtained from a TB calculation, which is used as the starting guess for the SCF calculation. In TB or PPP calculation, the density matrix is calculated using the diagonalization (Diag) or the O(M ) methods (Min). The O(M ) methods include the density matrix purifications and minimizations, combined with the numerical or radial truncations.

Set initial C 0 Fock matrix F = F [C n ] F C = C E Energies E n , Coefficients C n Diag. C n vs. C n+1 C ∞ Convergence n = n + 1 No Yes (a) Set initial D 0 Fock matrix F = F [D n ] F D = DF Density matrix D n Min. D n vs. D n+1 D ∞ Convergence n = n + 1 No Yes (b)

Perturbed dense matrix calculation

In this Section we aim to determine the perturbed density matrices D (1) , D (2) and D (3) for the evaluation of the static optical properties. For that purpose, we have considered a benchmark of 2D π-conjugated systems. We will first compare the reliability of the various density matrix pertubation theory (DMPT) methods, along with the FFD the 4 different methods. The responses using the exact method (AO-CPSCF) and the absolute errors for the same responses with respect to the AO-CPSCF (∆(CG-CPSCF), ∆(TC2-CPSCF) and ∆(HPCP-CPSCF) ) for the three other methods are reported. In both Tables, the error is evaluated with respect to the system size (number of cells) at each perturbation order. The perturbed density matrices D (x) , D (xx) and D (xxx) are explicitly calculated using the protocol in Figure {4.3}, then the responses α xx , β xxx and γ xxxx are deduced from Eq. (4.14). The FFD was excluded from our investigations because of the SCF instabilities. These issues might be related to unconsistencies between the size of the system with respect to the strength of the electric field.

At first glance, the results of Tables {4.4} and {4.5} are quite surprising. For both sets of polymers, the numerical accuracy is dramatically reduced for the CG-CPSCF and TC2-CPSCF, as the perturbation order and the size are increased (blue color in the while the HPCP-CPSCF always conserves a remarkable numerical accuracy, whatever the order of perturbation. In order to understand these results, we have probed the density matrix during the SCF iterations, for each method at all the orders. Figures {4.6} and {4.7} display the trace and the idempotency characters of the density matrix, respectively, the latter being evaluated by Basically, the trace of the unperturbed density matrix must correspond to the number of occupied states. On the other hand, at any perturbation order, the traces of the perturbed density matrices must be zero. For the HPCP-CPSCF, from first to third order, the trace of the perturbed density matrices is always zero during the iterations. This should be compared to the CG-CPSCF and TC2-CPSCF approaches, where the trace oscillates, to eventually reach non-zero values. The same behaviour is also observed for the idempotency in Figure {4.7}. In view of these results, the conservation of the numerical accuracy especially at third order for HPCP-CPSCF, is likely due to the conservation of the trace of the unperturbed density matrix D at each iteration step. Even though a good preconditioning may prevent convergence instabilities, the fact that the TC2-CPSCF and CG-CPSCF do not enforce N -representability conditions, they can not avoid any departure from the physical requirements. However, it is worthwhile to note that, in terms of error percentage CG-CPSCF and TC2-CPSCF remain valid. For example, for the largest error (PPV+ of Set 1, with n=32), %error ∼ 0.1% which is insignificant. Using the two forms of matrix illustration already used in Section 2.4.2, the Figures {4.9} and {4.10}, display the profile of the density matrix obtained at the end of the SCF procedure. As the order is increasing, the number nnz of significant elements in the density matrix is also increasing, whereas their magnitude decreases. (3) with nnz = 67% Fig. 4.9 First form of illustration for the converged density matrices using the HPCP, from zero to third order. nnz is the number of non zero elements at 10 -3 . This example is for PPV+ of Set 1 with n = 28. D (1) = D (x) , D (2) = D (xx) , D (3) = D (xxx) (d) D (3) with nnz = 67% Fig. 4.10 Second form of illustration for the converged density matrices using the HPCP, from zero to third order. nnz is the number of non zero elements at 10 -3 . This example is for PPV+ of Set 1 with n = 28. D (1) = D (x) , D (2) = D (xx) , D (3) = D (xxx) .

∆ Idemp =∥ D (k) - k l=0 D (l) D (k-l) ∥ (4.

Applications to non-linear optical properties of π-conjugated systems

From Figures {4.7}, we note the number of iterations for both TC2-CPSCF and HPCP-CPSCF at third order is large compared at first and second orders. This result is explained by the norm of the error vector [cf. Eq. (3.69)] represented in Figure {4.8}. The value of ∼ 10 3 at the early steps of the 3rd-order calculation compared to ∼ 1 at first and second order demonstrates that we are beyond the domain of applicability of the D-DIIS [START_REF] Garza | Comparison of self-consistent field convergence acceleration techniques[END_REF][START_REF] Konstantin N Kudin | A black-box selfconsistent field convergence algorithm: One step closer[END_REF]. As a result, in our implementation, we found a convergence instability for TC2-CPSCF and HPCP-CPSCF at third order. However, as shown on Figure {4.11}, the number of iterations obtained for lower orders and for each of the density matrix perturbation methods, shows clearly that the perturbed projections are the most efficient in terms of SCF iterations. For each density matrix perturbation method, the calculation time is pictured along with its fit.

scales as a power of 4 compared to the cubic scaling of the CG-CPSCF, TC2-CPSCF and HPCP-CPSCF. This demonstrates that the density matrix methods are already more efficient than the diagonalization. In Figure {4.13} is presented the number of iterations obtained for different values of numerical threshold τ (for a fixed size of polymer). We observe that the perturbed projections TC2-CPSCF and HPCP-CPSCF lead to better performances compared to the CG-CPSCF. However, the CG-CPSCF is more stable with a regular number of iterations whatever is the value of τ . The calculation time and the convergence of the

Conclusions and outlooks

In this thesis, we have confirmed that solving electronic structure without relying on the resolution of the Schrödinger equation is a promising approach which can be extended to the perturbation theory. As a matter of fact, the alternative solution deriving from the Liouville-von Neumann equation is the kernel of this work. For single-determinant approximation, the one-particle density matrix is necessary and sufficient to access all the electronic properties of the system and allows to bypass the computational demanding task related to the eigenvalue problem resolution, ie. the diagonalization. In a first part, from the description of the general framework of the density matrix minimization and polynomial expansion, we have proposed a canonical density matrix purification which respects the N -representability constraints. We have emphasized that this purification method is self-consistent, in the sense that, it does not rely on heuristic adjustement of the polynomial during the iterative process. From numerical experiments, the purification polynomial has shown good performances with respect to the other schemes, although its efficiency degrades for the pathological cases. Furthermore, when combined with sparse-matrix algebra to reach the linear-scaling regime, this new purification method is the second most efficient in terms of CPU time. However, as the other density matrix purifications, our variant approach presents the same instability symptoms when using the radial truncation. As a solution to this issue, we have proposed to relax the radial truncation when close to critical points. It is important to specify that this solution is not fully satisfactory. The problem of the radial truncation to the density matrix purifications truly deserves more attention.

In a second part, assuming an orthogonal basis, a detailed development of the density matrix perturbation theory has been presented. One of the presented perturbative methods is related to our new purification variant, corresponding to a new canonical and non heuristic density matrix perturbation method. Comparisons with other perturbed purification and minimization methods have been performed. A detailed analysis of the results has revealed that our method is more robust with a remarkable numerical accuracy, despite its important number of matrix multiplications. The new canonical density matrix perturbation method developed in this thesis is very promising thanks to the explicit consideration of the N -representability properties. The different density-matrix based pertubation methods discussed in this manuscript have been implemented in a code based on the Pariser-Parr-Pople Hartree-Fock model. The most likely next step is to extend this work to a general code including explicit non-orthogonal basis sets where the filling factor significantly deviates from 1/2. The idea will be then to investigate in more details the performance of the method within the linear-scaling regime, especially the influence of the sparsity on the accuracy of the perturbed quantities. Application to dynamic response calculation through the resolution of the time-dependent Liouville-von Neumann equation can also be envisaged.

From our very recent works, we finally noted that the N -representability properties for the density matrix can also be applied to the energy functional minimization. We found that it is possible to enforce these properties during the minimization, that is, the trace is conserved and the density matrix eigenvalues lied in the range 0 to 1, by using a suitably preconditioned conjugate gradient algorithm. Solving this point would be benificial for curing deficiencies observed in energy functional based density matrix perturbation theory.

Derivative direct inversion of iterative subspace

1: ! Initialization 2: D (k) 0 , ∥ e (k) 0 ∥, n = 0 3: ! Iterations 4: while ∥ D n+1 -D n ∥> tolerance do 5: n = n + 1 6:
Build the density matrix D (k) n and Fock matrix F (k) n .

7:

Compute the error matrix e (k) n using Eq. (3.69).

8:

Store F (k)

n and e (k) n .

9:

After ten iterations around, start the DIIS/D-DIIS extrapolations:

10: if ∥ e (k) n ∥ ≲ 10 -3 ∥ e (k) 0 ∥ 11:
Keep m (6 to 8) latest error matrix e (k) n to assemble B (k) using Eq. (3.75).

12:

Resolve the c

Ω LNV (D n ) = Tr{ F (3D 2 n -2D 3 n )} 12: ∇Ω LNV (D n ) = 3(D n F + F D n ) -2(D 2 n F + D n F D n + F D 2 n ) 13: b n = -Tr{H n G n } 14: c n = Tr{3H 2 n F -2 H 2 n D n F + H n D n H n F + D n H 2 n F } 15: d n = -2Tr{H 3 n F } 16: Minimal root of (b n + 2c n λ n + 3d n λ 2 n ) = 0 17: D n+1 = D n + λ n H n 18: G n+1 = -∇Ω LNV (D n+1 ) 19: γ n =        G n+1 Gn GnGn : (FR) or (G n+1 -Gn)G n+1 GnGn : (PR) 20
µ n+1 = µ n + α XS N -Ñe 12: Fn = F n -µ n I 13: Ω XS (D n ) = Tr{ Fn (3D 2 n -2D 3 n )} 14: ∇Ω XS (D n ) = 3(D n Fn + Fn D n ) -2(D 2 n Fn + D n Fn D n + Fn D 2 n ) 15: b n = -Tr{H n G n } 16: c n = Tr{3H 2 n Fn -2 H 2 n D n Fn + H n D n H n Fn + D n H 2 n Fn } 17: d n = -2Tr{H 3 n F } 18:
Minimal root of (b n + 2c n λ n + 3d n λ 2 n ) = 0

19:

D n+1 = D n + λ n H n 20: G n+1 = -∇Ω XS (D n+1 ) 21: γ n =        G n+1 Gn GnGn : (FR) or (G n+1 -Gn)G n+1 GnGn : (PR)
22: c n = 0.5

H n+1 = G n+1 + γ n H n 23: δN = | Ñe -Tr{D n }| 24
µ n = Tr{2D n F D n -F (3D n -2D 2 n ) -(3D n -2D 2 n )F }/M 11: Ω MS (D n ) = Tr{F (3D 2 n -2D 3 n )} + µ n (Tr{D n } -N ) 12: ∇Ω MS (D n ) = 3(DF n + F D n ) -2(D 2 n F + D n F D n + F D 2 n ) + µ n I 13: b n = -Tr{H n G n } 14: c n = Tr{3H 2 n Fn -2 H 2 n D n Fn + H n D n H n Fn + D n H 2 n Fn } 15: d n = -2Tr{H 3 n F } 16: Minimal root of (b n + 2c n λ n + 3d n λ 2 n ) = 0 17: D n+1 = D n + λ n H n 18: G n+1 = -∇Ω MS (D n+1 ) 19: γ n =        G n+1 Gn GnGn : (FR) or (G n+1 -Gn)G n+1 GnGn : (PR) 20
D n+1 = [(1 + c n )D 2 n -D 3 n ] /c n 18: else 19: D n+1 = [(1 -2c n )D n + (1 + c n )D 2 n -D 3 n ] /(1 -c n ) 20
F (D n ) = D 2 n (4D n -3D 2 n ) 11: G(D n ) = D 2 n (1 -D n ) 2 12: if ∥ Tr{G(D n )} ∥< 10 -4 then 13:
γ n = 3.0 c n = 0.5

3: ϵ min , ϵ max ← F 4: μ = Tr{F }/M , θ = N/M , λ 1 = N M (ϵmax-μ) , λ 2 = M -N M (μ-ϵ min ) , 0 < α < 1 5: λ o = min {λ 1 , λ 2 }, λ q = max {λ 1 , λ 2 } 6: D min = λ o (μI -F ) + θI, D max = λ q (μI -F ) +
c n = Tr{D 2 n -D 3 n }/Tr{D n -D 2 n } 16:
end if 17: b: input which is the RHS.

D n+1 = (1 -2c n )D n + 2(1 + c n )D 2 n -2D 3
F (k) = F (k) D (k)

12:

APROD: external routine required by SYMMLQ and which supplies the matrixvector product Ax, so the LHS. where A is an N by N symmetric matrix and b is a given vector. The matrix A is not required to be positive definite. (If A is known to be definite, the method of conjugate gradients might be preferred, since it will require about the same number of iterations as SYMMLQ but slightly less work per iteration.) The matrix A is intended to be large and sparse. It is accessed by means of a subroutine call of the form call APROD(N, x, y) which must return the product y = Ax for any given vector x. More generally, SYMMLQ is designed to solve the system (A -SHIFT I N )x = b where SHIFT is a specified scalar value. If SHIFT and b are suitably chosen, the computed vector x may approximate an (unnormalized) eigenvector of A, as in the methods of inverse iteration and/or Rayleigh-quotient iteration. Again, the matrix (A -SHIFT I N ) need not be positive definite. The work per iteration is very slightly less if SHIFT = 0.

A further option is that of preconditioning, which may reduce the number of iterations required. If M = CC t is a positive definite matrix that is known to approximate (A -SHIFT I N ) in some sense, and if systems of the form M y = x can be solved efficiently, the parameters PRECON and MSOLVE may be used (see below). When PRECON = .true., SYMMLQ will implicitly solve the system of equations No iterations were performed.

1: ∥r∥ appears to be less than the value RTOL ∥ Ā∥ ∥x∥.

The solution in x should be acceptable. This means that the residual is as small as seems reasonable on this machine. x may not contain an acceptable solution.

5: The iteration limit was reached before any of the previous criteria were satisfied.

6:

The matrix defined by APROD does not appear to be symmetric.

For certain vectors y = A v and r = A y, the products y t y and r t v differ significantly.

7:

The matrix defined by MSOLVE does not appear to be symmetric.

For vectors satisfying M y = v and M r = y, the products y t y and r t v differ significantly.

8: An inner product of the form x t M -1 x was not positive, so the preconditioning matrix M does not appear to be positive definite.

If ISTOP ≥ 5, the final x may not be an YNORM is an estimate of ∥x∥.

The principle of finite differences at a single point. For example, f (x + ξ) = f (x) + ξf (1)

(x) + ξ 2 2! f (2) + ... = ∞ k=0 ξ k k! f (k) (x) f (x -ξ) = f (x) -ξf (1) (x) + ξ 2 2! f (2) + ... = ∞ k=0 (-1) k ξ k k! f (k) (x) (E.5a) (E.5b)
where f (k) denotes the kth order derivative of f . By Rearranging for instance Eq. (E.5a), such that

f (x + ξ) -f (x) ξ -f (1) (x) = ξ 2! f (2) (x) + ξ 2 3! f (3) (x) + ... Truncation Error (E.6)
one can observe that the FED in Eq. (E.2) corresponds to a Taylor series truncated after the second term. The rhs of Eq. (E.6) is the error in terminating the series and is referred to as the truncation error (TE). [START_REF] Gordon | Numerical solution of partial differential equations: finite difference methods[END_REF][START_REF] John | Finite difference schemes and partial differential equations[END_REF][START_REF] Randall | Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems[END_REF] The TE can be defined as the difference between the partial derivative and its finite difference representation. The finite difference representation described above evaluates the function to be derived in a 2 points approximation. [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF]. In order to maximize the numerical accuracy of a finite difference representation, one can evaluate the function with a higher number of points. That involves a higher order Taylor series. For example, let us first set In Eq. (E.9) and Eq. (E.10), the function is now evaluated in 3 points, while in Eq. (E.11) the function is now evaluated in 4 points. The relation (E.8) means that the Taylor series evaluated in many different points, are combined in some appropriate way so that one can generalize a relation between the d th order derivative function and its finite difference representation by [START_REF] Eberly | Derivative approximation by finite differences[END_REF] where p (> 0) is the integer order of the TE, selected as desired. c m are the coefficients of the finite difference represensation. In order to determine c m , one uses the Taylor series for f (x + mξ), which is

f (x + 2ξ) = ∞ k=0 (2ξ) k k! f (k) (x) f (x -2ξ) = ∞ k=0 (-1) k (2ξ) k k! f (k) (x) (E.
ξ d d! f (d) (x) =
f (x + mξ) = ∞ k=0 m k ξ k k! f (k) (x) (E.13)
Introducing Eq. (E.13) into Eq. (E.12) yields The Table {E.1} gives the number of coefficients c m , so the different points m in the finite difference representation from the derivative order d and the TE order p, for the differente approximation. In other words, Table {E.1} gives a relationship between the number of terms in the finite difference representation and the number of terms in the Taylor series. For this reason, in Eq. (E.15), the sum over k or the number of terms in the Taylor series is no more infinite. Note that in Table {E.1}, for CED approximation, d + p is necessarily an odd number, while p can be chosen to be even regardless of the parity of d.

ξ d d! f (d) (x) = mmx m=mmn c m ∞ k=0 m k ξ k k! f (k) (x) + O(ξ d+p ) = ∞ k=0 mmx m=mmn m k c m ξ k k! f (k) (x) + O(ξ d+p ) = d+p-1 k=0 mmx m=mmn m k c m ξ k k! f (k) (x) + O(ξ d+p ) (E.
In order to understand Eq. (E.16) which gives the coefficients c m of the finite difference representation, let us approximate for example f (3) (x) with a FED and a 1st order TE ie. O(ξ), so d = 3 and p = 1. Using Table {E.1} gives m mn = 0 and m mx = 3. The linear system given by Eq. (E.16) is then

       (0) 0 (1) 0 (2) 0 (3) 0 (0) 1 (1) 1 (2) 1 (3) 1 (0) 2 (1) 2 (2) 2 (3) 2 (0) 3 (1) 3 (2) 3 (3) 3               c 0 c 1 c 2 c 3        =        0 0 0 1        (E.
              c 0 c 1 c 2 c 3        =        0 0 0 1        (E.18)
This equation is easily resolved by hand and its solution is (c 0 , c 1 , c 2 , c 3 ) = (-1, 3, -3, 1)/6. In other words, using Eq. (E.12), f (3) (x) = -f (x) + 3f (x + ξ) -3f (x + 2ξ) + f (x + 3ξ)

ξ 3 + O(ξ) (E.19)
Let us now approximate f (3) (x) with a CED and error O(ξ 2 ). Proceeding in the same way, d = 3 and p = 2, gives m mx = -m mn = 2. The resulting linear system is

         
(-2) 0 (-1) 0 (0) 0 (1) 0 (2) 0 (-2) 1 (-1) 1 (0) 1 (1) 1 (2) 1 (-2) 2 (-1) 2 (0) 2 (1) 2 (2) 2 (-2) 3 (-1) 3 (0) 3 (1) 3 (2) 3 (-2) 4 (-1) 4 (0) 4 (1) 4 (2)

4                     c -2 c -1 c 0 c 1 c 2           =           0 0 0 1 0           (E.
                    c -2 c -1 c 0 c 1 c 2           =           0 0 0 1 0           (E.21)
and has solution (c -2 , c -1 , c 0 , c 1 , c 2 ) = (-1, 2, 0, -2, 1)/12. Finally, the expression for

f (3) (x) is f (3) (x) = -f (x -2ξ) + 2f (x -ξ) -2f (x + ξ) + f (x + 2ξ) 2ξ 3 + O(ξ 2 ) (E.22)

Finite differences of bivariate functions

The CED approximation is enough used for finite field difference methods applied to the response tensors [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF][START_REF] Lu | Semi-empirical calculations of the nonlinear optical properties of polycyclic aromatic compounds[END_REF][START_REF] Zakharov | Molecular nonlinear optical parameters of π-conjugated nonalternant hydrocarbons obtained in semiempirical local coupled-cluster theory[END_REF][START_REF] Li | Electron correlation effects on the nonlinear optical properties of conjugated polyenes[END_REF][START_REF] Sim | Electron correlation effects in hyperpolarizabilities of p-nitroaniline[END_REF][START_REF] Henry A Kurtz | Calculation of the nonlinear optical properties of molecules[END_REF]. For example, we have implemented the finite field differ-

The principle of finite differences ence of Ref. [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF] which uses a CED with 7 points. That involves m = {0, ±1, ±2, ±3}, hence m mx = -m mn = 3. We can now deduce the approximate expressions for the derivatives of the response tensors given in Eq.

                                c -3 c -2 c -1 c 0 c 1 c 2 c 3                 =                                 (E.23)
which approximates the first order derivative of the energy with respect to the electric field as

µ = ∂E( ⃗ E ) ∂E ≈ 1 60ξ (E +3 -9E +2 + 45E +1 -45E -1 + 9E -2 -E -3 ) + O(ξ 6 ) (E.
                                c -3 c -2 c -1 c 0 c 1 c 2 c 3                 =                                 (E.25)
which approximates the second order derivative of the energy with respect to the electric field as

α = ∂ 2 E( ⃗ E ) ∂E 2 ≈ 1 180ξ 2 (2E +3 -27E +2 + 270E +1 -490E 0 + 270E -1 -27E -2 + 2E -3 ) + O(ξ 5 ) (E.
                                c -3 c -2 c -1 c 0 c 1 c 2 c 3                 =                 0 0 0 1 0 0 0                 (E.27)
which approximates the third order derivative of the energy with respect to the electric field as

β = ∂ 3 E( ⃗ E ) ∂E 3 ≈ 1 8ξ 3 (-E +3 + 8E +2 -13E +1 +13E -1 -8E -2 + E -3 ) + O(ξ 4 ) (E.
                                c -3 c -2 c -1 c 0 c 1 c 2 c 3                 =                 0 0 0 0 1 0 0                 (E.29)
which approximates the fourth order derivative of the energy with respect to the electric field as

γ = ∂ 4 E( ⃗ E ) ∂E 4 ≈ 1 6ξ 4 (-E +3 + 12E +2 -39E +1 +56E 0 -39E -1 + 12E -2 -E -3 ) + O(ξ 3 ) (E.30)
In Eq. (E.24), Eq. (E.26), Eq. (E.28) and Eq. (E.30), ξ is the differentation step of the finite field difference representation for the energy derivative with respect to the electric field ⃗ E . The unit of ξ is that of the electric field strength, ξ ∼ 10 -3 au [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF]. And E ±i is

E ±i = E( ⃗ E ± iξ) (E.31)
The principle of finite differences the energy of the system calculated for the strength of the electric field equal to ⃗ E ± iξ with i = {0, 1, 2, 3}. On the other hand, it is important to emphasize that the finite field difference representation of Eq. (E.24), Eq. (E.26), Eq. (E.28) and Eq. (E.30) is for an univariate function. That implies the electric field changes only in one direction, ie. E ±i = E(E x ± iξ) in the x direction. In the case where the electric field changes in two directions x and y [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF], the finite difference representation basically requires two differentiation steps ξ x and ξ y , respectively. Supposing that ξ x and ξ y can be comparable, then we can find a step ξ so that: ξ x = iξ and ξ y = jξ. We may write the energy as

E ±i,±j = E(E x ± iξ, E y ± jξ) (E.32)
This energy is associated to Eq. ( 4.3) which the expression is

h λ µν = h µν + i⟨x⟩ξδ µν + j⟨y⟩ξδ µν (E.33)
where (i, j) = {0, ±1, ±2, ±3}. and ⟨x⟩ and ⟨y⟩ represent the position vector components along the x and y directions. We obtain an energy matrix such as In this matrix, a row(column) means that the field changes in x(y) direction while is fixed in y(x) direction. The fourth row(column) corresponds to the univariate case where the field exits only in x(y). As a result, a derivative of this energy with respect to E x (E y ) requires the seven energies along a row (column) of this matrix. The formulas in Eq. (E.24), Eq. (E.26), Eq. (E.28) and Eq. (E.30) using the energies of the 4 th row (column) give the diagonal components of the tensors along x (y), since the field is applied only in a single direction. For example for the component β yyy , we need to take 165 the formula of Eq. (E.28) with the energies

E i,j =                 E -3,-3 E -2,-3 E -1,-3 E 0,-3 E +1,-3 E +2,-3 E +3,-3 E -3,
                E 0,-3 E 0,-2 E 0,-1 E 0,0 E 0,+1 E 0,+2 E 0,+3                
While for the non-diagonal components of tensors, one has to apply some combinations of formulas (E.24), (E.26), (E.28) and (E.30), as the field is applied in several directions. The value of the field changes in a direction while it is fixed in the other directions. As an example, for

γ xxyy = ∂ 2 ∂E 2 x ∂ 2 E(E x , E y ) ∂E 2 y = ∂ 2 X(E x ) ∂E 2
x we first use Eq. (E.26) for each column, which leads to a row vector

X -3 X -2 X -1 X 0 X +1 X +2 X +3
For the elements of this vector, X +1 for instance means the E y E y -second derivative while the field at E x is fixed to i = +1, and so on. Then, by applying once more Eq. (E.26) to this vector, we finally obtain the non-diagonal component γ xxyy . The below algorithm outines the finite field difference method which evaluates the optical properties. This method is from Ref. [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF] and uses a seven points representation. The energy is calculated at different points. The resulting energies are then used in combinations of sums of energies that define the derivatives for the response tensors.

166

The principle of finite differences 1: ! Data: F , ξ, tolerance. end for 15: end for 16: E(i, j) = Ẽ∞ 17: E(i, j) Eq.(E.24),Eq.(E.26),Eq.(E.28),Eq.(E.30) ----------------------→ µ, α, β, γ Algorithm 12 Finite field difference method for µ, α, β and γ (optical properties).

1. 1

 1 Flow diagram of the self-consistent field processus. . . . . . . . . . . . . . 1.2 Flow diagram of the SCF scheme including the CDA. . . . . . . . . . . . 1.3 Flow diagram of the SCF scheme including the DIIS. . . . . . . . . . . . 1.4 SCF convergence profiles obtained for the carbon nanotube

2. 1 4 . 5

 145 Influence of the fermion temperature on the Fermi-Dirac distribution. . . 2.2 McWeeny purification polynomial P (x) = 3x 2 -2x 3 . . . . . . . . . . . . 2.3 Palser and Manolopoulos polynomials for different ajustement parameter c n . 2.4 Projection polynomials used for the trace-correcting density matrix purifications. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Projection polynomials used of the trace-resetting density matrix purification. 2.6 Scattered eigenvalues from a chunk of the set of test Hamiltonians. . . . 2.7 (a) Color maps displaying the average number of purifications (p) as the function of the filling factor (θ) and energy gap (∆ϵ gap ). Results obtained from the PMCP, HPCP, TRS4 and TC2 methods. (b) Energy convergence profiles with respect to the first 15 iterations for selected values of θ, and (c) the corresponding density matrix trace conservation profiles. . . . . . 2.8 Comparison of four density matrix purifications in terms of matrix multiplications (MMs) for varying filling factors, and two band gap values ∆ϵ gap = 0.5 and ∆ϵ gap = 10 -6 . Heavy lines represent averages over 256 random Hamiltonians and shaded areas are the corresponding standard deviations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Rolling of a graphene sheet to generate a carbon nanotube. . . . . . . . . xx List of figures 2.10 First form of illustration of the sparsity pattern of the density matrix truncated at τ = 10 -8 during the SCF iterations, following the sequence (a) to (d). Results obtained for the CNT (11,5). . . . . . . . . . . . . . . 2.11 Second form of illustration of the sparsity pattern of the density matrix truncated at τ = 10 -8 during the SCF iterations, following the sequence (a) to (d). Results obtained for the CNT (11,5). . . . . . . . . . . . . . . 2.12 Chart for the radial truncation scheme displaying the circle (blue solid line) of radius R c . The largest circle (red dashed line) is of radius equal to a/2 inscribed in the square unit cell. . . . . . . . . . . . . . . . . . . . . 2.13 Convergence of the LNV energy with respect to the exact value obtained from diagonalization (no truncation) as a function of the radial cutoff R c . (a) Convergence profile obtained for the zigzag CNT. (b) Convergence profile obtained for the chiral CNT. . . . . . . . . . . . . . . . . . . . . . 2.14 Convergence of energy during the density matrix (a) minimization (LNV), and (b) purification (HPCP) for the CNT (8,0). A cutoff radius of 15 Å have been used (cf. text for more details). . . . . . . . . . . . . . . . . . 2.15 First form of illustration for the progression of the density matrix truncated at R c = 10 Å during the iterations, throughout the four respective sequences (a), (b), (c) and (d). . . . . . . . . . . . . . . . . . . . . . . . . 2.16 Second form of illustration for the progression of the density matrix truncated at R c = 10 Å during the iterations, throughout the four respective sequences (a), (b), (c) and (d). . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Calculation time as a function of the number of atoms. Linear scaling regime is achieved using the radial truncation. (a) R c = 50 Å. (b) R c = 10 and 50 Å Results were obtained for the replicated CNT (11,5). . . . . . . 2.18 Illustration of the C coefficients matrix for the diagonalization and of the untruncated D density matrix for the minimizations and purifications methods. (a) and (b) are the first form of illustration, while (c) and (d) are the second form of illustration, for the CNT (11,5). . . . . . . . . . . 4.1 Outline of the implementation giving the steps of the unperturbed SCF PPP calculation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 SCF procedure using: (a) the diagonalization [Diag], and (b) the density matrix energy minimizations [Min]. . . . . . . . . . . . . . . . . . . . . . 4.3 CPSCF density matrix perturbation methods. (a) AO-CPSCF and CG-CPSCF, (b) TC2-CPSCF and HPCP-CPSCF. . . . . . . . . . . . . . . . . . . . . 4.4 Benchmark of molecules. . . . . . . . . . . . . . . . . . . . . . . . . . . . List of figures xxi Benchmark of polymers. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Trace of the density matrices during the SCF iterations for the HPCP-CPSCF,

4 . 13

 413 Histogram of number of SCF iterations with respect to the numerical threshold τ . Results are obtained for the four density matrix perturbation methods at first and second orders. The number of cells is fixed at 250 for TPA+ and 62 PPV+ of Set 1. . . . . . . . . . . . . . . . . . . . . . . . . xxii List of figures 4.14 In (a) and (b), representation of the total calculation time as a function of the number of atoms, for TPA+ (Set 1). This representation compares the AO-CPSCF, CG-CPSCF, TC2-CPSCF and HPCP-CPSCF. The density matrix methods are all truncated at τ = 10 -4 and τ = 10 -8 . In (c) and (d), convergence for (hyper)polarizability per atom as a function of number of atoms. The polarizability α xx and first hyperpolarizability β xxx are calculated using the HPCP-CPSCF truncated at τ = 10 -8 for the 4 polymers of Set 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Nomenclature Subscripts α, β, γ, ... indices for atomic orbitals i, j, k, ... indices for molecular orbitals Other Symbols [A, B] Commutator AB -BA {A, B} Anticommutator AB + BA Tr{X} Trace of a matrix X ∥X∥ Frobenuis norm of a matrix or a vector X

  (1.13c) and (1.40), we obtain for the one-particle density operator the following expression D := µ,ν |µ⟩ D µν ⟨ν| , with: D µν :=

  number of unoccupied states, such that: M = N + N . Conversely, one can also define the one-hole density matrix built from the set of unoccupied eigenstates. By analogy with Eq. (1.53), this yields to introduce D = C ŌC † (1.55) 1.3 Restricted Hartree-Fock energy 17

Fig. 1 . 1

 11 Fig. 1.1 Flow diagram of the SCF processus.

Fig. 1 . 2

 12 Fig. 1.2 Flow diagram of the SCF scheme including the CDA.

Fig. 1 . 3

 13 Fig. 1.3 Flow diagram of the SCF scheme including the DIIS.

Fig. 1 .

 1 Fig. 1.4 SCF convergence profiles obtained for the carbon nanotube (11,5) using a simple approach, the CDA and the DIIS optimization. In (a) is displayed the energy minimization during the iterative process. In (b) is represented the energy error during the iterative process using a logarithmic scale.

  4 and Eq. (1.46)].

5 TFig. 2 . 1

 521 Fig. 2.1 Influence of the fermion temperature on the Fermi-Dirac distribution. The inflection point of the distribution is located at (µ = 0, η = 0.5)

Fig. 2 . 2

 22 Fig. 2.2 McWeeny purification polynomial P (x) = 3x 2 -2x 3 . The initial (2 × 2) mixed state, with occupation numbers {η 1 , η 2 } n=0 = {0.3, 0.7}, is purified to eventually reach {η 1 , η 2 } n=4 = {0, 1}.

Fig. 2 . 3

 23 Fig. 2.3 Polynomials of the Palser and Manolopoulos canonical purification.

  . For m = 2, both functions merge into the McWeeny polynomial of Figure{2.2}. Note that the computational ressources, measured in terms of matrix multiplication (MM), are given by the value of m. For instance: for m = 1, each purification requires one MM, for m = 2, two MMs are necessary, and so on.

3 Fig. 2 . 4

 324 Fig. 2.4 Projection polynomials used for the trace-correcting density matrix purifications at different orders.

Fig. 2 . 5

 25 Fig. 2.5 Projection polynomials used of the trace-resetting density matrix purification for different adjustement parameters of γ ∈ [0, 6].

1

  Lionel A. Truflandier, 1,a) Rivo M. Dianzinga, 1 and David R. Bowler 2 Institut des Sciences Moléculaires, Université Bordeaux, CNRS UMR 5255, 351 cours de la Libération, 33405 Talence cedex, France 2 London Centre for Nanotechnology, UCL, 17-19 Gordon St., London WC1H 0AH, United Kingdom; Department of Physics and Astronomy, UCL, Gower St., London WC1E 6BT, United Kingdom; and International Centre for Materials Nanoarchitechtonics (MANA), National Institute for Materials Science (NIMS), 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan (Received 6 January 2016; accepted 22 February 2016; published online 3 March 2016)

FIG. 1 .

 1 FIG. 1. (a) Convergence of the McWeeny Lagrangian and density matrix eigenvalues during the course of the minimization using a test Hamiltonian and an occupation factor θ = 0.10. A grey scale is used to guide the eye during the processes of purification. Each curve is a plot of the function L McW (x;γ n ) computed at each iteration n. The red line corresponds to L McW (x;0) = Ω McW . (b) Convergence of L McW (green circles) and the trace conservation Tr{D n } -N (black dots). (c) Convergence of ∥∇L McW ∥ F (green circles) and ∥D n ∥ F -N (black dots).

  which guarantees Tr{D n+1 } = N and [H , D n+1 ] = 0, ∀n. Added to the preconditioning λ(D 0 ) ∈ [0, 1], the iterative process should approach the (one-particle) ground-state energy E = Tr{H D} variationally. The parameter c [Eq. (7b)] is recognized as the unstable fixed point introduced in Ref. 21, where c ∈ [0, 1]. As a result, the interval [-1 3 , 1 3 ] constitutes the stable variational domain of γ. The variation of the McWeeny Lagrangian function and the DM eigenvalues during the course of the minimization is presented in Fig. 1(a) for a test Hamiltonian with N = 10, M = 100, and a suitably conditioned initial guess (vide infra). The corresponding convergence profiles of L McW and ∥∇L McW ∥ (green circles) are reported on Figs. 1(b) and 1(c), respectively, along with the trace conservation Tr{D n } -N and the DM norm convergence ∥D n ∥ -N (black dots). We may notice first that for γ = 0 (or c = x m = 1 2 ), L McW simplifies to Ω McW . For intermediate states, γ ∈ [-1 3 , 0] ∪ [0, 1 3 ], the symmetry of Ω McW is lost and the shape of L McW (x, γ n ) drives the eigenvalues in the hole (left) or in the particle (right) well. From the grey scale in Fig. 1(a), we observe how γ n influences L McW (along the y-axis) at x p and the abscissa of the second stationary point x m which is free to move in [x p, x p ]. This yields to transform the hole well from a local (n = 0) to a global (n = 15) minimum (or conversely the particule well from a global to a local minimum). At the boundary values γ = {-1 3 , 1 3 }, x p and x m merged to a saddle point in such a way that only one global minimum left at x p . Notice that, for situations where γ [-1 3 , 1 3 ], the saddle point transforms to a maximum and runaway solutions may appear. Nevertheless, as long as D 0 is well conditioned, such kind of critical problem should not be encountered. Figs. 1(b) and 1(c) highlight the minimization mechanism: (i) from iterate n = 0 to 12; γ → 0 + , L McW follows the search direction and decreases monotonically. (ii) At iterate n = 13; γ ≃ 0, L McW is close to the target value but the gradient residual is nonzero. (iii) From n = 14 to 15; γ < 0, the search direction is inverted. (iv) At iterate n = 16, all the eigenvalues are trapped in their respective wells. (iii) From iterate n = 17 to 23, γ → 0 -, we are in the McWeeny regime [Eq. (4)] and L McW eventually reaches the global minimum.

FIG. 2 .

 2 FIG. 2. (a)Color maps displaying the average number of purifications ( p) as the function of the filling factor (θ) and energy gap (∆ϵ gap ). Results obtained from the PMCP and HPCP methods using the initial guess of Eqs. (2)-(11) and (2)-(14) (notated PMCP+ and HPCP+). Each pixel on the maps corresponds to an average over 32 test Hamiltonians. (b) Energy convergence profiles with respect to the first 15 iterations for selected values of θ. (c) Average number of purifications as a function of ln(1/∆ϵ gap ).

)

  On multiplying Eqs. (A1a) and (A3a) by (1 -c n ) and c n , respectively [or multiplying Eqs. (A1b) and (A3b) by c n and (1 -c n )], and adding, we obtain

  and C t C = I 5. Ẽ := shift{E | N, ∆ϵ gap }, such that: ∆ϵ gap = ϵ N +1 -ϵ N 6. H := C ẼC t In step 5, given a filling factor, θ = N/M , and a band gap, ∆ϵ gap , a shift operator is applied to the eigenvalues located at the middle of the band gap in order to verify: ∆ϵ gap = ϵ N +1 -ϵ N . The final dense Hamiltonian matrix is recovered in step 6. Given a fixed energy band gap of 1.0 au and eigenspectrum width of about 6 au, examples of eigenvalues distributions obtained from a chunk of test Hamiltonians are represented in the Figure {2.6}.

Fig. 2 . 6

 26 Fig. 2.6 Scatter plot of the eigenvalues (in blue) from a chunk of the test Hamiltonians. The red, yellow and green pixels correspond to ϵ N , ϵ N +1 and middle of the band gap, respectively. Each panel correponds to a set 32 randomized symmetric matrices for filling factor θ ∈ {0.01, 0.90}.

Fig. 2 .

 2 Fig. 2.7 (a) Color maps displaying the average number of purifications (p) as the function of the filling factor (θ) and energy gap (∆ϵ gap ). Results obtained from the PMCP, HPCP, TRS4 and TC2 methods. (b) Energy convergence profiles with respect to the first 15 iterations for selected values of θ, and (c) the corresponding density matrix trace conservation profiles.

Figure {2. 7 Fig. 2 . 8

 728 Figure {2.7(b,c)} sheds some light on the HPCP trace-preserving property regarding the strong fluctuation of the number of occupied states observed for the trace-correcting method, and to a lesser extent, the trace-resetting mechanism. It is noteworthy that

2. 4

 4 Applications to carbon nanotubes 59 (a) Illustration of the chiral vector for a the graphene sheet. T denotes the tube axis. ⃗ a 1 and ⃗ a 2 are the unit vectors of graphene in real space. (b) Armchair type (p, p). (c) Zigzag type (p, 0). (d) Chiral type (p, q).

Fig. 2 . 9

 29 Fig. 2.9 Rolling of a graphene sheet to generate a carbon nanotube.

Fig. 2 .

 2 Fig. 2.10 First form of illustration of the sparsity pattern of the density matrix truncated at τ = 10 -8 during the SCF iterations, following the sequence (a) to (d). Results obtained for the CNT (11,5).

Fig. 2 .

 2 Fig. 2.11 Second form of illustration of the sparsity pattern of the density matrix truncated at τ = 10 -8 during the SCF iterations, following the sequence (a) to (d). Results obtained for the CNT (11,5).

2. 4

 4 Applications to carbon nanotubes 65 density matrix during the iterative procedure where, according to the legend color bar, the largest matrix elements are located on the diagonal. Starting from a sparse initial guess, 2 we observe how the DM is gradually filled up from step (a) to (d), to eventually reach convergence at nnz = 89%. The alternative representation of Figure{2.11}, where the elements (in absolute value) are sorted in decreasing order (from the top to the bottom), emphasizes on how really sparse are the matrices we are dealing with, and the amplitude of variations of the nnz. Each colored map in Figure{2.11} presents an important decaying constrast from black (the largest elements) to grey (the smallest elements). From (a) to (b), one can see that the nnz region is progressively widening to reach half of the picture. Then, from (c) to (d), the amplitude of those elements -the intensity of the red color-is increasing.

Fig. 2 .

 2 Fig. 2.12 Chart for the radial truncation scheme displaying the circle (blue solid line) of radius R c . The largest circle (red dashed line) is of radius equal to a/2 inscribed in the square unit cell.

Fig. 2 .

 2 Fig. 2.13 Convergence of the LNV energy with respect to the exact value obtained from diagonalization (no truncation) as a function of the radial cutoff R c . (a) Convergence profile obtained for the zigzag CNT. (b) Convergence profile obtained for the chiral CNT.

Fig. 2 .

 2 Fig. 2.14 Convergence of energy during the density matrix (a) minimization (LNV), and (b) purification (HPCP) for the CNT (8,0). A cutoff radius of 15 Å have been used (cf. text for more details).

Fig. 2 .

 2 Fig. 2.15 First form of illustration for the progression of the density matrix truncated at R c = 10 Å during the iterations, throughout the four respective sequences (a), (b), (c) and (d).

Fig. 2 .

 2 Fig. 2.16 Second form of illustration for the progression of the density matrix truncated at R c = 10 Å during the iterations, throughout the four respective sequences (a), (b), (c) and (d).

Fig. 2 .Fig. 2 .

 22 Fig. 2.17 Calculation time as a function of the number of atoms. Linear scaling regime is achieved using the radial truncation. (a) R c = 50 Å. (b) R c = 10 and 50 Å Results were obtained for the replicated CNT (11,5).

3. 1 background 77 3. 1 ( 3 F

 17713 Theoretical Theoretical background Let us consider an unperturbed system of N occupied and N unoccupied molecular states expanded onto a linear combination of M atomic basis functions, such that N + N = M . Considering an orthonormalized set of basis functions [cf. Section 1.4], we recall the generalized constraints on the one-particle density matrix D for the occupied states and the one-hole density matrix D for the unoccupied eigenstates, are given as Idempotency: D 2 = D and D2 = D Trace conservation: Tr{D} = N and Tr{ D} = N Complementarity: D + D = I Orthogonality: D D = 0 Eqs. (3.1a) and (3.1b) are the N -representability conditions for a single Slater determinant. Along with Eq. (3.1), the ground state is guaranteed if the density matrix and the Fock matrix fulfill the SCF conditions F D -DF = 0

  9e) Then, by repeating the perturbation identification on the Liouville-von Neumann equation, that is by inserting Eqs. (3.7a) and (3.7b) into Eq. (3.6a), we obtain the perturbed SCF conditions:

  ov,ij (ϵ i -εj ) -F

  . (3.14) and (3.23) into (3.28), simplifies Eq. (3.27) to:

  14)] = ( DD (1) D)(DD (1) D) = D(D (1) DD (1) ) D [using Eq. (3.1a)] (3.40)

52 ) 3 . 4

 5234 Perturbed projection by trace-correcting purification 89 On introducing the perturbative expansion (3.7) in Eqs. (3.51) and (3.52), the unperturbed TCP equation transforms to

Fig. 4 . 1 4 . 2

 4142 Fig. 4.1 Outline of the implementation giving the steps of the unperturbed SCF PPP calculation.

Fig. 4 .

 4 Fig. 4.2 SCF procedure using: (a) the diagonalization [Diag], and (b) the density matrix energy minimizations [Min].

  Fig. 4.5 Benchmark of polymers.

Fig. 4 . 6

 46 Fig. 4.6 Trace of the density matrices during the SCF iterations for the HPCP-CPSCF, TC2-CPSCF and CG-CPSCF, at zero (D), first (D (1) = D (x) ), second (D (2) = D (xx) ) and third (D (3) = D (xxx) ) orders. Results obtained for PPV+ of Set 1 with n = 28.

17

 17

  )

Fig. 4 . 7

 47 Fig. 4.7 Idempotency of the density matrices during the SCF iterations for the HPCP-CPSCF, TC2-CPSCF and CG-CPSCF, at zero (D), first (D (1) = D (x) ), second (D (2) = D (xx) ) and third (D (3) = D (xxx) ) orders. Results obtained for PPV+ of Set 1 with n = 28.

Fig. 4 . 8

 48 Fig. 4.8 Frobenius norm of the error vector [cf. Eq. (3.69)] of the density matrices during the SCF iterations for the HPCP-CPSCF, TC2-CPSCF and CG-CPSCF, at zero (D), first (D (1) = D (x) ), second (D (2) = D (xx) ) and third (D (3) = D (xxx) ) orders. This example is for PPV+ of Set 1 with n = 28.

  Fig.4.9 First form of illustration for the converged density matrices using the HPCP, from zero to third order. nnz is the number of non zero elements at 10 -3 . This example is for PPV+ of Set 1 with n = 28. D(1) = D (x) , D(2) = D(xx) , D(3) = D(xxx) .

Fig. 4 .Fig. 4 .

 44 Fig. 4.11 Histogram of number of SCF iterations with respect to the size of the systems, for the four density matrix perturbation methods at first and second orders. The results are obtained for TPA+ and PPV+ of Set 1.

Fig. 4 .

 4 Fig. 4.13 Histogram of number of SCF iterations with respect to the numerical threshold τ . Results are obtained for the four density matrix perturbation methods at first and second orders. The number of cells is fixed at 250 for TPA+ and 62 PPV+ of Set 1.

Fig. 4 .

 4 Fig. 4.14 In (a) and (b), representation of the total calculation time as a function of the number of atoms, for TPA+ (Set 1). This representation compares the AO-CPSCF, CG-CPSCF, TC2-CPSCF and HPCP-CPSCF. The density matrix methods are all truncated at τ = 10 -4 and τ = 10 -8 . In (c) and (d), convergence for (hyper)polarizability per atom as a function of number of atoms. The polarizability α xx and first hyperpolarizability β xxx are calculated using the HPCP-CPSCF truncated at τ = 10 -8 for the 4 polymers of Set 1.

1 :

 1 ! Data: F , D, µ, tolerance,threshold. 2: ! Initialization 3: F = F -µI 4: D = D 0 5: G 0 = H 0 = -∇Ω LNV (D 0 ) 6: n = 0 7: ! Line search by CGFAM routine 8: while ∥ Dn+1 -Dn ∥> tolerance do 9: n = n + 1 10: Ñe = Tr{D n } 11:

27 :

 27 : H n+1 = G n+1 + γ n H n 21: δN = | Ñe -Tr{D n }| 22: ! Result: Converged density matrix 28: D ∞ = Dn+1 Algorithm 2 LNV minimization at constant µ (µ-LNVm) 1: ! Data: F , D, N , tolerance,threshold. 2: ! Initialization 3: 10 -3 < α XS < 10 -2 4: D = D 0 5: G 0 = H 0 = -∇Ω XS (D 0 ) 6: n = 0 7: ! Line search by CGFAM routine 8: while ∥ Dn+1 -Dn ∥> tolerance do 9: n = n + 1 10: Ñe = Tr{D n } 11:

29 : 9 :

 299 ! Result: Converged density matrix 30: D ∞ = Dn+1 Algorithm 3 Xu-Scuseria's modified LNV minimization (XS-LNVm) 142 LNV minimizations and conjugate gradient routine by Jorge Nocedal 1: ! Data: F , D, N , tolerance, threshold. 2: ! Initialization 3: D = D 0 4: G 0 = H 0 = -∇Ω MS (D 0 ) 5: n = 0 6: ! Line search by CGFAM routine 7: while ∥ Dn+1 -Dn ∥> tolerance do 8: n = n + 1 Ñe = Tr{D n } 10:

27 :

 27 : H n+1 = G n+1 + γ n H n 21: δN = | Ñe -Tr{D n }| 22: ! Result: Converged density matrix 28: D ∞ = Dn+1 Algorithm 4 Millam-Scuseria's modified LNV minimization (XS-LNVm) Purification algorithms 1: ! Data: F , M , N , tolerance. 2: ! Initialization 3: ϵ min , ϵ max ← F 4: μ = Tr{F }/M 5: α = min N ϵmax-μ , M -N μ-ϵ min 6: D 0 = α(μI -F ) + (N/M )I 7: n = 0 8: ! Density matrix purification 9: while ∥ D n+1 -D n ∥> tolerance do 10: n = n + 1 11: if ∥ Tr{D n -D 2 n } ∥< 10 -4 then 12:

: end if 21 :n = n + 1 9 :D 15 :

 21915 end while 22: ! Result: Converged density matrix 23: D ∞ = D n+1 Algorithm 6 Canonical purification (Cp)1: ! Data: F , N , tolerance, m(> 2). 2: ! Initialization 3: ϵ min , ϵ max ← F 4: D 0 = (1 -2β m )(ϵ max I -F )/(ϵ max -ϵ min ) + β m I 5: n = 0 6: ! Density matrix purification 7: while ∥ D n+1 -D n ∥> tolerance do 8: if Tr{D n } < N then 10: D n+1 = I -(I -D n ) m (I + mD n ) n+1 = D m n [I + m(I -D n )]! Result: Converged density matrix 16: D ∞ = D n+1 Algorithm 7 Generalized Trace Correcting purification (TCp) 1: ! Data: F , N , tolerance. 2: ! Initialization 3: ϵ min , ϵ max ← F 4: γ min = 0, γ max = 6 5: D 0 = (ϵ max I -F )/(ϵ max -ϵ min ) 6: n = 0 7: ! Density matrix purification 8: while ∥ D n+1 -D n ∥> tolerance do 9: n = n + 1 10:

14 : else 15 :if γ n > γ max then 18 : 2 n 19 :D 148 Purification algorithms 1 :

 1415182191481 γ n = (N -Tr{F (D n )}) /Tr{G(D n )} D n+1 = 2D n -D else if γ n < γ min then n+1 = F (D n ) + γ n G(D n ) 23:end if 24: end while 25: ! Result: Converged density matrix 26: D ∞ = D n+1 Algorithm 8 Trace Resetting purification (TRSp) ! Data: F , M , N , tolerance. 2: ! Initialization

θI 7 :

 7 D 0 = αD min + (1 -α)D max 8: n = 0 9: ! Density matrix purification 10: while ∥ D n+1 -D n ∥> tolerance do 11:n = n + 1 12: if ∥ Tr{D n -D 2n } ∥< 10 -4 then 13:

n 18 : end while 19 : 8 :n = n + 1 9 : 16 :

 18198916 ! Result: Converged density matrix 20: D ∞ = D n+1 Algorithm 9 Hole-particle canonical purification (HPCP) 1: ! Data: F , N , tolerance, R c . 2: ! Initialization 3: ϵ min , ϵ max ← F 4: D 0 = (ϵ max I -F )/(ϵ max -ϵ min ) 5: n = 0 6: ! Density matrix purification 7: while ∥ Dn+1 -Dn ∥> tolerance do Dn = FILTER( D n , R c ) ! Result: Converged density matrix 17: D ∞ = Dn+1 Algorithm 10 Trace correcting purification (TC2) using radial truncation 1: Compute D and F at 0th order 2: Assemble the density matrices from lower orders to construct b. 3: ! Call APROD routine which constructs LHS 4: function APROD(D (k) ,D,F ) 5:

6 :

 6 return F, D, D (k) + 2DF (k) D -D, F (k) 7: end function 8: Initialize D (k) 9: ! Call SYMMLQ, the Saunders routine which resolves Ax = b 10: function SYMMLQ(D (k) ,APROD,b,...) 11:

13 :D

 13 (k) : output 14: end function 15: ! Result: Converged density matrix 16: D (k) Algorithm 11 Resolution of D-CPSCF by Ax = b solver external routine which straight supplies the matrix-vector product Ax. SYMMLQ is the routine written by Michael Saunders. This routine can be found at: http://web.stanford.edu/group/SOL/software/symmlq/ SYMMLQ is designed to solve the system of linear equations Ax = b

P

  (A -SHIFT I N )P t x = P b, i.e. Āx = b where P = C -1 , Ā = P (A -SHIFT I N )P t , b = P b, and return the solution x = P t x. The associated residual is r = b -Āx = P (b -(A -SHIFT I N )x) = P r. EPS refers to the machine precision computed by SYMMLQ. subroutine SYMMLQ( N, B, R1, R2, V, W, X, Y, APROD, MSOLVE, CHECKA, GOODB, PRECON, SHIFT, NOUT, ITNLIM, RTOL, ISTOP, ITN, ANORM, ACOND, RNORM, YNORM ) Subroutine parameters: external: APROD, MSOLVE integer: N, NOUT, ITNLIM, ISTOP, ITN logical: CHECKA, GOODB, PRECON 152 D-CPSCF equation solver and routine by Michael Saunders double precision: SHIFT, RTOL, ANORM, ACOND, RNORM, YNORM, B(N), R1(N), R2(N), V(N), W(N), X(N), Y(N) N (input) = The dimension of the matrix A B(N) (input) = The right hand side vector bR1(N) (input) = Workspace R2(N) (input) = Workspace V(N) (input) = Workspace W(N) (input) = WorkspaceX(N) (output) = Returns the computed solution x Y(N) (input) = Workspace APROD (input) = The external subroutine defining the matrix A For a given vector x, the statement call APROD ( N, x, y ) must return the product y = A x without altering the vector x MSOLVE (input) = The optional external subroutine defining a preconditioning matrix M , which should approximate (A -SHIFT I N ) in some sense. M must be positive definite. For a given vector x, the statement call MSOLVE ( N, x, y ) must solve the linear system M y = x without altering the vector x. In general, M should be chosen so that Ā has clustered eigenvalues. For example, if A is positive definite, Ā would ideally be close to a multiple of I N . If A or (A -SHIFT I N ) is indefinite, Ā might be close to a multiple of I N . NOTE: The program calling SYMMLQ must declare APROD and MSOLVE to be external. CHECKA (input) = If CHECKA = .true., an extra call of APROD will be used to check if A is symmetric. Also, if PRECON = .true., an extra call of MSOLVE ISTOP (output) = An integer giving the reason for termination... -1: β = 0 in the Lanczos iteration; i.e. the second Lanczos vector is zero. This means the RHS is very special. If there is no preconditioner, b is an eigenvector of A. Otherwise (if PRECON is true), let M y = b. If SHIFT is zero, y is a solution of the generalized eigenvalue problem Ay = λM y, with λ = α from the Lanczos vectors. In general, (A -SHIFT I N )x = b has the solution x = (1/α)y where M y = b. 0: b = 0, so the exact solution is x = 0.

2 :

 2 ∥r∥ appears to be less than the value EPS ∥ Ā∥ ∥x∥.

3 :

 3 ∥ Ā∥ ∥x∥ exceeds ∥b∥/EPS, which should indicate that x has essentially converged to an eigenvector of A corresponding to the eigenvalue shift. 4: ACOND (see below) has exceeded 0.1/EPS, so the matrix Ā must be very ill-conditioned.

  acceptable solution. ITN (output) = The number of iterations performed. ANORM (output) = An estimate of the norm of the matrix operator Ā = P (A -SHIFT I N ) P t , where P = C -1 . ACOND (output) = An estimate of the condition of Ā above. This will usually be a substantial under-estimate of the true condition. RNORM (output) = An estimate of the norm of the final transformed residual vector, P (b -(A -SHIFT I N ) x). YNORM (output) = An estimate of the norm of x This is √ x t M x. If PRECON is false, P (b -(A -SHIFT I N ) x).

  (E.5a) -Eq. (E.7a) -8Eq. (E.5b) + Eq. (E.7b)] FED with a second order TE,f (1) (x) = -f (x + 2ξ) + 4f (x + ξ) -3f (x) 2ξ + O(ξ 2 ) (E.9)159to a BED with a second order TE,f (1) (x) = 3f (x) -4f (x -ξ) + f (x -2ξ) 2ξ + O(ξ 2 ) (E.10)and to a CED with a fourth order TE,f (1) (x) = -f (x + 2ξ) + 8f (x + ξ) -8f (x -ξ) + f (x -2ξ) 12ξ + O(ξ 4 ) (E.11)

  mmx m=mmn c m f (x + mξ) + O(ξ d+p ) (E.12)

160 The

 160 principle of finite differences From the last equation, the simplest way to determine the c m coefficients is to constrained c m to have the property of the Lagrange polynomials[91-94] such that mmx m=mmn m k c m =    0, 0 ≤ k ≤ d + p -1 and k ̸ = d 1, have a unique solution for c m . The relation (E.16) corresponds a set of (d + p) linear equations in (m mx -m mn + 1) unkowns. FED BED CED(d,p)m mn = 0 m mx = d + p -1 m mn = -(d + p -1) m mx = 0 m mn = -(d + p -1)/2 m mx = (d + p -1)/2 TableE.1 Indices (m mn ,m mx ) for c m given by (d,p) corresponding to the type of finite representation difference approximation.

17

 17

) 161 Assuming

 161

28 )

 28

2 : ! Initialization 3 :, 3 do 7 :, 3 do 8 : 9 :n = n + 1 10 :

 2378910 D = D 0 4: n = 0 5: ! Computation of seven energy points 6: for i = -3for j = -3while ∥ Dn+1 -Dn ∥> tolerance do Fn = F [D n ] + ξ (i⟨x⟩ + j⟨y⟩) I

Table 1 .

 1 1 Ohno parametrization. t µµ and t µν are the on-site and off-site energy terms (in eV). t µν is assigned with respect to a distance criteria R d .

	.1}. HF-PPP self-consistent field (SCF) calculation (vide infra) is
	initiated after a first tight-binding (TB) calculation. The TB matrix elements are the
	on-site energy t µµ , and the hopping term t µν which is assigned with respect to the C-C

  max , and, εmin < ϵ min . It can be demonstrated[START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF] that, from the Nrepresentable initial guess built from Eq. (2.19), recursive application of the polynomials of Eq. (2.17) maintains the N -representability conditions while converging monotonically to the ground state energy associated with the ground state idempotent density matrix. As outlined by Niklasson,[START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF] the problem with the trace-preserving PMCP is that it slowly converges at very low or high filling factor[START_REF] Anders Mn Niklasson | Expansion algorithm for the density matrix[END_REF][START_REF] Adam | Canonical purification of the density matrix in electronic-structure theory[END_REF].

	20a)
	(2.20b)
	such that: εmax > ϵ

Table 2 .

 2

1 Density matrix solvers and their features.

Table 2 .

 2 2 Carbon nanotubes investigated in this work.

			a( Å) # cells
	(8, 0)	32	4.27	33
	(17, 0)	68	4.26	15
	(12, 0)	48	4.27	21
	(5, 4)	244	33.31	3
	(15, 5)	260	15.37	3
	(11, 5)	268	20.15	3
	Metallic systems are challenging when using linear-scaling methods.[63] Consequently,
	we have only considered the zigzag and chiral carbon nanotubes described in Table {2.3}.
	These CNTs were generated from the nanotube structure generator website TUBGEN[146],
	using a carbon-carbon bond length of 1.421 Å. Semi-empirical approaches such as the
	PPP model are generally used to have a qualitative picture of electronic structure and
	related properties. Nevertheless, it would be interesting to compare it to the state-
	of-the-art Kohn-Sham DFT calculations in order to appreciate the robustness of the
	parametrization described in Section 1.4.2.
	The Table {2.3} reports the energy gaps of the CNTs using two different PPP
	parametrizations, along with DFT calculations. We note a good agreement between the

Table 2 .

 2 3 Energy gap for the π-π * (in eV) at the Γ point of the first Brillouin zone.CNTs # atoms a( Å) # cells Zhang a Ohno a B3LYP b Present work using the Ohno[START_REF] Ohno | Some remarks on the pariser-parr-pople method[END_REF] and Zhang[START_REF] Zhang | Excitation energy calculation of conjugated hydrocarbons: A new pariser-parr-pople model parameterization approaching caspt2 accuracy[END_REF] PPP parametrizations. b KS-DFT calculations using the B3LYP exchange-correlation functional reported from Ref.[START_REF] Matsuda | Definitive band gaps for single-wall carbon nanotubes[END_REF].

	(8,0)	32	4.27	33	1.2481 1.3238	1.283
	(17,0)	68	4.26	15	0.6939 0.6940	0.734
	(12,0)	48	4.27	21	0.0488 0.0413	0.041
	(5,4)	244	33.31	3	0.9192 0.8069	-
	(15,5)	260	15.37	3	0.5860 0.6836	0.66
	(11,5)	268	20.15	3	0.0481 0.0664	0.00

a

  , we need to re-formulate the diagonal components in terms of D. Then, if we start from Eqs. (3.14), (3.21a) and (3.21b), we can recast Eq. (3.23) as

Table 4 .

 4 2 Calculated π-polarizabilities, first and second π-hyperpolarizabilities of fulvenes in au.

			Polarizability		1st hyperpolarizability	2nd hyperpolarizability (×10 4 )
		αxx	αyy	⟨α⟩	βxxx	βxyy	⟨β⟩	γxxxx	γyyyy	γxxyy	⟨γ⟩
						Benzene
	Ref. [158]	38.91		25.94				0.158			0.084
	FFD *	38.9173	38.9173	25.9449	0	0	0	0.1583	0.1583	0.0528	0.0844
	FFD	38.8585	38.8585	25.9057	0	0	0	0.1576	0.1576	0.0525	0.0841
	AO-CPSCF	38.8585	38.8585	25.9057	0	0	0	0.1577	0.1576	0.0524	0.0842
	CG-CPSCF	38.8585	38.8585	25.9057	0	0	0	0.1577	0.1576	0.0525	0.0842
	TC2-CPSCF	38.8616	38.8616	25.9077	0	0	0	0.1576	0.1576	0.0525	0.0841
	HPCP-CPSCF	38.8585	38.8585	25.9057	0	0	0	0.1576	0.1576	0.0525	0.0841
						Naphthalene
	Ref. [158]	98.57		54.59				2.73			0.643
	FFD *	98.5811	65.1941	54.5917	0	0	0	2.7249	1.6631	-0.5855	0.6434
	FFD	98.5899	65.1123	54.5674	0	0	0	2.7468	1.6600	-0.5767	0.6507
	AO-CPSCF	98.5899	65.1123	54.5674	0	0	0	2.7468	1.6600	-0.5766	0.6507
	CG-CPSCF	98.5899	65.1123	54.5674	0	0	0	2.7469	1.6601	-0.5765	0.6508
	TC2-CPSCF	98.5899	65.1123	54.5674	0	0	0	2.7468	1.6600	-0.5766	0.6508
	HPCP-CPSCF	98.5899	65.1123	54.5674	0	0	0	2.7468	1.6600	-0.5767	0.6507
						Anthracene
	Ref. [158]	178.37		93.62				11.95			2.180
	FFD *	178.3843	102.4862	93.6235	0	0	0	11.9506	2.5780	-1.8144	2.1799
	FFD	178.5480	102.3152	93.6211	0	0	0	12.0748	2.6028	-1.7922	2.2186
	AO-CPSCF	178.5480	102.3152	93.6211	0	0	0	12.0748	2.6027	-1.7922	2.2185
	CG-CPSCF	178.5480	102.3152	93.6211	0	0	0	12.0747	2.6027	-1.7922	2.2186
	TC2-CPSCF	178.5480	102.3152	93.6211	0	0	0	12.0748	2.6028	-1.7922	2.2186
	HPCP-CPSCF	178.5480	102.3152	93.6211	0	0	0	12.0750	2.6028	-1.7922	2.2186
						Phenanthrene
	Ref. [158]	159.19		83.83				6.09			2.35
	FFD *	159.2031	92.2974	83.8335	0	0	0	6.0937	1.5063	2.0834	2.3533
	FFD	159.2354	92.3038	83.8464	0	0	0	6.2105	1.5112	2.0920	2.3811
	AO-CPSCF	159.2354	92.3038	83.8464	0	0	0	6.2106	1.5112	2.0921	2.3811
	CG-CPSCF	159.2354	92.3038	83.8464	0	0	0	6.2106	1.5112	2.0621	2.3812
	TC2-CPSCF	159.2354	92.3038	83.8464	0	0	0	6.2105	1.5112	2.0920	2.3812
	HPCP-CPSCF	159.2354	92.3038	83.8464	0	0	0	6.2105	1.5112	2.0920	2.3811
						Azulene
	FFD *	129.2114	70.1698	66.4604	700.9358	-4.9089	417.6161	-0.5129	0.0419	1.3646	0.4516
	FFD	129.3624	70.1411	66.5012	703.2813	-6.7156	417.9394	-0.5148	0.0365	1.3498	0.4443
	AO-CPSCF	129.3624	70.1411	66.5012	703.2816	-6.7156	417.9396	-0.5147	0.0365	1.3498	0.4443
	CG-CPSCF	129.3625	70.1420	66.5015	703.2816	-6.7156	417.9396	-0.5168	0.0367	1.3498	0.4238
	TC2-CPSCF	129.3624	70.1411	66.5012	703.2817	-6.7157	417.9396	-0.5148	0.0365	1.3498	0.4443
	HPCP-CPSCF	129.3625	70.1411	66.5012	703.2816	-6.7151	417.9399	-0.5148	0.0365	1.3498	0.4443
	Table 4.1 Calculated π-polarizabilities, first and second π-hyperpolarizabilities of aromatic
	hydrocarbons in au.

Table 4 .

 4 3 Calculated π-polarizabilities, first and second π-hyperpolarizabilities of fulvalenes in au.

table ) , 4

 )4

	.3 Perturbed dense matrix calculation					109
	Order	n cells	AO-CPSCF	∆(CG-CPSCF)	∆(HPCP-CPSCF)	∆(TC2-CPSCF)
			Set 1	Set 2	Set 1	Set 2	Set 1	Set 2	Set 1	Set 2
	TPA
	1 st , αxx
		18 26 34 42	1.40 × 2.07 × 2.73 × 3.39 ×	4.49 × 10 4 1.08 × 10 5 1.98 × 10 5 3.11 × 10 5	4.60 × 10 -2 6.65 × 10 -2 2.28 × 10 -2 2.86 × 10 -2	4.21 × 10 -3 6.51 × 10 -1 5.52 × 10 -4 5.38 × 10 -3	3.05 × 10 -6 4.56 × 10 -6 6.20 × 10 -6 7.82 × 10 -6	7.97 × 10 -9 2.81 × 10 -8 1.14 × 10 -6 2.00 × 10 -7	3.05 × 10 -6 4.56 × 10 -6 6.20 × 10 -6 7.82 × 10 -6	7.97 × 10 -4 2.79 × 10 -4 1.27 × 10 -5 2.02 × 10 -4
		50	4.05 ×	4.41 × 10 5	3.43 × 10 -2	1.66 × 10 -2	2.91 × 10 -8	4.94 × 10 -6	2.91 × 10 -8	4.94 × 10 -3
		114 122 130	8.69 × 9.36 × 1.00 ×	5.83 × 10 5 8.92 × 10 5 1.05 × 10 6	7.45 × 10 -2 8.03 × 10 -2 8.60 × 10 -2	5.92 × 10 -2 2.69 × 10 -1 4.88 × 10 -1	8.02 × 10 -5 1.00 × 10 -7 1.09 × 10 -7	4.93 × 10 -5 6.87 × 10 -7 1.19 × 10 -7	8.02 × 10 -5 1.00 × 10 -7 1.09 × 10 -7	4.40 × 10 -5 1.59 × 10 -7 4.59 × 10 -7
	3 rd , γxxxx
		18 26 34 42 50 114 122 130	7.84 × 1.21 × 1.65 × 2.09 × 2.52 × 7.32 × 7.75 × 8.19 ×	4.36 × 10 10 4.41 × 10 11 2.16 × 10 12 6.83 × 10 12 1.61 × 10 13 3.12 × 10 13 7.94 × 10 13 1.11 × 10 14	6.71 × 10 -2 1.05 × 10 -2 8.90 × 10 -1 4.53 × 10 -1 2.11 × 10 -1 4.17 5.86 1.74 × 10 1	5.46 × 10 -2 1.51 × 10 -2 2.16 × 10 -1 1.16 × 10 -1 2.65 9.08 1.93 × 10 1 2.47 × 10 1	7.77 × 10 -7 2.01 × 10 -8 6.39 × 10 -7 1.65 × 10 -6 4.17 × 10 -6 5.78 × 10 -5 9.98 × 10 -7 3.13 × 10 -6	4.71 × 10 -7 2.01 × 10 -8 3.64 × 10 -6 1.54 × 10 -7 9.78 × 10 -6 4.85 × 10 -6 7.14 × 10 -5 6.41 × 10 -6	6.79 × 10 -7 4.00 × 10 -7 4.99 × 10 -7 2.60 × 10 -6 2.00 × 10 -6 3.47 × 10 -3 5.06 × 10 -3 6.03 × 10 -3	2.97 × 10 -4 9.76 × 10 -4 1.00 × 10 -2 3.20 × 10 -1 1.20 4.85 3.07 × 10 1 6.00 × 10 1
	TPA+
	1 st , αxx
		18 26 34 42 50 114 122 130	1.62 × 2.28 × 2.95 × 3.61 × 4.27 × 7.58 × 8.25 × 8.91 ×	7.50 × 10 4 1.90 × 10 5 3.56 × 10 5 5.38 × 10 5 7.10 × 10 5 1.48 × 10 6 1.64 × 10 6 1.81 × 10 6	3.36 × 10 -2 1.35 × 10 -2 1.94 × 10 -2 2.52 × 10 -2 3.09 × 10 -2 5.97 × 10 -2 6.54 × 10 -2 7.12 × 10 -2	1.68 × 10 -1 1.41 × 10 -1 6.21 1.04 × 10 1 1.55 × 10 1 1.62 × 10 1 1.46 × 10 1 1.25 × 10 1	6.09 × 10 -7 7.13 × 10 -7 8.63 × 10 -7 9.89 × 10 -7 3.39 × 10 -6 1.10 × 10 -5 1.57 × 10 -6 1.98 × 10 -6	6.82 × 10 -7 2.40 × 10 -6 3.21 × 10 -7 2.14 × 10 -6 3.36 × 10 -5 1.63 × 10 -7 1.32 × 10 -7 3.49 × 10 -6	6.09 × 10 -7 7.12 × 10 -7 8.64 × 10 -7 9.89 × 10 -7 3.39 × 10 -6 1.10 × 10 -5 1.58 × 10 -6 1.99 × 10 -6	2.25 × 10 -6 2.38 × 10 -3 1.80 × 10 -6 2.14 × 10 -2 2.49 × 10 -5 8.51 × 10 -2 6.38 × 10 -2 2.00
	2 nd , βxxx
		18 26 34 42 50 114 122 130	2.946 × 10 4 2.958 × 10 4 2.966 × 10 4 2.972 × 10 4 2.976 × 10 4 2.986 × 10 4 2.987 × 10 4 2.988 × 10 4	-4.80 × 10 7 -3.89 × 10 8 -1.56 × 10 9 -3.60 × 10 9 -5.61 × 10 9 -7.96 × 10 9 -8.55 × 10 9 -8.96 × 10 9	2.61 1.23 1.25 1.26 1.27 1.28 1.28 1.28	1.63 4.94 7.65 1.11 2.82 × 10 1 3.41 × 10 1 5.74 × 10 1 3.08 × 10 1	9.03 × 10 -6 3.58 × 10 -7 3.32 × 10 -6 8.90 × 10 -6 5.18 × 10 -5 5.21 × 10 -5 3.19 × 10 -6 1.95 × 10 -5	7.68 × 10 -6 6.52 × 10 -6 1.55 × 10 -7 3.43 × 10 -6 4.00 × 10 -5 8.54 × 10 -7 7.54 × 10 -7 4.29 × 10 -6	9.05 × 10 -6 3.51 × 10 -7 3.30 × 10 -6 8.92 × 10 -6 5.18 × 10 -4 5.21 × 10 -4 2.76 × 10 -6 1.94 × 10 -4	7.62 × 10 -4 6.52 1.41 × 10 -3 3.43 × 10 2 3.98 × 10 -1 3.64 × 10 -3 4.54 × 10 1 9.42 × 10 1
	3 rd , γxxxx
		18 26 34 42 50 114 122 130	1.41 × 1.84 × 2.27 × 2.71 × 3.14 × 3.58 × 4.01 × 4.88 ×	-1.75 × 10 10 1.25 × 10 12 1.69 × 10 13 7.12 × 10 13 1.49 × 10 14 3.10 × 10 14 3.45 × 10 14 3.81 × 10 14	8.84 × 10 -2 3.05 × 10 -1 4.17 × 10 -1 2.93 5.27 6.31 × 10 -1 9.62 × 10 -1 1.23	5.46 × 10 -1 1.51 × 10 -1 2.16 1.16 × 10 1 2.65 9.08 × 10 2 1.93 × 10 2 2.47 × 10 2	7.24 × 10 -6 1.66 × 10 -5 4.54 × 10 -5 2.38 × 10 -7 9.54 × 10 -6 2.71 × 10 -5 3.44 × 10 -6 1.33 × 10 -6	1.34 × 10 -7 7.48 × 10 -5 6.21 × 10 -6 9.04 × 10 -6 7.77 × 10 -7 3.22 × 10 -8 4.01 × 10 -6 7.87 × 10 -5	3.31 × 10 -3 3.68 × 10 -5 6.01 × 10 -5 2.44 × 10 -4 1.14 × 10 -3 5.91 × 10 -3 5.46 × 10 -3 4.83 × 10 -2	6.90 × 10 -3 1.89 × 10 -1 1.30 3.15 × 10 1 7.10 × 10 1 1.78 × 10 1 2.11 × 10 2 7.01 × 10 2

Table 4 .

 4 4 Numerical accuracy ∆ with respect to the AO-CPSCF for the CG-CPSCF, TC2-CPSCF and HPCP-CPSCF, at each perturbation order and for increasing molecular size. Results are obtained for TPA and TPA+.

	Order	n cells	AO-CPSCF	∆(CG-CPSCF)	∆(HPCP-CPSCF)	∆(TC2-CPSCF)
			Set 1	Set 2	Set 1	Set 2	Set 1	Set 2	Set 1	Set 2
	PPV
	1 st , αxx
		4 6 8 10	5.59 × 10 2 8.32 × 10 2 1.10 × 10 3 1.37 × 10 3	1.18 × 10 3 1.80 × 10 3 2.42 × 10 3 3.05 × 10 3	2.48 × 10 -6 3.88 × 10 -6 9.15 × 10 -7 1.47 × 10 -6	5.96 × 10 -7 1.57 × 10 -6 2.68 × 10 -6 3.90 × 10 -6	7.55 × 10 -8 9.58 × 10 -8 6.89 × 10 -7 7.53 × 10 -7	6.09 × 10 -8 4.86 × 10 -8 6.69 × 10 -8 1.06 × 10 -7	7.56 × 10 -8 9.58 × 10 -8 6.89 × 10 -7 7.53 × 10 -7	6.08 × 10 -8 4.86 × 10 -8 6.71 × 10 -8 1.06 × 10 -7
		12	1.65 × 10 3	3.67 × 10 3	1.99 × 10 -6	5.22 × 10 -6	8.80 × 10 -7	1.78 × 10 -7	8.80 × 10 -7	1.77 × 10 -7
		28 30 32	3.56 × 10 3 3.83 × 10 3 4.10 × 10 3	9.29 × 10 3 9.92 × 10 3 1.11 × 10 4	6.53 × 10 -6 7.18 × 10 -6 6.95 × 10 -6	1.93 × 10 -5 2.11 × 10 -5 2.46 × 10 -5	1.16 × 10 -8 1.24 × 10 -8 1.34 × 10 -8	1.82 × 10 -8 1.34 × 10 -8 5.80 × 10 -9	1.16 × 10 -8 1.25 × 10 -8 1.34 × 10 -8	1.81 × 10 -8 1.33 × 10 -8 6.00 × 10 -9
	3 rd , γxxxx
		4 6 8 10	1.32 × 10 6 2.10 × 10 6 2.89 × 10 6 3.69 × 10 6	1.17 × 10 7 2.05 × 10 7 2.94 × 10 7 3.83 × 10 7	7.01 × 10 -4 2.31 × 10 -4 2.16 × 10 -3 9.48 × 10 -3	4.15 × 10 -5 1.52 × 10 -4 7.04 × 10 -3 6.21 × 10 -3	9.05 × 10 -8 3.30 × 10 -8 5.18 × 10 -7 1.93 × 10 -8	8.74 × 10 -9 6.01 × 10 -8 2.61 × 10 -7 5.03 × 10 -7	8.00 × 10 -8 2.06 × 10 -6 1.76 × 10 -5 7.78 × 10 -5	1.99 × 10 -7 1.40 × 10 -6 1.97 × 10 -5 2.01 × 10 -5
		12	4.48 × 10 6	4.72 × 10 7	5.04 × 10 -2	9.81 × 10 -2	5.21 × 10 -6	9.10 × 10 -6	2.49 × 10 -4	7.29 × 10 -5
		28 30 32	1.08 × 10 7 1.16 × 10 7 1.23 × 10 7	1.36 × 10 8 1.54 × 10 8 1.63 × 10 8	4.42 × 10 -1 6.90 × 10 -1 8.09 × 10 -1	6.23 × 10 -1 5.78 × 10 -1 8.46 × 10 -1	1.94 × 10 -6 3.19 × 10 -7 3.74 × 10 -7	1.17 × 10 -7 7.78 × 10 -7 8.36 × 10 -7	3.42 × 10 -2 4.99 × 10 -2 7.09 × 10 -2	3.72 × 10 -2 7.30 × 10 -2 9.89 × 10 -2
	PPV+
	1 st , αxx
		4 6 8 10	7.43 × 10 2 1.01 × 10 3 1.28 × 10 3 1.56 × 10 3	1.55 × 10 3 2.17 × 10 3 2.79 × 10 3 3.41 × 10 3	3.73 × 10 -3 1.03 × 10 -3 1.40 × 10 -3 3.86 × 10 -3	2.88 × 10 -3 1.06 × 10 -2 3.05 × 10 -7 1.68 × 10 -7	1.48 × 10 -8 5.90 × 10 -7 7.53 × 10 -7 8.33 × 10 -7	4.97 × 10 -7 2.62 × 10 -7 2.61 × 10 -7 4.21 × 10 -8	6.10 × 10 -11 7.18 × 10 -9 2.28 × 10 -9 6.42 × 10 -8	5.09 × 10 -7 2.63 × 10 -7 2.64 × 10 -7 3.80 × 10 -8
		12	1.83 × 10 3	4.04 × 10 3	6.37 × 10 -3	5.13 × 10 -8	2.90 × 10 -6	5.52 × 10 -8	3.93 × 10 -6	5.88 × 10 -8
		28 30 32	3.19 × 10 3 3.47 × 10 3 3.74 × 10 3	9.66 × 10 3 1.02 × 10 4 1.09 × 10 4	1.93 × 10 -2 2.20 × 10 -2 2.47 × 10 -2	6.19 × 10 -6 4.91 × 10 -6 5.54 × 10 -6	6.52 × 10 -9 1.12 × 10 -6 6.57 × 10 -9	2.10 × 10 -6 2.91 × 10 -7 4.06 × 10 -7	2.02 × 10 -7 1.26 × 10 -6 1.87 × 10 -7	2.10 × 10 -6 2.88 × 10 -7 4.05 × 10 -7
	2 nd , βxxx
		4 6 8 10	-1.097 × 10 4 -1.099 × 10 4 -1.102 × 10 4 -1.103 × 10 4	-5.970 × 10 4 -5.981 × 10 4 -6.014 × 10 4 -6.040 × 10 4	1.60 1.57 1.58 1.58	7.66 × 10 -2 1.85 × 10 -1 9.57 × 10 -4 1.92 × 10 -3	6.83 × 10 -8 7.20 × 10 -7 5.70 × 10 -8 4.94 × 10 -6	1.16 × 10 -5 1.05 × 10 -5 2.13 × 10 -7 1.17 × 10 -6	1.21 × 10 -6 2.88 × 10 -5 1.86 × 10 -4 6.33 × 10 -4	4.97 × 10 -4 1.05 × 10 -5 2.21 × 10 -7 1.15 × 10 -6
		12	-1.105 × 10 4	-6.058 × 10 4	1.58	6.89 × 10 -4	4.31 × 10 -5	1.09 × 10 -8	2.28 × 10 -3	5.07 × 10 -8
		28 30 32	-1.108 × 10 4 -1.1093 × 10 4 -1.1096 × 10 4	-6.105 × 10 4 -6.110 × 10 4 -6.115 × 10 4	1.56 1.55 1.53	1.23 × 10 -3 1.10 × 10 -3 1.42 × 10 -3	9.72 × 10 -7 1.47 × 10 -5 1.11 × 10 -6	3.56 × 10 -7 6.21 × 10 -7 6.67 × 10 -7	2.21 × 10 -2 3.17 × 10 -2 8.71 × 10 -2	2.86 × 10 -7 4.51 × 10 -7 8.24 × 10 -7
	3 rd , γxxxx
		4 6 8 10 12 28 30 32	3.18 × 10 6 3.94 × 10 6 4.72 × 10 6 5.50 × 10 6 6.29 × 10 6 1.02 × 10 7 1.10 × 10 7 1.18 × 10 7	3.10 × 10 7 3.83 × 10 7 4.69 × 10 7 5.56 × 10 7 6.45 × 10 7 1.09 × 10 8 1.17 × 10 8 1.53 × 10 8	6.64 × 10 -1 9.07 5.15 6.01 × 10 1 2.10 × 10 2 3.08 × 10 3 6.42 × 10 4 7.90 × 10 4	5.20 × 10 -2 7.76 × 10 -2 2.26 × 10 -1 3.04 × 10 -1 9.48 7.41 6.64 × 10 -1 1.42 × 10 -1	9.05 × 10 -8 3.30 × 10 -7 5.8 × 10 -6 3.33 × 10 -6 1.01 × 10 -6 2.74 × 10 -8 5.21 × 10 -6 4.24 × 10 -6	8.03 × 10 -8 2.34 × 10 -7 6.00 × 10 -7 5.74 × 10 -6 1.49 × 10 -5 3.08 × 10 -8 8.25 × 10 -7 9.41 × 10 -7	5.46 × 10 -2 1.77 1.51 × 10 1 6.01 × 10 1 2.16 × 10 2 9.08 × 10 3 1.42 × 10 4 1.01 × 10 4	2.05 × 10 1 9.59 × 10 -6 3.14 × 10 -4 2.60 × 10 -3 1.14 × 10 -2 9.84 × 10 -4 9.95 × 10 -4 3.21 × 10 -3

Table 4 .

 4 5 Numerical accuracy ∆ with respect to the AO-CPSCF for the CG-CPSCF, TC2-CPSCF and HPCP-CPSCF, at each perturbation order and for increasing molecular size. Results are obtained for PPV and PPV+.

  -2 E -2,-2 E -1,-2 E 0,-2 E +1,-2 E +2,-2 E +3,-2 E -3,-1 E -2,-1 E -1,-1 E 0,-1 E +1,-1 E +2,-1 E +3,-1 E -3,0 E -2,0 E -1,0 E 0,0 E +1,0 E +2,0 E +3,0 E -3,+1 E -2,+1 E -1,+1 E 0,+1 E +1,+1 E +2,+1 E +3,+1 E -3,+2 E -2,+2 E -1,+2 E 0,+2 E +1,+2 E +2,+2 E +3,+2 E -3,+3 E -2,+3 E -1,+3 E 0,+3 E +1,+3 E +2,+3 E +3,+3

	
	              	(E.34)

Mamy Rivo DianzingaDecember 2016

For for the ensemble of N particles, there exist N ! possible permutations p, which constitute the elements of the (symmetric) group of permutation S N . As a result, we can show that: A = A † and A

= √ N !A .

In atomic units: ℏ = m e = e = 1, where ℏ, m e and e are the reduced Planck constant, the electron mass and elementary charge, respectively.

Greek indices refer to atomic orbitals, whereas Roman indices refer to molecular orbitals.

The permutation operator interchanges the coordinate of two electrons. When applied to the right of the atomic orbitals product, this gives: P 12 ϕ µ (r 2 )ϕ ν (r 1 ) = ϕ µ (r 1 )ϕ ν (r 2 ).

Note that: (µµ|ηη)δ µν δ ηκ = (νν|ηη)δ µν δ ηκ = (νν|κκ)δ µν δ ηκ = (µµ|κκ)δ µν δ ηκ

This is also clearly apparent from Eq. (1.3).

At first sight, passing from Eq. (1.83) to Eq. (1.84) might not be that straightforward. Given a functional F , such that F : R M ×M → R, the variation of F (X) with respect to X is formally given by: δF (X)/δX ≡ ∇F (X) = f (X) † , where f is the scalar derivative of F .[2,[START_REF] Brandt Petersen | The matrix cookbook[END_REF]

For a given symmetric (Hermitian) non-degenerate Fock martrix F ∈ R M ×M (∈ C M ×M ), there always exists a similarity transformation, such that:X t F X (X † F X) = diag{ϵ 1 ϵ 2 • • • ϵ M }, where {ϵ i } Mi=1 are the (real) eigenvalues of F , and the transformation matrix X is orthogonal (unitary), ie.X t X = X -1 X = I (X † X = X -1 X = I).From the definition (1.52) of the coefficient matrix C, it is obvious that X ≡ C.

Indeed, the Roothaan-Hall equation[START_REF] Carel | New developments in molecular orbital theory[END_REF][START_REF] George | The molecular orbital theory of chemical valency. viii. a method of calculating ionization potentials[END_REF] corresponds to a generalized eigenvalue problem F C = SCE deriving from the HF equations, expressed in a non-orthogonal basis. Nevertheless, we will retain this naming convention.

We recall that the first term within the braces is the expression to be minimized, and the second term is the constraint ; ∥ • ∥ stands for the Frobenius norm, and ( • ) for the scalar product.

Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP:[START_REF] Matsuda | Definitive band gaps for single-wall carbon nanotubes[END_REF].210.56.150 On: Wed, 08 Jun 2016 11:45:29

R. McWeeny, Proc. R. Soc. A

235, 496 (1956); 237, 355 (1956); 241, 239 (1957). 2 R. McWeeny, Rev. Mod. Phys. 32, 335 (1960). 3 X.-P. Li, R. W. Nunes, and D. Vanderbilt, Phys. Rev. B 47, 10891 (1993). 4 A. D. Daniels, J. M. Millam, and G. E. Scuseria, J. Chem. Phys. 107, 425 (1997).

Notice that ∀D, such that D 2 = D, then ∥D∥ F =  Tr{D}.[START_REF] Challacombe | Linear scaling computation of the fock matrix[END_REF] In that case, it can be shown that λ(D + 0 ) ∈ [-1 2 , 3 2 ].Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 147.210.56.150 On: Wed, 08 Jun 2016 11:45:29

This parameter correponds to the keyword "tolerance" introduced in the algorithms of the Appendices B and C.

In this work the sparsity of the initial guess is dictated by the sparsity of the semi-empirical ZDO Fock matrix [cf.Section 1.4].

Latter in this Chapter the dependence of the perturbed Fock matrices over the perturbed density matrices will be considered as implicit.

Obviously, the projection of Eq. (3.43) onto the subspace H occ-virt or H virt-occ leads to the second order AO-CPSCF equations(3.24).

Acknowledgements

This thesis would not be achieved without the help of some persons. For this reason, I can not fail to thank a number of people. First of all, I deeply thank Dr. Lionel Truflandier who was involved in an important part in this thesis. Dr. Lionel Truflandier was not only my instructor who supported me during these three years. But he was also a wise friend trying to teach me the professional habits necessary to carry on scientific work. I thank Pr. David Bowler, Pr. Anders Niklasson, Pr. Matt Challacombe, Dr. Jörg Kussmann and Pr. Christian Ochsenfeld for their interesting explanations of the linear scaling methods. Once more, I thank Pr. David Bowler and Pr. Anders Niklasson for agreeing to be members of the thesis defense committee.

I can not forget to thank Pr. Ionel Navon and Pr. Vladimir Ivanov who have taken some of their time to explain to me methodological details of their papers.

I also thank the Theoretical Chemistry & Modelling group of the Institut des Sciences Moléculaires which have accepted to invite Pr. Anders Niklasson. I also thank all the group members, especially Pr. Alain Fritsch and the computer engineer Philippe Aurel. I thank all the staff members of the Institut des Sciences Moléculaires, Professors and students, who came to my presentations.

I also thank Pr. Roland Hayn who has been my internship supervisor during my Master degree. Pr. Roland Hayn was my first professional relationship I had in France. I thank him for all he has done for me since I met him, for having helped me during my internship and for accepting to be in the thesis defense committee.

Of course, I thank my whole family, especially my father and my mother who both have always done everything for making me an educated person. Once more, I thank my parents who have always wanted me to obtain a Doctor of Philosphy degree. I deeply thank my parents, from the bottom of my heart, for calling me at least once a week to make sure I was doing well and to mentally and psychologically supported me in my work.

Declaration

I declare that this thesis is the presentation of the outcome of my original research work, after three years of hard labor. However, I clearly certify that contributions from other persons are involved. This dissertation contains equations, tables, figures, algorithms, appendices and a bibliography. This work is officially recognized under the supervision of Pr. Alain Fristch, but technically under the guidance of Dr. Lionel Truflandier at the Institut des Sciences Moléculaires of Bordeaux in France.

List of tables

Density matrix purifications and minimizations

Chapter 3

Density matrix perturbation theory

Set initial D

(1) 0

∞ )

n vs. D

n , D approach. [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF] In a second step, we will increase the size of the molecules and investigate the convergences.

Applications and comparison for small systems

The benchmark set of molecules is given in Figure {4.4}. All the carbon-carbon distances have been fixed to 1.4 Å. All the results presented in the following are in au. [START_REF] Ad Mclean | Computed ground-state properties of fh and ch[END_REF] A very tight convergence parameter of 10 -14 has been used for the all SCF calculations. The mean values of the (hyper)polarizability tensors, [START_REF] Williams | Finite field calculations of molecular polarizability and hyperpolarizabilities for organic π-electron systems[END_REF] have been calculated according to the standard definitions:

We have performed the FFD calculation using the same parameterization as in Ref. [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF], which is referred as FFD * in Table {4.1} to {4.3}. It is rather clear that our FFD * and results from Ref. [START_REF] Klimenko | Dipole polarizabilities and hyperpolarizabilities of the small conjugated systems in the π-electron coupled cluster theory[END_REF] are in good agreement. We have also performed a FFD calculation using the conventional Ohno's parameterization [cf. Table {1.1}]. On comparing the values obtained from the FFD and DMPT approaches using the same parameterization, we found that the methods are in perfect agreement demonstrating the reliability of our implementation. Note that the 1st hyperpolarizability β is expected to be zero for the centro-symmetric structures such as Phenanthrene compared to molecules without inversion symmetry such as Triaheptafulvalene, which is related to the rank of the tensor.

It is worth to note that in our implementation, no symmetry constraints were applied to the Fock matrix.

Methods efficiency for larger systems

In order to expand the comparison of the methods, we have investigated the convergence of (hyper)polarizabilities with respect to the system size. For this purpose, we have used dense density matrices, ie. without employing a truncation scheme. The models used in this section are the polymers presented in Figure {4.5}. Polymer A and B are the transpolyacetylene (TPA) and the polyphenylene vinylene (PPV), respectively. Polymer C (TPA+) and polymer D (PPV+), are structurally derived from polymer A and polymer B. TPA+ and PPV+ were added to the benchmark in order to bypass inversion symmetry, and obtain non zero first hyperpolarizabilities. We also specify that two sets of calculation were performed: (i) for polymers with carbon-carbon distances varying between 1.35 and 1.45 Å (Set 1), (ii) for polymers with fixed carbon-carbon Appendix A

Derivative direct inversion of iterative subspace

The Algorithm {1} outlines how we have implemented the DIIS/D-DIIS extrapolation for the calculations we have performed in this thesis work. In a SCF procedure, we start the DIIS/D-DIIS procedure only after ten iterations around, so that the error vector norm is about the thousandth of the initial error vector norm (step 10). First, this allows to get the solution closer to the convergence region as the norm of the error vector falls gradually. And secondly, this allows to use a sufficient number m of c (k) i optimization coefficients in order to get the averaged effective Fock matrix in the convergence domain. The chosen number m of c (k) i coefficients has to be reasonable (not too small, not too large) as indicated at step 11. The resolution of c (k) i coefficients at step 12 can be performed using a standard linear-equation solver such as the DGESV function from LAPACK library [START_REF] Anderson | LAPACK Users' guide[END_REF]. The step 14 requires the method to be used to compute the density matrix. As a result, one can compute one density matrix (unperturbed order), and even several density matrices at the same time (perturbed propjection). In the case of the perturbed projection, each order is defined by a density matrix, a Fock matrix, and the error vector requiring the lower orders density and Fock matrices.

Appendix B LNV minimizations and conjugate gradient routine by Jorge Nocedal

The algorithms in this appendix are for the three versions of LNV minimization. Actually, these algorithms present the key steps of the line search performed by the CGFAM routine. [START_REF] Charles | Global convergence properties of conjugate gradient methods for optimization[END_REF] The CGFAM routine is described below. After the line search, the density matrix is purified under the threshold parameter. In our calculations, threshold = 10 -2 .

In order to minimize functions, the present works uses the routine CGFAM, written by Jorge Nocedal. [START_REF] Charles | Global convergence properties of conjugate gradient methods for optimization[END_REF] This routine briefly described below is included in the CG+ code. CG+ is a conjugate gradient code for solving large scale, unconstrained, nonlinear optimization problems. CG+ implements three different versions of the conjugate gradient method: the Fletcher-Reeves method, the Polak-Ribiere method, and the positive Polak-Ribiere method (β always non-negative). A web-based server which solves unconstrained nonlinear optimization problems using this CG code can be found at: http://users.iems.northwestern.edu/~nocedal/CG+.html when the termination test is satisfied.

Appendix C

Purification algorithms

D-CPSCF equation solver and routine by Michael Saunders

Generalized D-CPSCF equations at k order:

In order to determine D (k) , we use a conjugate gradient solving Ax = b where

Ax and b are vectors while the terms of D-CPSCF equations are matrices. However, we can suppose the terms of D-CPSCF equations are reshaped in vectors. The resolution seems more technical. b is known since it involves the density matrices from lower orders.

x is D (k) . A does not need neither to be known nor explicitly extracted in some way. On the contrary, we straight need the matrix-vector product Ax which is LHS. Algorithm {11} outlines how we implement this resolution. We implement a routine (APROD) which constructs LHS. The APROD routine has only one variable, D (k) . D and F are like parameters since they are already calculated at 0th order (unperturbed order). SYMMLQ is the routine which performs the conjugate gradient where D (k) is the only variable which changes during the iterations. Of course, there are other parameters required in SYMMLQ routine in order to control the convergence. The feature of SYMMLQ routine is to solve Ax = b without explicitly requiring the matrix A. SYMMLQ only needs an will be used to check if M is symmetric.

GOODB (input) = Usually, GOODB should be .false.

Appendix E

The principle of finite differences

Finite differences of univariate functions

Let be a scalar function f of unidimensional variable x. The derivative of f , denoted here by f (1) , is commonly defined by

Since lim ξ→0 can not be computed, a discrete analogue is used instead,

where ξ (> 0) is a finite small step on the discret set of points x. The relation (E.2) is known as the forward Euler difference (FED) approximation [START_REF] Kendall E Atkinson | An introduction to numerical analysis[END_REF][START_REF] Uri | Computer methods for ordinary differential equations and differential-algebraic equations[END_REF][START_REF] Taras | Simple euler method and its modifications[END_REF] since it uses forward differencing. There exists also the backward Euler difference (BED) approximation:

and the centered Euler difference (CED) approximation:

The difference between these three approximations is given by their intrinsic error. Let us remind that a Taylor series is a representation of function, infinitely differentiable, by an infinite sum of terms, which are calculated from the values of the function derivatives