RAQET: Large-scale two-component relativistic quantum chemistry program package.

The Relativistic And Quantum Electronic Theory (RAQET) program is a new software package, which is designed for large-scale two-component relativistic quantum chemical (QC) calculations. The package includes several efficient schemes and algorithms for calculations involving large molecules which contain heavy elements in accurate relativistic formalisms. These calculations can be carried out in terms of the two-component relativistic Hamiltonian, wavefunction theory, density functional theory, core potential scheme, and evaluation of electron repulsion integrals. Furthermore, several techniques, which have frequently been used in non-relativistic QC calculations, have been customized for relativistic calculations. This article introduces the brief theories and capabilities of RAQET with several calculation examples.

Introduction

Both a qualitative and quantitative understanding of any elements in the periodic table is one of the most significant problems in quantum chemistry (QC). The solution is essential to theoretically promote the element strategy in Japan, whose concept is the understanding and the effective use of the roles of the key elements that determine the properties and functions of intelligent materials. One possible solution is relativistic QC (RQC), which is the combination of special relativity and QC. The fundamental treatment of RQC is based on four‐component (4c) Dirac theory, which utilizes a one‐particle Dirac operator as the Hamiltonian for large‐ and small‐components of the Dirac spinor. For many‐electron systems nonrelativistic (NR) Coulomb interaction1 and lower‐order quantum electrodynamics effects such as Breit2/Gaunt3 interaction, have been considered. Several sophisticated program packages for 4c relativistic calculations, such as RELCI,4 MCHF,5 MCDFGME,6 GRASP,7, 8, 9 and CI‐MBPT10 for atoms, MOLFDIR11 BERTHA,12, 13, 14 DIRAC,15 ReSpect,16 and BAGEL17 for atoms and molecules, have accurately predicted electronic energies, molecular geometries, and electronic/magnetic properties. However, their application to realistic systems containing hundreds or more atoms with heavy elements is still difficult due to the high computational cost.

An alternative approach to the 4c relativistic theory is a two‐component (2c) relativistic one, which relies solely upon information related to the electronic state or a large‐component of the Dirac spinor. Several approximated 2c relativistic theories have been proposed such as the Pauli approximation,18 Foldy–Wouthuysen (FW) transformation,19 Douglas–Kroll–Hess (DKH) transformation,20, 21, 22, 23, 24, 25, 26, 27 regular approximation (RA),28, 29, 30, 31 and relativistic scheme by eliminating small components (RESC).32 Furthermore, recent developments in the 2c treatment of one‐ and many‐electron molecular systems have provided a degree of accuracy comparable to the 4c treatment. These approaches are called Dirac‐exact 2c schemes, and include the normalized elimination of small components (NESC),33, 34, 35 infinite‐order DKH (IODKH),36, 37, 38, 39, 40 exact 2c (X2C),41, 42, 43 and infinite‐order 2c (IOTC) 44, 45 methods. The details and history of 2c relativistic methods are discussed in the review articles by Liu,46 Saue,47 and Peng and Reiher.48

Some one‐particle 2c relativistic Hamiltonians are now available in most QC program packages49, 50, 51, 52, 53, 54, 55 in spin‐free (SF) formalisms. The SF effect is important for describing bond contractions/elongations resulting from the relativistic shrinking/expansion of atomic orbitals.56 The reason for this easy implementation is that only the one‐electron NR integrals are replaced by 2c relativistic ones. However, a few QC program packages are able to treat spin‐dependent (SD) formalisms,15, 57, 58, 59 which include spin‐orbit interactions. For example, the open‐source NWCHEM package contains the X2C implementation with the SF and SD formalisms.60, 61 The SD formalism is important for the accurate description of molecular geometries62 and properties such as excitation energies,63 spectroscopic constants,64 and magnetic shielding constants65 in heavy‐element systems. The extension to include the SD formalism requires a considerable change of the NR program. For example, many variables such as orbital coefficients, one‐ and two‐electron integrals for the SD operators, and Fock matrices should be treated as complex numbers.

For the practical calculation of molecules containing heavy elements in the 2c relativistic framework, with accuracy close to the 4c relativistic one, further requirements must be satisfied in terms of accuracy and efficiency. For example, to maintain accuracy, the Hartree–Fock (HF) and electron correlation theories should be reformulated to describe the two‐electron SD interaction.42, 43, 45, 66, 67, 68 In addition, the 2c transformation changes the picture of molecular properties and interactions, which is termed the picture change (PC) effect.69, 70, 71, 72, 73, 74, 75, 76 This effect also contributes to density functional theory (DFT) calculations through the PC of electron density.31, 77 The 2c transformation should be quicker78, 79, 80, 81 and therefore more efficient. The basic theory and algorithms should be reconstructed for calculations involving molecules which contain heavy elements82, 83, 84 and for SD calculations.85, 86, 87, 88 Furthermore, large‐scale molecular theories with overall linear‐scaling techniques are required for calculations relating to large molecules.89

