Coupled Cluster and Quantum Monte-Carlo potential energy curves of the ground state of Be2 and Be2+ molecules

Accepted Manuscript Coupled Cluster and Quantum Monte-Carlo potential energy curves of the ground state of Be2 and Be2 + molecules Saeed Nasiri, Mansour Zahedi PII: DOI: Reference:

S2210-271X(17)30203-7 http://dx.doi.org/10.1016/j.comptc.2017.04.009 COMPTC 2482

To appear in:

Computational & Theoretical Chemistry

Received Date: Revised Date: Accepted Date:

6 February 2017 11 April 2017 19 April 2017

Please cite this article as: S. Nasiri, M. Zahedi, Coupled Cluster and Quantum Monte-Carlo potential energy curves of the ground state of Be2 and Be2 + molecules, Computational & Theoretical Chemistry (2017), doi: http:// dx.doi.org/10.1016/j.comptc.2017.04.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Coupled Cluster and Quantum Monte-Carlo potential energy curves of the ground state of Be2 and Be2+ molecules

Saeed Nasiri and Mansour Zahedi* Department of Chemistry, Faculty of Chemistry, Shahid Beheshti University G.C., P.O. Box 19839-63113, Evin, Tehran, Iran 19839

*

Corresponding Author

E-mail: [email protected] or [email protected]. Tel: 98-21-29902889; Fax: 98-21-22431661;

1

ABSTRACT Accurate potential energy curves (PECs) have been computed by the Diffusion Quantum MonteCarlo (DMC) method for the ground state of the Be2 and Be2+ molecules. Multi-determinant trial wave functions (TWF) have been used in the DMC computations. The TWFs have been computed by the Complete Active Space Self-consistent Field (CASSCF) method with the augcc-pVTZ basis set. For atomic species (Be and Be+) the full active space and for molecular species (Be2 and Be2+) the valence full active space have been used. Despite using proper active spaces, the calculated binding energy is not fully compatible with the available experimental value. The source of this discrepancy has been investigated for the present and previous works. Furthermore, the PECs are recomputed by the CCSD(T) level of theory (with all electrons correlated) with the aug-cc-pCVQZ, cc-pV5Z, cc-pCV5Z, cc-pV6Z and 6ZaPa-NR basis sets. In addition, the correlation energy, the dissociation energy ( ), the ionization potential (IP) and spectroscopic constants have also been computed by the obtained PECs. The Hartree-Fock (HF) energy has been computed by the aforementioned basis sets for all of the studied species. All of the obtained HF energies in this work are more negative than the previous computations, except for the Be atom. 1. Introduction Beryllium is one of the lightest atoms in the periodic table; however, this atom is a big challenge to the quantum chemistry. Many theoretical [1-8] and experimental investigations [915] have been done in the last decades to improve our knowledge about this atom. Moreover, despite this apparent simplicity of the electronic structure and a huge amount of investigations, there is no clear understanding even of the beryllium smallest cluster, namely the beryllium dimer. In addition, experimental and traditional theoretical studies are not fully compatible. The isolated beryllium atom has a closed-shell electronic structure (1s22s2) and all of the 2p orbitals are empty and have very close energy to the occupied 2s orbital [16]. Because of this near degeneracy and existing two electrons in the valence shell, the Be could form sp hybridized bonds (such as BeH2). Bondybey [17-19] showed that the Be2 is a stable molecule in the gas phase. However, it has very weak bond (De = 934.6 cm-1) with an unusual long bond length (2.45 Å) [12, 13, 20]. The De, similar to the other usual quantities in chemistry, has been computed from the total energy differences and computing the total energies for the whole bond length region is crucial to calculate an accurate data. The computed PEC by the restricted Hartree-Fock (RHF) method is purely repulsive and could not predict a stable molecule [21]. It seems that the formation of this molecule is the result of the very small part of the total energy which is called the correlation energy (E corr). The Ecorr is absent in the HF computations. The E corr is classified into two types, dynamical and nondynamical (or static), with very unclear boundary. In order to calculate both types of the 2

correlation, different approaches must be used. The static correlation is usually associated with the degeneracy of orbitals (e.g. 2s and 2p orbitals in the current study) and determinant expansion could calculate most part of this type of the correlation. The Complete Active Space Self-consistent Field (CASSCF) method (in the following the method is denoted as CAS(n,m) for the sake of simplicity which indicates a calculation incorporating an active space explicitly correlating n electrons in m orbitals) is a popular method in this case and has been used widely in studying the beryllium species [21-23]. However, the recently published paper [21] show that the accuracy of CAS computations are extremely dependent on the employed active space and basis set (cc-pV6Z was used in the latter paper); although larger active space could not give necessarily more accurate results [21]. In addition to the CAS method, different versions of CI method have also been used. For example, Helal et al. performed the valence-Full Configuration Interaction (vFCI) computation in the singlet and triplet states of the Be 2 molecule. Two FCI computations [5,24] performed on the Be atom with very large basis set (e.g. 13s12p11d11f10g9h8i7k6l5m5n4o3q2r in the latter reference) and the calculated total energies are very accurate. The CI and Multi-reference CI (MR-CI) computations have been performed by Koput [25] with correlation consistent basis sets and obtained accurate results. The dynamical correlation is related to the Coulomb interaction among the electrons. This type of the correlation has been calculated by a very large number of excited determinants of the reference wave function. The excited determinants have a very small portion of the dynamical correlation energy and computing 100% of this correlation demands extremely large number of excited determinants. Even for small molecules like Be2, practically it is impossible to perform such a huge computation with even moderate basis sets. In order to compute the Ecorr more accurately, in the last decades, a new type of wave function has been used which involves the electron-electron distance (r12 = |r1−r2|) explicitly. By the means of the Explicitly Correlated Functions (ECF), more than of 99.99% of the Ecorr have been calculated for small atoms and molecules such as Be, Li and LiH. For smaller systems with 2 or 3 electrons, the accuracy provided by the computations while employing ECF is more than other accurate methods such as the Diffusion Quantum Monte-Carlo (DMC) or the FCI. However, the ECF method may suffer from higher computational cost in a system with more electrons and performing such computations may be so difficult for systems with 5 or 6 electrons while the obtained accuracy may be similar or even less than that of the DMC results. To the extent of the author's knowledge, there are only a few computations for the Be [26-31] and Be+ [32-35] atoms with Gaussian ECF (ECG). In contrast to the ECF, the computational cost scales with N 3-4 (N is the number of electrons) for the DMC method and it is possible to perform a computation for larger systems such as the Be2 molecule. In this method, usually the Slater-Jastrow wave function (SJ) is employed as the trial wave function (TWF). The SJ wave function is a multiplication of a Slater-based wave function and the Jastrow Factor [36]. The Slater part is computed by both of the HF-based methods (CAS, CI, …) or the density functional theory (DFT). This part controls the accuracy of the DMC computation, because the nodal surface of the system described by this part if the Jastrow part is totally symmetric [37,38]. Moreover, the static correlation is imported 3

