J I Amaro-Estrada

A Scemama

M Caffarel

A Ramírez-Solís

On the stability of Be 3 : A benchmark complete active space self-consistent field + averaged quadratic coupled cluster study

The optimized geometries and binding energies for the linear and triangular isomers of the beryllium trimer have been obtained through benchmark multireference averaged quadratic coupled cluster (AQCC) calculations using very large complete active space SCF (CASSCF) references (12 active electrons in 13 and 14 orbitals). Geometries were optimized with the cc-pV5Z basis, while the binding energies (including counterpoise correction) were obtained with the significantly larger aug-cc-pV5Z basis set. The binding energies (27.3 and 16.3 kcal/mol for the equilateral and linear isomers, respectively) are larger than the previous full CI benchmark values, while the corresponding Be-Be equilibrium distances of 4.101 and 4.088 a.u. are smaller. In view of the near-size consistency character of the CASSCF + AQCC method, the fact that all 12 electrons are fully correlated, the active reference space includes 14 orbitals, and the very large basis set used here, we propose to consider these results as reference data for Be 3 . Using the electron pair localization function obtained at the CASSCF(12,15) level, it is clearly illustrated that the 2p orbitals lying in the molecular plane play a dominant role in the bonding pattern for the equilateral isomer.

I. INTRODUCTION

The structure and binding energies of the beryllium dimer and trimer have been widely studied, both theoretically and experimentally. The dimer is experimentally known to be very weakly bound and, since it contains only eight electrons, the accurate theoretical description of its ground (singlet) state potential surface has been the subject of many works since the early days of quantum chemistry. A large number of studies have shown that the bonding in the dimer is actually due to contributions of the hybridization of the 2s and 2p orbitals, and the non-dynamic correlation effects that arise are, therefore, crucial to properly define zeroth-order electronic wavefunctions. Many electronic structure methods ranging in complexity (perturbational up to fourth order, singleand multi-reference variational CI, density functional theory (DFT)-based, coupled cluster, and even full CI (FCI)) have been applied to study Be 2 (Refs. 3-5, 8-10, 17, 18, 21, 24-26) and Be 3 (Refs. 4, 5, 8, 10-13, 15, 20-24, 27, 28) using, also, many types of atomic basis sets. In the case of the Be dimer, a number of quantum Monte Carlo (QMC) approaches have also been used (see Refs. 30-36). Since Be exists in the bulk as a metallic solid, it is natural to raise the question of the evolution of the binding energy (BE) per atom as the number of atoms grows. For larger Be clusters, it has been found that, in spite of the fact that no d orbitals are involved in the bond shells, as mentioned by Junquera et al., [START_REF] Junquera | [END_REF] the presence of the unoccupied 2p orbitals, which are quasi-degenerate with a) Electronic mail: alex@uaem.mx. the 2s orbitals, allows the formation of Be clusters by this sp2 hybridization process. Actually, it has been shown that this hybridization grows with the number of Be atoms in the cluster. 5 Kaplan et al. 21 have argued that the binding in Be 3 has a mixed physical (van der Waals) and chemical (nonadditive exchange) nature, while the 2s-2p quasi-degeneracy leads to strong multiconfigurational character in the description of the wavefunctions for both dimer and trimer (see, for instance, Refs. 8, 12, 13, 18, 19, 24). Lee et al. 11 made a careful comparison of the electronic description at the MRCI and CCSD levels for Be n (n = 3-5) clusters. They concluded that the coupled cluster approach restricted to single and double excitations is not capable of quantitative accuracy in the description of the electronic structure for Be clusters, since this method underestimates the dissociation energies.

For the Be 2 van der Waals dimmer, it is well known, both theoretically 3-5, 8-10, 16-18, 21, 24-26, 30-36 and experimentally, 6,7,14,19,29 that the interaction energy is quite small, around 2.4 kcal/mol, while the Be-Be distance is rather large, around 2.4 Å. This rather small binding energy for the dimer is in contrast with the non-negligible binding energy per bond (around 13 kcal/mol) for the bulk. 11 So, it has been proposed that the existence of Be in the bulk can be explained through an increase of the BE for larger Be n clusters. Based on the previous idea, density functional techniques were applied to address beryllium clusters Be n (n = 2-8) by Beyer et al., 20 using a large 6-311++G(3df) basis set to study the structural and energetic properties as functions of cluster size and to determine their convergence towards metallic beryllium.

Unlike Be 2 , to the best of our knowledge, for the trimer there are no experimental values to compare with. However, there is a large number of theoretical works using many types of electronic structure methods. 4, 5, 8, 10-13, 15, 20-24, 27, 28 All these studies highlight the importance of an accurate treatment of both dynamic and non-dynamic correlation effects, coupled to the use of extended atomic basis sets including higher l-angular momentum orbitals. For Be 3 , both the geometry and the binding energy are known to be strongly dependent on such factors. Here, our goal is to present the results of new benchmark-type multireference variational calculations for the linear and triangular isomers of Be 3 using the largest augmented correlation-consistent Gaussian atomic basis sets defined for the beryllium atom. Section II deals with the theoretical method and the computational details, including the treatment of the basis set superposition error (BSSE) which is known to be important for the binding energy. In the same section, we recall the fundamentals of the electron pair localization function (EPLF) we will be using to study the binding pattern of the equilateral isomer. In Sec. III, we present the structural and energetic results compared with previous values. Finally, Sec. IV presents the conclusions and some perspectives.

II. METHOD AND COMPUTATIONAL DETAILS

A. Basis sets and non-dynamic/dynamic correlation treatments

Since we are interested in improving the most accurate values of the binding energy presently available for the Be trimer, some important aspects must be taken into account. First, in view of the strong dependence of the binding energy on the largest l-angular momentum atomic functions used, the highest-quality available basis sets have to be chosen. Thus, we have employed two of the most sophisticated atomic Gaussian basis sets existing for this atom, namely, the cc-pV5Z and aug-cc-pV5Z (AV5Z) basis sets of Dunning et al. [START_REF] Dunning | MOLPRO internal cc-pV5Z Be basis set from[END_REF]38 The latter correlation-consistent 15s9p5d4f3g2h basis set is contracted to [7s6p5d4f3g2h] and leads here to a molecular basis of 381 orbitals for 12 electrons in Be 3 . Second, as stressed by many authors, both non-dynamic and dynamic correlation effects are crucial to properly describe the atomic interactions in Be 3 (see, for instance, Refs. 8, 11, 12, 23, 24, 27, 28). In particular, the Be-Be bonds in Be 3 must be described using a zeroth-order wavefunction that includes, at least, the 2p orbitals that lie in the molecular plane. It is well known that the core-valence correlation effects are crucial to study this system; 8,17,24,27,28 therefore, we have used the complete active space self-consistent field (CASSCF) approach to build zeroth-order wavefunctions considering four active electrons per atom, thus including all 12 active electrons of the trimers. Since the 2p shells are nearly degenerate with the 2s one, the ideal chemically relevant active space would consist of five orbitals per atom (1s, 2s, 2p x , 2p y , and 2p z ), thus leading to reference wavefunctions built from CASSCF (12,15) calculations. However, with the present basis set, this ideal active orbital space yields zeroth-order wavefunctions containing more than 1.6 × 10 6 configuration state functions (CSFs) in the D 3h case; this reference space leads to computationally intractable MRCI calculations involving more than 9 × 10 10 uncontracted CSFs. Therefore, we have used a two-step procedure that will be explained below. The CASSCF wavefunctions were used as references to perform averaged quadratic coupled cluster (AQCC) (Ref. 39) calculations which, by construction, yield size-consistent energies. Here, it is crucial to highlight that this method has been shown to yield potential energy surfaces parallel to FCI ones with a much reduced computational cost. 17

B. Geometry optimizations at the CASSCF(12,12) + AQCC level

In order to perform the geometry optimizations, since the Be-Be bonding occurs mainly in the plane defined by the three atoms (here chosen as the xy plane), a reasonable choice is to exclude the three 2p z orbitals from the active space, thus obtaining CASSCF (12,12) references for both trimers. For the triangular molecule, the geometry optimizations have been performed under the C 2v point group (which includes the more symmetric D 3h group). The optimizations of the linear cluster (D ∞h structure) have been performed in the D 2h finite group. The geometries of both trimers were optimized using the automatic numerical algorithm at the AQCC level using the CASSCF (12,12) references and the cc-pV5Z basis set 44 with the MOLPRO-2006.1 code. [START_REF] Molpro | a package of ab initio programs designed[END_REF] Note that even with this basis set, which yields 273 molecular orbitals, the AQCC calculations include 4 024 119 522 (27 890 603) uncontracted (contracted) CSFs. Therefore, using as starting point, the optimized geometries reported by Junquera et al. [START_REF] Junquera | [END_REF] for the triangular cluster and from Vetere et al. 24 for the linear isomer, the CASSCF(12,12) + AQCC/cc-pV5Z optimizations required each more than a week of CPU time allocating up to 18 GB of RAM.