This article introduces a new 2c RQC program package, the Relativistic And Quantum Electronic Theory (RAQET), which satisfies the above requirements for accurate and practical relativistic calculations. This article is organized as follows: we begin with a brief overview of RAQET as well as an explanation of the general theory in the following section. Subsequent sections discuss the capabilities of the Hamiltonian, wavefunction theory (WFT), DFT, and integral evaluations together with a brief explanation of the theories. Basic run types and useful options are given then the following section provides the graphical user interface (GUI) for user‐friendly calculations. Finally, we discuss the performance such as the parallel capability for large systems before making some concluding remarks

Overview of RAQET

This section provides an overview of RAQET. In this program, the fundamental equation for many‐electron systems is solved in the 2c formalism, defined bywhere is the 2c wavefunction for the electronic state, E is an energy eigenvalue, and is the many‐electron 2c relativistic Hamiltonian given bywhere and are the electronic components of the one‐ and two‐electron Hamiltonians. This program mainly adopts the 2c IODKH/IODKH method40 for the Hamiltonian. This method approximates the block‐diagonalized many‐electron 4c Dirac Hamiltonian aswithandHere, is the one‐electron 4c Dirac Hamiltonian, and is the general expression of the two‐electron 4c operator. In RAQET, the simplest Coulomb interaction is adopted as G 4(i, j), i.e. G 4(i, j) = 1/r . In eq. (5), the many‐electron unitary transformation is approximated by the product of the one‐electron unitary transformation. The subsequent sections discuss the capabilities of the Hamiltonian , WFT for , and DFT in RAQET.

RAQET utilizes Gaussian‐type orbitals (GTOs) as basis functions. In heavy‐element systems, the molecular orbitals have high angular momentum such as d‐ and f‐orbitals and several nodes between the core and the valence regions. To describe these orbitals accurately, GTOs can be constructed from many primitive functions. Two types of GTOs are frequently adopted for molecular basis sets containing heavy elements: the segmented‐/general‐contractions (SC/GC),90 in which GTOs are treated as a linear combination of primitive functions with the different/same set of exponents. To calculate molecular integrals, such as the electron repulsion integrals (ERIs) and their derivatives (DERIs), which are the time‐consuming steps in QC calculations, the algorithm should be optimized for each type of contracted GTO. RAQET employs optimized algorithms for the evaluation of integrals for both SC and GC.83, 84, 91 The section on integral evaluation discusses the capabilities of the algorithms for ERIs and DERIs.

Hamiltonian

Several 2c relativistic Hamiltonians are available in RAQET, which are summarized in Table 1. For the one‐electron parts, of the 2c relativistic Hamiltonian in eq. (2), the first‐order DKH (DKH1) and IODKH have been implemented as well as the NR Hamiltonian. For the two‐electron parts, the NR Coulomb and IODKH transformed interactions were implemented.

Table 1

Capabilities of Hamiltonians

Keyword	h2+	g2++	Spin	Note
NR/NR	NR	NR	SF
DKH1/NR	DKH1	NR	SF/SD
IODKH/NR	IODKH	NR	SF/SD
LUT‐IODKH/NR	LUT‐IODKH	NR	SF/SD
IODKH/IODKH	IODKH	IODKH	SF/SD
LUT‐IODKH/LUT‐IODKH	LUT‐IODKH	LUT‐IODKH	SF	Default

In the 2c relativistic calculations, the unitary transformation, U in eq. (3), is determined by the matrix transformation method and the resolution of identity technique proposed by Hess.21, 22 The unitary transformation is a potential bottleneck because the transformations of one‐ and two‐electron Hamiltonians for the entire system scale as O(n 3) and O(n 5), respectively. Here, n denotes the number of primitive basis functions. To reduce this scaling, local unitary transformation (LUT) schemes for one‐78 and two‐electron81 Hamiltonians are implemented. The LUT scheme approximates the one‐electron unitary transformation, U(i) in eq. (5), as a block‐diagonal form of the subsystem contributions,where {A, B, …} denote subsystems that do not intersect. RAQET adopts an atomic partition as the subsystem. The LUT scheme has been designed based on the locality of relativistic interactions. For the one‐electron Hamiltonian, the relativistic effect is dominant for each atom and in interatomic interactions with nearest‐neighbor atoms. Other interactions behave like those in NR. The cutoff distance is introduced to separate the interatomic interactions, whose default value is set to 3.5 Å. For two‐electron interactions, the relativistic effect is dominant for each atom and the interatomic interactions can be approximated as the NR ones. Thus, the scaling for the transformations of one‐ and two‐electron Hamiltonians become O(n 1).