in the DMC method by the Slater part and the Jastrow part defines the dynamical correlation. The latter part involves electron-electron and electron-nucleus distances explicitly. For the mentioned reasons, a wise choice for the Slater part could improve the accuracy of the DMC computations significantly. For this purpose in this paper, the CAS method has been used to compute the Slater part of the TWF which is described in details in the next section. In addition to the aforementioned papers, there are many theoretical and experimental studies in the literature about the Be2 molecule which some of the most accurate computations are presented in the following. As it was mentioned, performing a FCI computation with a large basis set for beryllium dimer is difficult. For this purpose, the researchers have been trying to improve the method by using explicitly correlated basis sets in the calculations. Sims and Hagstrom [39,40] used the Hylleraas-Configuration-Interaction (Hy-CI) method with only (s, p, d, f) basis functions which multiplied by a power of inter-electronic distance rij; with this modification, accurate results have been computed with small basis sets. Also, a computation with the (r 12)-MR-CI method (explicitly correlated multi-reference configuration interaction) was performed and the obtained result had considerable accuracy [3]. In addition to the CI method, this modification has applied to the Coupled Cluster method; for example, Noga et al. [4] used the explicitly correlated coupled cluster method with Single and Double and triple excitations (CCSDT-R12) to compute the beryllium atom energy and could obtain about 99.9% of the correlation energy. Chan et al. [41] reported a very accurate computation with the density matrix renormalization group (DMRG), the initiator full configuration interaction quantum Monte Carlo (i-FCIQMC) and the explicitly-correlated CCSDT-R12 methods. In the aforementioned computations, moderate basis sets (cc-pCVnZ-F12, n= D, T and Q) have been employed and both of the calculated total energy and beryllium dimer binding energies are exceedingly accurate. Deible et al. [42] used various active spaces and basis sets (cc-pVQZ without g and f,g functions) to prepare the TWF for the DMC computations in calculating the total energies of the Be and Be 2 species and obtained very accurate total energies. In addition to the CASSCF method, other methods (HF, DFT, CI) have also been used to prepare the TWF in the latter paper. Despite calculating very accurate results for total energies in equilibrium bond length (Re) and dissociation limit, no PEC has been computed by the latter papers and calculated De’s are not as accurate as the previous methods. In addition to single particle orbitals, Bertini et al. [43,44] used an expansion of the explicitly correlated basis functions as the TWF for the Variational quantum Monte Carlo (VMC) method; despite using a small basis function (86 terms), the calculated energies are accurate (more than 97% of the Ecorr was calculated). In addition to the ab initio methods, the empirical potential energy functions have also been used to calculate the spectroscopic constants of the ground states and the obtained results have considerable accuracy.[6, 12-14, 20, 45] In contrast to the neutral beryllium dimer, the Be2+ has a relatively strong bond (De = 16 438 cm-1) [14]; however, it has not been studied as widely as the Be2. Recently, Li et al. [46] performed the CCSD(T) and the vFCI computations with the aug-cc-pVXZ and the aug-cc4

pCVXZ (X=3, 4, 5). The total energies at the equilibrium bond length for the neutral and cationic forms have been computed and the ionization potential has also been calculated. In addition, a series of computations, for calculating the PEC, with the MRDCI method (multi-reference singles and doubles configuration interaction) has been performed by various groups with different basis sets [32, 47, 48]. Also, the PEC has been computed by RCCSD(T)/cc-pV5Z level of theory by Merritt et. al. [13]. In addition, the CCSD(T) and vFCI computations with correlation consistent basis sets reported by Banerjee, which the total energies are extrapolated to complete basis set (CBS) limit [49]. In addition to theoretical studies, experimental measurements have also been reported by a few groups [13, 45]; the reported spectroscopic constants are used here to examine the computed results. In this work, the PECs of the beryllium dimer in the ground state in the neutral and cationic forms have been investigated. For this purpose, the DMC method has been employed; the ability of this method in calculating the E corr is the reason for this choice. To the extent of our present knowledge, there is no study with the DMC method of beryllium dimer PEC in the aforementioned forms. The obtained PECs have been used to calculate the dissociation energies and spectroscopic constants. In addition, the PECs of the studied system have been calculated by the CCSD(T) level of theory with the aug-cc-pCVQZ, 6ZaPa-NR, cc-pV5Z, cc-pCV5Z and ccpV6Z basis sets. The four former basis sets are obtained from the EMSL Basis Set Library [50, 51] and the two latter are obtained from Ref. [52]. In the next section, the details of the computations for calculating the PECs and the fitting process of the vibrational energy levels are discussed. In the Results and Discussion section, the obtained results for the aforementioned states are presented. 2. Computational details As it was mentioned in the introduction, a multi-determinant wave function is essential for describing the electronic structure of beryllium species and in this work the wave functions have been constructed using the CAS method with the aug-cc-pVTZ basis set. Usually, all of the computed configurations by the CAS method are not used in the TWF and the configurations with larger coefficients are employed to construct the TWF [41, 42]. For this purpose, it has been tried to choose an active space for the studied systems which makes possible to use more configurations in the TWF even with very small coefficients. However, because of near degeneracy of the 2s and 2p orbitals in beryllium species, it is crucial to use these orbitals in the active space. For atomic species (Be and Be+) all electrons and 1s, 2s, 2p orbitals are considered in the active spaces; while the CAS (4,5) and CAS(3,5) computations have been used for the Be and Be+ species respectively. For the aforementioned species, 1×10-7 threshold on coefficients on configurations has been applied in the CAS computations. However, another DMC computation with CAS(4,6) active space (3s orbital is added in the active space) has also been performed with the same threshold on coefficients. In the case of studied molecular species, full valence CAS computations have been employed to construct the TWF. Furthermore, valence electrons and 5

σg(2s), σu *(2s), σg(2pz), σu*(2pz), πu(2px, 2py), πg*(2px, 2py) orbitals have been used in the active

spaces of CAS(4,8) and CAS(3,8) computations (for the Be2 and Be2+ species respectively). In addition to the DMC computations, the singles and doubles coupled-cluster methods including perturbative triple excitations, CCSD(T), with all electrons correlated has been utilized for the ground state of the Be, Be+, Be2 and Be2+ species with the aug-cc-pCVQZ, 6ZaPa-NR, cc-pV5Z, cc-pCV5Z and cc-pV6Z basis sets. As it was mentioned in the introduction, Koput [25] has studied the ground state of Be2 molecule with CCSD(T) level of theory with correlation consistent basis sets which the employed basis sets are larger than the aforementioned basis set. Despite using larger basis set, the obtained results are not so accurate. In this work, the VMC [53] and the DMC[54] methods as implemented in the CASINO package [55] have been employed to calculate the total energies of the aforementioned systems. In the CASINO code, the importance sampling DMC method [56] is implemented to compute the energy of the many-body systems. In addition, the CAS computations have been used to compute the Slater part of the TWF; these computations have been performed by the Gaussian09 software package [57]. In the next step in order to construct a more accurate TWF, the Jastrow factor [36] has been added to the Slater part which has been calculated by the CAS computations. The form of the Jastrow factor that was used in the CASINO code was proposed by Drummond et al. [37]. The employed Jastrow factor in this work is completely the same as that used in previous work (see equation 3) [58]. After introducing the Jastrow factor, the resulting TWF is not an optimized wave function. The parameters in the Jastrow factor and the coefficients of the configurations have been optimized via minimization of the variance by the VMC method [53]. The best-optimized wave function was chosen as the TWF for the DMC computation; the DMC calculations have been continued until the statistical error reached 0.01 kcal ⁄mol (≈16 μha) with a small time step (τ = 0.001 a.u.). The calculated total energies by the DMC and the CCSD(T) method were employed to solve the radial Schrödinger equation by a computer program [59]. To evaluate spectroscopic constants, all of the computed vibrational energy levels as a function of vibration quantum number were fitted to the Dunham energy-level expression [60]. 3. Results and Discussion In this section the obtained results by the DMC and the CCSD(T) method are presented in two subsections. The first subsection includes the calculated total energies of atomic and molecular species by the DMC and the CCSD(T) methods. In addition, a comprehensive discussion about the obtained energies has also been presented. At the next part, the calculated PECs have been discussed in details. 3.1 Total energy

