Anthony Scemama

Investigating the volume maximizing the probability of finding ν electrons from variational Monte Carlo data

An algorithm is introduced for the search of a volume, in the threedimensional space, which maximizes the probability of finding να up electrons and ν β down electrons inside the volume, all the other electrons being outside of it. This search is performed after a Variational Monte Carlo sampling of the N -particle density generated by the wave function.

Introduction

Quantum Monte Carlo (QMC) methods have been found very successful in calculating ground-state energies of various molecular systems in recent years. Many applications were presented, including the study of large clusters [START_REF] Belomoin | Effects of surface termination on the band gap of ultrabright Si 29 nanoparticles: Experiments and computational models[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF] and the determination of reaction paths in organic chemistry [START_REF] Moriarty | The addition reaction of propargyl and acetylene : Pathways to cylic hydrocarbons. 2-nd Joint Meeting of the US Sections of the Combustion Institute[END_REF][START_REF] Akramine | Quantum Monte Carlo study of the CO interaction with a dimer model surface for Cr(110)[END_REF]. The precision obtained using these techniques is comparable if not better than this obtained with more traditional methods like coupled cluster using single and double excitations (CCSD), and complete active space self consistent field improved by second order perturbation theory (CAS-PT2). But accuracy in total energies is not always sufficient for a good understanding of the problem. The chemists need some other tools adapted to their needs, for example the analysis of the electron density.

Methods of partition of space based upon the topological analysis of a local function, like the electron density [START_REF] Bader | Atoms in Molecules -A Quantum Theory[END_REF] or Becke and Edgecombe's Electron Localization Function (ELF) [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF][START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF], separate space in different domains. Integrating the density in the subspaces gives the average electron population in each one of these domains, and generally one does not get an integral number of electrons. Savin [START_REF] Savin | Reviews of modern quantum chemistry: A celebration of the contributions of Robert G. Parr[END_REF], following the early work of Daudel's loge theory [START_REF] Aslangul | Aspects of the Localizability of Electrons in Atoms and Molecules: Loge Theory and Related Methods[END_REF], proposed a complete partition of space obtained by looking at the probabilities of finding a given number of electrons, ν, in each subspace. Then, Cancès et al. [START_REF] Cancès | How electrons guard the space: shape optimization with probability distribution[END_REF] adapted the level-set algorithm to investigate the volume which maximizes the probability of finding ν electrons inside the volume Ω and the N -ν other electrons of the system outside of it (inside Ω, the complement of Ω). However, electron correlation was not considered in these studies. A simple and efficient way to compute probabilities from correlated wave functions is the use of quantum Monte Carlo sampling. In this paper, an algorithm for the investigation of such a volume within the Variational Monte Carlo (VMC) framework is proposed.

The organization of the paper is as follows: In Section 2 the basic features of VMC needed to understand the present work are briefly presented. Section 3 is devoted to the definition and presentation of the calculation of the probabilities inside a volume. Section 4 introduces the volume search algorithm, and Section 5 presents some results obtained on very simple molecules. Finally, in Section 6 a summary of the main results of this work is presented.

Variational Monte Carlo

In the Variational Monte Carlo (VMC) method [START_REF] Ceperley | Monte Carlo simulation of a many-fermion study[END_REF] the N -particle probability density of a trial wave function, Ψ 2 T , is sampled by generating a series of electron configurations via the Metropolis algorithm [START_REF] Metropolis | Equations of state calculations by fast computing machine[END_REF]. The configurations are defined as the set of the 3N -electronic coordinates: R = ( r 1 , ..., r N ).