C. CASSCF + AQCC/aug-cc-pV5Z calculations

Once the optimal geometries were determined at the CASSCF (12,12) + averaged coupled pair functional (ACPF)/cc-pV5Z level, we proceeded to perform single point calculations using the larger aug-cc-pV5Z basis set considering the largest possible active reference space that our computational resources allowed, i.e., 12 electrons in 13 and 14 active MOs. In the latter case, this means that only one of the molecular orbitals involving the {2p z (1), 2p z (2), 2p z (3)} set of AOs was excluded from the optimal reference space, thus having the lowest natural occupation number in the CASSCF (12,15) calculation. Although there are several variational CI-based studies for Be 3 , 8,[START_REF] Junquera | [END_REF]24,27 we note that the largest FCI calculations done to date are those reported in 2009 by Vetere et al., 24 where the dimension of their CI spaces was 344 545 308 Slater determinants with the 5s3p2d atomic natural orbital (ANO) basis set (with 72 molecular orbitals for Be 3 ), but these authors kept the six 1s core electrons frozen and only studied the linear trimer. The final variational space of the present AQCC/AV5Z calculation using the CASSCF (12,14) As has been noted before, both the dimer 17 and the trimer, 8,[START_REF] Junquera | [END_REF]24,27 are strongly multiconfigurational due to the fact that the nearly degenerate 2p orbitals play an important role in the bonding pattern of these species. As an example of this, Table II presents the composition of the optimized CASSCF(12,15)/AV5Z wavefunction of the D 3h trimer, where a rather large number of CSFs appear with CI coefficients larger than 0.05 in absolute value; in particular, note the very low 0.86*0.86 = 0.79 weight of the HF configuration.

D. The binding energies and inclusion of the counterpoise correction

It is well known that in order to obtain accurate binding energies for Be n clusters, the counterpoise (CP) correction 41 has to be taken into account, since the basis set superposition error plays a very important role. 17 Therefore, we shall report BE obtained with and without the CP corrections (with respect to 3Be) for both trimers. The CP correction can be automatically calculated using the Molpro algorithm; however, for the present case, there is a technical limitation that arises due to the high symmetry of these Be 3 isomers. With the MOLPRO program, it is not possible to make the CP correction for symmetry-equivalent ghost atoms, so we had to reduce the symmetry from C 2v to use the C s point group by using slightly modified geometries involving changes of 10 -4 a.u. in the optimized x-y coordinates of both isomers. For methodological coherence, the CP-corrected binding energies were calculated using the CASSCF(4,13/14) + AQCC CP-corrected atomic reference energies and the energies of the optimized trimers at the CASSCF(12,13/14) + AQCC/AV5Z level, with the usual definition, BE(D 3h ) = E(D 3h ) -3 E(Be) for the equilateral isomer and BE(linear) = E(linear) -2 E(Be*) -E(Be c ), where Be c and Be* are the central and end atoms of the linear isomer, respectively. For instance, the CP-corrected atomic energies are -14.6476853 and -14.6476855 a.u. at the CASSCF(4,13) + AQCC and CASSCF(4,14) + AQCC levels using the optimized D 3h geometry for the ghosts. For comparison, the atomic CP-uncorrected AQCC/AV5Z energy is -14.646173 a.u. using the 1s,2s,2p x ,2p y ,2p z orbitals as active in the CASSCF(4,5) atomic reference; this energy was used to obtain the CP-uncorrected BE for the trimers at the CASSCF(12,13) + AQCC and CASSCF (12,14) + AQCC levels.

E. The electron pair localization function

Let us briefly present the electron pair localization function. The EPLF is a local scalar function defined in the ordinary 3D space, bounded above and below, which focuses essentially on the localization of electron pairs. It is a good descriptive tool for chemical bonds, since pairs of electrons play a central role in our everyday interpretation of chemical structure and reactivity (Lewis model, VSEPR). The original framework proposed to calculate such a localization function is that of quantum Monte Carlo approaches, 42

EPLF(r) = d σ σ (r) -d σ σ (r) d σ σ (r) + d σ σ (r) , (1)
where

d σ σ (r) = ψ| N i=1 δ(r -r i ) min j =i;σ i =σ j r ij |ψ (2)
is the average of the shortest distance between two electrons with the same spin, one "reference" electron being located at r, and

d σ σ (r) = ψ| N i=1 δ(r -r i ) min j ;σ i =σ j r ij |ψ (3)
is the average of the shortest distance between two electrons with opposite spins, one "reference"electron being located at r. By definition, the EPLF takes its values within the interval [-1,1]. It gives a local indicator of electron pairing as follows. In regions of space where electrons are unpaired, the average shortest distance between spin-like and spin-unlike electrons are similar, and the EPLF goes to zero. When antiparallel electrons are paired, EPLF takes positive values and when parallel electrons are paired, EPLF takes negative values. We have used this approach with great success to study the bonding in tetraoxygen 43 and the complex intermolecular interaction of a lithium atom with a thiophene ring. 44 Quite recently, an approximation of Eqs. ( 2) and ( 3) has been proposed in order to allow the analytical evaluation of the EPLF for the standard forms of wavefunctions employed in computational chemistry without resorting to QMC simulations. 45 It has been shown that the modified EPLF leads to nearly identical shapes when compared to the original EPLF with the advantage of having much faster calculations while avoiding the presence of statistical noise inherent to any Monte Carlo method. The modified formula can be implemented with HF-, CASSCF-, or DFT-based representations by employing any type of basis set. The new EPLF form is based on the following modified definition of the average of the shortest distances:

d σ σ (r) = - 1 γ (r) ln ψ| N i=1 δ(r-r i ) N j =i;σ i =σ j e -γ (r)|r i -r j | 2 |ψ d σ σ (r) = - 1 γ (r) ln ψ| N i=1 δ(r -r i ) N j ;σ i =σ j e -γ (r)|r i -r j | 2 |ψ (4)
where γ (r) is an r-dependent exponent directly related to the local electronic density (see details in Ref. 45). In order to study the bonding pattern of the lowest energy isomer of Be 3 , we have calculated the EPLF using the most dominant part of our best reference wavefunction, namely, the function issued from the CASSCF(12,15) calculation. The construction of the EPLF used all the determinants (around 30 000) whose CI coefficients are larger than 10 -5 . A discussion of the 3D plots of the EPLF will be presented below.

III. RESULTS AND DISCUSSION

A. Structures of both Be 3 isomers

We start by noting the dramatic effect on the reduction of the Be-Be distances when the 1s core electrons are included in the reference space as active; for instance, the optimal Be-Be distance of the linear isomer goes from 4.16 to 4.09 a.u., while for the triangular structure it is reduced from 4.17 to 4.10 a.u. at the CASSCF (12,12) + AQCC/cc-pV5Z level. Table III presents the optimized geometries for the linear and triangular trimers along with previous results at different levels of theory, most of which kept the 1s core elec- trons frozen. For the latter case, we note that even though the optimizations were done within the more general C 2v group, the optimal structure turned out to be an equilateral triangle (D 3h point group), as found before using other methods, like MP2, MP4, MRCI, and FCI. We can compare the present equilibrium distances with previous results. Sudhakar et al. 15 report an equilibrium MP2/6-3ll+G* bond length of 4.23 a.u. for the D 3h isomer. Using a larger basis set, Harrison and Handy 8 obtained a Be-Be distance of 4.22 a.u. through the MRCI/[7s4p2d] method. Junquera et al. [START_REF] Junquera | [END_REF] using FCI with the [3s2p1d] ANO basis set obtained an equilibrium distance of 4.29 a.u. for the equilateral isomer. Rendell et al. 12 reported much shorter bond lengths of 4.22 a.u.

[CCSD(T)] and 4.20 a.u. at the MRCI level, but note that their basis set (5s3p2d1f) already includes one f function. The same authors have stressed that both basis set improvements and electron correlation (coupled-cluster, MRCI) shorten the bond length and increase the binding energy. The shortest Be-Be distance reported for this isomer, 4.08 a.u., was obtained by Beyer et al. 20 with a DFT-based approach, using the B3LYP/6-311++G(3df) method. Lee et al. 11 obtained R e = 4.199 a.u. by MRCI with a [5s3p2d1f] basis set and their MRCI calculations used a threshold of 0.025 for the CI coefficient of the reference configurations. More recently we note the active space CCSDt and MRCI calculations of Ref. 27, performed in bases cc-pVTZ and cc-pVQZ, respectively. The active-space CCSDt calculations were an extension of an earlier study 28 and reproduce quite well the FCI results of Junquera et al.; [START_REF] Junquera | [END_REF] in particular, they report a Be-Be distance of 4.20 a.u for the equilateral isomer. However, in comparison with our 5-active orbital per atom procedure, we stress that they kept the 1s core orbitals frozen, thus neglecting all the core-core and core-valence correlation effects.

