David Amblard

Xavier Blase

Ivan Duchemin
email: ivan.duchemin@cea.fr

Many-body GW

come

I. INTRODUCTION

The description of the electronic properties of a quantum subsystem embedded in a polarizable, or dielectric, environment (a molecular interface, a solvant, an electrode, etc.) remains a central issue in many fields pertaining to solid-state physics, chemistry or biology. Starting from the historical image charge models for electronic distributions close to a metallic surface or within a dielectric cavity, 1,[START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] the need to describe the response of a polarizable environment to a charged (photoemission) or neutral (optical) excitation in a specific subsystem is still triggering significant developments to combine accuracy with numerical efficiency. In particular, the stabilization of an added hole or electron by the induced electronic rearrangements in a surrounding polarizable environment can be as large as several electronvolts. This so-called polarization energy, together with the additional effects of the electrostatic environment and wavefunction delocalization, strongly renormalize the electronic properties. In many situations, the environment is a complex, potentially infinite, system that cannot be fully described at the same quantum level as the subsystem of interest.

Conceptually close to the historical models of image charges, the polarizable continuum model (PCM) [START_REF] Miertus | [END_REF]4 considers a quantum subsystem located in a cavity carved into a medium described by an homogeneous macroscopic dielectric constant. As a more expensive alternative, discrete polarizable models, where atoms are described as polarizable centers, allow for a more realistic description of screening inhomogeneities at short range in response to an electronic excitation in the quantum subsystem. [5][6][7][8][9] Polarization energies converging slowly with environment size, thousands of polarizable centers may be needed in order to enter a regime where long-range extrapolation can be achieved on the basis of the calculated values. This comes as a challenge to fully ab initio approaches, triggering in practice the description of the environment at a semiclassical empirical level. In such semi-empirical discrete approaches, labeled generically QM/MM, or QM/MMpol to emphasize the polarizable nature of the environment, atoms are provided with effective polarizabilities that reproduce the correct molecular polarizability tensor and/or the macroscopic dielectric tensor of the material. 10,11 Concerning the quantum mechanical formalism used to describe electronic excitations in the central subsystem, a specific class of many-body perturbation theories, the GW [12][13][14][15][16] and Bethe-Salpeter equation (BSE) [17][18][19][20][21] formalisms for the study of charged and neutral electronic excitations, have been recently combined with polarizable models of environment, both at the continuum [22][23][24] and discrete [25][26][27][28] levels. Indeed, while recent studies demonstrated that GW calculations with cubic or even lower scaling could be achieved, [29][30][31][32][33][34][35][36][37][38][39][40] the slow convergence of electrostatic and dielectric (screening) effects with respect to system size forbids a brute force treatment of complex environments within such approaches.

The cost associated with the building of the irreducible electronic susceptibility is usually the bottleneck in time-dependent density-functional theory (TD-DFT) and GW calculations. As a cure, fragments approximations can dramatically reduce such a cost by neglecting wavefunction overlaps between weakly interacting subsystems. These fragment or subsystem approaches have been recently implemented at the GW and BSE levels in the case of systems presenting weakly interacting subunits, including molecular systems, [41][42][43][44][45][46] interfaces, [47][48][49][50] 2D materials, 51,52 but also in the less obvious case of covalent 2D systems. 53 It remains that obtaining the screened Coulomb potential W from the irreducible susceptibility, requiring a matrix inversion, restricts the number of fragments that can be dealt with at the fully ab initio level. As such, the largest fragment-based GW calculations were obtained for a system containing about 800 benzene molecules (4800 non-H atoms) in the case of molecular systems, an already remarkable achievement. 43 Similarly, we recently used a fragment GW approach to partition a multilayer h-BN system in up to 259 h-BN fragments containing 66 non-hydrogen atoms each. 53 The limiting factor to the simulation of larger systems was then the inversion of the Dyson equation to obtain W.

In the present study, we introduce and assess a fully ab initio scheme for embedded GW calculations with hundreds of thousand atoms in the environment. Besides adopting the fragment approximation, we search for an efficient low-rank representation of the susceptibility matrix associated with each fragment, projecting them on-the-fly onto a minimal polarization basis preserving the dipolar, quadrupolar, etc. fragments polarizability tensors. As a result, the size of the Dyson equation for W is dramatically reduced, while preserving the accuracy for the polarization energies at the meV level. We explore the trade-off between accuracy and efficiency in the case of a fullerene crystal, both in the bulk, surface, subsurface and few-layers-slab limits.

II. THEORY

We very briefly outline the GW formalism, directing the reader to thorough reviews and books for a more detailed account on the subject. [54][START_REF] Farid | Electron Correlation in the Solid State -Chapter 3[END_REF][START_REF] Onida | [END_REF][START_REF] Martin | Interacting Electrons: Theory and Computational Approaches[END_REF][START_REF] Ping | [END_REF][59] We further introduce our embedding scheme associated with the definition of the environmental screening, or reaction field. Finally, we describe our fitting scheme allowing to dramatically reduce the size of the fragments dielectric matrix expressed in an effective polarization basis that preserves short to long-range screening effects.

A. The GW formalism

Departing from the use of the electronic density in DFT, the GW formalism takes as a central variable the time-ordered one-body Green's function built from input {ε n , ϕ n } Kohn-Sham eigenstates, namely:

G(r, r ′ ; ω) = n ϕ n (r) ϕ * n (r ′ ) ω -ε n + iη × sgn(ε n -E F ) ( 1
)
where η is a positive infinitesimal and E F the Fermi energy. Relying on perturbation theory to low order in the screened Coulomb potential W, the energy-dependent exchange-correlation self-energy Σ(r, r ′ ; E) can be approximated by the GW operator under the form:

Σ GW (r, r ′ ; E) = i 2π +∞ -∞ dω e iηω G(r, r ′ ; E + ω) W(r, r ′ ; ω) (2)
with v the bare Coulomb potential and W the dynamically screened Coulomb potential built within the random phase approximation (RPA):

W(r, r ′ ; ω) = v(r, r ′ )
+ dr 1 dr 2 v(r, r 1 ) χ 0 (r 1 , r 2 ; ω) W(r 2 , r ′ ; ω).

Such an equation adopts a self-consistent Dyson-like form that needs to be inverted once the independent-electron susceptibility has been built from Kohn-Sham one-body eigenstates:

χ 0 (r, r ′ ; ω) = m,n ( f m -f n ) ϕ * m (r) ϕ n (r) ϕ m (r ′ ) ϕ * n (r ′ ) ω -(ε n -ε m ) + iη × sgn(ε n -ε m ) (4)
with { f m/n } level occupation numbers. The cost of calculating this independent-electron susceptibility grows as O(N 4 ) with respect to the number of electrons N in the system. Equivalently, dropping the space variables for compactness, the Dyson equation can be formulated as:

W(ω) = v + v χ(ω) v, (5)
with χ the RPA interacting susceptibility:

χ(ω) = χ 0 (ω) + χ 0 (ω) v χ(ω). (6)
FIG. 1. Symbolic representation of the fragment approximation. The standard fragment approximation on the independent-electron susceptibility (the X (I) 0 blocks in Eq. 9) can be strictly reformulated by constructing first the block diagonal matrix made of the gas phase interacting susceptibility of the isolated fragments (the

X (I)
g in Eq. 10). It is such an object that we represent here. The resulting fully interacting susceptibility is non-diagonal due to the coupling of the blocks by the off-diagonal V (I J) Coulomb matrix elements (see Eq.

10)

.

The knowledge of Σ GW allows to correct the Kohn-Sham eigenvalues, replacing the density-based exchange-correlation potential v XC by the GW self-energy:

ε GW n = ε KS n + ⟨ϕ n | Σ GW (ε GW n ) -v XC |ϕ n ⟩
with ε GW n the so-called quasiparticle energies at which the self-energy operator must be calculated.

B. Fragment approximation

Our implementation 60,61 of the GW approach adopts a resolution-of-the-identity (RI) formalism [62][63][64] where the density and its variations are expressed over a Gaussian auxiliary basis set {P}. The auxiliary basis functions must thus approximate the space generated by the products of molecular orbitals (MO) {ϕ n }:

ϕ n (r) ϕ m (r) RI ≃ P F P (ϕ n ϕ m ) P(r)
leading in particular to fitting coefficients F P (ϕ n ϕ m ) = 0 when the ϕ n and ϕ m molecular orbitals are non-overlapping. Within this representation, Eqn. 4 rewrites:

χ 0 (r, r ′ ; ω) RI ≃ P,Q X 0 (P, Q ; ω) P(r) Q(r ′ ) (7)
with coefficients

X 0 (P, Q ; ω) = m,n ( f m -f n ) F P (ϕ * m ϕ n ) F Q (ϕ m ϕ * n ) ω -(ε n -ε m ) + iη × sgn(ε n -ε m ) . (8)
In the fragment approximation, the full system auxiliary basis is simply the union of the subsystem basis sets, while each subsystem density can be expressed in its own corresponding basis.

In such a case, the analysis of Eqn. 8 indicates that the joint contribution of two non overlapping subsystems to the representation of the independent-electron susceptibility χ 0 (r, r ′ ; ω) should be zero. Within that limit, the RI-fitted independent-electron susceptibility matrix is thus block diagonal, with blocks corresponding to the constituting subsystems gas phase (isolated) irreducible susceptibilities.

Labelling X (I) 0 the block of fit coefficients corresponding to the independent-electron susceptibility of subsystem (I), the Dyson equation for the full system interacting susceptibility coefficients matrix X(ω) (eqn. 6) reads:

X(ω) -1 = I X (I) 0 (ω) -1 - I,J V (IJ) (9)
with V (IJ) the block corresponding to the Coulomb interactions between auxiliary basis elements of fragments (I) and (J). This equation can also be conveniently rewritten:

X(ω) -1 = I X (I) g (ω) -1 - I J V (IJ) (10)
with X (I) g (ω) the block of coefficients for the isolated (gas phase) interacting susceptibility of fragment (I). In this latter formulation, only the off-diagonal Coulomb interactions V (I J) that accounts for inter-fragments coupling are considered (see Fig. 1). Similar equations can be found in the framework of subsystem TD-DFT. 65,66 In the fragment approximation, the cost of calculating all fragments X (I) 0 and X (I) g blocks scales linearly with respect to the number of fragments. This is a considerable saving and it is now the inversion of the Dyson equation (Eqns. 3 or 10), with cubic scaling with respect to the total number of fragments, that becomes the bottleneck in the limit of a very large number of subsystems. Due to the V (IJ) inter-blocks Coulomb matrix elements, the full interacting susceptibility X(ω), and resulting screened Coulomb potential: W(ω) = V + VX(ω)V, are no longer diagonal by blocks. As such, the inversion of the Dyson equation for X(ω), rewriting Eq. 10 as:

X(ω) = X g (ω) + X g (ω)V off X(ω), (11)
where X g (ω) = I X (I) g (ω) and V off = I J

V (IJ) , becomes a huge effort in terms of computation time and memory in the limit of large systems. It is such a problem that we address below by optimally reducing the size of each X (I) g (ω) susceptibility representation.

C. Constrained reduction of the fragment susceptibilities

A typical calculation that expands MOs over a triple-zeta def2-TZVP basis set 67 involves the corresponding optimized auxiliary def2-TZVP-RI basis set 68 that is composed of 95 orbitals for e.g.

B, C, N, O atoms. In this situation, an environment containing of the order of 10 5 atoms will result in susceptibility matrices of the order of 10 7 in size to be dealt with in the Dyson equation.

For the fragment of interest, for which we want to actually perform a GW correction, we preserve the full auxiliary basis optimized for the corresponding MO basis sets. However, concerning the fragments in the environment, we emphasize that we are mainly interested in their contribution to the reaction field, namely to the induced dipoles, quadrupoles, etc., developed as a response to an electronic excitation in the central subsystem. As such, the full details of the susceptibility in the auxiliary {P} basis may not be necessary.

In order to reduce the computational effort of the Dyson equation, we therefore look for an efficient and compact way to represent the interacting susceptibility of the fragments in the environment. We seek a lower-rank approximation to the gas phase interacting susceptibilities χ (I) g (ω):

χ (I) g (r, r ′ ; ω) RI ≃ P,Q X (I) g (P, Q ; ω) P(r) Q(r ′ ) ( 12
) MODEL ≃ γ,γ ′ X (I) g (γ, γ ′ ; ω) γ(r) γ ′ (r ′ ) ( 13
)
where we use a small basis sets {γ}, which could be for example a minimal Gaussian (sp 3 ) 4-orbitals basis per atom in order to mimic the induced charges-and-dipoles models developed in QM/MM techniques. 10 For the isolated fragment (I), the resulting errors in the interacting susceptibility and its corresponding contribution to the screening, or reaction, field respectively read:

∆χ (I) g (r, r ′ ; ω) = γ,γ ′ X (I) g (γ, γ ′ ; ω) γ(r) γ ′ (r ′ ) - P,Q X (I) g (P, Q ; ω) P(r) Q(r ′ ) (14)
and

∆v (I) screen (r, r ′ ; ω) = dr 1 dr 2 v(r, r 1 ) ∆χ (I) g (r 1 , r 2 ; ω) v(r 2 , r ′ ). (15)
Once the model "polarization" basis {γ} is fixed, the associated X (I) g set of coefficients can thus be simply deduced by minimizing the error in the screening field ∆v (I) screen through a set {t} of test functions:

X (I) g (ω) = argmin X (I) g (γ,γ ′ ;ω) t,t ′ ⟨t |∆v (I) screen (ω)| t ′ ⟩ 2 . ( 16
)
The salient features of this equation is that, since we measure the difference in the reaction field on a set of test functions, we are free to focus on the specific components of the screening field vχv that we want to preserve. Details about the resolution of this equation are given in Appendix A.

A naive choice for the test functions {t} could be the set of auxiliary functions {P} used to construct the reference χ (I) g response function. As shown below, this strategy is rather inefficient. The reason for this failure is that we do not intend to use the model susceptibility to perform a GW calculation on fragment (I) itself. Instead, we need it to build the fully interacting screening potential vχv through the Dyson equation (Eq. 10) that couples (I J) fragments, focusing in the end on the central fragment (I = 0 in Fig. 2) on which we perform the GW calculation. As such, priority should be given to the interactions between the different model fragments, starting from neighboring fragments up to the long range interactions dominated by low order momenta of their polarizability tensors.

A very simple yet successful strategy consists in keeping the test functions localized on the atoms of the fragment (I) for which we seek the model susceptibility, but building the test set {t} out of very diffuse orbitals that will sample the surrounding fragments. The basic idea is that using such diffuse test functions allows to "reach out" for the effect of vχ (I) v on neighboring molecules.

The test set can then be completed with the auxiliary basis {P} associated with fragment (I), but down-weighted, in order to keep the emphasis on the diffuse orbitals during the minimization process. More details about this test basis can be found in Appendix B.

Simultaneously, the preservation of long range interactions can be guaranteed by enforcing low order Cartesian momenta of the susceptibility through Lagrange multipliers:

x m y n z p ∆χ (I) g (ω) x m ′ y n ′ z p ′ = 0. ( 17
)
FIG. 2. Sketch of the susceptibility blocks compression associated with the fragments in the environment (shaded C 60 ). The susceptibility of the central fragment of interest (fragment I=0), on which will be performed the embedded GW calculation, is not reduced.

In the following, we label l max the maximum order m + n + p enforced for a specific model susceptibility. For example, l max = 1 corresponds to the preservation of the fragment neutral monopole, as well as the dipolar polarizability tensor. Imposing such a constraint, along with the use of diffuse test functions, ensures that the reaction field will be well reproduced not only in the vicinity of fragment (I), but in the long-range as well. Imposing higher order polarizability tensors momenta (e.g. l max = 2), can be achieved but with a mild impact as discussed below.

D. Minimal effective polarizability basis

So far the choice of the {γ} has been left arbitrary. Contrary to the induced charge-and-dipole models used in semi-empirical QM/MM techniques, we follow here a more automated route. The polarization basis can be obtained as the result of a generalized minimization process, that is we include the {γ} in the minimization process:

argmin {γ}         min X (I) g (γ,γ ′ ;ω) t,t ′ ⟨t |∆v (I) screen (ω)| t ′ ⟩ 2         (18)
where the only input choice is now the number N P of γ polarization vectors. In practice, the γ functions are expressed in the {P} auxiliary basis set associated with fragment (I), namely

γ(r) = P C γP P(r) , (19)
so that the {C γP } coefficients are now the minimization variables. Formally, setting N p to the dimension of the original auxiliary basis leads to a model susceptibility X g strictly equal to the reference X g . Such an equality holds independently of the choice of the test basis or the constraints.

We tackle this somewhat complex minimization problem by iterating over two distinct steps: i) inner optimization of the X (I) g (ω) matrix elements at fixed {γ} (see equation 16), which is solved exactly through linear algebra; ii) outer optimization of the {C γP } coefficients at fixed X (I) g (ω) matrix elements which is done using gradient descent techniques. More details about this last step can be found in Appendix C.

In the simplest case where test functions {t} only span the auxiliary {P} set, and in the absence of constraint, the Eckart-Young-Mirsky theorem states that the N P optimal {γ} polarization vectors are similar to those defined in Refs. 69,70, namely the leading eigenvectors of the so-called symmetrized susceptibility. The present minimization formulation allows further flexibility with the introduction of test functions and constraints, allowing on-the-fly design of model dielectric functions for specific purposes, emphasizing short-to-long-range or on-site accuracy.

We conclude this Section by emphasizing again that the operations described above (calculations of the reference X (I) g (ω) and X (I) g (ω) matrices, SVD decomposition of related operators, etc.) are performed on isolated fragments, leading to a computational cost that is linear in the number of distinct fragments. In turn, the number of operations related to inverting the Dyson equation for the total screened Coulomb potential, involving interactions between all fragments, is dramatically reduced through reduction of the associated prefactor. As shown below, the optimal polarization basis can be made typically 10 2 times smaller than the original auxiliary basis set, preserving the polarization energy in the meV range, leading to a reduction of the order of 10 6 of the cost associated with obtaining an accurate W operator on the central fragment (I=0).

E. Technical details

The present subsystem approach with minimal representation of the fragments electronic susceptibility has been implemented in the beDeft (beyondDFT) package. 60,61 Input Kohn-Sham eigenstates are generated at the def2-TZVP PBE0 71,72 level with the Orca package. 73,74 We adopt the corresponding def2-TZVP-RI auxiliary basis sets associated with the Coulomb-fitting resolution-ofthe-identity (RI-V) approach. 62,75 The molecular geometries for C 60 and the pentacene are obtained at the def2-TZVP PBE0 level. The face-centered cubic (fcc) C 60 dense phase is constructed taking experimental lattice parameters 76 (a = 14.17 Å), neglecting orientational disorder. The C 60 surface we consider is the (111) surface.

Even though the present scheme allows to compute reaction fields at finite (imaginary) frequen-cies, following the analytic continuation approach to the dynamical GW self-energy implemented in beDeft, 60 we calculate here the polarization energies at the static Coulomb-Hole plus Screened-Exchange (COHSEX) level, 12 with:

Σ SEX (r, r ′ ) = - occp n ϕ n (r) ϕ * n (r ′ ) W(r, r ′ ; ω = 0) (20)
Σ COH (r, r ′ ) = 1 2 δ(r -r ′ ) W(r, r ′ ; ω = 0) -v(r, r ′ ) ( 21
)
with the screened-exchange term involving a summation over occupied (occp) levels only. Following previous studies, 26,27,[41][42][43][44] our polarization energy P n for a given energy level is taken to be the difference between the static COHSEX energy level in the presence of a polarizable environment and its analog in the gas phase, namely:

P n = ε GW e n -ε GW g n ≃ ε COHSEX e n -ε COHSEX g n =P ∆COHSEX n
where the index (e) and (g) in GW e/g and COHSEX e/g stand for embedded (e) and gas (g) phases.

The static COHSEX formalism can be understood as an approximation to the GW operator in the limit where the dielectric response is taken to be energy-independent, as indicated by the (ω = 0) restriction in Eqs. 20 and 21. Since only the static susceptibility needs to be calculated, it is one order of magnitude cheaper than a full GW approach requiring the calculation of nω ≃ 12 independentelectron χ 0 (iω) susceptibilities along the imaginary axis in our analytic continuation scheme.

The accuracy of using such a static limit for the polarisation energies induced by a neighboring dielectric medium was first analyzed in the case of molecules close to a metallic substrate. 47 Such an approximation, that considers only the low-frequency limit of the environment electronic dielectric response, is consistent with standard PCM implementations where the macroscopic dielectric constant originating from electronic degrees of freedom is taken to be the optical one in the low frequency limit (square of the refractive index). Similarly, in standard QM/MM implementations, the semi-empirical atomic polarizabilities are designed to reproduce the fragment electronic polarizability in the static limit. 10 Extension to dynamical reaction fields will be discussed in subsequent studies.

From the knowledge of such a polarization energy P n , calculated at the static ∆COHSEX level, the absolute quasiparticle energy can be obtained as:

ε GW e n = ε GW g n + P ∆COHSEX n .
In the fragment approximation, and in the absence of wavefunction hybridization, such a value yields the energy of the corresponding band center. 77 In particular, one can recover the experimental peakto-peak gap in the dense phase, namely the difference of energy between the highest-occupied and lowest-unoccupied molecular orbital (HOMO/LUMO) band centers. In the following, when needed, the gas phase ε GW g n quasiparticle energy levels will be calculated at the partially self-consistent evGW@PBE0 level, that has been shown 78,79 to be more accurate than non-self-consistent calculations, unless an optimally tuned functional is used for the starting DFT Kohn-Sham calculation. 78,80 Finally, we only need to perform an explicit GW correction for the fragment of interest (I = 0).

In other words, while the interacting susceptibility matrix X(ω = 0) of Eq. 10 is defined for the full system, the screened Coulomb potential matrix W = V + VXV entering Eqs. 20 and 21 is only computed explicitly for the corrected fragment. This enables us to save both on memory footprint and CPU time aspects.

III. RESULTS

A. Validation

We start by looking at the evolution of the static dipolar polarizability tensor for a given fullerene, in the gas phase, obtained with the model susceptibility matrix as a function of N P (the number of polarization vectors we keep). Namely, we compute

[α] i j = -dr dr ′ r i         N P γ,γ ′ γ(r) X g (γ, γ ′ ; ω = 0)γ ′ (r ′ )         r ′ j . ( 22
)
This tensor is a key quantity for the long-range screening effects originating from a given fragment and represents thus a direct measure of the accuracy of the fitted susceptibility. First, we explore the strategy where the test functions {t} are taken to span the auxiliary basis {P} located on the fragment (a fullerene) for which we build the model susceptibility. Relative errors on the (Frobenius) norm of the dipolar polarizability tensor, with respect to a reference calculation using the full auxiliary basis (5700 vectors), are represented in Fig. 3 (blue dots). As expected, this error decreases as the number of polarization vectors increases. For N P = 240, the relative error is of the order of 10%, that is still rather large. This number N P = 240 corresponds to a typical minimal sp 3 basis per atom of the kind used in semi-empirical induced charges-and-dipoles polarizable models.

We now perform the same exercise but adding to the test functions {t} a set of atom-centered for larger values of N p .

The quality of the polarizability tensor obtained with the low-rank susceptibility insures that longrange interactions will be accurately reproduced. We now focus on nearest-neighbor interactions.

We study in particular the HOMO/LUMO energy gap associated with a fullerene (in red in Fig. 4 Inset) surrounded by its first shell of 12 nearest-neighbors (in blue).

In We now study the effect of reducing the size N P of the polarization basis on the 12 surrounding C 60 . As previously, we start by using test functions taken only in the span of the auxiliary basis {P} of the fragment (a fullerene) whose model susceptibility is fitted. The results are provided on Fig. 4 (blue dots). As expected, the error on the central C 60 HOMO-LUMO gap, as compared to the reference calculation, decreases with the number of polarization vectors. For N P = 240, the error is of the order of 100 meV, allowing to have a qualitative result but representing still an error of the order of 10% with respect to the targeted polarization energy. As discussed above, this error converges to zero with N p equal to the size of the original atom-centered Gaussian auxiliary basis set (5700 vectors), but the convergence is too slow.

Similarly to the previous study of the polarizability, we now add diffuse functions in the test basis, with ζ = 0.01 bohr -2 . The related evolution of the error on the polarization energy for the gap is represented in Fig. 4 (orange dots). Clearly, the addition of diffuse functions, allowing to test the quality of the model V X (I) g V reaction field in the vicinity of molecule (I), dramatically accelerates the convergence of the polarization energy with respect to the size of the model susceptibility matrix. For N P = 180, namely ≃ 3% of the original susceptibility matrix size (5700) for one C 60 , the error is now of the order of 1 meV, reaching quantitative accuracy.

We further add the constraints (equation 17), with and without diffuse functions, to enforce the exact dipolar polarizability tensor with l max = 1 (up triangles in Fig. 4), or up to second order moments with l max = 2 (crosses in Fig. 4). In all cases, the constraints improve the accuracy, in particular in the small N P limit, even though their impact is not as important as adding diffuse test functions. Such a behaviour can be understood by looking, e.g., at Fig. 3 for ζ = 0.01 bohr -2 and N P = 120. The dipolar polarizability is already quite well reproduced so that the constraint leads to a small improvement. Fig. 4 reveals that the constraint l max = 2 improves very slightly the error on the gap, in comparison to the constraint l max = 1. When diffuse function are added to the test basis, the differences between l max = 1 (orange up triangles) and l max = 2 (oranges crosses) are less than 0.1 meV for N p ≥ 120. Since imposing the constraint comes at no cost, we keep l max = 1 in the forthcoming calculations.

We provide in the Supplementary Material (Fig. S1) the same study but for the HOMO and LUMO energies. Very similar results are obtained as compared to the gap study. In particular, including diffuse orbitals in the test set, together with the l max = 1, 2 constraints, the errors on the HOMO and LUMO are of the order of 3 to 4 meV, namely just slightly larger than the sub-meV error obtained on the gap.

Very similar results can also be obtained by considering Kohn-Sham basis sets, and associated auxiliary basis sets, containing diffuse orbitals. As a matter of fact, it is well documented that the description of response properties, such as the polarizability, converges faster when diffuse channels are included for the description of molecular orbitals and associated density. 81,82 For sake of consistence with our def2-TZVP/def2-TZVP-RI calculations, we adopt the def2-TZVPD/def2-TZVPD basis sets 82 that add diffuse channels for better response properties. The resulting data are presented in the Supplementary Material (Fig. S2). Due to the addition of diffuse orbitals in the Kohn-Sham basis, the closing of the gap increases by about 2%, from ∆Gap= -0.98 eV to ∆Gap= -1.00 eV, indicating a slightly enhanced screening, consistent with an enhanced molecular polarizability. The evolution of the error with the number N p of polarization vectors is also very similar to what we obtained in Fig. 4, demonstrating that our scheme (diffuse test functions and constraints) is universal, that is rather independent of the chosen input Kohn-Sham basis set. The test provided above for the polarization energy originating from the first-nearest-neighbors is the most stringent test. For fragments located farther away, the dipolar component of the reaction field, that we strictly impose, becomes more and more dominant. This is illustrated in Fig. 5 where we study the HOMO-LUMO gap of a C 60 surrounded now by its two nearest-neighbor shells (see Inset Fig. 5). The size of this cluster amounts to 55 fullerenes. When all fullerenes are described at their full def2-TZVP/def2-TZVP-RI level (in the fragment approximation), the gap of the central fullerene closes by 1.25 eV from the gas phase to the 55-C 60 cluster geometry. We study the effect of reducing the size N p of the polarization basis used to describe the susceptibility of the 42 surrounding C 60 in the second shell, keeping the full auxiliary basis to describe the central fullerene and its first-nearest neighbors shell. As such, we mainly focus on the error induced by the fitting process on fragments located at middle to long-range of the central subsystem of interest.

Consistently with the results obtained for the first shell of neighbors (Fig. 4), these calculations confirm that the addition of diffuse orbitals in the test set dramatically helps in reducing the error below the meV with a small number of N P polarization vectors (compare orange and blue data in Fig. 5). Further, as compared to the first-nearest neighbors case, a smaller number N P of polarization vectors is needed to go below the meV error when the constraint on the dipolar polarizability (l max = 1) is imposed. This is the signature that in the long-range, the dipolar response dominates the screening potential, or reaction field. The combination of diffuse test orbitals with ζ = 0.01 bohr -2 with the constraint l max = 1 leads to an error of the order of 0.1 meV for N P = 60, namely one polarization vector per atom. This is a dramatic reduction of the size of the polarization basis needed to describe the susceptibility blocks entering the Dyson equation.

