Implementation of Picture Change Corrected Density Functional Theory Based on Infinite-Order Two-Component Method to GAMESS Program

Abstract

We implemented the picture change correction method for two-electron Coulomb interaction and density operator, which was based on the infinite-order two-component method with the local unitary transformation, to the GAMESS program. Numerical assessments for molecules containing heavy element confirmed the accuracies and efficiencies of the implementation. Furthermore, the comparison with several types of treatments indicated that whole picture change corrections of one- and two-electron operators and the density operator are necessary for accurate two-component density functional theory calculations.

1 Introduction

Relativistic effect is essential to elucidate molecular properties or chemical reactions in heavy-element compounds. Two-component (2c) relativistic methods are one of the effective schemes to describe relativistic effects in quantum chemical calculations. In our group, accurate and efficient relativistic methods based on the infinite-order two-component (IOTC) method [1] with the local unitary transformation (LUT) [2,3,4], have been developed. Especially, we implemented the one-electron (1e) LUT-IOTC method for energy [2] and analytical energy gradient [5] calculations into the public version of GAMESS program package [6, 7]. The GAMESS documentation explains that the LUT-IOTC method is the fastest and most numerically reliable scalar relativistic method in the program.

In the 2c relativistic methods, the electronic Hamiltonian is extracted by the unitary transformation of the four-component (4c) Dirac Hamiltonian. At the same time, the pictures of any operators become different from those in 4c, which is termed the picture change (PC). Previous studies examined the PC effects in two-electron (2e) interaction [3] and density operator for the density functional theory (DFT) [8]. This letter describes the implementations of these PC-corrected operators at the spin-free LUT-IOTC level to the GAMESS program.

2 Theory and implementation

This section briefly explains the theoretical aspects of PC to evaluate expectation values for arbitrary operators. In the IOTC method, 1e Dirac Hamiltonian h^D is completely block-diagonalized into electronic (h⁺) and positronic (h⁻) Hamiltonians by the unitary transformation U, which is written as    .    (1)   

Here, the relationship between 4c and 2c wave functions Ψ₄ and Ψ₂ is given by    .    (2)   

Using equation (2), the 2c expression of the expectation value for an arbitrary 4c operator X is written as    ,    (3)    where [U^†XU]₁₁ represents the (1,1) (electronic) component of the U^†XU matrix. Thus, the picture changes from X to U^†XU [9]. For example, DFT uses the density operator (X = δ) to obtain electron density. Furthermore, in order to obtain the results close to the 4c method, PC correction in 2e Coulomb interaction (X = 1/r_ij) is also required [10].

The unitary transformation for the entire system becomes a time-consuming step. The number of unitary transformations in DFT is the same as that of grid points. Here, the computational cost in each unitary transformation scales O(N³) with respect to the number of primitive basis functions N. In addition, the scaling for 2e Coulomb interaction is O(N⁵). To avoid the problem, we have developed the LUT scheme, which is based on the locality of the relativistic effect. The LUT scheme approximates U for the entire system as the direct summation of the contributions in subsystems like,    ,    (4)    where U^A, U^B, … are the subsystems, which are normally set to atoms. In addition, the approximated unitary transformation is applied to the restricted region near each nucleus. Thus, the computational scalings become O(N) for one- and two-electron operators and O(N²) for density operator.

Table 1 shows the major capabilities in GAMESS, which involves the present codes. First, we implemented LUT-IOTC for 2e Coulomb interaction, which is termed LUT-IOTC/IOTC. The method was connected with several wave function theories such as the second-order Møller–Plesset (MP2), configuration interaction (CI), and the coupled cluster (CC) as well as DFT. Geometry optimization even at the LUT-IOTC/IOTC level can be performed without any additional coding because the analytical energy gradient requires the LUT-IOTC transformation only for 1e atomic interaction term, which has been implemented so far. For DFT, PC-corrected density operator (δ⁺) based on LUT-IOTC can be also treated.

Table 1. Major capabilities of LUT-IOTC/IOTC.

Energy	SCF, MP2, CI, CC, DFT (w/ δ+)
Gradient	SCF, MP2, CI, CC, DFT (w/ δ+)
Others	Divide-and-conquer, semi-numerical Hessian

3 Results and discussion

Figure 1 shows the deviations in energy (ΔE) from the 4c method at the BLYP/Sapporo (-DKH3)-DZP-2012 level for diatomic molecules (CuH, Cu₂, AgH, Ag₂, AuH, and Au₂). For Hamiltonians, nonrelativistic (NR/C), LUT-IOTC for 1e term (LUT-IOTC/C), LUT-IOTC for 1e and 2e terms (LUT-IOTC/IOTC) were adopted. The effect of PC-corrected density operator (w/ δ⁺) by the LUT-IOTC method was also examined. NR/C gives large deviations from 4c. On the other hand, LUT-IOTC/IOTC w/ δ⁺ gives good agreement with 4c in all molecules. The largest deviation was 0.007 hartree in Au₂. In addition, the effects of PC-corrected 2e Coulomb term in BLYP between LUT-IOTC/C and LUT-IOTC/IOTC, and of PC-corrected density operator between LUT-IOTC/IOTC and LUT-IOTC/IOTC w/ δ⁺ are essential to reproduce the 4c results. The deviations are several hartrees.