(

In other words, a configuration R may be viewed as a "snapshot" of the molecule showing the instantaneous positions of all the electrons. Ψ T can be of any kind (Hartree-Fock, Kohn-Sham, Multi-Configuration. . . ), and is generally improved by taking account of the electron correlation via a Jastrow factor [START_REF] Gaskell | The Collective Treatment of a Fermi Gas: II[END_REF][START_REF] Garcia-Lekue | An inhomogeneous and anisotropic Jastrow function for non-uniform many-electron systems[END_REF]:

ψ T ( R) = D ↑ ( R)D ↓ ( R) exp [ α i,j J(r iα , r jα , r ij )] (2)
where the sum over α denotes a sum over the nuclei, i,j a sum over the pair of electrons, and D σ (σ =↑ or ↓) are determinants made of one-particle space-orbitals. The exponential term is the Jastrow factor, which is introduced to reproduce the electron-electron cusp condition of the exact wave function and, also, to incorporate some explicit coupling between electron-nucleus and electron-electron coordinates.

Different expressions for the Jastrow part have been presented in the literature. Here, the following form was chosen [START_REF] Assaraf | Computing forces with quantum monte carlo[END_REF]:

J(r iα , r jα , r ij ) = s(x ij ) + p (α) (x iα ) + c 1 x 2 iα x 2 jα + c 2 (x 2 iα + x 2 jα )x 2 ij + c 3 x 2 ij (3) with x ij = r ij 1 + b σ r ij x iα = r iα 1 + b α r iα s(x) = s 1 x + s 2 x 2 + s 3 x 3 + s 4 x 4 p (α) (x) = p (α) 1 x + p (α) 2 x 2 + p (α) 3 x 3 + p (α) 4 x 4 ,
b σ can take two different values depending on the spin of the pairs of electrons considered. In this latter expression the quantities {b σ , b α , c i , s i , p (α) i } play the role of parameters which need to be optimized. A standard approach consists in searching for parameters minimizing the fluctuations of the local energy defined as

E L ( R) ≡ HΨ T ( R)/Ψ T ( R). (4)
This criterion is based on the fact that if Ψ T is the exact wave function, it is an eigenfunction of the Hamiltonian operator. In that case, the local energy is the associated eigenvalue and the fluctuations entirely vanish. Accordingly, small fluctuations are associated with "good" trial wave functions [START_REF] Umrigar | Optimized trial wave functions for quantum Monte Carlo calculations[END_REF].

3 Calculation of the probabilities

Probabilities in VMC

Cancès et al. [START_REF] Cancès | How electrons guard the space: shape optimization with probability distribution[END_REF] proposed a partition of space obtained by trying to find the volume Ω which maximizes the probability of finding ν electrons inside the volume and the N -ν remaining electrons outside of it (inside Ω, the complement of Ω).

P Ω (ν) = N ! ν! (N -ν)! Ω d1d2 . . . dν Ω d(ν + 1) . . . dN |Ψ| 2 (5)
From the data of the VMC simulation, the computation of the P Ω (ν) is straightforward. At every step of the VMC simulation, the coordinates of the electrons (the current configuration) are stored in a file. At the end of the simulation this file contains all the history of the random walk. As a posttreatment, this data is used in the following way. For all the configurations one counts how many electrons are inside Ω. The probability of finding ν electrons is obtained by calculating the ratio N ν /N conf , where N ν is the number of times ν electrons were inside Ω all the other electrons being inside Ω, and N conf is the total number of configurations. Figure 1: Probability distributions on a sphere with radius r centered on the nucleus in the first excited state of 1 He. P Ω (0), P Ω (1) and P Ω (2) are given at the Hartree-Fock level, and in presence of a Jastrow factor.

In the first singlet excited state of the Helium atom (1s 1 2s 1 ), P Ω (0), P Ω (1) and P Ω (2) were computed for Ω being a sphere centered on the nucleus with a variable radius r. The probabilities were computed from a VMC simulation using an Unrestricted Hartree-Fock (UHF) trial wave function improved or not by a Jastrow factor, and the values of P Ω (0), P Ω (1) and P Ω (2) are plotted on figure 1 as a function of the radius r of the sphere. As the two electrons are in different atomic shells, the 1s electron is closer to the nucleus than the 2s electron. Hence, as in the loge theory [START_REF] Aslangul | Aspects of the Localizability of Electrons in Atoms and Molecules: Loge Theory and Related Methods[END_REF], there exists a value for the radius r which maximizes the probability of finding one electron inside Ω and one electron outside of Ω (r = 1.4 a.u.). The following section presents an algorithm to perform the search of the volume corresponding to this maximum.

On figure 1 one can also notice the influence of the dynamic electron correlation introduced by the Jastrow factor, which tends to push away the electrons from each other. In all the region where the probability of having zero or one electron is predominant (0.0 ≤ r < 1.4 a.u.), the curves with or without the Jastrow factor are super imposable. For larger spheres (r ≥ 1.4 a.u.) when the probability of finding two electrons rises the two curves become more and more separated.

Then a similar computation for the aluminum atom (13 electrons) was performed. The radius of the sphere was varied from r = 0.0 to r = 4.0 Å, and the P Ω (ν) corresponding to the sphere were computed. Figure 2 shows that the usual shell structure is recovered. Maxima with a high value (∼ 0.8) are obtained for P Ω (2) and P Ω [START_REF] Cancès | How electrons guard the space: shape optimization with probability distribution[END_REF], which correspond respectively to the shells n = 1 and n = 2.

P( 8) P( 9) P( 10) P( 11) P( 12) P( 5) P( 4) P( 7) P( 6)

Figure 2: Probability distributions inside a sphere with radius r centered on the nucleus of the Aluminum atom.

Refinement of the definition

In order to give a better description of the chemically important regions of molecules, it is necessary to take account of the spin of the electrons in the computations of the probabilities. The volume maximizing the probability of finding two anti-parallel electrons is a better description of a chemical bond than the volume maximizing the probability of finding two electrons whatever their spins. Let P Ω (ν, m) be the probability of finding ν electrons inside Ω, where m = 2S + 1 and S is the sum over the spins of the ν electrons. This way, finding two anti-parallel electrons with all the other electrons outside of Ω contributes to P Ω (2, 1), and finding two up electrons or two down electrons, all the other electrons outside of Ω, both contribute to P Ω (2, 3).

Volume search algorithm

Space is represented as a regular three-dimensional regular boolean grid X. 1 is assigned to the elementary cells of space corresponding to Ω, and 0 is assigned to the elementary cells of space corresponding to Ω. To perform the search, a trial volume is needed. It can be, for example, a sphere placed in a reasonable position, like between two atoms. It can also be an atomic basin within Bader's Atoms In Molecules (AIM) theory [START_REF] Bader | Atoms in Molecules -A Quantum Theory[END_REF], or a valence basin obtained from the topological analysis of the Electron Localization Function [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF][START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF].

The volume is moved and distorted in the elementary cells in such a way that a particular function f (Ω) is maximized. This function is linked to the maximization of P Ω (ν, m).

The maximization function

The following function is maximized during the search process of P (ν, m), for a N -electron system, with M possible values of m:

f (Ω) = P Ω (ν, m) × N i = 0 i = ν M j = 1 j = m [1 -P Ω (i, j)] (6)
When f (Ω) is maximized, P Ω (ν, m) is maximized. At the same time, all the other probabilities are minimized. This way a volume is obtained, where the variance of the average population inside Ω is smaller than that when only P Ω (ν, m) is maximized. Hence, the algorithm converges faster to a higher value of P Ω (ν, m).

Distortions of the volume

In this section, one iteration consists in going through all the cells belonging to the surface separating Ω and Ω. f k is the value of the maximization function at iteration k with volume Ω k , and f 0 is its initial value with the trial volume Ω 0 . First, the cells which belong to the surface separating Ω k and Ωk are identified. These are the cells which have a value equal to 1, and at least one neighbor with a value equal to 0.

Then, for each one of these cells, a contraction of Ω is attempted, by removing the current cell of the boundary (0 is assigned to this cell). If the value of f computed with the contracted volume is greater than f k , the removal of the cell is favorable to the maximization of P Ω (ν, m). The move is accepted, f k is set equal to f and the next cell is treated.

If the contraction of the volume is not favorable to the maximization of f , an expansion of the volume is attempted by assigning the value 1 to all of the first neighbors of current cell. If the value of f computed with the expanded volume is greater than f k the expansion of the volume is favorable to the maximization of P Ω (ν, m). The move is accepted, f k is set equal to f and the next cell is treated.

If neither the contraction nor the expansion of the volume were favorable to maximize f , the volume is not moved at that point, and the next cell is processed.

Once all the cells of the boundary have been moved, a convergence test is realized by checking if f k -f k-1 is below a fixed threshold. If convergence is not reached, the list of cells constituting the boundary between Ω k and Ωk is updated and all the cells of the boundary are treated again.

In this algorithm, note that first the volume is contracted, and only if the move is not favorable the volume is expanded. This sequence leads to the smallest volume when the function f is very flat. Moreover, the cells of the surface are sorted as a function of their distance to the centroid of the volume. This way the symmetry of the trial volume is preserved if it corresponds to a symmetry imposed by the wave function. One can notice that the execution time is proportional to the number of cells constituting the surface.

Acceleration of the process

As the computation time is proportional to the number of points of the boundary between Ω and Ω and to the number of electron configurations ("snapshots" of the positions of the electrons), the process is accelerated in the following way. During the first computation of the probabilities, for each elementary cell the list of configurations which have at least one electron in this particular cell is associated to the cell. The computation of the probabilities is refreshed depending only on the configurations which have at least one electron in the deleted or added cells.

If one cell b has to be removed from Ω, for each element of the list of configurations which have at least one electron in b, the associated values of ν and m are computed and the corresponding value in the array P containing the probabilities is decreased. The values of the probabilities are now independent on the cell b and this cell can be deleted without changing any value of P . The contributions of the configurations associated to the cell b are added by computing again the values of ν and m for each one of these configurations in the new volume.

If one cell has to be added to Ω, the contributions of the configurations associated to the cell are removed from P , the cell is added to the volume, then the contributions due to the configurations associated to the cell are computed and added to P .

The gain of time is remarkable. For an N electron system in a 80 × 80 × 80 grid with 10 6 configurations, there are in average 2N electrons per cell. When a cell is added or deleted, the computation is performed in average over 2N configurations instead of all the 10 6 configurations. As the electron density is not uniform, the computation of the probabilities is faster in the regions of low density than in the regions of high density, but the regions of high density are enclosed in very few elementary cells.

Regularization

When the computation has converged to the optimal volume, the surface is not smooth because of the statistical noise coming from the QMC simulation. To erase the irregularities of the surface, the well known Gaussian or average filters can be used. In that case a blurred region is obtained between Ω and Ω, with real values in the [0.0, 1.0] interval. Then it is necessary to find which isosurface is the best for the separation of the volumes.

Another possibility is the use of the median filter which operates as follows: Each value of the raw data is considered as being at the center of a larger box containing several elementary cells of the 3D grid (typically 3×3×3 or 5×5×5). The list of values within this box are sorted, and the middle value of the list is assigned to the central elementary cell. After the use of this filter, the numerical values of the grid are values which belong to the same set as the initial values: For a boolean grid, the application of this filter does not introduce intermediate real values between 0 and 1. This technique was used to present the results.

Some applications

The inter atomic equilibrium distances were computed with the Gaussian98 [START_REF] Frisch | Gaussian 98[END_REF] series of programs, using the B3LYP density functional [START_REF] Becke | Density-Functional Thermochemistry III: the Role of Exact Exchange[END_REF][START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF] and Dunning's correlation consistent double-ζ Gaussian basis set [START_REF] Dunning | Gaussian basis sets in correlated molecular calculations I[END_REF] (cc-pVDZ): 1.6037 Å in LiH and 1.1044 Å in N 2 . The trial wave functions used for the VMC simulations were composed of the Hartree-Fock determinant built on the cc-PVDZ Gaussian basis set, and a Jastrow factor optimized both for the given molecule and the determinantal part of the wave function.

For each molecule, a set of 800 independent walkers was chosen to achieve a 10,000 step random walk using the Variational Monte Carlo method (VMC). Every 10 steps1 , for each walker the coordinates of the electrons were saved, so we end up with 800,000 quasi-independent electronic configurations. Symmetry operations with respect to the xy, yz, and xz planes let us multiply the number of configurations according to the symmetry of the wave function. For homonuclear diatomic molecules, 6,400,000 configurations were obtained. The errors were estimated by block averaging: The full set of configurations was divided into ten different smaller sets of configurations (blocks), and the probabilities were calculated in the final volume using each block of data. The usual formula was used to compute the variance of the probabilities, and thus the statistical error was estimated. The maximum value of this error was 4.10 -3 , so we report safely the values of the probabilities with two significant digits.

All the volumes were rendered using the Molekel visualization program [START_REF] Flkiger | MOLEKEL 4.0[END_REF].

H 2 towards dissociation

It is well known that the Restricted Hartree-Fock (RHF) method does not describe correctly the dissociation of the dihydrogen molecule. The constraint that the two electrons have to occupy the same space function imposes the weights of the ionic (H + +H -↔H -+H + ) and covalent (H-H) contributions to be equal to 50%, even towards dissociation. This does not correspond to the physical situation as the dihydrogen molecule is known to dissociate homolytically. Therefore the H 2 molecule was studied towards dissociation, with an inter-atomic distance of 4 Å, by sampling with VMC a RHF wave function, and a RHF wave function improved by a Jastrow factor. If space is separated into two symmetric subspaces A and B, each containing one proton (H A and H B ), three different situations arise. The first one describes the homolytic dissociation of H 2 with one electron in subspace A and one electron in subspace B. In the second situation, there are two electrons in subspace A and zero in subspace B, and in the third situation there are two electrons in subspace B and zero in subspace A. The last two situations correspond to an ionic dissociation of H 2 .

From these configurations one can compute the probabilities of finding zero, one or two electrons in one subspace. With the RHF wave function, one obtains a probability of 0.50 to find one electron, a probability of 0.25 to find zero electron and a probability of 0.25 to find two electrons in a subspace. This corresponds to a half-resonant and half-covalent character. When electron correlation is included in the calculation via the Jastrow factor, the probability of finding one electron in a half-space gets close to one (0.96 was obtained). This describes the high localization of one single electron in each subspace, and thus the homolytic dissociation of H 2 .

LiH

H Li Ω (2,1)
Ω(2,1)

1.6037 Angstroms A Hartree-Fock trial wave function of LiH was sampled with VMC. The volumes which maximize the probability of finding two anti-parallel electrons around the lithium atom and around the hydrogen atom are represented in figure 3. A volume surrounding the Lithium atom is obtained, and the rest of the space corresponds to the hydrogen atom. Indeed, looking for the volume which maximizes P Ω (2, 1) around the hydrogen atom and the volume which maximizes P Ω (2, 1) around the lithium atom consists in maximizing the same function. Hence there is a virtually unique separation between these two volumes.

P Ω(2,1),Li (2, 1) = P Ω(2,1),H (2, 1) = 0.98. These two probabilities are very close to one, which may be interpreted as a high localization of the anti-parallel electron pairs since the variance of the population in each volume is very low. The values of the probabilities illustrate the strongly ionic interaction between the Li + ion and the H -ion.

N 2

A Hartree-Fock trial wave function, improved by a Jastrow factor was sampled with VMC. In the nitrogen molecule, one can describe the lone pairs and the core regions as volumes which maximize P (2, 1) (denoted respectively Ω(2, 1),lone pair and Ω(2, 1),core), and the triple bond as the volume which maximizes P (6, 1) in the central region of the molecule (denoted Ω(6, 1),bond). These volumes have been searched, and are displayed on Volumes which maximize the probability of finding two, four or six electrons between the nuclei (Ω(2, 1),bond, Ω(4, 1),bond, Ω(6, 1),bond), volumes corresponding to the cores (Ω(2, 1),core) and volumes corresponding to the lone pairs (Ω(2, 1),lone pair).

In the triple bond region, as expected, Ω(2, 1) ⊂ Ω(4, 1) ⊂ Ω(6, 1). From the values of table 1, the core electrons are, as expected very localized. Hence, we interpret the high value of P Ω(6,1),bond [START_REF] Silvi | Classification of chemical bonds based on topological analysis of electron localization functions[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF] in the triple bond as one electron (whatever its spin) coming from a lone pair. This is confirmed by the important value of P Ω(2,1),lonepair [START_REF] Belomoin | Effects of surface termination on the band gap of ultrabright Si 29 nanoparticles: Experiments and computational models[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF] in each lone pair. Identically, the high value of P Ω(6,1),bond [START_REF] Bader | Atoms in Molecules -A Quantum Theory[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF] in the triple bond region is due to one electron escaping to a lone pair region, which is also confirmed by the high value of P Ω(2,1),lonepair [START_REF] Moriarty | The addition reaction of propargyl and acetylene : Pathways to cylic hydrocarbons. 2-nd Joint Meeting of the US Sections of the Combustion Institute[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF].

Both processes occur simultaneously, and one can describe four major cases. In the first case, one electron of spin σ goes from a lone pair to the triple bond, and one electron of spin σ goes from the triple bond to the same lone pair, nothing has changed. If the latter electron goes to the other lone pair, this corresponds to the the second case where an electron is transferred from one lone pair to the other. This contributes to P Ω(2,1),lonepair (1, 2) and P Ω(2,1),lonepair [START_REF] Moriarty | The addition reaction of propargyl and acetylene : Pathways to cylic hydrocarbons. 2-nd Joint Meeting of the US Sections of the Combustion Institute[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF]. In Table 1: Probability distributions in the N 2 molecule. P Ω(ν,m),X (ν , m ) denotes the probability of finding ν electrons with spin multiplicity m in the volume Ω(ν, m), X. Ω(ν, m), X is the volume which maximizes the probability of finding ν electrons with spin multiplicity m , in the region X. X can be the bond regions, a lone pair region or a core region. Ω(2, 1),bond Ω(4, 1),bond Ω(6, 1),bond P Ω(2,1),bond (1, 2) = 0.28 P Ω(4,1),bond (3, 2) = 0.22 P Ω(6,1),bond (5, 2) = 0.20 P Ω(2,1),bond (2, 1) = 0.32 P Ω(4,1),bond (4, 1) = 0.23 P Ω(6,1),bond (6, 1) = 0.21 P Ω(2,1),bond (3, 2) = 0.19 P Ω(4,1),bond (5, 2) = 0.19

P Ω(6,1),bond (7, 2) = 0.20 P Ω(4,1),bond (4, 3) = 0.12 P Ω(6,1),bond (6, 3) = 0.11 Ω(2, 1),core Ω(2, 1),lone pair P Ω(2,1),core (1, 2) = 0.08 P Ω(2,1),lonepair (1, 2) = 0.25 P Ω(2,1),core (2, 1) = 0.85 P Ω(2,1),lonepair (2, 1) = 0.39 P Ω(2,1),core (3, 2) = 0.08 P Ω(2,1),lonepair (3, 2) = 0.20 P Ω(2,1),lonepair (2, 3) = 0.07 the third case, one electron of spin σ goes from a lone pair to the triple bond, and one electron of spin σ goes from the triple bond to the same lone pair. This contributes to P Ω(2,1),lonepair [START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF][START_REF] Moriarty | The addition reaction of propargyl and acetylene : Pathways to cylic hydrocarbons. 2-nd Joint Meeting of the US Sections of the Combustion Institute[END_REF] in the lone pair and to P Ω(6,1),bond [START_REF] Becke | A simple measure of electron localization in atomic and molecular systems[END_REF][START_REF] Moriarty | The addition reaction of propargyl and acetylene : Pathways to cylic hydrocarbons. 2-nd Joint Meeting of the US Sections of the Combustion Institute[END_REF] in the triple bond. Now, if the latter electron goes to the other lone pair, this fourth case contributes to P Ω(6,1),bond (6, 3), P Ω(2,1),lonepair (1, 2) and P Ω(2,1),lonepair [START_REF] Moriarty | The addition reaction of propargyl and acetylene : Pathways to cylic hydrocarbons. 2-nd Joint Meeting of the US Sections of the Combustion Institute[END_REF][START_REF] Filippi | Quantum Monte Carlo calculations of H2 dissociation on Si(001)[END_REF]. As a conclusion, from the calculated probabilities one can observe the delocalization of two electrons (parallel or anti-parallel) between the lone pairs and the triple bond.

Conclusions

In this paper an implementation of a recent electron localization method is proposed, based upon the probability distributions of the electrons in volumes of space. The volumes can overlap, as opposed to localization domains obtained from the topological analysis of local functions, and this aspect is closer to the quantum mechanical point of view of the system where the electrons are delocalized in space. The exposed method can describe bonds, single electrons and lone pairs, but one can also estimate the delocalization of the electrons. Further work will focus on the analysis of the probabilities in the regions of overlap of these volumes, a better understanding of electron delocalization, and the study of larger systems.

The presented results show that the Variational Monte Carlo method owns very rich information which can be easily used to interpret the behavior of the Nparticle density even in highly correlated computations. As every wave function can be sampled with the Metropolis algorithm, the method presented here can be applied to every wave function (open shell, MCSCF, plane waves, excited states, solids,. . . ). The wave function Ψ T as well as the N -particle density |Ψ T | 2 are unchanged under a unitary transformation of the basis functions. The Nparticle density |Ψ T | 2 , as well as the wave function Ψ T , is unchanged under a unitary transformation of the basis functions. As the computed probabilities depend only on the N -particle density generated by the wave function, this algorithm can be used after the sampling of a trial wave function expressed in localized orbitals, widely used in linear scaling algorithms. Moreover, it is important to point that the calculation of the probabilities can be generalized to Fixed-Node Diffusion Monte Carlo, using the second order approximation to correct the mixed estimators of the probabilities.

Figure 3 :

 3 Figure 3: Slice of the surface separating the volumes which maximize the probability of finding two opposite-spin electrons around Li and around H in LiH.

Figure 4 :

 4 Figure4: Slice of the N 2 molecule: Volumes which maximize the probability of finding two, four or six electrons between the nuclei (Ω(2, 1),bond, Ω(4, 1),bond, Ω(6, 1),bond), volumes corresponding to the cores (Ω(2, 1),core) and volumes corresponding to the lone pairs (Ω(2, 1),lone pair).

This number was chosen for convenience. Taking account of the correlation time would not have changed the values of the averages, but only their error bars.

Acknowledgments

The author gratefully acknowledges Andreas Savin and Michel Caffarel for interesting discussions. Numerical computations have been carried out using QMC-MOL, a Quantum Monte Carlo program for the electronic structure of molecules developed in the laboratory (http://www.lct.jussieu.fr) and using the computers of the Centre de Calcul Recherche et Réseau Jussieu (CCR, Paris 6).