For the linear isomer, there are far fewer studies; however, we point out the recent work by Vetere et al. 24 using the FCI method to address linear beryllium chains. For the linear Be 3 cluster, they obtained an equilibrium distance of 4.22 a.u. with the 5s3p2d ANO basis. However, note here that with the 4s3p1d basis set, the FCI distance is significantly longer (4.28 a.u.), thus revealing the very strong dependence of the equilibrium distance on the basis set quality. The newest highly correlated results come from Klos et al., 27 where they obtain a Be-Be distance of 4.15 a.u. for the linear isomer at the CCSDt/cc-pVTZ level but, again, neglecting all the core-core and core-valence correlation effects. The present CASSCF (12,12)+AQCC/cc-pV5Z optimized geometry yields R e = 4.088 a.u. and it was obtained with a significantly larger basis set including much higher l-angular momentum functions up to 3f2g1h. The crucial factor here is that our reference CASSCF wavefunctions were built considering the six 1s core electrons as active, while in the calculations of Refs. 23, 27, and 28 these were kept frozen, i.e., the three lowest MO are doubly occupied in all the configurations of their FCI and CCSDt expansions.

B. Binding energies

The abundant literature on Be 3 shows that the binding energy is extremely dependent upon the type of correlation treatment and the quality of the basis sets employed. Let us first begin with the triangular isomer. Sudhakar and Lammertsma 15 reported a detailed analysis with the 6-31G* and 6-31+G* basis sets at various levels of electronic correlation using coupled cluster (CC) and Möller-Plesset theories; their dissociation energies were 31.4 kcal/mol at MP2/6-31G* and 18.7 kcal/mol at MP4/6-31G* which illustrate a large overestimation of correlation effects with MP2, while their 21.6 kcal/mol value at MP4/6-311+G* shows the dependence on basis set quality. Rendell et al. 12 report a D e of 22.5 kcal/mol at the MRCI/[5321] level and 20.6 kcal/mol after correction for harmonic zero-point vibrational energy contributions. They found a much smaller binding energy of 11.3 kcal/mol at the CCSD level with the same basis set, but a value of 20.4 kcal/mol when the triple excitations are included at the CCSD(T) level. Watts et al. 13 reported 16.2, 11.3, and 13.9 kcal/mol values using fourth-order many-body perturbation theory (MBPT)(4), CCSD(T), and MRCI respectively, all using a 4s2p1d ANO basis set. However, their MRCI calculations used the CASSCF reference that arises with the 2s and 2p orbitals as active, keeping the three 1s orbitals frozen, thereby excluding all the core-valence correlation effects. Harrison and Handy 8 calculated a MRCI value of 19 kcal/mol with a 7s4p2d basis set, but they estimated a true binding energy of 24 ± 2 kcal/mol by including extrapolation of higher excitations and basis set extensions. Junquera et al., [START_REF] Junquera | [END_REF] through full CI calculations using the 3s2p1d ANO basis set, reported a value of 17.2 kcal/mol; however, they also kept the three 1s orbitals frozen (doubly occupied) in all the determinants of their FCI wavefunction. The newest estimations of the binding energy for Be 3 come from the active space CCSDt/cc-pVTZ calculations of Klos et al. 27 where they report BE of 21.2 and 12.63 kcal/mol for the equilateral and linear isomers, respectively; however, we emphasize that they also kept the 1s core orbitals frozen, thereby neglecting the core-core and core-valence correlation effects. Our CASSCF (12,14)+AQCC/AV5Z binding energy is 29.5 kcal/mol but, as expected, the BSSE is non-negligible (around 10% of the BE) so that when the CP correction is applied, a value of 27.3 kcal/mol is obtained. This new benchmark value is larger than all previous estimates, but it remains close to the upper limit of the prediction made by Harrison and Handy 8 twenty-five years ago concerning the exact binding energy.

For the linear cluster, Rendell et al. 12 report an MRCI dissociation energy of 8.8 kcal/mol with the [5s3p2d1f] basis set. Vetere et al. 24 performed systematic FCI calculations using increasingly larger ANO basis sets. They found that the smaller basis sets lead to dissociative potential curves and that only the largest basis sets yield stable linear Be 3 clusters. Their FCI/[5s3p2d] dissociation energy is 11.4 kcal/mol 46 without taking into account the basis set superposition error, but we recall that their FCI calculations kept the 1s orbitals frozen.

Our CASSCF (12,14)+AQCC/AV5Z dissociation energy is 17.9 kcal/mol and, even when the CP correction is applied, it remains slightly larger (16.3 kcal/mol) than previous estimates.

C. Analysis of the bonding pattern in the equilateral isomer through the EPLF

Since the EPLF is a scalar field in ordinary space, there are two easy ways to visualize it. The first one is to plot the points in space where the EPLF takes on a specified value (called the threshold in what follows). Of course the shape and extension of the resulting isosurface will change with varying thresholds. Figures 1 and2 show the surfaces of the EPLF obtained with two isosurface values. In both isosurfaces, the three 1s 2 couples of spin-paired electrons appear as spheres centered around each nucleus. However, it is quite interesting to note that the first isosurface (at 0.011) also reveals a larger "ring" which completely surrounds the nuclear skeleton and has three radial arms between each couple of Be atoms. If the threshold is slightly increased to 0.0113, the isosurface varies, the surrounding ring disappears, and the regions between the Be atoms become more dominant.

IV. CONCLUSIONS AND PERSPECTIVES

In this work, we have reported the results of benchmark type multireference variational CI calculations. The optimal geometries and binding energies for the linear and triangular isomers of the beryllium trimer have been obtained with very large basis sets and considering the six 1s core electrons as active. The calculations have been done using a twostep procedure. First, starting from the best previously reported geometries, we have reoptimized the geometries of both isomers at the AQCC level using very large complete active space SCF references with 12 active electrons in 12 orbitals with the cc-pV5Z basis set. Then, using these optimized FIG. 1. EPLF of the equilateral Be 3 isomer. EPLF calculated with all the determinants in the CASSCF (12,15) wavefunction whose CI coefficients are larger than 10 -5 . Isosurface value of 0.011. geometries, we have performed single-point AQCC/aug-cc-pV5Z calculations using even larger complete active space SCF references with 12 active electrons in 13 and 14 orbitals. The binding energies were counterpoise (CP) corrected to account for the basis set superposition errors. The present CPcorrected binding energies, 27.3 and 16.3 kcal/mol for the equilateral and linear isomers, are larger than the previous FCI benchmark values. The equilibrium distances of 4.101 a.u. for the equilateral isomer and 4.088 a.u. for the linear cluster are slightly shorter than the previous best estimates. These effects are mainly due to the inclusion of the six 1s core electrons as active in our AQCC benchmark calculations.

To visualize the nature of the bonding for the D 3h isomer, we have calculated the EPLF function and the plots clearly show the crucial role played by the 2p orbitals lying in the molecular plane.

Since the AQCC method has been shown to yield parallel potential curves to the FCI ones for Be clusters and since FIG. 2. EPLF of the equilateral Be 3 isomer. EPLF calculated with all the determinants in the CASSCF (12,15) wavefunction whose CI coefficients are larger than 10 -5 . Isosurface value of 0.0113.

we have included the six 1s core electrons as active (which leads to the inclusion of up to 4.7 × 10 10 CSFs), we believe that the geometries and energies presented here are the most accurate to date and, thus, should be used as reference data. The MRCI calculations based on the ideal chemically relevant reference space (namely, all 12 electrons in 15 active orbitals) using the very large AV5Z basis set are still unfortunately too large to be feasible with present-day computational resources. Finally, let us note that work is under way to obtain the best estimate of the binding energy for Be 3 through fixed-node diffusion quantum Monte Carlo calculations using the present CASSCF (12,15) reference as trial wavefunction at the optimized geometry with the cc-pV5Z basis set.

TABLE I .

 I Details of the CASSCF and ACPF calculations for the 1 A 1 state of triangular Be 3 with the aug-cc-pV5Z basis set. The numbers in parentheses correspond to the A 1 , B 1 , B 2 , and A 2 irreps of the C 2v group, respectively. Note that CI energy convergence to 10 -6 a.u. took 10 CPU days allocating 32 GB of RAM for the D 3h isomer. TableIpresents the details of the CASSCF and AQCC variational spaces for this isomer.

		CASSCF(12,13)	CASSCF(12,14)
	Active orbitals	13 (7 5 1 0)	14 (7 5 2 0)
	Slater determinants	740 720	2 262 121
	CSFs in reference	198 598	545 857
	Averaged quadratic coupled cluster calculations with CASSCF references
	Internal configs	198 598	545 857
	Contracted configs	94 763 087	240 820 790
	Uncontracted configs	19 875 254 627	46 667 364 646
	cluster. For the linear isomer, the CASSCF(12,14) includes
	277 697 CSFs, which generate variational AQCC spaces
	with 23 268 073 773 (120 696 813) uncontracted (contracted)
	CSFs.