6

Table 1 shows the computed total energies by the DMC method for the Be and Be + species in this paper along with some of the accurate results of previous works. In addition to the calculated energies by the previous methods, total energies of all studied systems have been computed by the CCSD(T) level of theory with aug-cc-pCVQZ, cc-pV5Z, cc-pCV5Z, 6ZaPa-NR and cc-pV6Z basis sets. Moreover, the Ecorr has also been computed for all of the studied systems in this and previous works; this quantity may reveal the accuracy of the different methods more clearly. For this purpose, the calculated HF energies in complete basis set limit (HFlimit) by Koga et al. [2] are presented in Table 1 which a numerical method was used in the computations. In addition, the HF energy of the Be with large basis set (cc-pV7Z) [25] is also represented in the table. The former value is more negative than the latter value by 0.13 µha. So the calculated energy by the numerical method is used as HFlimit in the Ecorr calculations for the Be atom. For the Be+ cation, there is only the HF energy which is calculated by the numerical method [2]; for this purpose, in addition to the numerical method, the HF energy has also been computed by the cc-pCV5Z, ccpV6Z and 6ZPA-NR basis sets. In contrast to the neutral form, for the cationic form the calculated energy by the numerical method is not the most negative energy and the obtained value by cc-pV6Z basis set is at least 63 µha more negative than the former energy. Therefore, the computed value by the cc-pV6Z basis set is used in the calculations of the E corr (see Table 1). As it is clear from the table, the ECG method could calculate all of the Ecorr accurately (100%). However, the calculated values by the DMC method are accurate as the former method; although, the results of the latter method are dependent on the employed TWFs in the calculations. For example, 99.90 and 88.40 percent of the Ecorr have been calculated by the DMC method with different TWFs. CAS(4,16)/cc-pVQZ-g (without g function) [42] and CAS(2,4)/ccpVTZ [61] methods have been used in the aforementioned papers to compute the TWF respectively. By using more an accurate TWF, the computed energies by this method become more accurate. As it was mentioned previously, the CAS(4,5) and CAS(3,5)/aug-cc-pVTZ computations are employed to calculate the TWF for the Be and Be + species respectively. The calculated total energies and the Ecorr for both of the studied systems are not as accurate as the previous DMC computations; however, 99.41 and 99.96% of the E corr have been calculated for the Be and Be+ atoms respectively. The CAS(4,6) computation has also been used to obtain the TWF; with the new TWF 99.63% of the E corr has been computed for the Be atom. Besides the DMC method, the CCSD(T) method with different and large basis sets has been employed to calculate the total energies of the Be and Be+ moieties and are shown in the table. The latter method has calculated 56-99 percent of the Ecorr; despite using larger basis sets (cc-pV6Z or 6ZaPa-NR), the most negative values for the Be and Be+ is calculated by cc-pCV5Z basis sets. This improvement could be because of correlation effects of the 1s core electrons of the Be species which is described more accurately in the core-valence correlation consistent basis sets. Also, Koput calculated the total energies at the CCSD(T)/aug-cc-pCVnZ (n=Q–7) level of the theory [25] and obtained more negative values for total energies of the studied systems (see Table 1). Additionally, Sharma et al. [41] have performed the DMRG with different basis set which the calculated energies are as accurate as the aforementioned methods (see Table 1). 7

In contrast to the atomic species, the HFlimit energies have not been reported for the Be2 and Be2+ moieties. In this work the HF energies for the molecular species have been computed by the cc-pCV5Z, cc-pV6Z and 6ZPA-NR basis sets; the most negative energy has been used in the calculation of the E corr for each species (see Table2). For the molecular systems there are no computations with ECF method; also, there are no FCI computations with large basis sets. However, FCI computations with a moderate basis set are shown in Table 2. Additionally, Sharma et al. [41] have performed the i-FCIQMC and the DMRG computations for the Be2 molecule with moderate basis sets (cc-pCVnZ-F12 n=D, T and Q) and the binding energy has been calculated. The total energy of the Be2 for the latter method is only available in the paper and is represented in the table which is the most negative reported total energy. In addition to the FCI computations, few computations with the vFCI method [21,23,49] with relatively large basis sets have been shown in the table; however, the obtained total energies of the studied molecular species are not as accurate as the DMC or even the CCSD(T) methods. Various DMC computations are performed for the Be2 molecule; Deible et al. performed a series of the DMC computation with different active spaces and basis sets [42]. The CAS(4,8), CAS(4,16) and CI computations with the cc-pVQZ basis set have been used (QZ abbreviation was used in the paper) for constructing the TWF. In the paper, the basis set is used in two different versions, at the first version f and g functions and at the second version only the g function are omitted. The first version [5s4p3d] is smaller than the employed basis set in this work (aug-cc-pVTZ or [5s4p3d2f]); with the same active spaces, the calculated energy in the present work for the Be2 molecule is more negative than the previous paper. This comparison makes clear that the effect of the basis set size is as important as the active space size. In addition, it seems that for the Be2 molecule sufficing just to the employed basis set used for the Be atom is not appropriate. For smaller systems, namely the Be atom, only the size of the active space affects the calculated energy by the DMC method and omitting the g or f function could not change the accuracy of obtained energies so much. In the aforementioned paper, despite using smaller basis set, due to employing larger active space (CAS(4,8)) the obtained energy for the Be atom is more negative than the present result. However, the second version of the basis set (5s4p3d2f) is similar to the employed basis set here and this basis set has only been used with CAS(4,16) active space; due to the larger active space, the calculated energies of the Be and Be 2 are more negative than the energies in the present work. Despite obtaining more accurate energies with larger active spaces for both of the Be and Be2 by Deible, the obtained De is less accurate than the present work (see Table 3). In this work, the CAS(4,5) and CAS(4,8) computations have been used for the Be and Be2 species as the Slater part of the TWF respectively and the calculated energies by the DMC method for both of the species are less accurate than the DMC computation with CAS(4,16) for both species; however, the calculated De by these energies (997.08 cm-1) is more accurate than the latter paper (845 cm-1). In order to reveal the source of this discrepancy and the effect of the basis set and the active space in the calculated energies by the DMC method and the De, the absolute error (ΔE = Ecal – EExp.) of the computed energies for the Be2 and Be species in the present work along with the most accurate values are shown in Figure 1. In the last paper (Ref. 8

[42]), various methods (CAS(4,8), CAS(4,16) and CI) have been used; as it is clear from the figure, by using larger basis sets or active spaces the absolute error in the total energy of the Be atom has not been changed so much (about 30 cm-1). However, this change for the Be2 molecule is up to 300 cm-1; despite this improvement, the error in the total energy of the Be 2 molecule is larger than the two Be atoms (see Figure 1). Furthermore, the difference between the absolute error of the largest CAS computations (CAS(4,16)) in the two Be atoms and the Be2 molecule is about 100 cm-1 which this difference directly enters in the calculation of the De. Besides, larger active spaces or even CI computations could not lower the error of the total energy of the Be2 molecule as the two Be atoms in the DMC computations. In fact, the discrepancy arises from the different accuracy in the describing the nodal surface of the Be and Be2 species. For example, in the CAS(4,8) computation for the Be atom 1s, 2s, 2pz,x,y , 3s, 3pz,x/y orbitals are employed in the active space; however, in the active space of the Be2, the molecular orbitals (MOs) derived predominantly only from the 2s and 2p atomic orbitals (AOs) are imported. It seems that much larger active spaces (larger than CAS(4,16)) are required for the Be2 molecule to obtain the same accuracy for both of the two Be atoms and the Be2 molecule. In this work, it has been tried to calculate the total energies both of the Be and Be2 with similar accuracy; for this purpose, in the active spaces of both species similar orbitals are used (for the Be atom 1s, 2s, 2p and for the Be2 molecule σg(2s), σu *(2s), σg(2pz), σu*(2pz), πu(2px, 2py), πg*(2px, 2py)). The absolute errors in the total energy of the two Be (2Be) atoms and the Be2 molecule are 242.7 and 180.1 cm-1 respectively which are larger than the other DMC methods; however, the difference between these errors is 62.6 cm-1 which is the absolute error in the computed De in this work. In contrast to the previous computations, in this work, the absolute error in the calculation of the energy of the two Be atoms are larger than the Be 2 molecule. In order to improve the computations, in the next step 3s orbital is also added to the active space (CAS(4,6)) and the DMC computation is performed again for the Be atom. This modification makes the total energy 87.79 cm-1 more negative and the new value for the De is 904.018 cm-1 which is only smaller than the experimental measurement by 30.5 cm-1 value. Table 3 shows the absolute error of the De for all of the present and previous methods. As it is evident from the Table 3 and Figure 1, the most accurate De is computed by the DMRG method; this method has computed the energy of the Be atom with the same error as the other methods. Although, the error in the energy of the Be2 molecule is much smaller than the other computations; however, the error in the calculated energies is not zero. Despite these errors in the calculated total energies of two Be atoms and the Be2 molecule, the De has been obtained with high accuracy. This accuracy in the De is due to the same accuracy in calculating of the total energies of two Be atoms and the Be2 molecule (65.4 and 68.9 cm-1 respectively). In addition to the aforementioned methods, the energies of the studied systems have been recomputed by CCSD(T) with the augcc-pCVQZ, cc-pV5Z, cc-pCV5Z, 6ZaPa-NR and cc-pV6Z basis sets. The calculated energies are not as accurate as the previous results; however, up to 98.47% of the Ecorr has been computed which is obtained by the cc-pCV5Z basis set.