The one‐ and two‐electron parts in the 2c relativistic Hamiltonians of eq. (2) are divided into the SF and SD terms respectively asFor division, the Dirac relation is applied to operators in the transformed Hamiltonian, defined bywhere is the vector of three Pauli spin matrices, is the momentum operator, X is the arbitrary operator, and i is the imaginary unit. In RAQET, all 2c relativistic Hamiltonians for one‐electron interactions are implemented based on both the SF and SD formulations. For two‐electron interactions, the SF formulations of IODKH and LUT‐IODKH for are implemented. For the SD term of IODKH and LUT‐IODKH for , further investigations will be discussed in future. From the viewpoints of accuracy and efficiency, LUT‐IODKH/LUT‐IODKH is the default option in RAQET for calculations relating to heavy‐element systems.

Wavefunction

Hartree–Fock

There are various HF scheme modules in RAQET which are summarized in Table 2. For SF calculations in the 2c relativistic framework, the restricted HF (RHF), unrestricted HF (UHF), and open‐shell HF (ROHF) are available. In ROHF, we can select various parameters as keywords92, 93, 94, 95, 96, 97, 98, 99 which are listed, as well as the parameter values, in the Supporting Information.

Table 2

Capabilities of HF schemes

Spin	Keyword	Restriction of symmetrya
SF	RHF	S 2, S z, K, K 0
	UHF	S z, K 0
	ROHF	S 2, S z, K, K 0
SD	KRHF	K
	KUHF	K (partially)
	KROHF	K
	GHF	–

S 2, Sz, K, and K 0 mean square and z‐component of spin, time‐reversal, and complex conjugation symmetries, respectively.

For SD calculations, the Kramers‐restricted HF (KRHF),86 Kramers‐unrestricted HF (KUHF),87 Kramers‐restricted open‐shell HF (KROHF),88 and generalized HF (GHF) 85 are available. These schemes use spinors , which are expressed as a linear combination of ‐ and ‐spin basis functions,where C is the complex orbital coefficient matrix, and is the basis function. The orbitals represented in eq. (9) require the treatment of complex numbers to describe the SD contributions, such as the spin‐orbit interaction, in the 2c relativistic Hamiltonian. In addition, the dimension of the basis function space increases to twice its original size to represent the ‐ and ‐spin mixing. Due to these differences in the SF formulations, the whole procedures in KRHF, KROHF, KUHF, and GHF have been reconstructed for the SD calculations.

The working equation in GHF is derived using eq. (9) without any symmetry constraints, which, given in block form, iswhere , , and are the Fock matrix, the overlap matrix, and orbital energies, respectively. KRHF preserves only time‐reversal symmetry, and the working equation is defined bywhere the presuperscript Q denotes a quaternion matrix, is one of the quaternion units. KRHF works as a generalization of the spin restriction and induces twofold degeneracy for the orbitals. Thus, KRHF is regarded as the relativistic counterpart of RHF. The difference between KRHF and RHF is that the Fock matrix, orbital coefficients, and density matrix are defined in quaternionic space rather than real space. Using the characteristics of KRHF, KUHF, and KROHF with quaternionic matrices, the algorithms become similar to those for RHF, UHF, and ROHF, respectively. The working equation for KUHF is written aswhere the superscripts and denote the pseudo‐alpha and pseudo‐beta spins, respectively. The working equation for KROHF is defined bywhere is Roothaan's effective Fock operator extended to quaternions,The subscripts C, O, and V show closed‐, open‐, and virtual‐shell MOs, respectively. The same manner and parameters can be used in eq. (14) as for ROHF.

KROHF preserves the time‐reversal symmetry while KUHF partially breaks this symmetry for open‐shell calculations. KRHF, KUHF, and KROHF have the advantage that they can perform with smaller degrees of freedom than GHF. These constraints are effective for rapid SCF convergence.87, 88

Post‐Hartree–Fock

Electron correlation schemes are implemented in RAQET and they are summarized in Table 3. In SF calculations, the Møller–Plesset perturbation theory (MPPT), including MP2, MP2.5, and MP3, and the coupled‐cluster (CC) theory, including LCCD, CCD, LCCSD, CCSD, CCSD[T], CCSD(T), CCSDT, CCSDTQ, are available. These were generated using the tensor contraction engine,100, 101, 102 a computerized symbolic algebra system. All of these schemes can be combined with the divide‐and‐conquer (DC) technique,89, 103, 104, 105, 106, 107, 108 which is one of the fragmentation‐based linear‐scaling electron correlation approaches.