As an additional validation, we plot in Fig. 6 the screening potential, or reaction field, v screen (r, r 0 ) = [vχv](r, r 0 ) associated with an elementary positive source point-charge located in r 0 , in the vicinity of a single pentacene molecule. Namely, for a test charge located in r 0 , we plot v screen (r, r 0 ) as a function of r in the pentacene plane. In the case of a single fragment molecule, the reaction field reduces to vχ (I) v, with (I) the index of that molecule. The reference reaction field is provided in Fig. 6(a) while the error associated with the model susceptibility, for a fixed N P = 70

B. The C 60 crystal and surface environments

Beyond the small cluster models, we now study the evolution of the gap of a fullerene from the gas phase to a C 60 face-centered-cubic crystal (fcc) environment. We thus want to calculate the closing of the gap by screening effects in the limit of an infinite environment. Such a quantity, labeled ∆Gap below, is also coined the polarization energy. We will consider the cases of bulk C 60 (fullerite) and further of a C 60 at the (111) surface and sub-surface. Experimental photoemission experiments are very much surface sensitive for organic systems, with a limited penetration depth of the input photons or electrons, so that comparison with the surface location is more appropriate. The bulk limit is obtained by immersing a C 60 molecule in a sphere of fullerenes with increasing radius.

The surface and subsurface limits are obtained by using half-a-sphere of polarizable fullerenes in the environment. In the absence of wavefunction delocalization (or band dispersion) in the fragment approach, we focus on the peak-to-peak gap, namely the gap between the center of the HOMO and LUMO bands. We emphasize that the absence of permanent ground-state dipole, quadrupole, etc., in fullerene molecules, precludes the influence of any electrostatic crystal field in the ground-state.

In the present case of fullerene crystals, the reference X (I) g susceptibility can be constructed for a single fullerene and the resulting fitted X (I) g matrix can be "copied" to form the model susceptibility block associated with each fullerene in the environment. Even though rotational disorder was not explored in this study, rotating the X (I) g matrix, to follow the rotation of a given fullerene, can be easily implemented. Beyond rotations, the effect on the polarization energy of changes in the susceptibility matrix associated with slight atomic distortions around some average equilibrium geometry, is expected to be small but may be explored in future studies.

The calculations are performed with the parameters described above, namely keeping the full auxiliary def2-TZVP-RI basis set to describe the susceptibility of the C 60 of interest for which we perform our embedded GW calculation. The same full auxiliary basis set is used for its first-nearest neighbors. For the rest of the environment, namely the second-nearest neighbors and those located farther away, we keep N P = 60 optimized polarization vectors for each fullerene. Diffuse test functions with ζ = 0.01 bohr -2 are adopted with the l max = 1 constraint. We focus on the gap closing (∆Gap) from the gas to the dense phase. Such a polarization energy, originating from the screening by the environment, is described at the ∆COHSEX level as emphasized above.

We plot in Fig. 7 the evolution of the polarization energy, or ∆Gap, associated with the peak-topeak gap for a bulk C 60 (orange dots), a surface C 60 (blue dots) and a C 60 at sub-surface (red dots), as a function of 1/N 1/3 C 60 , where N C 60 represents the number of fullerenes retained in the sphere or half-sphere we use for the environment. Analytic derivations show indeed that such a polarization energy converges slowly with a 1/N 1/3 C 60 behavior in the asymptotic regime. This asymptotic behavior is confirmed numerically in Fig. 7 by the straight dashed-line fits, one for each type of system, going through the calculated energies in the large N C 60 limit.

In the asymptotic infinite bulk size limit, the gap of the central C 60 is closing by 1.69 eV (orange dashed line) with respect to the gas phase. This can be compared to the polarization energy of the biggest studied system. Our calculations are performed for spheres containing up to 3367 C 60 , representing 202 020 carbon atoms. The gap of such a system closes by 1.58 eV, which represents a difference of 0.11 eV with respect to the extrapolated infinite size value. This highlights the difficulty to capture the polarization energy with an accuracy within the 0.1 eV threshold when limited size environments are considered.

1 -1/3 2 -1/3 5 -1/3 16 -1/3 125 -1/3 ∞ -1/3 N -1/3 C 60
In order to compare to the G 0 W 0 @LDA 3.0 eV peak-to-peak gap in a fully periodic bulk calculation performed by Shirley and Louie in their pioneering study, 83 we compute the def2-TZVP G 0 W 0 @LDA gap for an isolated fullerene. Subtracting the 1.69 eV polarization energy in the bulk limit to our gas phase 4.46 eV G 0 W 0 @LDA HOMO-LUMO gap, we end up with a 2.8 eV peak-to-peak gap, in fair agreement with the periodic G 0 W 0 @LDA value.

We also compute the asymptotic infinite limit for the gap of one C 60 located at the (111) surface (see Inset Fig. 7(b)). The blue dashed line gives an asymptotic closing of the gap amounting to 1.40 eV, which is in good agreement with the experimental values of 1.1 eV, 84 1.2 eV 85,86 or 1.4 eV. 87,88 Such values were obtained by taking the provided experimental peak-to-peak gaps subtracted to the experimental 4.9 eV gap value for a fullerene in gas phase. 89 Alternatively, using the analytic continuation approach to the dynamical GW self-energy implemented in beDeft, 60 we compute the def2-TZVP evGW@PBE0 gap value for one C 60 in the gas phase. In the partially self-consistent evGW approach, the HOMO and LUMO energy levels are converged at the 0.1 meV level at the 7 th cycle. Taking the computed 5.1 eV gas phase gap value, and adding the calculated polarization energy at the surface, one obtains a surface peak-to-peak gap of 3.7 eV, within the 3.5-3.8 eV experimental range. The presence of a metallic substrate when performing photo-emission experiments, potentially enhancing the screening in the limit of few C 60 layers, and alternatively the limited screening from the fullerene crystal in the few layer limits, may explain variations between experimental values. On the theoretical side, the influence of the fragment approximation, together with treating screening effects in the static COHSEX limit, remains to be studied.

We further compute the polarization energy for one C 60 at the sub-surface (see Inset Fig. 7(c)).

The red dashed fit provides an asymptotic infinite size polarization energy of 1.63 eV. This value, closer to the bulk limit (difference of 0.06 eV) than the surface limit (difference of 0.23 eV), tends to show a rapid convergence of the polarization energy with respect to the depth of the considered fullerene. Namely, a fullerene in the subsurface presents properties already close to bulk case. We finally conclude this study by making a connection to periodic slab calculations, namely the traditional approach for the study of surfaces with periodic boundary conditions. Along that line, we converge the polarization energy through the addition of C 60 infinite sublayers, rather than by increasing the radius R of a half-sphere of polarizable molecules. Such a representation can be obtained by creating stacks of disks with increasing lateral radius R 2D , as represented in the Inset Fig. 8(a), extrapolating the polarization energy to infinity for a given number (n) of layers with an asymptotic scaling low in 1/R 2 2D (see Ref. 53). Taking as a reference an infinite size monolayer, as illustrated in the Inset Fig. 8(b), the polarization energy scales as (1/n). We emphasize that for (n=5), only 85% of the polarization energy is captured. Summing the 0.85 eV polarization energy of one C 60 embedded in an infinite size monolayer, and the asymptotic polarization energy of 0.54 eV coming from an infinite number of layers with respect to a single monolayer, we find a total polarization energy of 1.39 eV. This value is nearly identical to the 1.40 eV value obtained with the "half-sphere" surface approach.

C. Discussion on CPU and memory requirements

The calculation for the biggest studied system made of 3367 C 60 , amounting to 202 020 atoms as discussed in section III B, required around 8000 total CPU hours distributed on 720 cores.

Alternatively, this represent a typical wall-time (time to completion of the run) of approximately 11 hours. Such a calculation was performed using the Irene SKL partition (Intel Skylake 8168 processors, with a base frequency of 2.70 GHz) of the IRENE supercomputer from GENCI-IDRIS. Only 2.5 terabytes of memory were used to study this system, making possible to run such calculation on smaller size computer clusters. Such a small memory footprint was made possible thanks to the small dimension (60) of the polarization basis for each fragment in the environment. Let's stress out here that describing each fullerene by its full def2-TZVP-RI basis, of dimension 5700, would lead to a memory footprint of 2.6 petabytes for each related matrices, with similar requirements for the total Coulomb potential V or the total susceptibility matrix X.

It would have been impossible to store any of these big matrices on the 310 terabytes available on the supercomputer used for this study. The present numbers are certainly indicative and may change depending on the chosen parameters, but they illustrate how efficient QM/QM (GW/RPA) calculations can be when using compressed susceptibility blocks for the environment.

IV. CONCLUSION

We have presented a fully ab initio QM/QM embedded GW calculation with a polarizable environment containing up to 200 000 atoms, with total typical CPU timings below 10 000 hours at the def2-TZVP level. Such calculations are made possible by adopting a fragment, or sub-system, approximation. Further, the susceptibility matrices associated with the fragments in the environment are reduced on-the-fly to a very low-rank representation, with block dimensions equivalent to the number of atoms in the fragment. This approach yields to representations even more compact than standard semi-empirical approaches based on polarizable atoms described by onsite {p x , p y , p z } local induced dipoles, namely 1 degree of freedom per site as compared to the 3 degrees of polarizable atoms. Such a reduction allows inverting the Dyson equation for the screened Coulomb potential W with very limited CPU and memory requirements.

The present scheme is certainly far from being optimal. The choice of the same localization parameter (ζ) for all (s,p,d,f,g) diffuse channels in the ensemble of test functions can e.g. lead to forthcoming improvements. Our main point was however to show that very simple choices could already dramatically help in constructing model susceptibility operators both extremely compact and accurate in reproducing medium-to-long-range reaction fields.