9

In contrast to the neutral form, the cationic form of the beryllium dimer has not been paid much attention to as the former. For this species there is no DMC or other accurate computations and the most accurate computation has been performed by the CCSD(T) level of theory with the aug-cc-pCV5Z basis set [46]. Similar to the previous species, the total energy of the Be 2+ molecule is also computed by the DMC and the CCSD(T) methods and are shown in Table 2. Also, the calculated energies by Koput [25] with the CCSD(T) method but larger basis sets are also shown in the table. It is clear that larger basis sets give more negative energies. Similar to the Be2, despite using larger basis sets (cc-pV6Z or 6ZaPa-NR) the most negative value for the Be2+ species is calculated by the cc-pCV5Z basis sets. In the case of the Be2+, the calculated total energy by the DMC method is more negative than the other energies. In the literature there is no estimation of the exact total energy for the Be2+ molecule; so the CE% has not been computed. Despite importance of the E corr in computational studies, this quantity is not a measurable property and the accuracy of the computations becomes clear in comparison with experimental data. In this work, total energies of neutral and cationic forms of the beryllium atom and dimer are calculated; these energies have been used to compute the ionization potential (IP). Table 4 shows theoretical and experimental ionization potentials (IP) of the Be and Be2 in this work and some of the previous works. For the Be atom the most accurate IP belongs to the ECG wave function; however, the difference between the calculated value in this work by the DMC method and the experimental data is only 8.6 meV. In addition to the DMC results, IPs of the studied systems are calculated by CCSD(T) with different basis sets. In general, calculated IPs by the latter method is more accurate than the former. This accuracy in calculating the IPs by the CCSD(T) method may result from the similar ability of the method in computing the total energy of the species in neutral and cationic forms. 3.2 The PECs In the last step, the PECs of the Be2 and Be2+ species have been calculated with the DMC and the CCSD(T) methods. To the extent of the authors’ knowledge, no PEC has been reported with the DMC method for the aforementioned species. The calculated curves by the aforementioned methods are shown in Figure 2. The computed PECs in this work are shown in relative energies (in mha); the relative energies have been computed by the following equation: (1) In addition, Gdanitz [3] computed the PEC for the ground state of the Be2 molecule with r12MRCI method and obtained an accurate curve (De = 903 cm-1). Also, Koput [25] used the CCSD(T) and the CCSD(T)-F12 methods with very large basis sets (aug-cc-pCVnZ and n=Q, 5, 6, 7 and cc-pCVQZ-F12); however the obtained results, namely the De, are not accurate enough and even extrapolating to the complete basis set limit could not improve the accuracy. For example, 623, 671, 688 and 696 cm-1 values have been computed by the aug-cc-pCVnZ n = Q, 5, 6, 7 basis sets for the De. Although, more accurate total energy for both of the Be and Be 2 species 10

have been computed by these basis sets. The last three aforementioned basis sets have been optimized only by minimizing the total energy of the Be atom; it seems that the optimization process is not enough for obtaining an efficient basis set to calculate the PEC. In this work, 657.85 and 827.94 cm-1 have been obtained by the cc-pCV5Z and cc-pV6Z basis sets meaning that more accurate data have been calculated with smaller basis set. The employed basis set in this work have been obtained from the Ref. [52]. Despite smaller size of the cc-pV6Z basis set, it could predict the De more accurately. Also, Koput [25] has reported a curve which is calculated by the CCSD(T) method and corrected by the vFCI method and has higher accuracy with respect to the available experimental values. The PECs of the aforementioned papers are also shown in Figure 2. It is clear from the figure that the computed PECs by the DMC, CCSD(T)+FCI and r12-MRCI methods are very close to one another; however, the obtained curve by the CCSD(T)+FCI computation is more accurate than the others. As it was mentioned in the previous section, the CAS(4,8) computation could not describe the Be 2 molecule in the Re accurately. It seems that such a deficiency in describing the system is clearer in the distance between 3 and 5 Å with respect to the curve of the last two aforementioned methods and a larger active space is required for more accurate description. Among the computed PECs by the CCSD(T) method, the cc-pV5Z basis set could not predict the PEC even qualitatively; although the cc-pV6Z basis set predicts the PEC more accurately. The two employed core-valence correlation consistent basis sets (aug-cc-pCVQZ and cc-pCV5Z) give very close curves; while the larger basis set gives more negative energies as expected. Also, the 6ZPa-Na basis set is as large as the former basis set; however, the calculated PEC is not compatible with the other computed PECs. In summary, the computed PECs by the CCSD(T) method show that larger basis set is crucial to calculate an accurate PEC; in addition to the size of the basis sets, basis functions with higher angular momentum is also essential. Also, Table 5 shows other spectroscopic constants which have been computed in the present and previous works. For the Be2+ cation, the latest experimental values for the dissociation energy is 16 438 cm-1 [45]; however smaller value has been measured two years earlier by Merritt et al. (16 072 cm-1) [13]. In this part, the former value has been chosen and the calculated values are compared to this value. In this work, by the means of the DMC method 15 652.509cm-1 value is obtained. Moreover, the De has also been computed by the CCSD(T) method and shown in Table 6 along with the other spectroscopic data in the present and previous works. The best prediction belongs to the cc-pV6Z basis set (16 412.412 cm-1) which the obtained value is more accurate than the DMC result. Furthermore, the predicted value for the Re by the DMC and the CCSD(T) methods is very close to the experimental value (2.211 Å [45]). Also, Figure 3 shows the calculated curves by the aforementioned methods for the Be2+ cation. In comparison with the other computational methods, the obtained spectroscopic constants in this work are more accurate with respect to the reported experimental values. 4.