Table 3

Capabilities of post‐HF schemes with and without the DC technique

Correlation method	Keyword	Spin
		w/o DC	w/ DC
MP	MP2	SF/SD	SF
	MP2.5	SF	SF
	MP3	SF	SF
CC	LCCD	SF	SF
	CCD	SF	SF
	LCCSD	SF	SF
	CCSD	SF	SF
	CCSD[T]	SF	SF
	CCSD(T)	SF	SF
	CCSDT	SF	SF
	CCSDTQ	SF	SF

MP2 is available in SD calculations and the SD formalism of MP2, which is termed generalized MP2 (GMP2), is reconstructed using the molecular orbitals of eq. (9) in GHF.68 The working equation is defined aswithHere, {i, j} and {a, b} are the occupied and virtual orbitals, respectively while {, } and {, } are the spin functions for electrons 1 and 2, respectively. For the SD term of the two‐electron 2c relativistic Hamiltonian in eq. (7), arbitrary combinations of spins are required for the evaluation of operators such as , , and , where X is an arbitrary two‐electron operator.

Core potential method

RAQET contains three core potential methods based on the valence‐only Huzinaga–Cantu equation:109 the model core potential (MCP),110, 111 ab initio model potential (AIMP),112, 113 and frozen core potential (FCP) methods.114, 115 These schemes reduce the number of electrons treated explicitly by replacing the core electrons with the potential. Thus, an efficient HF calculation is realized. These core potential methods are implemented in the SF formulations at the HF level of the theory. MCP and AIMP include the model potentials for the core region, whose parameters were determined as to reproduce reference relativistic calculations such as the one‐electron third‐order DKH Hamiltonian with the two‐electron NR Coulomb interactions. However, in FCP, the core potentials have been constructed from the frozen atomic core orbitals instead of atomic model potentials. Any 2c relativistic Hamiltonians, including (LUT‐)IODKH/NR and (LUT‐)IODKH/(LUT‐)IODKH, are seamlessly combined with FCP. In addition, if required, the core electrons can be self‐consistently relaxed using frozen molecular valence potentials following the valence SCF calculation.115 This relaxation is essential to evaluate the molecular properties related to core orbitals such as core level shifts.

Density Functional Theory

Ground and excited state calculation

DFT calculations, based on the SF and SD formulations for the ground state, are executable in RAQET. Table 4 summarizes the available exchange–correlation (XC) functionals. The Slater–Dirac exchange,116, 117 Vosko–Wilk–Nusair (VWN) correlation,118 and Perdew–Wang correlation119 are implemented as a local density approximation (LDA) functional. As for the generalized gradient approximation (GGA), the Becke's exchange proposed in 1988 (Becke88),120 Lee–Yang–Parr (LYP) correlation,121 Perdew–Burke–Ernzerhof (PBE) exchange–correlation,122 and revised PBE exchange123 functionals have been implemented. The 1998 Van Voorhis–Scuseria meta‐GGA (VS98) 124 and 2006 Minnesota local meta‐GGA (M06‐L) 125 functional are also included. The hybrid functionals with various degrees of HF exchange,126, 127, 128, 129, 130, 131, 132 as well as the double‐hybrid functionals combined with the HF exchange and the MP2 correlation,133, 134, 135, 136 are available. Some of the functionals have been implemented on the basis of Fortran code repositories.137, 138 The local response dispersion (LRD) method139, 140, 141 is included as a dispersion correction method. For excited states, the time‐dependent DFT (TDDFT) is implemented at the SF level to evaluate excited state energies, electronic transition moments, and oscillator strengths. TDDFT calculations based on the SD formalism will be available after the report on the numerical accuracy for the case of the highly accurate 2c Hamiltonian.

Table 4

Capabilities of XC functionals for DFT calculations

Type	Keyword	Ratio of HF exchange	Ratio of MP2 correlation
LDA	SVWN	0%	0%
	SPW92	0%	0%
GGA	BLYP	0%	0%
	PBE	0%	0%
	revPBE	0%	0%
Meta‐GGA	VS98	0%	0%
	M06‐L	0%	0%
Hybrid	B3LYP	20%	0%
	BHHLYP	50%	0%
	PBE0	25%	0%
	M05	28%	0%
	M05‐2X	56%	0%
	M06	27%	0%
	M06‐2X	54%	0%
	M06‐HF	100%	0%
Double‐hybrid	B2PLYP	53%	27%
	B2GPPLYP	65%	36%
	PBE0‐DH	50%	12.5%
	PBE0‐2	79.3701%	50%

Picture change in DFT