Figure 1.

 Deviations of BLYP energies in 2c from 4c.

Table 2 compares the central processing unit (CPU) times for several 2c DFT calculations of tris-(2-phenylpyridinato) iridium (III) (Ir (ppy)₃), as illustrated in Figure 2, at the PBE0/Sapporo (-DKH3)-DZP-2012 level. For comparison, the CPU times for the nonrelativistic treatment, i.e., NR/C, are also given. In addition with the total CPU times, the components for time-consuming steps, namely, evaluation and transformation of 2e integrals (TEI) and the first SCF cycle (SCF), are shown. For measurement of the CPU times, one CPU core of Intel® Xeon® Gold 5122/3.60 GHz was used. The numbers of total SCF cycles were 39 and 31 for the nonrelativistic and 2c relativistic calculations, respectively. The DFT calculations with the fine Lebedev grid (96 radial and 302 angular points) were performed in 15 cycles for NR/C, 14 cycles for LUT-IOTC/IOTC w/ δ⁺, and 13 cycles for others. In LUT-IOTC/C, the CPU times of TEI and SCF are similar to those in NR/C. In LUT-IOTC/IOTC and LUT-IOTC/IOTC w/ δ⁺, the CPU times of TEI are approximately 16% larger than the others because of the use of local transformation. In LUT-IOTC/C w/ δ⁺ and LUT-IOTC/IOTC w/ δ⁺, the CPU times of SCF are larger than those in the others. This is because the transformation for δ operator is performed in all grid points. 4c calculation was difficult in the present computing environment. In consequence, the most accurate LUT-IOTC/IOTC w/ δ⁺ calculation can be accomplished within three times of CPU times of nonrelativistic treatment.

Table 2. CPU time [second] in Ir (ppy)₃.

	TEI	SCF	Total
NR/C	1175.17	359.35	12392.61
LUT-IOTC/C	1177.74	345.89	9994.87
LUT-IOTC/C w/ δ +	1195.72	1933.32	32360.91
LUT-IOTC/IOTC	1368.86	341.39	10025.63
LUT-IOTC/IOTC w/ δ +	1392.46	1949.37	34332.76

Figure 2.

 Ir (ppy)₃.

4 Concluding remarks

We have newly implemented the PC corrections for 2e Coulomb interaction and density operators in GAMESS program package, which enable to perform accurate and efficient 2c relativistic calculations. This letter describes theoretical aspects and performances of the implemented methods. In particular, the numerical result for the whole PC corrections of one- and two-electron operators and the density operator has been clarified to give the results close to those of 4c treatment. Furthermore, its computational costs are slightly larger than but comparable with those of the nonrelativistic calculations Therefore, we highly recommend GAMESS users to select the present method for any all-electron calculations without any hesitation after the present codes are open to the public.

Acknowledgments

Parts of the calculations were performed at the Research Center for Computational Sciences (RCCS), Okazaki Research Facilities, National Institutes of Natural Sciences (NINS). Author J.S. is grateful for the PRESTO program, "Advanced Materials Informatics through Comprehensive Integration among Theoretical, Experimental, Computational, and Data-Centric Sciences," sponsored by the Japan Science and Technology Agency (JST).

Reference

• [1]   M. Barysz, A. J. Sadlej, J. Chem. Phys., 116, 2696 (2002). doi:10.1063/1.1436462
• [2]   J. Seino, H. Nakai, J. Chem. Phys., 136, 244102 (2012). doi:10.1063/1.4729463 PMID:22755560
• [3]   J. Seino, H. Nakai, J. Chem. Phys., 137, 144101 (2012). doi:10.1063/1.4757263 PMID:23061833
• [4]   J. Seino, H. Nakai, J. Chem. Phys., 139, 034109 (2013). doi:10.1063/1.4813595 PMID:23883012
• [5]   Y. Nakajima, J. Seino, H. Nakai, J. Chem. Phys., 139, 244107 (2013). doi:10.1063/1.4850638 PMID:24387357
• [6]   M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem., 14, 1347 (1993). doi:10.1002/jcc.540141112
• [7]   Y. Nakajima, J. Seino, M. W. Schmidt, H. Nakai, J. Comput. Chem. Jpn., 15, 68 (2016). doi:10.2477/jccj.2016-0029
• [8]   Y. Ikabata, T. Oyama, M. Hayami, J. Seino, H. Nakai, J. Chem. Phys., 150, 164104 (2019). doi:10.1063/1.5090523 PMID:31042880
• [9]   J. Seino, W. Uesugi, M. Hada, J. Chem. Phys., 132, 164108 (2010). doi:10.1063/1.3397070 PMID:20441259
• [10]   J. Seino, M. Hada, Chem. Phys. Lett., 461, 327 (2008). doi:10.1016/j.cplett.2008.07.009