Conclusions

11

Potential energy curves have been computed by means of the DMC and the CCSD(T) methods for the ground state of the Be2 and Be2+ molecules. The trial wave functions for performing the DMC computations are used as a multiplication of multi-determinants wave function and the Jastrow factor. Also, the CCSD(T) method has been performed with the aug-cc-pCVQZ, 6ZaPaNR, cc-pV5Z, cc-pCV5Z and cc-pV6Z basis sets. The calculated dissociation energy in the present and previous works has not been fully compatible with experimental values. For this purpose, the origin of this problem has been investigated and revealed that this discrepancy arises from the different accuracy of the calculated energies of the Be and Be2 species. In this work, by obtaining the total energies by the same accuracy, the dissociation energy has been computed with only 30 cm-1 difference relative to the experimental value. In addition, the ionization potential and correlation energy are calculated; the obtained values well corroborate with the available experimental or previous theoretical values. By the aforementioned basis sets, the HF energies for all of the studied systems have been computed; except for the Be atom, all of the computed values are more negative than the previous reported values. The most negative values for the HF energy have been used in computing the correlation energy. In contrast to the neutral form, the calculated quantities, namely De or IP, of the cationic form (Be2+) are more accurate with respect to the experimental data. Furthermore, the vibrational energy levels and spectroscopic constants have been computed by the calculated PECs by the DMC and the CCSD(T) method. Acknowledgments We are grateful to Prof. S. W. Ng for providing us with the Gaussian software suite of programs and hardware (machine time) facilities; and, for providing us with the opportunity to access some new features of the Gaussian products. Our thanks also go to Dr K. Doblhoff-Dier for her invaluable technical help. The authors further acknowledge financial support from the Research Council of Shahid Beheshti University. Technical support of the Chemistry Computational Center at Shahid Beheshti University is also gratefully acknowledged.

12

Figure 1. Absolute error in the energies of the two Be atoms (2×Be), Be2 molecule and the dissociation energy (De) with respect to the experimental values. In this work, DMC(4,8) has been used for the Be2 molecule and for the Be atom CAS(4,5) and CAS(4,6) have been used. Also DMRG[41], DMC/CAS(4,8)/QZ-fg, DMC/CAS(4,16)/QZ-fg, DMC/CAS(4,16)/QZ-g, DMC/IC/QZ-g[42] and DMC/CAS(4,8)[62] have been shown.

13

Figure 2. Potential energy curves of the ground state of the Be2 molecule. The CCSD(T)+FCI and r12-MRCI curves are obtained from Refs. [25] and [3] respectively.

14

Figure 3. Potential energy curves of the ground state of the Be2+ molecule.

15

Table 1. Computed total (E), correlation (Ecorr) energies and percentage of the obtained correlation energy (CE%) of the Be and Be+ species. Be+

Be Ecorr (mha)

Eh HF/cc-pCV5Z HF/cc-pV6Z HF/6ZPA HF HF/cc-pV7Z CI CI-R12 r12-CCSDT h DMRG i LRDMC j DMC/CAS(2,4) k l DMC/CAS(4,8)/QZ-fg DMC/CAS(4,8)/QZ-fgl DMC/CAS(4,8)/QZ-gl DMC/CI/QZ-gl ECG

HYL VMC DMC/CAS(4,8) y CCSD(T)/ . aug-cc-pCV5Z z aug-cc-pCV5Z aa aug-cc-pCV6Z aa aug-cc-pCV7Z aa aug-cc-pCVQZ cc-pV5Z cc-pCV5Z cc-pV6Z 6ZPA

-14.573 012 1 -14.573 020 8 -14.573 020 12 b -14.573 023 17 c -14.573 023 04 -14.619 121 88 d -14.667 3560 e -14.667 3523 f -14.667 356 411g -14.667 261 -14.667207 -14.657 5(1) -14.656 45 -14.667 23 -14.667 27 -14.667 25 -14.66727 -14.667 356 486m -14.667 356 463o -14.657 9 q -14.665 1(2) -14.667 27

x

-14.666 511 j -14.666 511 -14.666 883 -14.667 050 -14.664 951 17 -14.645 729 01 -14.666 506 3 -14.654 648 8 -14.649 803 40 -14.666 803 ab ±7.05×10-6 DMC -14.667 015 ac ±8.17×10-6 Exp. -14.667 36 ad a used as the HF energy in the Ecorr computation. e f Ref. [24]. Ref.[5]. j k Ref. [63]. Ref. [61]. o p Ref. [65]. Ref. [66]. y z Ref. [62]. Ref.. [46] ac The CAS(4,6) has been used for the Be atom.

Ecorr CE(%)

Ecorr (mha)

Ecorr CE(%)

-14.324 763 176 4 n -14.324 763 175 7 o -14.324 763 18 p -14.324 763 176 790 150 r

47.305 47.305 47.305

100.0 100.0 100.0

47.305

100.0

-14.324 034 j

46.576

98.46

-14.323 143 07 -14.304 153 58 -14.324 032 2 -14.312 341 9 -14.308 095 97

45.685 26.695 46.574 34.884 30.638

96.58 56.43 98.45 73.74 64.77

-14.324 742 ±4.51×10-6

47.284

99.96

Eh -14.277 439 1 a -14.277 458 3 -14.277 395 49 -14.277 394 81

46.099 94.333 94.329 94.333 94.238 94.184 84.477 83.427 94.207 94.247 94.227

48.87 100.0 99.99 100.0 99.90 99.84 89.55 88.44 99.86 99.90 99.88

94.333 94.333

100.0 100.0

84.877

89.97

92.077 94.247

97.60 99.90

93.488 93.488 93.860 94.027 91.928 72.706 93.483 81.626 76.780

99.10 99.10 99.49 99.67 97.45 77.07 99.10 86.53 81.39

93.780

99.41

93.992

99.63

94.337 100.0 b Ref.[2]. g Ref. [40]. l Ref. [42]. q Ref. [7]. aa Ref. [25]. ad Ref. [8].

16

c

d Ref. [25]. Ref. [23]. i Ref. [4]. Ref. [41]. m n Ref. [64]. Ref. [31]. r x Ref. [33]. Ref. [44]. ab The CAS(4,5) has been used for the Be atom. h

Table 2 . Computed total (E), correlation (Ecorr) energies and percentage of the obtained correlation energy (CE%) of the Be2 and Be2+ molecules. Be2+

Be2 Ecorr (mha)