Any calculation of the expectation values for molecular properties should consider the PC correction to reproduce the results of the 4c relativistic calculations. In DFT calculations with the 2c relativistic framework, we should consider the PC problem in the electron density ,77 which originates from the unitary transformation of the 4c wavefunction from to whereThe electron density at position r is defined for both the 4c and 2c formalisms bywhere is the delta function, and 1 n is the n × n identity matrix. The subscript elec denotes the electronic component of the unitary transformed delta function. This PC correction is now available for ground state DFT calculations. Furthermore, the LUT scheme will realize an efficient PC correction process for eq. (18) following the confirmation of the locality of the PC correction in an earlier article.77

Integral Evaluation

Efficient algorithms to evaluate two‐electron integrals are implemented in RAQET, and are summarized in Table 5. The Gaussian quadrature technique,142, 143 one of the most effective algorithms for ERIs composed of GTOs with high angular momentum and short contraction, has been adopted for Gauss–Rys. This algorithm is useful for the two‐electron IODKH transformation because the transformation is performed in the primitive basis function space. The Pople–Hehre and McMurchie–Davidson (PHMD) algorithm144, 145, 146 utilizes the expression of a linear combination of Hermite Gaussians in rotated space, which can be useful for calculations involving molecules which contain only light elements. PHMD accepts integrals for s, p, and d functions. The source code for PHMD was taken from an integral library in open‐source software for massively parallel QC calculations called Scalable Molecular Analysis Solver for High‐performance computing systems (SMASH) which was developed by Ishimura.147 The accompanying coordinate expansion and transferred recurrence relation (ACE–TRR) algorithm employs the ACE formula148 with the ACE‐vertical recurrence relations (RRs) and transfer relations.84 ACE provides the most effective expression for GTOs with long contractions because it is designed to minimize the floating operation count, according to the contraction length, by optimal coordinate expansion. The RR scheme reduces the operation count and improves memory usage for ERIs with high angular momentum GTOs. Furthermore, ACE–TRR is optimized for the SC and GC basis sets. This algorithm is effective for heavy‐element systems, whose basis sets contain long contractions and functions with high angular momentum. RAQET implements a hybrid of Gauss–Rys, PHMD, and ACE–TRR, and switches to the optimal algorithm according to the angular momentum, contraction length, and contraction type by default. Each algorithm is also designed for two orbital shapes of the Cartesian coordinates and spherical harmonics. Derivative two‐electron integral evaluations, such as DERIs, Gauss–Rys, ACE–TRR, and the hybrid algorithm are available up to second‐order derivatives in both SC and GC formulations.

Table 5

Capabilities of ERI evaluations

Keyword	Contraction typea	Orbital shape[b]	Note
GaussRys	SC	XYZ/SH	Used in 2‐electron IODKH transformation
PHMD	SC	XYZ/SH	For s, p, and d functions
ACE‐TRR	SC/GC	XYZ/SH	For all functions
Hybrid	SC/GC	XYZ/SH	For all functions (default)

SC and GC mean the optimized algorithms for segmented‐contraction and general‐contraction, respectively. [b] XYZ and SH mean Cartesian and spherical harmonics, respectively.

Basic Run Types

Three fundamental run types are now available in RAQET: single‐point energy, geometry optimization, and molecular property. Single‐point energy is available for all of the Hamiltonians, wavefunctions, and DFTs mentioned in the previous section. The details of the geometry optimization and property calculations are explained in subsequent sections.

Geometry optimization

For geometry optimization, the energy gradients of Hamiltonians, wavefunctions, and DFT formulations are required to evaluate the forces at the nuclei. The analytical gradients of NR, SF‐IODKH/NR, SF‐LUT‐IODKH/NR, and SD‐LUT‐IODKH/NR, as well as the numerical gradients for all of the Hamiltonians, are available in RAQET. The LUT analytical scheme is the default option for the calculations because it requires considerably less computational time than conventional 2c relativistic Hamiltonians149, 150, 151 with similar accuracy. For wavefunctions and DFT, the analytical gradients of RHF, RDFT, UHF, UDFT, RMP2, and GHF have been implemented. RAQET utilizes the redundant internal coordinates152 (default) and Cartesian coordinates for geometry optimization. Furthermore, rational functional optimization,153 which uses calculated forces, was adopted to change structures. The source codes of the Cartesian coordinates, redundant internal coordinates, and rational functional optimization were taken from SMASH.

Analysis and molecular property