TABLE II .

 II List of CSF in second quantized notation with CI coefficients larger than 0.05 in the CASSCF(12,15) wavefunction of Be 3 (D 3h ) at the equilibrium geometry.

		C 2v Irreps
	A 1	B 1	B 2	A 2	CI coefficient
	2222000	22000	000	0	0.8650686
	2222200	20000	000	0	-0.1771835
	2220020	22000	000	0	-0.1265278
	2222000	20200	000	0	-0.1015432
	2220000	22200	000	0	-0.0951988
	2220200	22000	000	0	-0.0948863
	222bb00	2aa00	000	0	-0.0938290
	222aa00	2bb00	000	0	-0.0938290
	2220ba0	22000	000	0	0.0904947
	2220ab0	22000	000	0	-0.0904947
	2220000	22000	200	0	-0.0887784
	2222ba0	20000	000	0	-0.0883776
	2222ab0	20000	000	0	0.0883776
	222ab00	2ba00	000	0	0.0863489
	222ba00	2ab00	000	0	0.0863489
	2222020	20000	000	0	-0.0762465
	2202000	22000	200	0	-0.0728661
	2220220	20000	000	0	0.0593958
	2220200	20200	000	0	0.0568272
	22ba000	22000	ab0	0	-0.0552713
	22ab000	22000	ba0	0	-0.0552713
	222a0b0	2ab00	000	0	-0.0545349
	222b0a0	2ba00	000	0	-0.0545349

TABLE III .

 III Be 3 optimized AQCC/cc-pV5Z geometries (a.u.) and AQCC/aug-cc-pV5Z energies (a.u.). Binding energies in kcal/mol. BE(CP) means counterpoise corrected values. Energies are reported using the CASSCF(12,13) and CASSCF(12,14) reference spaces (see text). Previous values in parentheses. BE are calculated using the CASSCF(4,5) + AQCC atomic energy; for details on the CP corrections see the text.

		CASSCF(12,13)	CASSCF(12,14)
	D 3h	R e = 4.101 (4.20 b , 4.21 c , 4.29 d , 4.22 e , 4.38 f , 4.08 g , 4.20 h )
	E(CASSCF) -43.825969	-43.841233
	E(AQCC)	-43.984548	-43.985563 (-43.882330 b )
	BE	28.88	29.51
	BE(CP)	26.67	27.29 (22.5 b , 20.4 c , 17.2 d , 19 e , 13.9 f , 21.2 h )
	Linear	R e = 4.088 (4.220 a , 4.26 b , 4.16 h )
	E(CASSCF) -43.832138	-43.838818
	E(AQCC)	-43.966717	-43. 967047 (-43.874507 a , -43.868969 b )
	BE	17.69	17.90
	BE(CP)	16.10	16.32 (11.4 a , 8.84 b , 12.63 h )
	a From Ref. 24, FCI with [5s3p2d] ANO basis.

b From Ref. 12, MRCI/[5s3p2d1f]. c From Ref. 12, CCSD(T)/[5s3p2d1f]. d From Ref. 23, FCI with [3s2p1d] ANO basis set. e From Ref. 8, MRCI/[7s4p2d]. f From Ref. 13 MRCI/[4s2p1d]. g From Ref. 20 B3LYP/6-311++G(3df). h From Ref. 27, active space CCSDt/cc-pVTZ.

Downloaded 03 Jun 2013 to 130.120.228.84. This article is copyrighted as indicated in the abstract. Reuse of AIP content is subject to the terms at: http://jcp.aip.org/about/rights_and_permissions

This value was obtained using the energy of the cluster reported at 4.25 a.u. (see TableIIIof Ref. 24), while their optimal Be-Be distance is 4.22 a.u., thus yielding a lower bound value.Downloaded 03 Jun 2013 to 130.120.228.84. This article is copyrighted as indicated in the abstract. Reuse of AIP content is subject to the terms at: http://jcp.aip.org/about/rights_and_permissions

ACKNOWLEDGMENTS A.R.-S. thanks support from CONACYT Project No. 130931. The present work summarizes the M.Sc. thesis of J.I.A.-E. who acknowledges a CONACYT graduate scholarship. We thank unlimited CPU time from the UAEM supercomputing center and from CALMIP-Université de Toulouse.


Anthony Scemama

Emmanuel Giner

An efficient implementation of Slater-Condon rules

come

An efficient implementation of Slater-Condon rules

In this work we consider wave functions Ψ expressed as linear combinations of Slater determinants D of orthonormal spin-orbitals φ(r):

Ψ = i c i D i (1)
Using the Slater-Condon rules, [START_REF] Slater | The theory of complex spectra[END_REF][START_REF] Condon | The theory of complex spectra[END_REF] the matrix elements of any one-body (O 1 ) or two-body (O 2 ) operator expressed in the determinant space have simple expressions involving one-and two-electron integrals in the spin-orbital space. The diagonal elements are given by:

D|O 1 |D = i∈D φ i |O 1 |φ i (2)
D|O 2 |D = 1 2 (i,j)∈D φ i φ j |O 2 |φ i φ j - φ i φ j |O 2 |φ j φ i
For two determinants which differ only by the substitution of spin-orbital i with spin-orbital j: For two determinants which differ by two spinorbitals:

D|O 1 |D j i = φ i |O 1 |φ j (3)
D|O 1 |D jl ik = 0 (4)
D|O 2 |D jl ik = φ i φ k |O 2 |φ j φ l -φ i φ k |O 2 |φ l φ j
All other matrix elements involving determinants with more than two substitutions are zero. An efficient implementation of those rules requires:

1. to find the number of spin-orbital substitutions between two determinants 2. to find which spin-orbitals are involved in the substitution 3. to compute the phase factor if a reordering of the spin-orbitals has occured This paper proposes an efficient implementation of those three points by using some specific bit manipulation instructions at the CPU level.

Algorithm

In this section, we use the convention that the least significant bit of binary integers is the rightmost bit.

As the position number of a bit in an integer is the exponent for the corresponding bit weight in base 2, the bit positions are numbered from the right to the left starting at position 0. To be consistent with this convention, we also represent the arrays of 64-bit integers from right to left, starting at position zero. Following with the usual notations, the spin-orbitals start with index one.

Binary representation of the determinants

The molecular spin-orbitals in the determinants are ordered by spin: the α spin-orbitals are placed before the β spin-orbitals. Each determinant is represented as a pair of bit-strings: one bit-string corresponding to the α spin-orbital occupations, and one bit-string for the β spin-orbital occupations. When the i-th orbital is occupied by an electron with spin σ in the determinant, the bit at position (i -1) of the σ bitstring is set to one, otherwise it is set to zero.

The pair of bit-strings is encoded in a 2dimensional array of 64-bit integers. The first dimension contains N int elements and starts at position zero. N int is the minimum number of 64-bit integers needed to encode the bit-strings:

N int = N MOs /64 + 1 (5)
where N MOs is the total number of molecular spinorbitals with spin α or β (we assume this number to be the same for both spins). The second index of the array corresponds to the α or β spin. Hence, determinant D k is represented by an array of N int 64-bit integers

I k i,σ , i ∈ [0, N int -1], σ ∈ {1, 2}.

Finding the number of substitutions

We propose an algorithm to count the number of substitutions between two determinants D 1 and D 2 (algorithm 1). This number is equivalent to the degree of excitation d of the operator Td which transforms

D 1 into D 2 (D 2 = Td D 1
). The degree of excitation is equivalent to the number of holes created in D 1 or to the number of particles created in D 2 . In algorithm 1, the total number of substitutions is calculated as the sum of the number of holes created in each 64-bit integer of D 1 .

On line 4, (I 1 i,σ xor I 2 i,σ ) returns a 64-bit integer with bits set to one where the bits differ between I 1 i,σ and I 2 i,σ . Those correspond to the positions of the holes and the particles. The popcnt function returns the number of non-zero bits in the resulting integer. At line 5, two d contains the sum of the number of holes and particles, so the excitation degree d is half of two d.