As future directions, the scheme described here above can be merged with subsystem-DFT techniques [90][91][92] used to improve the input Kohn-Sham eigenstates, allowing to exploit fragment approximations while accounting for the frozen density of the neighboring molecules. In particular, the electrostatic field generated by the environment in the ground-state, of crucial importance in organic media composed of molecules with a permanent dipole, quadrupole, etc. 10,11,[93][94][95] can be further accounted for at the DFT level.

Further, the block-diagonal form of the susceptibility can also be improved within a cluster expansion technique, 43 calculating the gas phase susceptibility of pairs of interacting neighboring fragments, building a susceptibility matrix that is tridiagonal by blocks rather than strictly diagonal.

Within the scope of the present study, we tentatively assess the quality of the fragment approximation in the Supplementary Material, section III, where we compare, in the limit of a small number of fullerenes, the impact of neglecting the overlap of molecular orbitals between fragments. Allowing the delocalization of molecular orbitals over up to 4 fullerenes in the environment, we find that the fragment approximation does not induce errors larger than the percent on the calculated polarization energy. While we cannot study much larger systems at the many-body level without the fragment approximation, the 4-fullerenes environment induces a polarization correction to the central C 60 gap that represents already ∼30% of the true infinite surface limit.

Finally, the GW formalism allows to define a dynamical dielectric response, going beyond the standard low-frequency limit in the optical range of common PCM or QM/MM approaches. The importance of such dynamical corrections needs to be assessed, calculating the polarization energies within the full GW formalism rather than its static COHSEX limit, generalizing the recent merging of the GW formalism with a dynamical PCM. 24

Optimization problem for A 2 = 0

To find the other coefficients {C γP }, we use a greedy strategy. We consider only the first and second terms of the left-hand side of the equation A3, so namely A 2 = 0, and write the factorized form of the general symmetric matrix A 1 = S∆S † , with ∆ a diagonal matrix. Independently of ∆, the required coefficients {C γP } defining the optimal directions are only determined by the S matrix. Using the N cstr first directions we have found, given Z cstr = Z -Γ † 0 U R C 0 U † R Γ 0 , and m = max(0, N -N cstr ), we thus search for solutions S [START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] opt of the low rank sub-problem

S (2) opt = argmin S,∆ / rank(S)≤m Z cstr -Γ † 0 K 0 S∆S † K † 0 Γ 0 .
Defining J = K † 0 Γ 0 , and using its SVD J = U J Σ J V † J , this problem is equivalent to

S (2) opt = argmin S,∆ / rank(S)≤m V † J Z cstr V J -Σ J U † J S∆S † U J Σ J . (C1) Using the SVD of V † J Z cstr V J = U 1 Σ 1 V † 1
, and defining Σ (m) 1 the diagonal matrix made of the first m singular values of Σ 1 , the Eckart-Young-Mirsky theorem leads to a solution S [START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] opt of Eq. C1 such that S [START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] opt = Σ

(m) 1 V † 1 Σ -1 J U † J † , ∆ = V † 1 U 1 .
The required coefficients {C γP } are given by the columns of K 0 S [START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] opt . At the end of the process, the coordinates of the N polarization vectors {γ} in the auxiliary basis, corresponding to the {C γP } coefficients, are given by the columns of U R and of K 0 S [START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] opt . In the specific case where no constraints is enforced, R 0 is not defined, nor U R 0 , and K 0 is the identity matrix. So the basis {γ} is given by the columns of S [START_REF] Jackson | Electrostatics of Macroscopic Media[END_REF] opt , with Z cstr = Z and J = Γ 0 . We emphasize that such a method can be easily adapted to another choice of the first guess {γ 0 }, different of the auxiliary basis.

I. ERROR INDUCED BY THE FIT METHOD FOR THE HOMO AND LUMO ABSOLUTE ENERGIES

In section III.A. of the main manuscript (Fig. 4), we study the error, as a function of the size N p of the polarization basis per C 60 , associated with the gap of one C 60 embedded in a sphere of 12 neighbors. We focus now on the related HOMO and LUMO energies (Fig. S1). To allow a log-scale representation, only the absolute value of the errors are plotted. Negative errors are indicated by markers surrounded by black bold edges. Taking N p = 180 polarization vectors, the value at which the error on the gap becomes lower than 1 meV, we find that the error on the absolute HOMO and LUMO energies is only slightly larger, at the 3-4 meV scale. As such, not only energy differences, but also absolute energy levels, are described quite accurately with very few optimal polarization vectors to describe the environment susceptibility.

II. IMPACT OF AUGMENTING THE KOHN-SHAM AND AUXILIARY BASIS SETS

In the main manuscript, the def2-TZVP/def2-TZVP-RI Kohn-Sham and auxiliary basis sets have been used. We explore the impact of using basis sets containing more diffuse channels known to provide improved polarizability values for a given basis size. S1,S2 To remain consistent with the def2-triple-zeta basis sets used in the main manuscript, we now adopt the def2-TZVPD/def2-TZVPD-RI basis sets that includes additional diffuse channels designed for the calculation of response properties. S2 We reproduce in Fig. S2

III. VALIDITY OF THE FRAGMENT APPROACH

In this section, we compare calculations performed using the fragment approximation to reference calculations where all fullerenes in the environment are treated as a single subsystem, allowing wavefunction delocalization over all molecules in the environment. We again focus on the closing of the HOMO-LUMO gap of a "central" C 60 (in red, Fig. S3) induced by the screening originating from the environment (blue molecules, Fig. S3). Due to the cost of calculating the susceptibility beyond the fragment approximation, the number of fullerenes in the environment is limited to four. To compare data on the same footing, we leave the central (red) fullerene as an isolated fragment to prevent delocalization of the associated HOMO and LUMO. What we want to explore is thus the impact of fragmenting the environment. The results are provided in Table S1. The impact of the fragment approximation, namely not allowing the delocalization of the molecular orbitals over all fullerenes in the environment, increases from 0.5 meV to 1.5 mev, remaining below the percent error. Even though limited to rather small systems, such an exploration provides a better understanding of the quality of the fragment approximation. Note that the ∆Gap associated with the 4 C 60 system (-0.45 eV) represents ∼32% of the infinite surface limit (∆Gap = -1.4 eV; see Fig. 7 main manuscript, blue symbols).

- 2 2 FIG. 3 .

 223 FIG. 3. Evolution of the relative error [α] fit -[α] ref / [α] ref (in percentage) of a fullerene RPA dipolar polarizability tensor [α] (see Eq. 22) as a function of the number N P of polarization vectors (maximum number 5700). || • || corresponds to the Frobenius norm. [α] ref is computed with the susceptibility described in the full auxiliary basis (with 5700 vectors), while [α] fit is calculated using only N p polarization vectors. Results for different choices of test functions are plotted.

2 FIG. 4 .

 24 FIG. 4. Error on the central C 60 [red, Inset] gap as a function of the number N P of polarization vectors per C 60 in the first-neighbors shell [blue, Inset]. Energies on the ordinates are in meV and log-scale. Results for test functions with and without diffuse orbitals, with and without constraints (see section II C), are shown. Diffuse functions use ζ = 0.01 bohr -2 . The reference ∆Gap = -0.98 eV corresponds to the gap reduction from the gas phase to the 13-C 60 cluster when all fragments are described by their full basis (N p = 5700).

2 FIG. 5 .

 25 FIG. 5. Error on the gap for a fullerene surrounded by its two first-shells of neighbors. The susceptibility of the central (in red, Inset) and 12 first-nearest-neighbors (in blue, Inset) C 60 are described by the full auxiliary basis (5700 orbitals), while the susceptibility for each of the 42 C 60 in the second shell of neighbors (in grey, Inset) is described by N P polarization vectors. Energies on the ordinates are in meV and log-scale. Results for test functions with and without diffuse orbitals, with and without constraints (see section II C), are shown. Diffuse functions use ζ = 0.01 bohr -2 (see text). The value ∆Gap = -1.25 eV corresponds to the reference gap reduction for the central C 60 from the gas phase to this 55-fullerenes cluster configuration when all fragments are described with their full basis sets.

FIG. 6 .

 6 FIG. 6. In plane [vχv](r, r 0 ; ω = 0) static reaction field generated by a pentacene molecule in response to a positive unit source charge (grey dot indicated by +e) in r 0 . In (a), the full reaction potential. In (b), (c) and (d) the error ∆v screen (r, r 0 ; ω = 0) (see Eq. 15) with respect to the full reaction potential upon substituting the full χ by its low-rank approximation (with 70 polarization vectors). In (b), low-rank χ model obtained without diffuse functions in the {t} test set; in (c) adding diffuse orbitals to the {t} test set, and in (d) adding further the constraint on the dipolar molecular polarizability. Positions are in angstroms, and units of the reaction field are in meV and in log scale. Equipotentials at (±1, ±10, ±100) meV in figure (a), and at (±0.1, ±1, ±10) meV in figures (b), (c) and (d) are represented by white lines. The value of the reference and of the error associated with the reaction field at the position of the source is indicated.

FIG. 7 .

 7 FIG. 7. Evolution of the gap closing ∆Gap, at the static ∆COHSEX level, between a gas phase fullerene and N C 60 systems, as a function of 1/N 1/3 C 60 . The orange/blue/red dots show calculations for one C 60 embedded at the center of a sphere [Inset (a)], at the surface of a half-sphere [Inset (b)], and at the subsurface of a half-sphere [Inset (c)], respectively. Dashed lines represent [a/N 1/3 C 60 + b] fits.