RAQET can perform several analyses and evaluations of molecular properties such as the populations, electric, and magnetic properties in the WFT framework. The Mulliken population analyses are available to obtain charges at the NR and SF levels. In addition, atomic spin vectors can be evaluated in GHF based on the Mulliken‐type atomic division at the NR, SF, and SD levels. For the electric properties, first‐ and second‐order analytical derivatives of the electronic energy with respect to the electric field were implemented to evaluate the electric dipole moment and static polarizability, which are formulated based on both the NR and IODKH/NR Hamiltonians at the SF and SD levels. For the magnetic properties, the first‐ and second‐order analytical derivatives with respect to the external magnetic field and/or the nuclear magnetic moment were implemented to evaluate the magnetic dipole moment and the nuclear magnetic shielding constant. These schemes are constructed in the IODKH/NR schemes at the SD level as well as the NR one. To perform gauge‐origin independent calculations in the NR framework, RAQET employs the gauge‐including/independent atomic orbital (GIAO).154, 155 In the 2c relativistic framework, RAQET employs the gauge factorization before unitary transformation (GF‐UT) formulation, which has been derived for gauge‐origin independent calculations by introducing the gauge factor and 2c unitary transformation into the molecular orbitals.156 The SD relativistic effect, which can include the spin‐dipolar and the Fermi‐contact contributions, is essential for molecules containing heavy elements,157 particularly for calculation of the nuclear magnetic shielding constant. The second‐order properties with respect to both the electric and magnetic fields require analytical derivatives of the density matrices. RAQET utilizes the coupled perturbed HF framework based on RHF, UHF, and GHF, which is solved by an iterative procedure.158, 159

Useful Options

Basis sets

RAQET has a basis set database for the easy preparation of input files. This contains various basis sets optimized for both NR and 2c relativistic calculations. For example, the NR basis sets contain STO‐types such as STO‐3G and STO‐6G, Pople‐types such as 6‐31G and 6‐311G(d), Dunning‐types such as cc‐pVDZ and aug‐cc‐pVTZ, and Sapporo‐types such as Sapporo‐DZP and Sapporo‐TZP+d. The relativistic basis sets contain cc‐pVDZ‐DK, Sapporo‐DKH3‐DZP, SARC‐DKH, DKH3‐Gen‐TK/NOSec‐V‐TZP, and ANO‐RCC. Furthermore, users can manually add new basis sets to the database. The details of the available basis sets are summarized in the Supporting Information.

Initial guess

Users can choose five types of initial guesses for the SCF calculations. Table 6 summarizes the techniques implemented in RAQET. The keyword Hcore produces the MO coefficients and density matrix from the diagonalization of the one‐electron Hamiltonian matrix. This guess can represent the relativistic effects of the initial core orbitals at the SF and SD levels without any two‐electron interactions. The keyword Huckel, which is the default option in RAQET, provides the MO coefficients and density matrix from the extended Hückel calculation.160 RAQET adopts the energy parameters reoptimized from the atomic calculations using the SF‐IODKH/NR Hamiltonian and minimal basis sets optimized by the DKH3 Hamiltonian161 at the ROHF level. This gives a good initial description of the chemical bonds at the SF level. The keyword Atomic provides the density matrix constructed from atomic density matrices. This represents the relativistic effects for core orbitals together with intra‐atomic two‐electron interactions at the SF and SD levels. The keyword Small provides the MO coefficients by projecting those generated from a SCF calculation using smaller basis sets. This guess gives good initial core and valence orbitals for any system at the SF and SD levels. Finally, the keyword MORead uses the MO coefficients written in an input file.

Table 6

Capabilities of initial guess techniques

Keyword	Note
Hcore	Use MO information obtained from diagonalization of one‐electron Hamiltonian matrix
Huckel	Use MO information obtained with extended Hückel method (default)
Atomic	Use the sum of atomic densities
Small	Use MO information by projecting MO obtained from SCF calculations with smaller basis sets
MORead	Read MO information from input file

SCF convergence

RAQET contains several convergence techniques for SCF calculations, which are summarized in Table 7. All of these techniques are available for KRHF, KROHF, KUHF, and GHF using complex numbers as well as RHF, ROHF, and UHF using real numbers. Now, users can use all of the convergence techniques for all SCF calculations at the NR, SF, and SD levels. The keyword sDamp utilizes the static damping technique, which gives a mixture of the density matrices from previous and present iterations with a weighting factor. This technique gives stable results, though the rate of convergence is slow. The zero value of the weighting factor gives a result which corresponds to the fixed‐point algorithm, which in turn corresponds to the SCF calculation without any acceleration techniques. The keywords C1‐DIIS and C2‐DIIS are Pulay's techniques of direct inversion in the iterative subspace (DIIS),162, 163 which are extrapolation schemes using a linear combination of the density or Fock matrices from previous SCF iterations. The extrapolation coefficients are optimized by minimizing the DIIS error vector. These techniques show rapid convergence when the initial guess is good. The keyword EDIIS uses an interpolation scheme based on the DIIS technique. The interpolation coefficients are optimized by minimizing the HF energy functional, which contains a linear combination of the density matrices.164 This technique is effective when the initial guess is poor, although the convergence rate is slower close to the minimum than DIIS. A combination of C2‐DIIS and EDIIS (DIIS + EDIIS) 165, 166, 167 is also available, and is the default option in RAQET. This technique has the advantages of both DIIS and EDIIS and provides rapid convergence even if the initial guess is poor. The keyword SOSCF uses second‐order SCF orbital optimization. This technique utilizes an exponential parameterization to obtain the orbital rotation angles, and requires the calculation of first and approximated second derivatives of the energy with respect to the orbital coordinates. This provides rapid convergence when the initial guess is good as well as DIIS.168, 169 The keyword FON adopts the fractional occupation number based on the Fermi distribution function.170 This technique can be useful for calculations involving degenerate systems such as metal clusters.