The Hamming weight is defined as the number of non-zero bits in a binary integer.The fast calculation of Hamming weights is crucial in various domains of computer science such as error-correcing codes [START_REF] Hamming | Error detecting and error correcting codes[END_REF] or cryptography [START_REF] Hilewitz | Comparing fast implementations of bit permutation instructions[END_REF]. Therefore, the computation of Hamming weights has appeared in the hardware of processors in 2008 via the the popcnt instruction introduced with the SSE 4.2 instruction set. This instruction has a 3-cycle latency and a 1-cycle throughput independently of the number of bits set to one (here, independently of the number of electrons), as opposed to Wegner's algorithm [START_REF] Wegner | A technique for counting ones in a binary computer[END_REF] that repeatedly finds and clears the last nonzero bit. The popcnt instruction may be generated by Fortran compilers via the intrinsic popcnt function.

Identifying the substituted spinorbitals

Algorithm 2 creates the list of spin-orbital indices containing the holes of the excitation from D 1 to D 2 . At line 4, H is is set to a 64-bit integer with ones at the positions of the holes. The loop starting at line 5 translates the positions of those bits to spinorbital indices as follows: when H = 0, the index of the rightmost bit of H set to one is equal to the number of trailing zeros of the integer. This number can be obtained by the x86 64 bsf (bit scan forward) instruction with a latency of 3 cycles and a 1-cycle throughput, and may be generated by the Fortran trailz intrinsic function. At line 7, the spin-orbital

for i ← 0 to N int -1 do 4 H ← I 1 i,σ xor I 2 i,σ
and I 1 i,σ ; index is calculated. At line 8, the rightmost bit set to one is cleared in H.

The list of particles is obtained in a similar way with algorithm 3.

Computing the phase

In our representation, the spin-orbitals are always ordered by increasing index. Therefore, a reordering may occur during the spin-orbital substitution, involving a possible change of the phase.

As no more than two substitutions between determinants D 1 and D 2 give a non-zero matrix element, we only consider in algorithm 4 single and double substitutions. The phase is calculated as -1 Nperm , where N perm is the number permutations necessary to bring the spin-orbitals on which the holes are made to the positions of the particles. This number is equal to the number of occupied spin-orbitals between these two positions.

We create a bit mask to extract the occupied spinorbitals placed between the hole and the particle, and we count them using the popcnt instruction. We have to consider that the hole and the particle may or may not not belong to the same 64-bit integer.

On lines 6 and 7, we identify the highest and lowest spin-orbitals involved in the excitation to delimitate the range of the bit mask. Then, we find to which 64-bit integers they belong and what are their bit positions in the integers (lines 8-11). The loop in lines 12-13 sets to one all the bits of the mask con- tained in the integers between the integer containing the lowest orbital (included) and the integer containing highest orbital (excluded). Line 14 sets all the m rightmost bits of the integer to one and all the other bits to zero. At line 15, the n + 1 rightmost bits of the integer containing to the lowest orbital are set to zero. At this point, the bit mask is defined on integers mask[j ] to mask[k ] (if the substitution occurs on the same 64-bit integer, j = k). mask[j ] has zeros on the leftmost bits and mask[k ] has zeros on the rightmost bits. We can now extract the spin-orbitals placed between the hole and the particle by applying the mask and computing the Hamming weight of the result (line 17). For a double excitation, if the realization of the first excitation introduces a new orbital between the hole and the particle of the second excitation (crossing of the two excitations), an additional permutation is needed, as done on lines 18-19. This if statement assumes that the Holes and Particles arrays are sorted.

Optimized Implementation

In this section, we present our implementation in the Fortran language. In Fortran, arrays start by default with index one as opposed to the convention chosen in the algorithms.

A Fortran implementation of algorithm 1 is given in Figure 1. This function was compiled with the Intel A static analysis of the executable was performed using the MAQAO tool [START_REF] Djoudi | MAQAO: Modular assembler quality Analyzer and Optimizer for Itanium 2[END_REF] : if all the data fit into the L1 cache, an iteration of the loop takes in average 3.5 CPU cycles on an Intel Sandy Bridge CPU core. A dynamic analysis revealed 10.5 cycles for calling the function, performing the first statement and exiting the function, and 4.1 CPU cycles per loop iteration. Therefore, to obtain the best performance this function will need to be inlined by the compiler, eventually using compiler directives or inter-procedural optimization flags.

For the identification of the substitutions, only four cases are possible (figure 2) depending of the degree of excitation : no substitution, one substitution, two substitutions or more than two substitutions. If the determinants are the same, the subroutine exits with a degree of excitation of zero. For degrees of exci- tation higher than two, the subroutine returns a degrees of excitation equal to -1. For the two remaining cases, a particular subroutine is written for each case. The cases were ordered from the most probable to the least probable to optimize the branch prediction. In output, the indices of the spin-orbitals are given in the array exc as follows:

• The last index of exc is the spin (1 for α and 2 for β)

• The second index of exc is 1 for holes and 2 for particles

• The element at index 0 of the first dimension of exc gives the total number of holes or particles of spin α or β

• The first index of exc refers to the orbital index of the hole or particle

The subroutine for single excitations is given in figure 3. The particle and hole are searched simultaneously. The ishift variable (line 19) allows to replace the integer multiplication at line 7 of algorithm 2 by an integer addition, which is faster. Line 38 contains a bit shift instruction where all the bits are shifted 6 places to the right. This is equivalent to doing the integer division of high by 64. Line 39 computes high (mod 64) using a bit mask. The compiler may recognize that those last two optimizations are possible, but it can not be aware that high is always positive. Therefore, it will generate additional instructions to handle negative integers. As we know that high is always positive, we can do better than the compiler. The test at line 35 is true when both the hole and the particle have been found. This allows to compute the phase factor and exit the subroutine as soon at the single excitation is obtained. If the hole and particle belong to the same integer (line 42), the bit mask is created and applied on the fly to only one 64-bit integer. Otherwise, the bit mask is created and applied on the two extreme integers j and k. The integers between j and k, if there are any, don't need to have a bit mask applied since the mask would have all bits set to one. Finally, line 52 calculates -1 Nperm using a memory access depending on the parity of N perm : phase dble is an array of two double precision values (line 10).

For double excitations, the subroutine is similar to the previous one. The nexc variable counts how many holes and particles have been found, in order to exit the loop (line 45) as soon as the double excitation is found. The calculation of the phase is the same as the case of single excitations (lines 48-67), but in the case of a double excitation of the same spin, orbital crossings can occur (lines 68-76).

Benchmarks

All the benchmarks were realized on a quad-core Intel Xeon CPU E3-1220 @ 3.10GHz (8 MiB cache). The benchmark program is single-threaded, and was run using a single CPU core at the maximum turbo frequency of 3.4 GHz. The numbers of CPU cycles 1 subroutine get_single_excitation(det1,det2,exc,phase,Nint) 2 implicit none 3 integer, intent(in) :: Nint 4 integer*8, intent(in) :: det1(Nint,2) 5 integer*8, intent(in) :: det2(Nint,2) 6 integer, intent(out) :: exc(0:2,2,2) 7 double precision, intent(out) :: phase 8 integer :: tz, l, ispin, ishift, nperm, i, j, k, m, n, high, low 9 integer*8 :: hole, particle, tmp 10 double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /) 11 12 exc(0,1,1) = 0 13 exc(0,2,1) = 0 14 exc(0,1,2) = 0 15 exc(0,2,2) = 0 16 do ispin = 1,2 1 subroutine get_double_excitation(det1,det2,exc,phase,Nint) 2 implicit none 3 integer, intent(in) :: Nint 4 integer*8, intent(in) :: det1(Nint,2), det2(Nint,2) 5 integer, intent(out) :: exc(0:2,2,2) 6 double precision, intent(out) :: phase 7 integer :: l, ispin, idx_hole, idx_particle, ishift 8 integer :: i,j,k,m,n,high, low,a,b,c,d,nperm,tz,nexc 9 integer*8 :: hole, particle, tmp 10 double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /) 11 exc(0,1,1) = 0 12 exc(0,2,1) = 0 13 exc(0,1,2) = 0 14 exc(0,2,2) = 0 15 nexc=0 16 nperm=0 17 do ispin = cases d ∈ {0, 1}, which are a very small fraction of the the number of cycles is very close to this of n excitations. Note that the calculation of the density matrix only requires N (N +1)/2 determinant comparisons, and this explains why the CPU time is smaller than for the n excitations benchmark.

All the source files needed to reproduce the benchmark presented in this section are available at https: //github.com/scemama/slater_condon.

Summary

We have presented an efficient implementation of Slater-Condon rules by taking advantage of instructions recently introduced in x86 64 processors. The use of these instructions allow to gain a factor larger than 6 with respect to their software implementation. 1: CPU time (seconds) and average number of CPU cycles measured for the n excitations function, the get excitation subroutine and the calculation of the one-electron density matrix in the molecular orbital basis. get excitation was also called using selected pairs of determinants such that the degree of excitation (d) was zero, one, two or higher.

As a result, the computation of the degree of excitation between two determinants can be performed in the order of 10 CPU cycles in a set of 128 molecular orbitals, independently of the number of electrons.