Eh HF/cc-pV5Z HF/cc-pCV5Z HF/cc-pV6Z a HF/6ZPA CI DMRG d e LRDMC f DMC/CAS(4,8)/QZ-fg DMC/CAS(4,8)/QZ-fg f DMC/CAS(4,8)/QZ-g f f DMC/CI/QZ-g g DMC/CAS(4,8) c CCSD(T)/aug-cc-pCV5Z CCSD(T)/aug-cc-pCVQZ h CCSD(T)/aug-cc-pCV5Z h CCSD(T)/aug-cc-pCV6Z h CCSD(T)/aug-cc-pCV7Z h CCSD(T)/aug-cc-pCVQZ CCSD(T)/cc-pV5Z CCSD(T)/cc-pCV5Z CCSD(T)/cc-pV6Z CCSD(T)/6ZPA

-29.146 020 40 -29.146 024 2 -29.146 041 4 -29.146 040 2 -29.241 305b -29.338 657 -29.322 95 -29.337 07(3) -29.338 32 -29.338 38(3) -29.338 48(2) -29.337 36 -29.336 080 -29. 334 196 -29.336 080 -29.336 901 -29.337 272 -29.328 423 -29.291 245 -29.336 010 -29.313 070 -29.304 368 -29.338 149 DMC ±1.00×10-5 Exp. i -29.338 97 a used as the HF energy in the Ecorr computation. d e Ref. [41]. Ref. [63]. i Ref. [67].

Ecorr CE(%)

95.264 192.616 176.909 191.029 192.279 192.339 192.439 191.319 190.039 188.155 190.039 190.860 191.231 182.382 145.204 189.969 167.029 158.327

49.38 99.84 91.70 99.02 99.66 99.69 99.75 99.17 98.50 97.53 98.50 98.93 99.12 94.53 75.26 98.47 86.58 82.06

192.108

99.57

192.929 Ref. [1]. f Ref.[42] . b

17

Eh

Ecorr (mha)

-28.921 390 -28.921 400 -28.921 425 -28.921 365 -28.970 500 c

49.074

-29.064 807

143.381

-29.061 010 -29.024 450 -29.062 996 -29.040 516 -28.921 358 -29.064 897 ±1.10×10-5

139.584 103.024 141.570 119.090 143.472

100.00 c the aug-cc-pCV5Z is used Ref. [46]. g h Ref.[62] . Ref. [25].

Table 3. The De and absolute error in the De. All values are in cm-1 unit. method De De(cal.) - De(Exp.) Ref. FCI 844.800 -89.8 [23] MR-CI 922.900 -11.7 [20] DMRG 930.572 -4.0 [41] i-FCIQMC 923.988 -10.6 [41] CCSD(T)-F12b 689.150 -245.5 [41] CCSD(T)+FCI/CBS 938 3.4 [68] CCSD(T) + FCI 935.1 0.5 [25] DMC/CAS(4,8)-fg 573 -361.6 [42] DMC/CAS(4,16)-fg 819 -115.6 [42] DMC/CAS(4,16)-g 845 -89.6 [42] DMC/CI 873 -61.6 [42] DMC/CAS(4,8) 618.64 -315.88 [62] LRDMC/JAS-HF 1153.4 218.8 [63] CCSD(T)/aug-cc-pCV5Z 671 -263.6 [46] CCSD(T)/aug-cc-pCV5Z 671.153 -263.4 [25] CCSD(T)/aug-cc-pCV6Z 688.053 -246.5 [25] CCSD(T)/aug-cc-pCV7Z 696.174 -238.4 [25] CCSD(T)/ aug-cc-pCVQZ 630.924 -303.7 CCSD(T)/cc-pCV5Z 657.853 -276.7 CCSD(T)/cc-pV6Z 827.946 -106.7 CCSD(T)/6ZPA 1044.963 110.4 997.08 62.5 DMC a 904.02 -30.6 DMC b Exp. 934.6 0.0 [20] a CAS(4,5) and CAS(4,8) have been used as the TWF for the Be and Be2 species respectively. b CAS(4,6) and CAS(4,8) have been used as the TWF for the Be and Be2 species respectively.

18

Table 4. The Ionization potential (in eV) of the Be and Be2 species. method Be Be2 CI a 7.369 CCSD(T)/aug-cc-pCV5Z b 9.319 7.382 ECG 9.322 c 9.322 d e HYL 9.065 CCSD(T)/aug-cc-pCVQZ 9.301 7.276 CCSD(T)/cc-pV5Z 9.295 7.260 CCSD(T)/cc-pCV5Z 9.319 7.429 CCSD(T)/cc-pV6Z 9.314 7.416 CCSD(T)/6ZPA 9.298 10.422 9.308 f DMC 7.491 9.314 g h Exp. 9.322 63 7.418 i a

b

c

E(Be2) and E(Be2+) are obtained from Refs. [1] and [46] respectively. E(Be2) and E(Be2+) are obtained from Refs. [64] and [31] respectively. e E(Be2) and E(Be2+) are obtained from Refs. [7] and [33] respectively. g The CAS(4,6) has been used for the Be atom.

d

19

E(Be2) and E(Be2+) are obtained from Ref. [46]. E(Be2) and E(Be2+) are obtained from Ref. [65]. f The CAS(4,5) has been used for the Be atom. h i Ref. [69]. Ref. [13].

Table 5. Calculated spectroscopic constants for the ground state of the Be2 molecule. method

Re (Å)

FCI a CCSD(T)/aug-cc-pCVQZ b CCSD(T)/aug-cc-pCV5Z b CCSD(T)/aug-cc-pCV6Z b CCSD(T)/aug-cc-pCV7Z b

2.47 2.4732 2.4675 2.4645 2.4636

CCSD(T)/aug-cc-pCVQZ

2.45

CCSD(T)/cc-pCV5Z

2.45

CCSD(T)/cc-pV6Z

2.43

CCSD(T)/6ZPA

2.55

DMC

2.42 c

Exp. a

Ref. [23].

2.438 Ref. [25].

b

20

ωe (cm-1) 219.9 186.1 192.1 194.5 195.8 300 ±10.5 280 ±5.5 240 ±11.3 229.7 ±1.25 245 ±3.9 222.6

ωexe (cm-1) 38.6

Be (cm-1) 0.5877 0.5913 0.5930 0.5938

67 ±8.5 50.2 ±2.5 16.58 ±5.91 30.18 ±3.6 45.73 ±1.2 c

0.6975 0.7831 ±0.01 0.8633 ±0.055 0.6422 0.7074

0.609 Ref. [20].

Table 6. Calculated spectroscopic constants for the ground state of the Be2+.molecule. method FCI a CCSD(T)/cc-pV5Z MRDCI

Re (Å)

De (cm-1)

2.223 2.211 2.236 b 2.249 c

16 172 16 684 16 212 15 646 15652.509 ± 4.89

DMC

2.23

CCSD(T)/aug-cc-pCVQZ

2.22

16 272.156

CCSD(T)/cc-pCV5Z

2.23

16 184.704

CCSD(T)/cc-pV5Z

2.26

16 639.754

CCSD(T)/cc-pV6Z

2.22

16 412.412

CCSD(T)/6ZaPa-NR

2.30

16 817.090

ωe (cm-1)

ωexe (cm-1)

Be (cm-1)

513

4.2

528.0

4.8

518.7 ± 0.6 518.6 ± 0.30 522.3 ± 0.30 521.9 ± 0.30 466.2 ± 0.4

3.596 ± 0.061 4.166 ± 0.034 3.777 ± 0.038 4.072 ± 0.036 3.58 ± 0.04