Table 7

Capabilities of SCF convergence techniques

Keyword	Note
sDamp	Static damping technique for density matrix
C1‐DIIS	DIIS extrapolation of Fock matrix using error vector from previous Fock matrix
C2‐DIIS	DIIS extrapolation of Fock matrix using error vector from SCF condition
EDIIS	Energy DIIS interpolation of Fock matrix
SOSCF	Second‐order SCF orbital optimization
FON	Use fractional occupation number based on the Fermi distribution function
EDIIS+DIIS	Hybrid of EDIIS and C2‐DIIS (default)

Graphical User Interface

RAQET utilizes sophisticated molecular modeling and visualization software, Winmostar.171 This software contains modules for easy creation of an input file and a visualization of the calculated results. Figure 1 shows an example of how an input file is created. We select the computational conditions such as the run type, Hamiltonian, wavefunction, basis sets, and so on from list and combo boxes. After the calculations in RAQET, Winmostar can use the output files to visualize the molecular structures and properties. As an example, Figure 2 shows the visualization of the atomic spin vectors for a Rb3 molecule. The vectors were obtained using the GHF method with the IODKH/NR Hamiltonian at the SD level, which indicates the non‐collinear spin states.

Figure 1

Input example of RAQET program on GUI. [Color figure can be viewed at wileyonlinelibrary.com]

Figure 2

Visualization of atomic spin expectation vectors for Rb3 on GUI. [Color figure can be viewed at wileyonlinelibrary.com]

Parallel Capability

RAQET partially supports a parallel computational environment using the hybrid parallelization technique with OpenMP/MPI modules. OpenMP parallelization is compatible with single‐point energy calculations, geometry optimizations, and molecular property calculations while MPI parallelization is only compatible with single‐point energy calculations. Users can choose a sequential, OpenMP/MPI, or hybrid process according to their computational environment. Figure 3 shows the parallelization efficiencies of OpenMP in a single‐point energy calculation for an Ir(ppy)3 molecule and Au10 cluster at the SF‐IODKH/NR level in UHF. The calculations were performed using an Intel(R) Xeon(R) CPU E7‐4870 with 512GB memory. The horizontal and vertical axes show the number of cores used and the effective cores, respectively. The number of effective cores was evaluated as t(n core)/t(2 core), where t is the wall time. The curve in RAQET is close to the ideal efficiency, which is indicated by the diagonal line. When the number of cores is less than 40, the parallel efficiency is greater than 99%. When the number is 80 cores, the efficiency is greater than 85%.

Figure 3

Parallel efficiency for calculations of single‐point energies for the Au10 cluster and Ir(ppy)3 complex. [Color figure can be viewed at wileyonlinelibrary.com]

Calculation Example

This section demonstrates the performance of RAQET. First, the computational time was measured in a single‐point energy calculation for the Au20 cluster. The calculation was performed using NR, IODKH/NR, LUT‐IODKH/NR, and LUT‐IODKH/LUT‐IODKH Hamiltonians at the SF level in the RHF framework. For other options, the DIIS + EDIIS technique for the SCF convergence, MO coefficients obtained from the extended Hückel calculation for the initial guess, and DKH3‐Gen‐TK/NOSec‐V‐TZP for the basis set, were employed. The calculations were performed on 32 cores of an Intel(R) Xeon(R) CPU E5‐2683 v4 with 64GB memory. Table 8 shows the wall time for each calculation step, the number of SCF iterations, and the converged total energies. Here, the SCF step includes evaluations of the NR two‐electron integrals (TEIs) because the direct SCF technique was adopted. For LUT‐IODKH/LUT‐IODKH, the relativistic two‐electron integrals for each atom were stored in disk storage before the SCF step. In NR/NR, two steps are required: one‐electron integrals (OEIs) and SCF. In (LUT‐)IODKH/NR, the transformation step for the OEIs based on (LUT‐)IODKH is included. In (LUT‐)IODKH/LUT‐IODKH, the transformation step for TEI based on LUT‐IODKH is also included. These relativistic frameworks also require additional OEIs and TEIs for these transformations. From the results, we see that the wall time for IODKH/NR is larger than NR/NR. This is due to the transformation of the OEIs. However, the wall time in LUT‐IODKH/NR and LUT‐IODKH/LUT‐IODKH are similar to that for NR/NR. Thus, the LUT scheme with the relativistic IODKH transformation for both OEIs and TEIs realizes accurate calculations for large molecules with little additional computational time compared with NR calculations.