Obtaining the list of holes and particles involved in a single or double excitation can be obtained in the order of 50-90 cycles, also independently of the number of electrons. For comparison, the latency of a double precision floating point division is typically 20-25 cycles, and a random read in memory is 250-300 CPU cycles. Therefore, the presented implementation of Slater-Condon rules will significantly accelerate determinant-driven calculations where the twoelectron integrals have to be fetched using random memory accesses. As a practical example, the oneelectron density matrix built from 10 000 determinants of a water molecule in the cc-pVTZ basis set was computed in 0.2 seconds on a single CPU core.

D|O 2 |D j i = k∈D φ i φ k |O 2

 2 |φ j φ k -φ i φ k |O 2 |φ k φ j * Laboratoire de Chimie et Physique Quantiques, CNRS-IRSAMC, Université de Toulouse, France.

1 arXiv: 1311 Algorithm 1 :

 113111 .6244v1 [physics.comp-ph] 25 Nov 2013 Compute the degree of excitation between D 1 and D 2 . Function n excitations(I 1 , I 2 ) Input: I 1 , I 2 : lists of integers representing determinants D 1 and D 2 . Output: d: Degree of excitation. 1 two d ← 0; 2 for σ ∈ {α, β} do 3 for i ← 0 to N int -1 do 4 two d ← two d+ popcnt (I 1 i,σ xor I 2 i,σ );

Algorithm 2 : 2 2 k

 222 Obtain the list of orbital indices corresponding to holes in the excitation from D 1 to D Function get holes(I 1 , I 2 ); Input: I 1 , I 2 : lists of integers representing determinants D 1 and D 2 . Output: Holes: List of positions of the holes. 1 for σ ∈ {α, β} do ← 0;

3

 3

5 while H = 0 do 6 position 8 H 9 k ← k + 1 ;

 56891 ← trailing zeros(H); 7 Holes[k, σ] ← 1 + 64 × i + position; ← bit clear(H, position); 10 return Holes;

Algorithm 3 : 4 P 5 while P = 0 do 6 position 8 P

 34568 Obtain the list of orbital indices corresponding to particles in the excitation from D 1 to D 2 Function get particles(I 1 , I 2 ) Input: I 1 , I 2 : lists of integers representing determinants D 1 and D 2 . Output: Particles: List of positions of the particles. 1 for σ ∈ {α, β} do 2 k ← 0; 3 for i ← 0 to N int -1 do ← I 1 i,σ xor I 2 i,σ and I 2 i,σ ; ← trailing zeros(P); 7 Particles[k, σ] ← 1 + 64 × i + position; ← bit clear(P, position); 9 k ← k + 1; 10 return Particles;

Algorithm 4 :|O|D 2 Function 8 k 9 m

 4289 Compute the phase factor of D 1 GetPhase(Holes, Particles) Input: Holes and Particles obtained with alorithms 2 and 3. Output: phase ∈ {-1, 1}. 1 Requires: n excitations(I 1 , I 2 ) ∈ {1, 2}. Holes and Particles are sorted. 2 nperm ← 0; 3 for σ ∈ {α, β} do 4 n σ ← Number of excitations of spin σ; 5 for i ← 0 to n σ -1 do 6 high ← max(Particles[i, σ], Holes[i, σ]); 7 low ← min(Particles[i, σ], Holes[i, σ]); ← high/64 ; ← high (mod 64); 10 j ← low/64 ; 11 n ← low (mod 64); 12 for l ← j to k -1 do 13 mask[l] ←not(0); 14 mask[k] ← 2 m -1; 15 mask[j] ← mask[j] and (not(2 n+1 )+1) 16 for l ← j to k do 17 nperm ← nperm+popcnt(I j,σ 1 and mask[l]); 18 if (n σ = 2) and ( Holes[2, σ] < Particles[1, σ] or Holes[1, σ] > Particles[2, σ] ) then 19 nperm ← nperm + 1; 20 return -1 nperm ; Fortran compiler 14.0.0 with the options -xAVX -O2.

Figure 1 :

 1 Figure 1: Fortran implementation of algorithm 1.

Figure 2 :

 2 Figure 2: Fortran subroutine for finding the holes and particles involved in a matrix element.

Figure 3 :

 3 Figure 3: Fortran subroutine for finding the holes and particles involved in a single excitation.

Figure 5 :

 5 Figure 5: Fortran subroutine for the calculation of the one-electron density matrix in the molecular orbital basis.

d ← two d/2;

return d;

were obtained by polling the hardware time stamp counter rdtscp using the following C function: double rdtscp_(void) { unsigned long long a, d; __asm__ volatile ("rdtscp" : "=a" (a), "=d" (d)); return (double)((d<<32) + a); } Two systems were benchmarked. Both systems are a set of 10 000 determinants obtained with the CIPSI algorithm presented in ref [START_REF] Giner | Using perturbatively selected configuration interaction in quantum monte carlo calculations[END_REF]. The first system is a water molecule in the cc-pVTZ basis set [START_REF] Thom | Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen[END_REF], in which the determinants are made of 5 α-and 5 β-electrons in 105 molecular orbitals (N int = 2). The second system is a Copper atom in the cc-pVDZ [START_REF] Balabanov | Systematically convergent basis sets for transition metals. i. all-electron correlation consistent basis sets for the 3d elements sczn[END_REF] basis set, in which the determinants are made of 15 α-and 14 β-electrons in 49 molecular orbitals (N int = 1). The benchmark consists in comparing each determinant with all the derminants (10 8 determinant comparisons). These determinant comparisons are central in determinant driven calculations, such as the calculation of the Hamiltonian matrix in the determinant basis set. As an example of a practical application, we benchmark the calculation of the oneelectron density matrix on the molecular orbital basis using the subroutine given in figure 5.

The programs were compiled with the Intel Fortran Compiler version 14.0.0. The compiling options included inter-procedural optimization to let the compiler inline functions via the -ipo option, and the instruction sets were specified using the -xAVX or the -xSSE2 options. Let us recall that the AVX instruction set includes all the instructions introduced with SSE4.2. Using AVX instructions, the popcnt function is executed through its hardware implementation, as opposed to SSE2 in which popcnt is executed through its software implementation.

In table 1 we report measures of the CPU time and of the average number of CPU cycles of all the subroutines presented in this paper for the water molecule and the Copper atom.

The measures of n excitations and get excitation are made using all the possible pairs of determinants. Most of the time the degree of excitation d is greater than 2, so the average number of cycles has a large weight on the d > 2 case, and this explains why the average number of cycles is much lower than in the d = 1 and d = 2 cases.

As the computation of the one-electron density matrix only requires to find the spin-orbitals for the


Anya Rothenbuhler
email: anya.rothenbuhler@free.fr

Chantal Lotton

Delphine Fradin

Three Common Variants of LEP, NPY1R and GPR54 Show No Association with Age at Menarche

Context: Twin and family studies indicate a significant heritability of pubertal timing and more specifically of age at menarche (AAM). Objective: Test the association of AAM with common variants of three candidate genes suspected to have a prominent role in reproductive physiology: leptin ( LEP ), neuropeptide Y receptor 1 ( NPY1R ) and GPR54 . Design and Methods: We selected the -2459 LEP, the rs7687423 NPY1R and the rs350132 GPR54 variants as the more common coding or regulatory variants (minor allelic frequency 1 0.10) in these gene regions. To avoid stratification problems that can impair association studies, we used the Q-TDT method based on allele transmission to evaluate the relationship of these variants with AAM in 245 healthy women from 107 families of European ancestry. Results: We found no association of AAM with any of the studied variants. Conclusions: Keeping in mind that common variants do not recapitulate the whole genetic variation in a given gene region, this study indicates that the studied LEP, NPY1R and GPR54 variants do not have a major influence upon pubertal timing in Caucasian women. Effects of these genetic loci on age at menarche can definitively be excluded only through determination of extended haplotypes in a larger cohort.

Introduction

The timing of menarche is influenced both by environmental and by genetic factors, each group of factors accounting for approximately half the AAM variation in developed countries [START_REF] Loesch | Genetic correlates of menarcheal age: a multivariate twin study[END_REF] . Genetic determinants of AAM remain unknown, but several association studies have started to unravel positional candidates across associated genomic regions [START_REF] Sedlmeyer | Determination of sequence variation and haplotype structure for the gonadotropin-releasing hormone (GnRH) and GnRH receptor genes: investigation of role in pubertal timing[END_REF][START_REF] Stavrou | Combined estrogen receptor alpha and estrogen receptor beta genotypes influence the age of menarche[END_REF][START_REF] Guo | Polymorphisms of estrogenbiosynthesis genes CYP 17 and CYP19 may influence age at menarche: a genetic association study in Caucasian females[END_REF][START_REF] Lai | Narod SA: CYP gene polymorphisms and early menarche[END_REF][START_REF] Kadlubar | The CYP3A4 * 1B variant is related to the onset of puberty, a known risk factor for the development of breast cancer[END_REF][START_REF] Gorai | Estrogen-metabolizing gene polymorphisms, but not estrogen receptor-alpha gene polymorphisms, are associated with the onset of menarche in healthy postmenopausal Japanese women[END_REF][START_REF] Xita | Association of SHBG gene polymorphism with menarche[END_REF][START_REF] Worda | The influence of Nos3 polymorphisms on age at menarche and natural menopause[END_REF][START_REF] Rothenbuhler | Weight-adjusted genome scan analysis for mapping quantitative trait loci for menarchal age[END_REF] .