0.8059 ± 0.003 0.7579 0.7485 0.756 5 0.754 2 0.6976

d

Exp. a

Ref. [49]

2.211 b Ref. [47].

16 072 16 438 e c Ref. [32].

21

524.88 4.44 d Ref. [13].

e

Ref. [14].

Reference. [1] S. Evangelisti, G.L. Bendazzoli, L. Gagliardi, Full configuration interaction calculations on Be2, Chem. Phys., 185 (1994) 47-56. [2] T. Koga, S. Watanabe, K. Kanayama, R. Yasuda, A.J. Thakkar, Improved Roothaan–Hartree–Fock wave functions for atoms and ions with N≤54, The Journal of Chemical Physics, 103 (1995) 3000-3005. [3] R.J. Gdanitz, Accurately solving the electronic Schrödinger equation of atoms and molecules using explicitly correlated (r12-)MR-CI.: The ground state of beryllium dimer (Be2), Chem. Phys. Lett., 312 (1999) 578-584. [4] J. Noga, D. Tunega, W. Klopper, W. Kutzelnigg, The performance of the explicitly correlated coupled cluster method. I. The four‐electron systems Be, Li−, and LiH, J. Chem. Phys., 103 (1995) 309-320. [5] O. Jitrik, C.F. Bunge, Atomic configuration interaction and studies of He, Li, Be, and Ne ground states, Phys. Rev. A, 56 (1997) 2614. [6] V. Špirko, Potential energy curve of Be2 in its ground electronic state, J. Mol. Spectrosc., 235 (2006) 268-270. [7] R.F. Gentner, E.A. Burke, Calculation of the S 1 State of the Beryllium Atom in Hylleraas Coordinates, Phys. Rev., 176 (1968) 63. [8] E. Lindroth, H. Persson, S. Salomonson, A.M. Mårtensson-Pendrill, Corrections to the beryllium ground-state energy, Phys. Rev. A, 45 (1992) 1493. [9] A.-M. Mårtensson-Pendrill, S.A. Alexander, L. Adamowicz, N. Oliphant, J. Olsen, P. Öster, H.M. Quiney, S. Salomonson, D. Sundholm, Beryllium atom reinvestigated: A comparison between theory and experiment, Phys. Rev. A, 43 (1991) 3355. [10] L.A. Kaledin, A.L. Kaledin, M.C. Heaven, V.E. Bondybey, Electronic structure of Be2: theoretical and experimental results1, J. Mol. Struct., 461–462 (1999) 177-186. [11] M.C. Heaven, J.M. Merritt, V.E. Bondybey, Bonding in beryllium clusters, Annu. Rev. Phys. Chem., 62 (2011) 375-393. [12] J.M. Merritt, V.E. Bondybey, M.C. Heaven, Beryllium dimer—caught in the act of bonding, Science, 324 (2009) 1548-1551. [13] J.M. Merritt, A.L. Kaledin, V.E. Bondybey, M.C. Heaven, The ionization energy of Be 2, and spectroscopic characterization of the (1) 3Σ+u,(2) 3Πg, and (3) 3Πg states, Phys. Chem. Chem. Phys., 10 (2008) 4006-4013. [14] I.O. Antonov, B.J. Barker, V.E. Bondybey, M.C. Heaven, Spectroscopic characterization of Be, J. Chem. Phys., 133 (2010) 1. [15] P.F. Bernath, Extracting potentials from spectra, Science, 324 (2009) 1526-1527. [16] S. Salomonson, H. Warston, I. Lindgren, Many-body calculations of the electron affinity for Ca and Sr, Phys. Rev. Lett., 76 (1996) 3092. [17] V.E. Bondybey, Electronic structure and bonding of Be 2, Chem. Phys. Lett., 109 (1984) 436-441. [18] V.E. Bondybey, Laser-induced fluorescence and bonding of metal dimers, Science, 227 (1985) 125131. [19] V.E. Bondybey, J.H. English, Laser vaporization of beryllium: Gas phase spectrum and molecular potential of Be2, J. Chem. Phys., 80 (1984) 568-570. [20] K. Patkowski, V. Špirko, K. Szalewicz, On the elusive twelfth vibrational state of beryllium dimer, Science, 326 (2009) 1382-1384. [21] M. El Khatib, G.L. Bendazzoli, S. Evangelisti, W. Helal, T. Leininger, L. Tenti, C. Angeli, Beryllium dimer: a bond based on non-dynamical correlation, J. Phys. Chem. A, 118 (2014) 6664-6673. [22] M.R.A. Blomberg, P.E.M. Siegbahn, Beryllium dimer, a critical test case of MBPT and CI methods, Int. J. Quantum Chem, 14 (1978) 583-592. [23] W. Helal, S. Evangelisti, T. Leininger, A. Monari, A FCI benchmark on beryllium dimer: The lowest singlet and triplet states, Chem. Phys. Lett., 568–569 (2013) 49-54. [24] C.F. Bunge, Configuration interaction benchmark for Be ground state, Theor. Chem. Acc., 126 (2010) 139-150. 22