FIG. 8 .

 8 FIG. 8. Gap evolution for one surface C 60 [in red, Inset (a)], as a function of the inverse number n of C 60 layers. Here the reference is the gap of a fullerene in an infinite C 60 monolayer [see Inset (b)]. The infinite size layer(s) limit is obtained by extrapolating disk(s) to their infinite R 2D radius. The dashed line represents the [a/n + b] fit.

  FIG. S1: Same as Fig. 4 (main manuscript) but for the ∆LUMO and ∆HOMO energies, where the ∆ means the evolution between the gas phase and the first-shell-of-neighbors environment [see Inset in (a)]. Data with black bold edges indicate negative values.

  the equivalent of the Fig.4of the main manuscript, showing very similar results. In the presence of diffuse orbitals (ζ = 0.01 bohr -2 ) in the test set, together with imposing the l max = 1, 2 constraints, the error reaches the meV error for N p larger than ∼120. The converged reduction of the gap ∆Gap = -1.0 eV is slightly larger than the ∆Gap = -0.98 eV obtained with the def2-TZVP/def2-TZVP-RI basis sets, indicating a slight enhancement of the polarizability in the presence of diffuse channels.

2 FIG. S2 :

 2S2 FIG. S2: Error on the central C 60 [red, Inset] gap as a function of the number N P of polarization vectors per C 60 in the first-neighbors shell [blue, Inset]. Energies on the ordinates are in meV and log-scale. Results for test functions with and without diffuse orbitals, with and without constraints (see section II.C. of the main manuscript), are shown. The reference ∆Gap = -1.00 eV corresponds to the gap reduction from the gas phase to the 13-C 60 cluster when all fragments are described by their full auxiliary basis def2-TZVPD-RI (N p = 6900).

- 1 FIG

 1 FIG. S3: Representations of the systems with (a) 3 C 60 , (b) 4 C 60 and (c) 5 C 60 used to test the fragment approach. In red, the fullerene for which we calculate the difference ∆Gap between the gap phase and the systems represented here. In (d) we represent the delocalized (HOMO -1) molecular orbital associated with the 4 C 60 in the environment of system (c) beyond the fragment approximation.

TABLE S1 :

 S1 Variation ∆Gap, from the gas phase to the 3, 4, and 5 C 60 systems, for the central C 60 (in red in Fig.S3) allowing wavefunction delocalization over all C 60 in the environment. The error induced by fragmenting the environment is provided in absolute value and in percentage.Calculations performed at the def2-TZVP/def2-TZVP-RI level (see details in the main text).

	System	∆Gap (meV)	Error (meV)	Relative error (%)
	3 C 60	-245	0.5	0.2
	4 C 60	-350	2.2	0.6
	5 C 60	-447	1.5	0.3

SUPPLEMENTARY MATERIAL

See the Supplementary Material for a study of (I) the error induced by the fit method on the HOMO and LUMO energies, (II) the impact of adding diffuse functions in the Kohn-Sham and associated auxiliary basis sets, and (III) a test of the fragment approach validity.

ACKNOWLEDGMENTS

DA is indebted to ENS Paris-Saclay for his PhD fellowship. This work was performed using HPC resources from GENCI-IDRIS (Grant 2021-A0110910016). XB and ID acknowledge support from the French Agence Nationale de la Recherche (ANR) under contract ANR-20-CE29-0005.

DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available within the article and its associated Supplementary Material.

Appendix A: Computation of the model susceptibility

In this Appendix, we compute the model susceptibility X (I) g (ω) solution of the equation 16. From now on and for compactness, we dropped the exponent (I) because we focus only on one fragment, and the frequency index ω. Results can be computed separately for each required frequency.

Keeping the same notations as the ones used in section II, we denote B the Coulomb matrix between the auxiliary basis and the test basis, such that its coefficients are (Q||t), where (•||•) indicates a Coulomb integral. We define Γ the Coulomb matrix between the polarization basis and the test basis, with coefficients (γ||t), and E (respectively R) the overlap matrix between the auxiliary basis set (respectively the polarization basis) and the constraint basis, with coefficients ⟨Q|x m y n z p ⟩ (respectively ⟨γ|x m y n z p ⟩). The equation 16 can be rewritten

with || • || the Frobenius norm, under constraints like the equation 17, which can also be rewritten

To find this N × N matrix X g , we build feasible solutions, namely matrices which satisfy all constraints enforced by the equation A2. Then, among all such feasible solutions, we compute the optimal one, which minimizes Eq. A1.

For compactness, we define

, and assuming that the rank of R is equal to the number of constraints N cstr , we search for solutions such that

Inverting U R requires to consider the nullspace of U † R , of which we write an orthonormal basis as the columns of the N × (N -N cstr ) matrix K such that U † R K = 0. Under these considerations, X has the form

So far, there are some redundancy in the M 1 and M 2 terms which can be lifted by projection on U and K supplementary subspaces. At the end, defining

R , and using the fact that both O and Z are symmetric, we find that

where A 1 and A 2 are computed via the equation A1.

To solve this equation, let write M ⊕ the pseudo-inverse of the matrix M. Let P be the matrix such that P = K † Γ ⊕ K † , let I N be the identity matrix of size N × N and let L be the matrix such

Using the formula of X of Eq. A3 in Eq. A1, and keeping in mind that U and K span supplementary subspaces, the derivation condition with respect to A 1 and A 2 leads to

If no constraints is enforced, we have R = L = 0, K = I N and P = Γ ⊕ . The formula simplifies thus into

Appendix B: Definition of the test basis

We detail here how the test basis {t} is set-up. It is in particular designed to reproduce as best as possible the effects of the susceptibility of reference at short to middle range. To do so, we rely on test sets that sample uniformly the Coulomb interaction: starting from a set of functions {t 0 }, we orthonormalize this set with respect to the Coulomb norm, leading to the test basis {t 1 }. Namely, writing V 0 the Coulomb matrix such that its coefficients are equal to (t 0 ||t ′ 0 ), where (•||•) indicates a Coulomb integral, the i th vector's coordinates of {t 1 }, in the basis set {t 0 }, is given by the i th column of V -1/2 0 . This basis set is such that (t 1 ∥t ′ 1 ) = δ t 1 ,t ′ 1 . We create two of such sets: i) one spanning the auxiliary basis, namely using the notations of the section II {t 0 } = {P}, and ii) another test basis {t d 1 } based on the set of diffuse orbitals {t d 0 }. In this study, we used for {t d 0 } a set of atom-centered (s,p,d,f,g) diffuse Gaussian orbitals with, for sake of simplicity, the same e -ζr 2 radial part.

The final test basis {t} is the direct sum of {t 1 } and {t d 1 }. The first set is down-weighted with respect to the second one, here by a factor 1/50, to strengthen the influence of the diffuse functions during the fitting process.

Appendix C: Optimization of the polarization basis

In this Appendix, we explain how we compute the polarization basis {γ} (see notation of section II), in which is computed X g . More exactly, we develop our method to compute the {C γP } coefficient set of Eq. 19. Starting from the equation A3, these coefficients were computed through two sub-problems.

Optimization for N = N cstr

Keeping the same notation as in Appendix A, we started by studying the specific problem of N = N cstr , with {γ 0 } = {P}. To be more explicit, we re-index by 0 all matrices which depends directly on {γ 0 }. C 0 being invertible, with rank(C 0 ) = N cstr , there is a matrix S such that A 1 = SC 0 S † , and A 2 = C 0 S † . We search for S such that S (1) opt = argmin

The solution S (1) opt of this non-linear optimization problem is computed by a gradient descent algorithm.

We keep U R = U R 0 + K 0 S (1) opt as a subset of the researched optimal directions {γ}. More exactly, the i th column of U R corresponds to the coefficients {C γP } of the i th vector of the required basis set {γ}.

Supplementary Material : Many-body GW calculations with very large scale polarizable environments made affordable: a fully ab initio QM/QM approach David Amblard, 1 Xavier Blase, 1 and Ivan Duchemin 2 1) Univ. Grenoble Alpes, CNRS, Inst NEEL, F-38042 Grenoble, France 2) Univ. Grenoble Alpes, CEA, IRIG-MEM-L_Sim, 38054 Grenoble, France a) a) Electronic mail: ivan.duchemin@cea.fr


Michel Vasquez

Djamal Habet

Complete and Incomplete Algorithms for the Queen Graph Coloring Problem

The queen graph coloring problem consists in cover ing a n x n chessboard with n 2 queens, so that two queens of the same color cannot attack each other. When the size, n, of the chessboard is a multiple of 2 or 3, it is bard to color the queen graph with only n colors. We have developed an exact algorithm which is able to salve exhaustively this problem for dimensions up to n = 12 and find one solution for n = 14 in one week of computing time. The 454 solutions of Queens_l2 2 show horizontal and vertical symmetries in the color repartition on the chessboard. From this observation, we design a new exact, but incomplete, a!M gorithm which leads us to color Queens_ n 2 problems with n colors for n = 15, 16, 18, 20, 21, 22, 24, 28 and 32 in less than 24 hours of computing time by the exploitation ofsymmetries and other geometM rie properties.

INTRODUCTION

Given an x n chessboard, a queen graph is a graph with n 2 vertices, each of them corresponding to a square of the board. Two vertices are connected by an edge if the corresponding squares are in the same row, column, or diagonals. This corresponds to the rules for moving the queen in a chess game.

In this paper, we consider the following problem: are n colors sufficient to place n 2 queens on the chessboard so that there is no clash between two queens of same color? If so, n is the chromatic number of Queens_ n 2 -noted x n -sin ce the maximum clique num M ber is n. Indeed, the rows, the columns and the 2 main diagonals constitute the 2n + 2 maximum cliques ofthis graph.