Table 8

Wall time in minutes and the total energy in Hartree for the Au20 cluster using NR, IODKH/NR, LUT‐IODKH/NR, and LUT‐IODKH/LUT‐IODKH Hamiltonians at the SF level in the RHF framework

Step	NR	IODKH/NR	LUT‐IODKH/NR	LUT‐IODKH/LUT‐IODKH
Wall time [min]
OEI	0.06	0.17	0.06	0.06
Transformation of OEI	–	5.42	0.04	0.05
Transformation of TEI	–	–	–	5.99
SCF	2260.46	1542.95	1465.85	1722.98
SCF per iter.	27.91	35.07	36.65	31.91
Total	2260.69	1548.71	1466.12	1729.25
Number of iteration	81	44	40	54
Total energy [Hartree]	−327,952.2329	−380,224.3549	−380,224.3549	−380,378.6063

Second, the computational time of the electron correlation schemes were evaluated for single‐point energy calculations. Three systems were adopted: zigzag‐structured (HAt) molecules (n = 5, 10, 20, and 30), Hg2+ in 100 water molecules, and bis‐(b‐carotene)Pt10. The calculations were performed using the LUT‐IODKH/NR Hamiltonian at the SF level in the RMP2 framework with and without the DC technique. The convergence was set to DIIS + EDIIS, the initial guess was Hückel, and the basis sets were Sapporo‐DZP. These calculations were performed on 16 cores of an Intel(R) Xeon(R) CPU E5‐2690 with 128GB memory. Table 9 shows the wall time for evaluations of the correlation energy. In larger molecules, the DC technique effectively reduces the computational time and enables calculations to be performed in a practical time. In particular for (HAt), linear scaling was achieved with respect to n.

Table 9

Wall time in minutes of MP2 calculations for systems containing heavy elements with and without the DC technique

System	Wall time [min]
	w/ DC	w/o DC
(HAt)5	91.9	69.9
(HAt)10	217.8	558.9
(HAt)20	462.0	1743.9
(HAt)30	675.5	–
Hg2+ with 100 H2O	1520.3	–
Bis‐(b‐carotene) Pt10	1842.5	–

Summary and Future Development

The RAQET program has been developed for QC calculations of systems containing heavy elements based on the 2c relativistic framework. In this article, we have introduced the main features of RAQET in terms of the 2c relativistic Hamiltonian, WFT, DFT, core potential method, two‐electron integral techniques, molecular properties, initial guesses, SCF convergence techniques, GUI, and hybrid parallelization. Several techniques were customized to allow efficient calculations for heavy elements and the relativistic formulations including the SF and SD Hamiltonians. Higher‐order electron correlation calculations provided from the tensor contraction engine with the DC technique are also available in the relativistic framework. RAQET has been registered in MateriApps,172 which is Japanese portal site of materials science simulation for computational materials science researchers, theoreticians, experimentalists, and computer scientists. Furthermore, RAQET is partially open to the public. We hope that RAQET will become a useful tool to push the frontiers of the element strategy.

We hope to make RAQET a more sophisticated package by adding extensions designed in the 2c relativistic framework. For the full SD treatment of energy and higher‐order electron correlation calculations, the SD formulation of two‐electron interactions will be implemented. In addition to the grand state calculation, excited state calculation schemes in both WFT and DFT are under development. For analyses of thermodynamics and reactions, we can use frequency calculations, transition‐state optimizations, and reaction path searches. The solvent effect will be also required for accurate analysis and prediction. For the further high‐performance computation, especially on supercomputer systems, massive‐parallelization will be adopted.

Acknowledgment

Some of these calculations were performed at the Research Center for Computational Science (RCCS), the Okazaki Research Facilities, and the National Institutes of Natural Sciences (NINS). M. H. is grateful to the Japan Society for the Promotion of Science (JSPS) for providing a research fellow. We acknowledge Dr. Ishimura for providing parts of libraries and advising on the designation of the RAQET program. We are also grateful to X‐ability Co. Ltd. for developing the GUI module of RAQET.

Supporting Information

Click here for additional data file.