Implementation of Laplace Transformed MP2 for Periodic Systems With Numerical Atomic Orbitals

We present an implementation of the canonical and Laplace-transformed formulation of the second-order Møller–Plesset perturbation theory under periodic boundary conditions using numerical atomic orbitals. To validate our approach, we show that our results of the Laplace-transformed MP2 correlation correction for the total energy and the band gap are in excellent agreement with the results of the canonical MP2 formulation. We have calculated the binding energy curve for the stacked trans-polyacetylene at the Hartree–Fock + MP2 level as a preliminary application.


INTRODUCTION
The second-order Møller-Plesset perturbation theory (MP2) is a post-Hartree-Fock approach to take the electron correlation effect into account. Although it is very simple in form, it can capture around 90% of the correlation energy (Bartlett and Stanton, 2007); so the MP2 method is still of high interest in the quantum chemistry (Schütz et al., 1999;Kobayashi and Nakai, 2006;Bartlett and Stanton, 2007) and solid-state physics communities (Suhai, 1983(Suhai, , 1992Sun and Bartlett, 1996;Pisani et al., 2008;Marsman et al., 2009;Schäfer et al., 2018).
However, the O(N 5 ) calculation scaling of the original (canonical) MP2 method has limited the application of the MP2 method in large systems. A series of algorithms have been proposed to speed up the calculations, such as local MP2 method (Saebø and Pulay, 1993;Pisani et al., 2005Pisani et al., , 2008Maschio, 2011), Lapace-transformed MP2 method (Häser and Almlöf, 1992;Häser, 1993;Ayala and Scuseria, 1999;Ayala et al., 2001;Schäfer et al., 2018), or resolution of the identity (RI) MP2 method (Katouda and Nagase, 2010;Ren et al., 2012). The local MP2 method proposed by Pulay (1983) and Saebø and Pulay (1993) has been efficiently implemented (Schütz et al., 1999) in the MOLPRO code for molecules, then the periodic version of the local MP2 method has been implemented (Pisani et al., 2005(Pisani et al., , 2008Maschio, 2011) in the CRYSCOR code and in the CP2K code (Usvyat et al., 2018) for extended systems. Since the spatially localized orbitals or Wannier functions are adopted, the computational scaling of the local MP2 method is O(N). The Laplace-transformed MP2 method is originally proposed by Häser and Almlöf (1992) and Häser (1993), and have been implemented for both the molecule (Ayala and Scuseria, 1999) and extended systems (Ayala et al., 2001) in the GAUSSIAN suite of programs. The localized atomic orbitals have been employed and the computational scaling is also O(N). The Laplace-transformed MP2 method has been combined with the resolution of identity (RI) technique to further improve the computational efficiency (Izmaylov and Scuseria, 2008). Further rigorous integral screening scheme (Lambrecht et al., 2005) has been introduced on top of the Laplace-transformed MP2 to perform the calculations for a system comprising 1,000 atoms (Doser et al., 2008). Recently, the Laplace-transformed MP2 method has also been implemented (Schäfer et al., 2018) in VASP using stochastic orbitals.
So far, most of the implementations of the MP2 are adopting the Gaussian-type orbital (GTO) as the basis set. However, in the calculation of the periodic system, too diffused GTO with a long tail will increase the number of cells in the auxiliary supercell, and therefore the computational cost will increase. Compared with GTO, the numerical atomic orbital (NAO) is strictly localized, which could naturally leads to lower order scaling of computational time vs. system size. Here in this work, we have implemented the canonical MP2 and Laplace-transformed MP2 for the extended systems using NAO, and the results obtained by these two approaches are consistent. Furthermore, we have investigated the MP2 correlation correction to the band structure with both the canonical and Laplace-transformed formulation; our implementation has been validated by comparing the MP2 correlation correction of the total energy and the band gap to the literature values.
The remainder of this paper is organized as follows. The fundamental theoretical framework and the implementation details for the canonical and Laplace transformed MP2 are presented in section 2. The benchmark calculations are presented in section 3. In section 4, we summarize our main achievement and highlight the possible future research direction related to this work.

Numerical Atomic Orbitals
The numerical atomic orbital is defined by a product of a numerical radial function and a spherical harmonic χ Ilmn (r) = ϕ Iln (r)Y lm (r) . (1) By solving the one-dimension radial Schrödinger equation we can get the radial part of the numerical atomic orbital ϕ Iln (r), where V(r) denotes the electrostatic potential for orbital ϕ Iln (r), and V cut ensures a smooth decay of each radial function, which is strictly zero outside a confining radius r cut . In order to perform the Hartree-Fock and MP2 calculation, the electron repulsion integrals (ERIs) are needed: we use NAO2GTO scheme described to calculate them as shown in the following section.

The NAO2GTO Scheme to Calculate ERIs
In the NAO2GTO scheme, we fit the NAO with GTOs, then we calculate the ERIs analytically; in this way, the strict cutoff of the atomic orbitals is satisfied with NAO and the construction of Hartree-Fock exchange (HFX) matrix can scale linearly with the system sizes (Shang et al., 2011). Since the angular part of the NAOs is spherical harmonic, while the GTOs are Cartesian harmonic function, a transformation between the Cartesian and spherical harmonic functions is performed within the NAO2GTO scheme.

Canonical MP2 Formulation
In extended systems, the normalized crystal orbital ψ i (k, r) is a linear combination of Bloch functions φ µ (k, r): in which N is the number of cells in extended systems, µ is the index of the atomic orbitals, i refers to the crystal orbital index, R denotes the cells in the extended systems (auxiliary supercell), refers to the atomic orbital whose center is displaced from the cell R by r µ , and C µ,i (k) are the coefficients of the crystal orbitals.
The MP2 correlation correction for the total energy of the unit cell is in which we use labeling i,j for occupied orbitals and a,b for unoccupied orbitals. I/J refer to the composite index (i,k i )/(j,k j ), V k is the volume of the Brillouin zone, and ǫ i (k i ) is the Hartree-Fock eigenvalue for the eigenstate ψ i (k i ). It should be noted that by using the identity ( R exp ik · R = Nδ k,0 ) derived with the Born-von Karman periodic boundary condition, we can remove one dimension integration over where T is the translation operator. The formalism of summation over 3-fold k points and over 4-fold k points (Equation 6) give the same results. Similarly, the MP2 correlation correction (ǫ g (k g ) (2) ) to the Hartree-Fock eigenstate ψ g (k g ) can be written as Frontiers in Chemistry | www.frontiersin.org When ψ g (k g ) is the occupied orbital at the valance band maximum (VBM), we can see U(g) < 0, V(g) > 0 and |U(g)| < |V(g)|, then the MP2 renormalization of the VBM is positive and will move the VBM orbital upward. When ψ g (k g ) is the unoccupied orbital at the conduction band minimum (CBM), we have |U(g)| > |V(g)|, so the MP2 renormaliztion of the CBM is negative, and will move the CBM downward. In total, the MP2 renormalization of the band gap is negative, and the MP2 band gap is smaller than the Hartree-Fock band gap.

Laplace-Transformed MP2 Formulation
The Laplace transform is defined as: which can be used to remove the denominator in the canonical MP2 formulation: The integration in Equation (10) can either be done by using a least square fitting method (Häser and Almlöf, 1992;Häser, 1993) or by using a Jacobian transform (Ayala and Scuseria, 1999;Kobayashi and Nakai, 2006) in which the Jacobian transform is Then the final integration ( 1 0 f (r)dr) in Equation (12) in evaluated with Romberg quadrature method, which uses refinements of the extended trapezoidal rule to reduce error in definite integrals.
In this way, the E mp2 correlation correction energy can be written as a new integration form: where µ,ν,σ ,λ and the following λ,δ,τ ,κ refer to the indexes of the atomic orbitals. (χ 0 µ χ R σ σ |χ R ν ν χ R λ λ ) is the electron repulsion integrals defined as and The 4-fold k points are treated independently within X In this way, the 4-fold integration over k-points in Equation (6) can be reduced to 1-dimensional k-points integral as shown in Equations (17) and (18). Furthermore, the locality of the atomic basis function can be adopted for the calculation of T R µ R λ R ν R σ µ,λ,ν,σ and electron repulsion integrals, and the total computational scaling could be O(N · N k ) if the distant screening between these ERIs are applied. Here in this work, such distance screening has not been used, so our implementation results in a O(N 2 · N k ) scaling. It is worth noting that in order to keep the exponential value (e (ǫ i )t ) in Equation (17)/Equation (18) to be smaller than unity, we have inserted the Fermi energy level into the exponential factor (e (ǫ i −ǫ f )t ) in order to make the calculation to be more stable.
Based on Equations (7) and (10), we have the Laplacetransformed MP2 correlation correction (ǫ g (k g ) (2) ) for the eigenstate: Similarly, in order to keep the exponential value to be smaller than unity and avoid computational divergence, we inserted the VBM/CBM value into the exponential factor when calculated the MP2 reformulation of the VBM/CBM. The canonical and Laplace-transformed MP2 methods described above have been implemented in the Order-N performance HONPAS code (Qin et al., 2014).

RESULTS
In order to validate our implementation, we perform benchmark calculations for 1-dimensional systems. We use norm-conserving pseudopotentials generated with the Troullier-Martins scheme to represent the interaction between core ion and valence electrons. The single-zeta (SZ), double-zeta (DZ), and doublezeta polarized (DZP) basis sets are generated using SIESTA. Then the NAOs are fitted with GTOs to perform the Hartree-Fock calculation as discussed in Shang et al. (2011).
First, we use a 1-dimensional hydrogen chain as an example to make the comparison between the results of canonical MP2 and those of Laplace-transformed MP2. The lattice parameter for the 1-dimensional H chain is set to 2.6 Å, and the H-H bond length is set to 1.346 Å. The SZ basis set is adopted, so that there are only two atomic orbitals in the unit cell. The Brillouin zone is sampled by 1 × 1 × 6 k-points. The unit cell is a 20 × 20 × 2.6 Å box, and the real-space integration mesh is set to be 100 Ry. In the Laplace-transformed MP2 method, the Romberg method is adopted to perform the final integration. The accuracy of integration results depends on the order of Romberg integration (n), since the results of the Romberg integration are obtained in a recursive manner, R(n, j) = R(n, j − 1) + R(n,j−1)−R(n−1,j−1) 4 j−1 −1 . As shown in Table 1, when the order of Romberg integration increases from 3 to 8, the MP2 correlation correction for the unit cell energy (E mp2 ) as well as the MP2 correlation correction for the band structure (ǫ (2) VBM , ǫ (2) CBM , ǫ (2) gap ) are converged to the results of canonical MP2. When using the media precision parameter (n = 5) in Laplace-transformed MP2, we get a absolute/relative error of 4 × 10 −6 eV/0.0007% for the correlation of the unit cell energy (E mp2 ), and we get a absolute/relative error of 3 × 10 −5 eV/0.02% for the MP2 correlation correction for the band gap (ǫ (2) gap ). Overall, we find an excellent agreement between the Laplace-transformed MP2 and the canonical MP2 benchmark results.
We also examine the relative error between the results of Laplace-transformed MP2 and those of canonical MP2 with respect to the basis set size (SZ, DZ, DZP), as shown in Table 2 with ethylene molecule as an example. Again, we find an excellent agreement between the Laplace-transformed MP2 and the canonical MP2 benchmark results.  Here, the MP2 correlation correction for the unit cell energy (E mp2 ) and the MP2 correlation correction for the band structure (ǫ gap ) have been examined with the above two approaches. Here, we use a 1-dimensional hydrogen chain as an example.
Second, we perform the Laplace-transformed MP2 calculation for the 1D polymer trans-polyacetylene as shown in Table 3. The order of Romberg integration is set to be n = 5. The SZ basis set is adopted in our calculation. The Brillouin zone is sampled by 1 × 1 × 30 k-points. The real-space integration mesh is set to be 200 Ry. We compare our calculated MP2 correlation correction for the total energy per unit cell with the one obtained in Sun and Bartlett (1996). The G3 geometry parameters as listed in Sun and Bartlett (1996) are adopted to keep the geometry to be the same for comparison. We get an absolute/relative error of 0.26 eV/8% for the correlation correction of the unit cell energy (E mp2 ). The difference comes from the usage of the different basis set, since in Sun and Bartlett (1996), the STO-3G basis set is used, whose shape is different from the SZ basis that we are using. For a similar reason, when using the same G6 geometry parameter of transpolyacetylene (Sun and Bartlett, 1996), we get an absolute/relative error of 0.09 eV/7% for the correlation correction of the band gap (ǫ (2) gap ) when compared with Sun's result. We then investigate the performance and scaling of our implementation, and we show timings for the transpolyacetylene molecules with variable number of atoms in Figure 1. We find a linear scaling for the calculation of ERIs and an O(N 2 ) scaling for the calculation of the Laplace-transformed MP2. This is not too surprising, since we can see from Equation (14) that there are two loops over the ERIs for the calculation of Laplace-transformed MP2, so we get the O(N 2 ) scaling.
Finally, we show the calculated binding-energy curves as functions of the distance between two trans-polyacetylene chains with the PBC-MP2 method. Although the MP2 theory gives overestimation of the dispersion interaction energy (Tkatchenko et al., 2009), it is still a superior starting point for the dispersion correction compared to Hartree-Fock and semi-local density functional theory (DFT). As shown in Figure 2, the MP2 method results in a binding states. On the contrary, Hartree-Fock and PBE functional fail to identify any binding between the two chains. We can see in Figure 2 that the energy profile   (Sun and Bartlett, 1996).
Trans-polyacetylene Sun and Bartlett (1996)   calculated with Hartree-Fock and PBE functional shows a repulsive behavior as the two chains are brought closer together.

CONCLUSIONS
We have implemented the canonical and Laplace-transformed algorithms to calculate the MP2 correlation correction for the total energy and the band gap of periodic systems in HONPAS code with numerical atomic orbitals. The results obtained by the canonical MP2 and Laplace-transformed MP2 are consistent with each other. We have also validated the implementation by comparing the results with the literature data. We have studied the binding-energy curves for the two stacked transpolyacetylene chains, which shows the MP2 method can well describe the correlation energy and the long-range van der Waals interactions. Future work will address the application of the Laplace-transformed MP2 method to 3-dimensional periodic systems in the HONPAS code.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
JY designed and directed the project. HS developed the theoretical formalism, carried out the implementation, and performed the numerical simulations. Both JY and HS contributed to the final version of the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
HS acknowledges Georg Kresse for inspiring discussions. HS thanks the Tianhe-2 Supercomputer Center for computational resources.