Another approach is a direct test of association of AAM with suspected candidate gene polymorphisms. In this respect, the current study tested if the individual variability of AAM is associated with major variants of LEP , NPY1R or GPR54 , genes that are known for their participation to the hypothalamic pathways that regulate the gonadotropic function. In each of these three gene regions, we selected a variant based on two criteria: the variant has to be common in Caucasians (mean allele frequency 1 0.10) and should possibly have a functional effect.

Studies in animal models and in humans support a role of leptin in the physiology of sexual maturation [START_REF] Ahima | Leptin accelerates the onset of puberty in normal female mice[END_REF][START_REF] Popovic | Leptin, nutrition and reproduction: new insights[END_REF][START_REF] Cervero | The role of the leptin in reproduction[END_REF] . Mutations in the LEP [START_REF] Strobel | A leptin missense mutation associated with hypogonadism and morbid obesity[END_REF] and LEPR [START_REF] Clement | A mutation in the human leptin receptor gene causes obesity and pituitary dysfunction[END_REF] genes are associated with hypogonadotropic hypogonadism. The LEP gene is located in 7q31.3, spans 20 kb and contains 3 exons separated by two introns. The -2548 SNP selected for the current study is frequent, is located in the promoter region and was suggested to be regulatory [START_REF] Stunff | A common promoter variant of the leptin gene is associated with changes in the relationship between serum leptin and fat mass in obese girls[END_REF][START_REF] Hoffstedt | A polymorphism in the leptin promoter region (-2548 G/A) influences gene expression and adipose tissue secretion of leptin[END_REF] . Our hypothesis was that through the control of circulat-Rothenbuhler/Lotton/Fradin ing leptin concentration, the LEP gene variant could modulate the variability of pubertal timing independently of fat mass.

Although there are no data implicating NPYR in humans, animal models suggest that NPY1R mediates the hypothalamic effects of NPY, an orexigenic neurotransmitter which influences gonadotropic function [START_REF] El | Neuropeptide Y: A hypothalamic brake restraining the onset of puberty in primates[END_REF] and interferes with leptin effects in NPY neurons [START_REF] Pralong | The neuropeptide Y Y1 receptor regulates leptin-mediated control of energy homeostasis and reproductive functions[END_REF][START_REF] Baskin | Leptin receptor longform splice-variant protein expression in neuron cell bodies of the brain and co-localization with neuropeptide Y mRNA in the arcuate nucleus[END_REF] . The NPY1R gene, located in 4q31.3-q32, is subject to alternative splicing, which makes the functional location of the common variant that we selected difficult to infer (intron, exon 1, or promoter of the NPY1R gene). NPY1R and NPY5R genes are transcribed in opposite directions from a common promoter region suggesting that their transcriptional regulation is common [START_REF] Herzog | Overlapping gene structure of the human neuropeptide Y receptor subtypes Y1 and Y5 suggests coordinate transcriptional regulation[END_REF] . If the studied variant of the NPY1R gene is a functional SNP which influences the transcriptional activity, it could also regulate the expression of NPY5R. NPY5R is known to be implicated in the NPY control of GnRH secretion in prepubertal rats [START_REF] Lebrethon | In vitro stimulation of the prepubertal rat gonadotropin-releasing hormone pulse generator by leptin and neuropeptide Y through distinct mechanisms[END_REF] .

The third common non synonymous SNP that we tested (rs350132, His-Leu 364) is located in exon 5 of GPR54 , a gene that is a key regulator of gonadotropic function [START_REF] De Roux | Hypogonadotropic hypogonadism due to loss of function of the KiSS1-derived peptide receptor GPR54[END_REF][START_REF] Gottsch | A role for kisspeptins in the regulation of gonadotropin secretion in the mouse[END_REF][START_REF] Seminara | The GPR54 gene as a regulator of puberty[END_REF] . While hypogonadotropic hypogonadism has been associated with deletions encompassing exon 5 in GPR54 [START_REF] De Roux | Hypogonadotropic hypogonadism due to loss of function of the KiSS1-derived peptide receptor GPR54[END_REF] , there are no data, however, indicating that this SNP has a functional role in GPR54 signal transduction.

Materials and Methods

Patients

The study was approved by the local Institutional Review Board and informed written consent was obtained for each participant. The study population, based on a sister-pair approach, included 245 healthy women of recent European ancestry from 106 families; 79 pedigrees of two, 23 pedigrees of three, two pedigrees of four and two pedigrees of five sisters volunteered for the study. Growth, body mass index and puberty data were collected from growth charts ('Carnet de Santé') and medical interviews. Age at menarche was recalled from women using both status quo and direct recall methods validated in previous studies [START_REF] Must | Recall of early menstrual history and menarcheal body size: after 30 years, how well do women remember?[END_REF][START_REF] Koprowski | Ability of young women to recall past body size and age at menarche[END_REF] .

Genotyping

Genomic DNA was extracted from whole blood with the Gentra Extraction Kit. Amplification (35 cycles) was carried out in 50 l containing 200 ng of genomic DNA, 1 M of each primer, 1 ! PCR buffer (Invitrogen), 0.2 m M of each dNTP, 1.5 m M MgCl 2 and 1.25 units of Taq polymerase (Invitrogen). The following primers and annealing temperatures were used: LEP forward: TTTCCTGTAATTTTCCCGTGA, LEP reverse: AAAGCAAA-GACAGGCATAAAAA, annealing temperature 50 ° C; NPY1R forward: AAGGCTCTACAGTCTGTACC, NPY1R reverse: GTG-GAGCACACCTGTAGTCC, annealing temperature 58 ° C and GPR54 forward: GCCTTTCGTCTAACCACCTTC, GPR54 reverse: GGAGCCGCTCGGATTCCCAC, annealing temperature 62 ° C. PCR products, 242 bp for LEP , 432 bp for NPY1R and 524 bp for GPR54 , were digested respectively with CfO1 restriction enzyme (Roche) for LEP , Pst1 restriction enzyme (Biolabs) for NPY1R and HpyCH4V restriction enzyme (Biolabs) for GPR54 , then analyzed by gel electrophoresis. For the -2459 LEP variant, two fragments of 181 and 61 bp identified homozygote CC subjects; three fragments of 242, 181 and 61 bp identified CA heterozygotes; and a single 242-bp fragment identified AA homozygotes. For the rs7687423 NPY1R variant, two fragments of 256 and 176 bp identified homozygote TT subjects; three fragments of 432, 256 and 176 bp identified CT heterozygotes, and a single 432 bp fragment identified CC homozygotes. For the rs350132 GPR54 variant, four fragments of 160, 153, 137 and 74 bp identified AA homozygotes; five fragments of 297, 160, 153, 137 and 74 bp identified AT heterozygotes, and three fragments of 297, 153 and 74 bp identified TT homozygotes.

Statistical Analysis

AAM showed a continuous distribution in the cohort. The studied variants followed Hardy-Weinberg equilibrium. To test for both linkage and association with the variation of AAM we used the Q-TDT method implemented in Merlin [START_REF] Abecasis | Merlin -rapid analysis of dense genetic maps using sparse gene flow trees[END_REF] . The orthogonal model of Abecasis et al. was used in our analyses, where the total association is partitioned into orthogonal within and between family components ( ␤ b and ␤ w , respectively). The between family component is specific for each nuclear family and could be confounded by population stratification. However, the within family component is significant only in the presence of linkage disequilibrium caused by close linkage. Thus, a Q-TDT test of the significance of ␤ w based on allelic transmission is robust to population stratification. The null model was 'means = Mu + ␤ b ' and the full model was 'means = Mu + ␤ b + ␤ w '. The analysis was adjusted for BMI and BMI Z score at age at menarche for reasons explained in a previous analysis [START_REF] Rothenbuhler | Weight-adjusted genome scan analysis for mapping quantitative trait loci for menarchal age[END_REF] . Results are expressed as means 8 SD (standard deviation).

Results

Table 1 describes the genotypic and phenotypic characteristics of the 245 women. Mean age of participants was 34.4 8 13 years. Age at menarche averaged 13.1 8 1.5 years as reported in comparable Caucasian women modern populations.