A fi rst exact algorithm, working on independent sets, salves the queen graph problems for n up to 14. We have obtained completeness for n = 10 (no solution) and n = 12 (454 solutions), but only one solution for Queens_l4 2 in 168 hours of computing tune. When examining Queens_l2 2 solutions, we find symmetries, both vertical and horizontal, between vertices belonging to different independent sets. That means that we could generate vertices by pairs rather than one by one. This observation leads us to consider these symmetries, and some other geometric operations, in order to reduce the size of the search space. This approach is no more complete but we expect to solve larger instances under the hypothesis that such geometric properties exist.

The paper is structured as follows. After a brief overview of the literature, we describe the main characteristics of our search tree alM gorithm based on independent sets. Theo we present the four kinds of geometric operations implemented to improve our results. Finally,

BRlEF OVERVIEW

Although the graph coloring problem bas been the subject of intense research, applications on the QueensMn 2 problem have received a limited attention: Mehrotra and Trick [START_REF] Mehrotra | A Column Generation Approach for Graph Coloring[END_REF] use a column generation ap� proach to the independent set formulation of the graph coloring prob lem, devising an efficient algorithm to solve the maximum weighted independent set problem arising in the column generation process, and are able to salve problems up ton= 9. Caramia and Dell'Olmo [ 1] suggest a sophisticated algorithm based on the iterative coloring extension of a maximum clique. Extensive computational results are gîven where QueensM n 2 problems are solved up ton= 9.

On the other band, heuristic methods are also used to treat the Queens_n 2 problem. For instance, Kochenberger, Glover et al. [START_REF] Kochenberger | An Unconstrained Quadratic Binary Programming Approach to the Vertex Coloring Ptob lem[END_REF] transfonn this problem into an unconstrained quadratic binary probM lem, and solve it by the tabu search method. Other recent works using local search and solution combining are designed to tackle Queens_ n 2 problems up ton= 16 [START_REF] Chiarandini | An Application of lterated Local Search to Graph Coloring Problem[END_REF][START_REF] Hamiez | Coloration de graphes el planification de rencontres sp ortives: heuristiques, algorithmes et analyses[END_REF]. However, non-exact methM ods fail to prove that X n = n and only give an upper bound for the chromatic number.

NESTED ENUMERATIONS

Our algorithm to solve the Queens_n 2 problem is based on the enu meration of the independent sets (J S) of the que en graph. An I S is defined as a subset of vertices which are not linked by an edge.

Hence, all the vertices of the same J S can have the same color.

The n squares belonging to the first row of the chessboard are definitively colored with n different colors and do not take part in the combinatoric. lndeed, these squares correspond to a maximum clique of the graph. This classical technique bas already been exploited in [START_REF] Marino | Breaking the symmetry of the graph colour ing problem with genetic algorithms[END_REF] and avoids to explore many symmetric confi gurations.

Then, the algorithm iteratively constructs independent sets while checking that there are no two queens on the same square. lt is based on the following statements.

3,1 IS Enumeration and Assignment

No JS can contain more than n vertices because there are n disjoint cliques of size equal ton (for instance, the n rows of the chessboar � ). Accordingly, coloring n 2 vertices with n col ors leads to find n dis joint independent sets I1 .. • In with exactly n vertices each, among a set J Sn of the I S candidates with n vertices. J Sn is generated by a standard depth-first search.

This enumeration pro�edure is embedded in a branching algorithm which tries to assign up to n I S subject to the non overlapping con straint: Vi # j ) Ii n I 3 = 0 (only one queen by square). Since choosing one J S corresponds to coloring n vertices ( or squares), the depth of the search tree is lower than n.

At each stage of the search we identify the subsets IS j i of J S that cover each free square (j, i) of the chessboard. The uncolored square which can be covered by the smallest number of remaining I S corresponds to the branching node. The process backtracks as soon as there is no more J S to cover an uncolored square.

The J S are not constructed once and for ail because of the huge amount of memory required to store them as soon as n is growing up: for n = 15, we have 1475300 IS and for n = 16, we have 9609410 I S (taking into account the filtering technique described below).

Filtering

Reducing the search space while exploring it is a classical technique for solving Constraint Satisfaction Problem, CSP (see, for instance, the study of Sabin and Freuder {9}). Caramia and Dell'Olmo [START_REF] Caramia | Constraint Propagation in Graph Color ing[END_REF] have also applied constraint propagation to the Graph Coloring Problem.

In our approach, the implementation of a filtering process aims at decreasing the huge nurpber of J S to be selected by using the non overlapping constraint that constitutes a• dense constraint network. Moreover, we have found another implicit constraint, based on the cliques of the subMgraph of the uncolored vertices. Hence, our fi lterM ing procedure is two fold. At first, every time an I S Ii is selected, propagation on the non overlapping constraint is carried out. Secondly, after the i th JS asM signment, n -i other I S have to be chosen to form a solution. If, at this stage of the search, there is a clique of size n -i C • in the subgraph of the uncolored squares, then ail the remai�in;JS to be chosen must cover one vertex in this clique. This means that if the condition I3nC n -i = 0 holds for an i i E !Sn, then we can remove I 3 from the search space under the current node of the search tree. By construction, the I 3 cannot produce such a condition with squares belonging to the same row or column (there are no more than n rows and n columns in the chessboard, and each I 3 counts n vertices). This is not the case for the diagonals. For example, at the root of the search tree, we can delete the independent sets which do not cover one square of each of the 2 main diagonals. At the next node, we can consider these 2 main diagonals plus the 4 with n -1 squares, and so on.

. ""

' ' ' '" ' • ' ' ?I "" ' ' ' ' 1 h ' ,, ' Id ' ' Id ' ' ' ' ' ' I"
' ' ''"'"

, 1 h ' ,, " ' " 1 shows the number of J S generated when considering -or not -the filtering on cliques. It also shows the number of nodes proM duced at lev el 3 of the search tree and the elapsed CPU time, in sec onds: half of the search space is pruned when fi ltering on the cliques of the sub-graph constituted by the uncolored vertices.

FIRST RESULTS

Up to QueensMl l 2 the answer is instantaneous : there is no solution for Queens_l0 2 • Thus Xto � 11 since xu = 11 we deduce that Xw = 11. Indeed, the first 10 rows and the fi rst 10 columns of the Queens_l 1 2 solution constitute a 11 colors correct assignment for the 10 x 10 chessboard . Exploring the search tree for Queens_l2 Those geometric properties and other ones, reflecting the repartiM tion of the colors on the chessboard, are defined and exploited in the next section in order to go ahead with Queens_n 2 problems.

USING SYMMETRIES

Before describing each geometric property, let us define the grid maintaining the chessboard of dimension n. The left higher corner square bas the coordinates (0, 0) and the central square has the coorM dinates ( n ;-

, n

;-1 ). Note that, if n is even, the value of n;-1 remains fractional.

The general idea is to use some geometric operators to generate more than one independent set per branching node. A first application of this idea when n = 2p consists in generating J S by pairs and consequently in reducing the depth of the search tree by 2. Secondly, depending on the n value, other geometric operations are combined and implemented.

We go on the assumption that there are solutions presenting such properties. Theo the resulting enumeration algorithm does not guarM antee to find solutions for a given Queens_n 2 instance, Hence, the elaborated algorithm is exact but incomplete.

Horizontal symmetry for n = 2p