[25] J. Koput, The ground-state potential energy function of a beryllium dimer determined using the single-reference coupled-cluster approach, Phys. Chem. Chem. Phys., 13 (2011) 20311-20317. [26] J. Komasa, Dipole and quadrupole polarizabilities and shielding factors of beryllium from exponentially correlated Gaussian functions, Phys. Rev. A, 65 (2001) 012506. [27] J. Komasa, Lower bounds to the dynamic dipole polarizability of beryllium, Chem. Phys. Lett., 363 (2002) 307-312. [28] J. Komasa, J. Rychlewski, K. Jankowski, Benchmark energy calculations on Be-like atoms, Phys. Rev. A, 65 (2002) 042507. [29] K. Pachucki, J. Komasa, Excitation energy of Be 9, Phys. Rev. A, 73 (2006) 052502. [30] M. Stanke, J. Komasa, S. Bubin, L. Adamowicz, Five lowest S 1 states of the Be atom calculated with a finite-nuclear-mass approach and with relativistic and QED corrections, Phys. Rev. A, 80 (2009) 022514. [31] M. Stanke, J. Komasa, D. Kędziera, S. Bubin, L. Adamowicz, Three lowest S states of 9Be+ calculated with including nuclear motion and relativistic and QED corrections, Phys. Rev. A, 77 (2008) 062509. [32] B. Meng, P.J. Bruna, J.S. Wright, Stability of highly excited electronic states of Be+2, Mol. Phys., 79 (1993) 1305-1325. [33] M. Puchalski, D. Kȩdziera, K. Pachucki, Ground state of Li and Be+ using explicitly correlated functions, Phys. Rev. A, 80 (2009) 032521. [34] M. Puchalski, K. Pachucki, Relativistic, QED, and finite nuclear mass corrections for low-lying states of Li and Be+, Phys. Rev. A, 78 (2008) 052511. [35] L.-Y. Tang, Z.-C. Yan, T.-Y. Shi, J. Mitroy, Dynamic dipole polarizabilities of the Li atom and the Be+ ion, Phys. Rev. A, 81 (2010) 042521. [36] R. Jastrow, Many-Body Problem with Strong Forces, Phys. Rev., 98 (1955) 1479-1484. [37] N. Drummond, M. Towler, R. Needs, Jastrow correlation factor for atoms, molecules, and solids, Phys. Rev. B, 70 (2004) 235119. [38] N.D. Drummond, R.J. Needs, Variance-minimization scheme for optimizing Jastrow factors, Phys. Rev. B, 72 (2005) 085124. [39] J.S. Sims, S. Hagstrom, Combined Configuration-Interaction—Hylleraas-Type Wave-Function Study of the Ground State of the Beryllium Atom, Phys. Rev. A, 4 (1971) 908. [40] J.S. Sims, S.A. Hagstrom, Hylleraas-configuration-interaction study of the 1 S ground state of neutral beryllium, Phys. Rev. A, 83 (2011) 032518. [41] S. Sharma, T. Yanai, G.H. Booth, C.J. Umrigar, G.K.-L. Chan, Spectroscopic accuracy directly from quantum chemistry: Application to ground and excited states of beryllium dimer, J. Chem. Phys., 140 (2014) 104112. [42] M.J. Deible, M. Kessler, K.E. Gasperich, K.D. Jordan, Quantum Monte Carlo calculation of the binding energy of the beryllium dimer, J. Chem. Phys., 143 (2015) 084116. [43] L. Bertini, D. Bressanini, M. Mella, G. Morosi, Linear expansions of correlated functions: Variational Monte Carlo case study, Int. J. Quantum Chem, 74 (1999) 23-33. [44] L. Bertini, M. Mella, D. Bressanini, G. Morosi, Explicitly correlated trial wavefunctions in quantum Monte Carlo calculations of excited states of Be and Be, Journal of Physics B: Atomic, Molecular and Optical Physics, 34 (2001) 257. [45] I.O. Antonov, B.J. Barker, V.E. Bondybey, M.C. Heaven, Spectroscopic characterization of Be2+X Σ2u+ and the ionization energy of Be2, J. Chem. Phys., 133 (2010) 074309. [46] H. Li, H. Feng, W. Sun, Y. Zhang, Q. Fan, K.A. Peterson, Y. Xie, H.F. Schaefer Iii, The alkaline earth dimer cations (Be2+, Mg2+, Ca2+, Sr2+, and Ba2+). Coupled cluster and full configuration interaction studies†, Mol. Phys., 111 (2013) 2292-2298. [47] H. Hogreve, Theoretical calculations for low-lying adiabatic states of Be2+, Chem. Phys. Lett., 187 (1991) 479-486. [48] I. Fischer, V.E. Bondybey, P. Rosmus, H.J. Werner, Theoretical study of the electronic states of BeLi and Be2+, Chem. Phys., 151 (1991) 295-308. 23

[49] S. Banerjee, J.N. Byrd, R. Côté, H.H. Michels, J.A. Montgomery, Ab initio potential curves for the and states of: Existence of a double minimum, Chem. Phys. Lett., 496 (2010) 208-211. [50] K.L. Schuchardt, B.T. Didier, T. Elsethagen, L. Sun, V. Gurumoorthi, J. Chase, J. Li, T.L. Windus, Basis set exchange: a community database for computational sciences, J. Chem. Inf. Model., 47 (2007) 1045-1052. [51] D. Feller, The role of databases in support of computational chemistry calculations, J. Comput. Chem., 17 (1996) 1571-1586. [52] B.P. Prascher, D.E. Woon, K.A. Peterson, T.H. Dunning Jr, A.K. Wilson, Gaussian basis sets for use in correlated molecular calculations. VII. Valence, core-valence, and scalar relativistic basis sets for Li, Be, Na, and Mg, Theor. Chem. Acc., 128 (2011) 69-82. [53] H. Conroy, Molecular Schrödinger Equation. I. One‐Electron Solutions, J. Chem. Phys., 41 (1964) 1327-1331. [54] D.M. Ceperley, B. Alder, Ground state of the electron gas by a stochastic method, Phys. Rev. Lett., 45 (1980) 566. [55] R. Needs, M. Towler, N. Drummond, P.L. Ríos, Continuum variational and diffusion quantum Monte Carlo calculations, J. Phys. Condens. Matter, 22 (2010) 023201. [56] C. Umrigar, M. Nightingale, K. Runge, A diffusion Monte Carlo algorithm with very small time‐step errors, J. Chem. Phys., 99 (1993) 2865-2890. [57] G.W.T. M.J. Frisch, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R., V.G.Z. Cheeseman, J.A. Montgomery, R.E. Stratmann, J.C., S.D. Burant, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain,, J.T. O. Farkas, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C., C.A. Pomelli, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q., K.M. Cui, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B., J.C. Foresman, J.V. Ortiz, B.B. Stefanov, G. Liu, A. Liashenko,, I.K. P. Piskorz, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A., C.Y.P. Al-Laham, A. Nanayakkara, C. Gonzales, M. Challacombe,, B.G.J. P.M.W. Gill, W. Chen, M.W. Wong, J.L. Andres, M.R.E.S., J.A.P. Head-Gordon, Gaussian 98, Gaussian, Inc., Pittsburgh, PA,, 1998. [58] S. Nasiri, M. Zahedi, Accurate potential energy curves of Li2 and LiH: A Quantum Monte-Carlo (QMC) study, Chem. Phys. Lett., 634 (2015) 101-107. [59] R.J. Le Roy, LEVEL 8.0: A computer program for solving the radial Schrödinger equation for bound and quasibound levels, University of Waterloo Chemical Physics Research Report CP-663, (2007). [60] J.L. Dunham, The energy levels of a rotating vibrator, Phys. Rev., 41 (1932) 721. [61] J.A.W. Harkless, K.K. Irikura, Multi‐determinant trial functions in the determination of the dissociation energy of the beryllium dimer: Quantum Monte Carlo study, Int. J. Quantum Chem, 106 (2006) 2373-2378. [62] J. Toulouse, C.J. Umrigar, Full optimization of Jastrow–Slater wave functions with application to the first-row atoms and homonuclear diatomic molecules, J. Chem. Phys., 128 (2008) 174101. [63] M. Marchi, S. Azadi, M. Casula, S. Sorella, Resonating valence bond wave function with molecular orbitals: Application to first-row molecules, J. Chem. Phys., 131 (2009) 154116. [64] K. Strasburger, Excited S-symmetry states of positronic lithium and beryllium, J. Chem. Phys., 144 (2016) 144316. [65] M. Stanke, D. Kȩdziera, S. Bubin, L. Adamowicz, Ionization potential of 9Be calculated including nuclear motion and relativistic corrections, Phys. Rev. A, 75 (2007) 052510. [66] S. Bubin, O.V. Prezhdo, Excited States of Positronic Lithium and Beryllium, Phys. Rev. Lett., 111 (2013) 193401. [67] R.J. Needs, M.D. Towler, N.D. Drummond, P.L. Ríos, Continuum variational and diffusion quantum Monte Carlo calculations, J. Phys. Condens. Matter, 22 (2009) 023201. [68] K. Patkowski, R. Podeszwa, K. Szalewicz, Interactions in diatomic dimers involving closed-shell metals, J. Phys. Chem. A, 111 (2007) 12822-12838. [69] R.D. Johnson Iii, NIST Computational Chemistry Comparison and Benchmark Database (CCCBDB), NIST Standard Reference Database Number 101, Release 15b, August 2011. 24

   

Potential energy curves have been computed by the DMC and the CCSD(T) methods. Multi-configuration expansions and the Jastrow factor used as trial wave function. Dissociation energy (De) and other spectroscopic constants have been calculated. The origin of discrepancy between calculated and experimental De has been proposed.

25

Graphical Abstract