None of the studied polymorphisms of the LEP , NPY1R and GPR54 genes were found to be associated with AAM when each SNP was considered individually. When Q-TDT analyses were corrected for BMI and BMI Z score at age at menarche, results were still nonsignificant. All p values were above 0.2. The p values for the studied variants of LEP , NPY1R and GPR54 were, respectively, 0.9, 0.8 and 0.6 when nonadjusted; 0.2, 0.6 and 0.7 when adjusted to BMI at age at menarche, and 0.7 for the three variants when adjusted to BMI Z score at age at menarche.

Due to the limited size of the sample, our study was underpowered to evaluate potential effects of gene-gene interactions using the combinations of genotypes at the three studied loci. In an attempt to combine information at two loci, we could only compare 11 women with LEP CC and NPY CC genotypes in whom median AAM was 13.2 compared with 13.1 years in the rest of the sample (n = 193).

Discussion

Since AAM averaged 13.1 years with a continuous distribution of values in the studied population sample, we considered that it was representative of modern Caucasian women of Western Europe [START_REF] De Muinich Keizer | Trends in pubertal development in Europe[END_REF] .

Association studies are a priori a powerful means for identifying the common variants that underly complex traits [START_REF] Risch | The future of genetic studies of complex human diseases[END_REF][START_REF] Cardon | Association study designs for complex diseases[END_REF][START_REF] Tabor | Opinion: Candidate-gene approaches for studying complex genetic traits: practical considerations[END_REF] . Our study, however, has several major weaknesses. (1) One is the size of the studied sample which is insufficient for quantitative traits like AAM that are likely to involve polygenic determinants of limited effect size [START_REF] Farrall | Quantitative genetic variation: a post-modern view[END_REF] . (2) Another one is the limited number and arbitrary choice of genes that we tested (not based on positional information from genomewide studies) versus the multiple ones that could be involved in the control of puberty timing. [START_REF] Stavrou | Combined estrogen receptor alpha and estrogen receptor beta genotypes influence the age of menarche[END_REF] The study is based on the belief that the selected variants could be functional, which was proven for none of them, nor were they in linkage disequilibrium with known functional variants. (4) Our strategy of analysis is entirely based on the search for common variants and nothing indicates yet that common variants are the basis for the genotypic variation of evolutionary traits, like AAM.

Several studies that searched for associations of candidate genes with AAM in Caucasian women of North America or Europe, and in Japanese, 'Hispanic Americans' and Asian women have remained negative overall [START_REF] Sedlmeyer | Determination of sequence variation and haplotype structure for the gonadotropin-releasing hormone (GnRH) and GnRH receptor genes: investigation of role in pubertal timing[END_REF][START_REF] Stavrou | Combined estrogen receptor alpha and estrogen receptor beta genotypes influence the age of menarche[END_REF][START_REF] Guo | Polymorphisms of estrogenbiosynthesis genes CYP 17 and CYP19 may influence age at menarche: a genetic association study in Caucasian females[END_REF][START_REF] Lai | Narod SA: CYP gene polymorphisms and early menarche[END_REF][START_REF] Kadlubar | The CYP3A4 * 1B variant is related to the onset of puberty, a known risk factor for the development of breast cancer[END_REF][START_REF] Gorai | Estrogen-metabolizing gene polymorphisms, but not estrogen receptor-alpha gene polymorphisms, are associated with the onset of menarche in healthy postmenopausal Japanese women[END_REF][START_REF] Xita | Association of SHBG gene polymorphism with menarche[END_REF][START_REF] Worda | The influence of Nos3 polymorphisms on age at menarche and natural menopause[END_REF] .

In an extensive haplotype-based association study of GnRH1 and GnRHR in 125 white men and women from North America, three intronic SNPs of GnRHR (rs3822196, hCV1960913 and rs3796718) were weakly associated with late pubertal development (p ! 0.05). In 506 American women of European, Asian or African descent, another GnRHR SNP (hCV3145733) was marginally associated with late AAM and a rare GnRHR haplotype with early AAM. None of the studied GnRH1 variants showed any relationship with the phenotypes [START_REF] Sedlmeyer | Determination of sequence variation and haplotype structure for the gonadotropin-releasing hormone (GnRH) and GnRH receptor genes: investigation of role in pubertal timing[END_REF] .

SNPs in the estrogen ␣ receptor (A/C 351, T/C 397) or ␤ receptor (A/G 1730) were reported to be possibly associated with AAM in 145 adolescents of a closed rural community in Greece [START_REF] Stavrou | Combined estrogen receptor alpha and estrogen receptor beta genotypes influence the age of menarche[END_REF] , in which none of the studied genotypes, however, reached Hardy-Weinberg equilibrium.

Among 38 SNPs covering the full transcript length of CYP 17 and CYP19 genes, Guo et al. [START_REF] Guo | Polymorphisms of estrogenbiosynthesis genes CYP 17 and CYP19 may influence age at menarche: a genetic association study in Caucasian females[END_REF] found that five SNPs of the CYP19 gene were significantly associated with AAM in 1,048 white American females, while none of the coding CYP17 variants showed any association. However, a functional variant in the promoter of the CYP17 gene was associated with AAM in 317 Japanese women, and in 583 Caucasian, African and Asian women living in North America [START_REF] Lai | Narod SA: CYP gene polymorphisms and early menarche[END_REF][START_REF] Gorai | Estrogen-metabolizing gene polymorphisms, but not estrogen receptor-alpha gene polymorphisms, are associated with the onset of menarche in healthy postmenopausal Japanese women[END_REF] . The polymorphic variants of the CYP3A4 , CYP1B1 (exon 3 Leu to Val) and CYP1A2 (intron 1) genes showed no association either with AAM in 583 Caucasian, African or Asian women living in North America [START_REF] Lai | Narod SA: CYP gene polymorphisms and early menarche[END_REF] , while an association was found between the CYP3A4 * 1B functional promoter variant and early puberty in 137 African-American, Hispanic and Caucasian girls [START_REF] Kadlubar | The CYP3A4 * 1B variant is related to the onset of puberty, a known risk factor for the development of breast cancer[END_REF] . Rothenbuhler/Lotton/Fradin

The functional (TAAAA)n polymorphism located in the promoter of the sex hormone-binding globulin (SHBG) gene was found to have a genetic contribution to AAM. Late menarche was more frequently associated with long alleles ( 1 8 repeats) of this variant in 130 adolescents of a closed Greek community [START_REF] Xita | Association of SHBG gene polymorphism with menarche[END_REF] . The Glu298Asp missense variant and T-786C gene polymorphisms located in the promoter region of the endothelial nitric oxide synthase (NOS1) gene showed no association with AAM in 87 Caucasian women [START_REF] Worda | The influence of Nos3 polymorphisms on age at menarche and natural menopause[END_REF] . A recent study showed no association of specific LEP (3 CTTT repeat) and LEPR (Gln 1 Arg substitution, exon 6) polymorphisms with constitutional delay of growth and puberty [START_REF] Banerjee | Phenotypic variation in constitutional delay of growth and puberty: relationship to specific leptin and leptin receptor gene polymorphisms[END_REF] .

The current study examined the possibility that yet untested common variants of LEP , NPY1R and GPR54 genes could be major genetic predicting markers or determinants of pubertal timing in humans. The study was designed to detect only a relatively large effect size of the studied polymorphisms considered separately. The negativity of our findings does not exclude that genetic variation in the considered genes could influence AAM, because other SNP combinations at the LEP , NPY1R or GPR54 gene loci that we have not tested could be involved. In such case, however, the genomic variability would have to involve variants that are not in LD with the studied SNPs. Our hypothesis was tested with the robust Q-TDT approach which excludes potential stratification effects in the studied women cohort [START_REF] De Muinich Keizer | Trends in pubertal development in Europe[END_REF] .

Since the level of influence of a given allele on the individual variability of a continuous trait is a direct function of both frequency and magnitude of the effect of this variant, the current association analysis suggests that none of the studied variant is involved at a major level ( 1 10% of the trait variance) in the genetic determination of pubertal maturation in humans.

Table 1 .

 1 Age at menarche for the whole sample and for each genotypic group of the studied variants of LEP, NPY1R and GPR54

		Whole	LEP			NPY1R			GPR54
		sample	CC	CA	AA	TT	CT	CC	AA	AT	TT
	Genotype
	Number	245	67	145	33	95	117	32	150	81	8
	Frequency		0.27	0.59	0.14	0.39	0.48	0.13	0.63	0.34	0.03
	Age at menarche
	(years, mean8SD)	13.181.5	13.481.5	13.281.6	13.281.4	13.181.4	13.281.6	13.581.5	13.281.5	13.381.5	13.481.6
	Range	9418	9.5417	9418	10416	9416	9.5418	10418	9.5418	9417	12417