The vertex corresponding to the square (j, i) of the chessboard is fixed at the same time with the vertex corresponding to the square (j 1 , i') obtained by the linear transformation:

( j' _ n:_;-1 ) = ( -1 O ) X ( j _ n;-1 ) i ' -n ; l Ü 1 i -":_;-1
Accordingly, the pseudo code expressing the horizontal symmetry (j ', i') = H(j, i) is: The independent sets { a, b} , • • • , { u, v} of Queens_22 2 are enu merated simultaneously: the depth of the search tree is divided by 2. This geometric operator works also for n = 12, 14, 16 and n = 18. However, the computing time becomes very important (see section 6) when n is growing up and no solution is found for n = 20 after 168 hours (i.e. one week) of computing time.

{ J' =n-j-1 i' = i 111"1' • e ' km 0 • ' u V l r o n 1 i h f dm ' • i ' Ill V ' k m o t ' 0 n 1 • u Ill d i h f ' Il'! k ' . i l n 0 u ' 0 V V ,m_ s i h f 1 °, • d ' ' ' i IE! 0 1 • ' n V U ml h k o llli d f • . " 111\'1 • ' u i . km . t n 1 f 1 V d r iPI o h m 1 i r V i lo Ill! h ' . f d • Ill o i u • • k n 0 u ,11!'11 n h f 1 d a i i r C k •• mill l V O k n o h f d u t i Ill 0 • l!'I i • V C •• r m 1 i h d n u ' Ill r • f 1 k • o a " t vm C: . 2" i t f 1 Ill � q 0 V • d n m c h u • r i El k . ' Ill i 0 d h f k m u t a r ' V n 1 •• C 0 i 11!1 d m f 1 • ' V Ill 1 h • 0 • ; 11!'.1 u • t k • n C h Ill ; k d • e • r v m n u a t f . C 1 i:l!'I 2 "' i h u • n • C t I!! k 1 D, d r ... f V li!: i . f d • t • mll!1 u 1 0 i 1 ' • k v Ill n f • h C . V o Ill m r t • i • k d ' 1 f ; h • • n I!! p u • ' V 0 1 1!!11 n h ' i f • i d g mlll k :o u • t • l u o k ' i f n • H""' hm e 1 d 1 0 V • ' i V n f t 1 d " 11!1 r h • • Ill o ' k • . "' u ; 1 ' mi 10 • < 0 V . .., f u r d h • i n t k u • ! •• m • h i Ill 1 C d k l!lll i • f n t 0 r V n k f V • Ill! • d f ; • h 1 0 ' l 11!1 • u " l m

Composition of symmetries for n = 4p

In this case, the transformation is a composition of the vertical and horizontal symmetries. Hence, we consider -at the same time the vertices (j, i), (j ', i'), 0", i 11

) and (j 111 , i' 11 ), such as (j', i') = H(j, i) (horizontal symmetry) and {j 11 , i 11 ) = V{j, i) where the op erator V corresponds to the vertical symmetry transformation defined by the linear operation: ( j'n;-1 ) ( 1 i ' -n ; l Q and the corresponding pseudo code is:

{ j' =j ., ' z =n-i-1 J --2 - • n-1 ) , n-1 i --2-
Then we deduce the coordinates of the vertex corresponding to the square (j '11, i'") by combining those two operators: The independent sets {a,v, k, +}, •--,{b, &, f,*} of Queens_32 2 are enumerated simultaneously: the depth of the search tree is divided by 4. This geometric operator improves the computing times for n = 12, 16 and 20. lt also enables to solve Queens_n 2 for n = 24 and 28.

(j 11' 1 i'") = 1-{, 0 V (j, i),
These first geometric operators are not appl icable to Queenss, 2 when n is odd. We consider this case in the sections below.

Centrnl symmetry for n = 3p

The vertices corresponding to the squares (j, i) and (j', i') are as signed simultaneously, by a symmetry to the central square of the chessboard expressed by the linear operation:

( j' -";-; ) = ( -1 0 ) X ( i ' -n :_;-1 Ü -1 J -- 2 - • n-1 )
. n-1 i--2and the pseudo code generating (i' , j') from (i, j) is:

{ ., ' 1 J =n -J - i' =n-i-1
This operation is wrong for I S that cover the central square of the chessboard, because this square is an invariant for this symmetry. In fact, it is not necessary to assign this last J S. Indeed, when the n -l first I S are assigned, n squares (inc/uding the central one) still have to be colored. Each uncolored square belongs to exactly one row, column and diagonal. Let us prove this statement by contradiction: if there are two uncolored squares in a row (column) then only n -2 squares are colored on this row (column), which is in contradiction with the hypothesis that n -1 I S of size n are fixed. If there are two uncolored squares on a diagonal, then the process backtracks because of the filtering on the cliques constituted by the diagonals.

The independent sets { a, b}, • • • , { m, n} of Queens.15 2 are gen erated simultaneously: the depth of the search tree is divided by 2. The I S { o} covering the central square of the chessboard is built implîcitly. Even if this transformation allowed us to solve the Queens_n 2 for n = 15, no solution was found for n = 21 after 68 hours of comput ing tirne. In the next section, we will deal with this last case.

Composition of rotations for

n = 4p + 1
Here the geornetric operator is 'If /2 rotation around the central square of the chessboard and defined by the linear transfom1ation:

and the eorresponding pseudo code of the operation (j ', i 1

) = n(j, i. ) is: { •/ . 1 J =n-i- i' = j
Iwo other vertices, (j", i") and (j "' , i "' ), are also generated simul taneously so that (j",i") = R,(j',i') = 'R 2 (j, i) and (j 111 ,i'11) 'R(j The independent sets { a, b, c, d}, • , •, {q, r, s, t} of Queens_2 1 2 are constructed sirnultaneously: the depth of the search tree is di vided by 4. As for the Queens_l 5 2 problem, the independent set { u}, covering the central square, is built implicitly.

11 , i 11 ) = R-3 (j, i). li'! h l ln n • i m b t f j ' l1'. r . k . • ,Ill! t m d h'n f b k i ' . ' 0 ' n i l • ultJ'• h u t i d i ' n • f 1 0 "'" 1 m • 0 � ' ' k 1 d h ' k m f b u i . •• n1n t ' f j Sl. (t b i n d h 1 11 ' 0 • n k f • ' . "'i ' t "' ' j f o r ' .. d h i 11 t "' o b k p l • n ' 0 n"' S D :k il'l! .-t d h b 1 f " ' ..l • ' i 0 l n ' "' ' !< b o ' f 1 d h . m , :u f In .. ' � " ' ' l 0 i t k h Ill ' n " d i • k l . b 1 " • r o ip m 0 n i f 'i d i • h ' . b r f 0 j . ' k n uln i • • t m h t ,1 • ' f • k b ' h 110 m o n ' l < .. j d • ' i ' k b uln • ll'il f i ' l rn i . � . • ( o n h n t " • f b h t . d i ' . . " ' m k Q .. ; ' •• " h l d f b ' • Il! k i nlo o I P m t t'l 1 -i:,-j n j dm o r u k f b C . • t • h t'1 . 0 r • 1 ' ' t o h i IP ' .,,. .. j f b n k d m • k t • r n " ' C g t u d h • j '0 f b i 1 .. • rn' • 0 j n d • h • 1" t k b 1 r u f ' 11 u • i klo t mi Q C 1 2 1 l d h • f b 0 r lfl I a ' i ' 1 ' .. ,� l h ' d 0 k m Il> n i f Ill

RESULT SYNTHESIS

The experiments are carried out on a Xeon 3 Ghz machine. The mem ory required to solve the considered instances is smaller than 1 MB, even ifa large number of IS is generated. However, the Queens_22 2 problcm was distributed on 20 P4 2.4 Ghz machines. Indeed, a ma chine is affected to each possible position of the second vertex of the first independent set_ , which covers the square of left higher corner of the ch essboard.

73790

In Table 2, for each dimension n, we express the elapsed time in seconds to reach afirst solution by using (or not: No Op. column) the geometric properties described in the previous sections. The colurnn R, 2 corresponds to the central symmetry (R is the rr/2 rotation de scribed in section 5.4 ). The line "-" means that either no solution is found after 168 hours of CPU time, or the geornetric operator does not match with the dimension n of the chessboard.

Regarding the results, using the horizontal symmetries (7-l column) reduces the time needed to reach a solution frorn 107 sec. to 1 sec. for n = 12, and from 240952 sec. to 5 sec. for n = 14. The combination of both horizontal and vertical symrnetries {7-l& V column) reduces again the required time, from 243 to 1 sec. for the Queens_16 2 prob lern and pennits to solve Queens.n 2 for n = 20, 24, 28 and 32. How ever, in considering the three first columns, the required time to find a solution grows exponentially with the dimension of the Queens_n 2 problem. The resolution of larger problems (dimensions greater than 32) appears to be very difficult.

Finally, readers who are interested in other certificates are in vited to consult the web page www. lgi2p . ema . fr/~vasquez/ queen .htm

CONCLUSION

In this paper, we proposed a branching algorithm to solve the Queens.n 2 coloring problems. It embeds an independent sets enu meration by a depth search procedure, and it is also reinforced by an efficient filtering based on the cliques belonging to the uncolored vertices of the queen graph . This first algorithm achieved the com pleteness for Queens_lQ 2 and Queens.12 2 problems and proved that x14 = 14. The 454 Queens.12 2 solutions show horizontal and vertical sym metries between vertices belonging to pairs of independent sets. This observation gives us the intuition of designing an exact, but incom plete, algorithm around these geometric properties. This idea was primordial to improve the previous results. Indeed, the exploitation of the structural properties of the colors repartition on the chessboard ( symmetries and rotations) allowed us to deal with larger Q ueens_n 2 problems. Thus, we improved significantly our first results: we suc• ceeded to color up to 1024 vertices of the Q ueens.32 2 graph with 32 colors ( against 196 vertices of the Q ueens.14 2 graph).

This work gives eleven counterexamples to the Gardner conjecture ( 4 ) stating that "it is possible to color Queens.n 2 with II colors if and only if n is not divisible by either 2 or 3".

Our future work consists in investigating two directions in order to deal with large instances: firstly we will attempt to find other geomet rîc operators to generate simultaneously more vertices, and secondly we will try to find other important properties that should lead us to implement a more efficient fi ltering procedure.

1

  École des Mines d'Alès-Commissariat à \ 'Énergie Atomique, site EERIE, Parc Scientifique Georges Besse, 30035 Nîmes cedex 1, France. email: {Michel.Vasquez, Djamal.Habet}@ema.fr we summarize the computing time needed to color Queens M n 2 up to n = 32 , before concluding the paper.

Figure 1 .

 1 Figure 1. Chessboard configuration at node 4 for Queens_l2 2

  2 requires less than 1 hour CPU. The enumeration finds 454 solutions proving that x12 = 12, Hence, the algorithm achieved complete ness for Queens_l0 2 and Queens_l 2 2 • However, only one solution is found for Queens_14 2 after a week of computing time. This result is sufficient enough to prove that x14 = 14. But we are not able to solve instances larger than QueensM14 2 • When examining the 454 solutions obtained on Queens_l2 2 , we observe that 98 of them present only horizontal symmetry, 98 present only vertical symmetry (the vertices of different J S corne by pair), and 258 present both vertical and horizontal ones (the vertices of different J S corne four at a time ).

Figure 2 .

 2 Figure 2. Certificate for x22 = 22

Figure 3 .

 3 Figure 3. Certificate for X32 = 32

Figure 4 .

 4 Figure 4. Certifi cate for xis = 15

Figure 5 .

 5 Figure 5. Certificate for x:21 = 21

Table 1 .

 1 Comparative statistics

	n	IS	Without filtering nodes	soc.	IS	With filtering nodes	sec.
	10	724	2883	0	544	991	0
	11 2680	45232	0 1744	15262	0
	12 14200	7559966	163 9440	3834450	75
	13 73712 1273054614 115 146 52008 571824811 45389
	Table

Table 2 ,

 2 Computational time needed to obtain the first solutions.

	n	No Op.	1-l	Geometric Operations 1-l&V 'R_2	'R..&'R 2 &R.3
	12 14 15 16 18 20	107 240952	5 243 2171	4897
	21 22• 24 28		233404'	10 1316	30844
	32