Search

Density functional calculations on atoms are often used for determining accurate initial guesses as well as generating various types of pseudopotential approximations and efficient atomic-orbital basis sets for polyatomic calculations. To reach the best accuracy for these purposes, the atomic calculations should employ the same density functional as the polyatomic calculation. Atomic density functional calculations are typically carried out employing spherically symmetric densities, corresponding to the use of fractional orbital occupations. We have described their implementation for density functional approximations (DFAs) belonging to the local density approximation (LDA) and generalized gradient approximation (GGA) levels of theory as well as Hartree-Fock (HF) and range-separated exact exchange [S. Lehtola, Phys. Rev. A 2020, 101, 012516]. In this work, we describe the extension to meta-GGA functionals using the generalized Kohn-Sham scheme, in which the energy is minimized with respect to the orbitals, which in turn are expanded in the finite element formalism with high-order numerical basis functions. Furnished with the new implementation, we continue our recent work on the numerical well-behavedness of recent meta-GGA functionals [S. Lehtola and M. A. L. Marques, J. Chem. Phys. 2022, 157, 174114]. We pursue complete basis set (CBS) limit energies for recent density functionals, and find many to be ill-behaved for the Li and Na atoms. We report basis set truncation errors (BSTEs) of some commonly used Gaussian basis sets for these density functionals and find the BSTEs to be strongly functional dependent. We also discuss the importance of density thresholding in DFAs and find that all of the functionals studied in this work yield total energies converged to $0.1\ \mu E_{h}$ when densities smaller than $10^{-11}a_{0}^{-3}$ are screened out., Comment: 30 pages, 2 figures. Fixed some typos

We have recently described the implementation of atomic electronic structure calculations within the finite element method with numerical radial basis functions of the form $\chi_{\mu}(r)=r^{-1}B_{\mu}(r)$, where high-order Lagrange interpolating polynomials (LIPs) were used as the shape functions $B_{\mu}(r)$ [S. Lehtola, Int. J. Quantum Chem. 119, e25945 (2019)]. In this work, we discuss how $\chi_{\mu}(r)$ can be evaluated in a stable manner at small $r$ and also revisit the choice of the shape functions $B_{\mu}(r)$. Three kinds of shape functions are considered: in addition to the $\mathcal{C}^{0}$ continuous LIPs, we consider the analytical implementation of first-order Hermite interpolating polynomials (HIPs) that are $\mathcal{C}^{1}$ continuous, as well as numerical implementations of $n$-th order ($\mathcal{C}^{n}$ continuous) HIPs that are expressed in terms of an underlying high-order LIP basis. Furnished with the new implementation, we demonstrate that the first-order HIPs are reliable even with large numbers of nodes and that they also work with non-uniform element grids, affording even better results in atomic electronic structure calculations than LIPs with the same total number of basis functions. We demonstrate that discontinuities can be observed in the spin-$\sigma$ local kinetic energy $\tau_{\sigma}$ in small LIP basis sets, while HIP basis sets do not suffer from such issues; however, either set can be used to reach the complete basis set limit with smooth $\tau_{\sigma}$. Moreover, we discuss the implications of HIPs on calculations with meta-GGA functionals with a number of recent meta-GGA functionals, and find most Minnesota functionals to be ill-behaved. We also examine the potential usefulness of the explicit control over the derivative in HIPs for forming numerical atomic orbital basis sets, but find that confining potentials are still likely a better option., Comment: 26 pages, 9 figures. Fixed typos

Sun et al. [J. Chem. Phys. 144, 191101 (2016)] suggested that common density functional approximations (DFAs) should exhibit large energy errors for excited states as a necessary consequence of orbital nodality. Motivated by self-interaction corrected density functional calculations on many-electron systems, we continue their study with the exactly solvable $1s$, $2p$, and $3d$ states of 36 hydrogenic one-electron ions (H-Kr$^{35+}$) and demonstrate with self-consistent calculations that state-of-the-art DFAs indeed exhibit large errors for the $2p$ and $3d$ excited states. We consider 56 functionals at the local density approximation (LDA), generalized gradient approximation (GGA) as well as meta-GGA levels, also including several hybrid functionals like the recently proposed machine-learned DM21 local hybrid functional. The best non-hybrid functional for the $1s$ ground state is revTPSS. The $2p$ and $3d$ excited states are more difficult for DFAs as Sun et al. predicted, and LDA functionals turn out to yield the most systematic accuracy for these states amongst non-hybrid functionals. The best performance for the three states overall is observed with the BHandH global hybrid GGA functional, which contains 50% Hartree-Fock exchange and 50% LDA exchange. The performance of DM21 is found to be inconsistent, yielding good accuracy for some states and systems and poor accuracy for others. Based on these results, we recommend including a variety of one-electron cations in future training of machine-learned density functionals., 11 pages, 4 figures, 1 table

The recombination ("dimerization") of peroxyl radicals (RO2 center dot) is one of the pathways suggested in the literature for the formation of peroxides (ROOR', often referred to as dimers or accretion products in the literature) in the atmosphere. It is generally accepted that these dimers play a major role in the first steps of the formation of submicron aerosol particles. However, the precise reaction pathways and energetics of RO2 center dot + R'O-2 center dot reactions are still unknown. In this work, we have studied the formation of tetroxide intermediates (RO4R'): their formation from two peroxyl radicals and their decomposition to triplet molecular oxygen (O-3(2)) and a triplet pair of alkoxyl radicals (RO center dot). We demonstrate this mechanism for several atmospherically relevant primary and secondary peroxyl radicals. The potential energy surface corresponds to an overall singlet state. The subsequent reaction channels of the alkoxyl radicals include, but are not limited to, their dimerization into ROOR'. Our work considers the multiconfigurational character of the tetroxides and the intermediate phases of the reaction, leading to reliable mechanistic insights for the formation and decomposition of the tetroxides. Despite substantial uncertainties in the computed energetics, our results demonstrate that the barrier heights along the reaction path are invariably small for these systems. This suggests that the reaction mechanism, previously validated at a multireference level only for methyl peroxyl radicals, is a plausible pathway for the formation of aerosol-relevant larger peroxides in the atmosphere.

We have assessed the accuracy of the magnetic properties of a set of 51 density functional approximations, including both recently published and already established functionals. The accuracy assessment considers a series of 27 small molecules and is based on comparing the predicted magnetizabilities to literature reference values calculated using coupled-cluster theory with full singles and doubles and perturbative triples [CCSD(T)] employing large basis sets. The most accurate magnetizabilities, defined as the smallest mean absolute error, are obtained with the BHandHLYP functional. Three of the six studied Berkeley functionals and the three range-separated Florida functionals also yield accurate magnetizabilities. Also, some older functionals like CAM-B3LYP, KT1, BHLYP (BHandH), B3LYP, and PBE0 perform rather well. In contrast, unsatisfactory performance is generally obtained with Minnesota functionals, which are therefore not recommended for calculations of magnetically induced current density susceptibilities and related magnetic properties such as magnetizabilities and nuclear magnetic shieldings. We also demonstrate that magnetizabilities can be calculated by numerical integration of magnetizability density; we have implemented this approach as a new feature in the gauge-including magnetically induced current (GIMIC) method. Magnetizabilities can be calculated from magnetically induced current density susceptibilities within this approach even when analytical approaches for magnetizabilities as the second derivative of the energy have not been implemented. The magnetizability density can also be visualized, providing additional information that is not otherwise easily accessible on the spatial origin of magnetizabilities.

Automatic differentiation represents a paradigm shift in scientific programming, where evaluating both functions and their derivatives is required for most applications. By removing the need to explicitly derive expressions for gradients, development times can be shortened and calculations can be simplified. For these reasons, automatic differentiation has fueled the rapid growth of a variety of sophisticated machine learning techniques over the past decade, but is now also increasingly showing its value to support ab initio simulations of quantum systems and enhance computational quantum chemistry. Here, we present an open-source differentiable quantum chemistry simulation code and explore applications facilitated by automatic differentiation: (1) calculating molecular perturbation properties, (2) reoptimizing a basis set for hydrocarbons, (3) checking the stability of self-consistent field wave functions, and (4) predicting molecular properties via alchemical perturbations.

Density fitting (DF), also known as the resolution of the identity (RI), is a widely used technique in quantum chemical calculations with various types of atomic basis sets - Gaussian-type orbitals, Slater-type orbitals, as well as numerical atomic orbitals - to speed up density functional, Hartree-Fock, and post-Hartree-Fock calculations. Traditionally, custom auxiliary basis sets are hand-optimized for each orbital basis set; however, some automatic schemes have also been suggested. In this work, we propose a simple yet numerically stable automated scheme for forming auxiliary basis sets with the help of a pivoted Cholesky decomposition, which is applicable to any type of atomic basis function. We exemplify the scheme with proof-of-concept calculations with Gaussian basis sets and show that the proposed approach leads to negligible DF/RI errors in Hartree-Fock (HF) and second-order M{\o}ller-Plesset (MP2) total energies of the non-multireference part of the W4-17 test set when used with orbital basis sets of at least polarized triple-$\zeta$ quality. The results are promising for future applications employing Slater-type orbitals or numerical atomic orbitals, as well as schemes based on the use of local fitting approaches. Global fitting approaches can also be used, in which case the high-angular-momentum functions produced by the present scheme can be truncated to minimize the computational cost., Comment: 41 pages, 3 figures

Most computational studies in chemistry and materials science are based on the use of density functional theory. Although the exact density functional is unknown, several density functional approximations (DFAs) offer a good balance of affordable computational cost and semi-quantitative accuracy for applications. The development of DFAs still continues on many fronts, and several new DFAs aiming for improved accuracy are published every year. However, the numerical behavior of these DFAs is an often overlooked problem. In this work, we look at all 592 DFAs for three-dimensional systems available in Libxc 5.2.2 and examine the convergence of the density functional total energy based on tabulated atomic Hartree-Fock wave functions. We show that several recent DFAs, including the celebrated SCAN family of functionals, show impractically slow convergence with typically used numerical quadrature schemes, making these functionals unsuitable both for routine applications or high-precision studies, as thousands of radial quadrature points may be required to achieve sub-$\mu E_{h}$ accurate total energies for these unctionals, while standard quadrature grids like the SG-3 grid only contain $\mathcal{O}(100)$ radial quadrature points. These results are both a warning to users to lways check the sufficiency of the quadrature grid when adopting novel functionals, as well as a guideline to the theory community to develop better behaved density functionals., Comment: 16 pages, 6 figures

This article summarizes technical advances contained in the fifth major release of the Q-Chem quantum chemistry program package, covering developments since 2015. A comprehensive library of exchange–correlation functionals, along with a suite of correlated many-body methods, continues to be a hallmark of the Q-Chem software. The many-body methods include novel variants of both coupled-cluster and configuration-interaction approaches along with methods based on the algebraic diagrammatic construction and variational reduced density-matrix methods. Methods highlighted in Q-Chem 5 include a suite of tools for modeling core-level spectroscopy, methods for describing metastable resonances, methods for computing vibronic spectra, the nuclear–electronic orbital method, and several different energy decomposition analysis techniques. High-performance capabilities including multithreaded parallelism and support for calculations on graphics processing units are described. Q-Chem boasts a community of well over 100 active academic developers, and the continuing evolution of the software is supported by an “open teamware” model and an increasingly modular design. This article summarizes technical advances contained in the fifth major release of the Q-Chem quantum chemistry program package, covering developments since 2015. A comprehensive library of exchange-correlation functionals, along with a suite of correlated many-body methods, continues to be a hallmark of the Q-Chem software. The many-body methods include novel variants of both coupled-cluster and configuration-interaction approaches along with methods based on the algebraic diagrammatic construction and variational reduced density-matrix methods. Methods highlighted in Q-Chem 5 include a suite of tools for modeling core-level spectroscopy, methods for describing metastable resonances, methods for computing vibronic spectra, the nuclear-electronic orbital method, and several different energy decomposition analysis techniques. High-performance capabilities including multithreaded parallelism and support for calculations on graphics processing units are described. Q-Chem boasts a community of well over 100 active academic developers, and the continuing evolution of the software is supported by an "open teamware" model and an increasingly modular design.

Long in the making, computational chemistry for the masses [J. Chem. Educ. 1996, 73, 104] is finally here. We point out the existence of a variety of free and open source software (FOSS) packages for computational chemistry that offer a wide range of functionality all the way from approximate semiempirical calculations with tight-binding density functional theory to sophisticated ab initio wave function methods such as coupled-cluster theory, both for molecular and for solid-state systems. By their very definition, FOSS packages allow usage for whatever purpose by anyone, meaning they can also be used in industrial applications without limitation. Also, FOSS software has no limitations to redistribution in source or binary form, allowing their easy distribution and installation by third parties. Many FOSS scientific software packages are available as part of popular Linux distributions, and other package managers such as pip and conda. Combined with the remarkable increase in the power of personal devices—which rival that of the fastest supercomputers in the world of the 1990s—a decentralized model for teaching computational chemistry is now possible, enabling students to perform reasonable modeling on their own computing devices, in the bring your own device (BYOD) scheme. In addition to the programs’ use for various applications, open access to the programs’ source code also enables comprehensive teaching strategies, as actual algorithms’ implementations can be used in teaching. We discuss the availability and use of various FOSS quantum chemistry packages and demonstrate what kinds of calculations are feasible with these programs, assuming only extremely modest computational resources. Our examples confirm that FOSS software enables decentralized approaches to computational chemistry education within the BYOD scheme, affording a democratization of the science of computational chemistry as well.

We have assessed the accuracy for magnetic properties of a set of 51 density functional approximations, including both recently published as well as already established functionals. The accuracy assessment considers a series of 27 small molecules and is based on comparing the predicted magnetizabilities to literature reference values calculated using coupled cluster theory with full singles and doubles and perturbative triples [CCSD(T)] employing large basis sets. The most accurate magnetizabilities, defined as the smallest mean absolute error, were obtained with the BHandHLYP functional. Three of the six studied Berkeley functionals and the three range-separated Florida functionals also yield accurate magnetizabilities. Also some older functionals like CAM-B3LYP, KT1, BHLYP (BHandH), B3LYP and PBE0 perform rather well. In contrast, unsatisfactory performance was generally obtained with Minnesota functionals, which are therefore not recommended for calculations of magnetically induced current density susceptibilities, and related magnetic properties such as magnetizabilities and nuclear magnetic shieldings. We also demonstrate that magnetizabilities can be calculated by numerical integration of the magnetizability density; we have implemented this approach as a new feature in the gauge-including magnetically induced current method (GIMIC). Magnetizabilities can be calculated from magnetically induced current density susceptibilities within this approach even when analytical approaches for magnetizabilities as the second derivative of the energy have not been implemented. The magnetizability density can also be visualized, providing additional information that is not otherwise easily accessible on the spatial origin of the magnetizabilities., 22 pages, 4 figures + 25 pages of supporting information

We develop a methodology for calculating, analyzing and visualizing nuclear magnetic shielding densities, which are calculated from the current density via the Biot-Savart relation. Atomic contributions to nuclear magnetic shielding constants can be estimated within our framework with a Becke partitioning scheme. The new features have been implemented in the GIMIC program and are applied in this work to the study of the $^1$H and $^{13}$C nuclear magnetic shieldings in benzene (C$_6$H$_6$) and cyclobutadiene (C$_4$H$_4$). The new methodology allows a visual inspection of the spatial origins of the positive (shielding) and negative (deshielding) contributions to the nuclear magnetic shielding constant of a single nucleus, something which has not been hitherto easily accomplished. Analysis of the shielding densities shows that diatropic and paratropic current-density fluxes yield both shielding as well as deshielding contributions, as the shielding or deshielding is determined by the direction of the current-density flux with respect to the studied nucleus instead of the tropicity. Becke partitioning of the magnetic shieldings shows that the magnetic shielding contributions mainly originate from the studied atom and its nearest neighbors, confirming the localized character of nuclear magnetic shieldings., 18 pages, 8 figures, includes supporting information

The homogeneous electron gas (HEG) is a key ingredient in the construction of most exchange-correlation functionals of density-functional theory. Often, the energy of the HEG is parameterized as a function of its spin density $n$, leading to the local density approximation (LDA) for inhomogeneous systems. However, the connection between the electron density and kinetic energy density of the HEG can be used to generalize the LDA by evaluating it on a weighted geometric average of the local spin density and the spin density of a HEG that has the local kinetic energy density of the inhomogeneous system, with a mixing ratio $x$. This leads to a new family of functionals that we term meta-local density approximations (meta-LDAs), which are still exact for the HEG, which are derived only from properties of the HEG, and which form a new rung of Jacob's ladder of density functionals. The first functional of this ladder, the local $\tau$ approximation (LTA) of Ernzerhof and Scuseria that corresponds to $x=1$ is unfortunately not stable enough to be used in self-consistent field calculations, because it leads to divergent potentials as we show in this work. However, a geometric averaging of the LDA and LTA densities with smaller values of $x$ not only leads to numerical stability of the resulting functional, but also yields more accurate exchange energies in atomic calculations than the LDA, the LTA, or the tLDA functional ($x=1/4$) of Eich and Hellgren. We choose $x=0.50$ as it gives the best total energy in self-consistent exchange-only calculations for the argon atom. Atomization energy benchmarks confirm that the choice $x=0.50$ also yields improved energetics in combination with correlation functionals in molecules, almost eliminating the well-known overbinding of the LDA and reducing its error by two thirds., Comment: 9 pages, 5 figures including supporting information

Fermi--L\"owdin orbitals (FLO) are a special set of localized orbitals, which have become commonly used in combination with the Perdew--Zunger self-interaction correction (SIC) in the FLO-SIC method. The FLOs are obtained for a set of occupied orbitals by specifying a classical position for each electron. These positions are known as Fermi-orbital descriptors (FODs), and they have a clear relation to chemical bonding. In this study, we show how FLOs and FODs can be used to initialize, interpret and justify SIC solutions in a common chemical picture, both within FLO-SIC and in traditional variational SIC, and to locate distinct local minima in either of these approaches. We demonstrate that FLOs based on Lewis' theory lead to symmetry breaking for benzene -- the electron density is found to break symmetry already at the symmetric molecular structure -- while ones from Linnett's double-quartet theory reproduce symmetric electron densities and molecular geometries. Introducing a benchmark set of 16 planar, cyclic molecules, we show that using Lewis' theory as the starting point can lead to artifactual dipole moments of up to 1 Debye, while Linnett SIC dipole moments are in better agreement with experimental values. We suggest using the dipole moment as a diagnostic of symmetry breaking in SIC and monitoring it in all SIC calculations. We show that Linnett structures can often be seen as superpositions of Lewis structures and propose Linnett structures as a simple way to describe aromatic systems in SIC with reduced symmetry breaking. The role of hovering FODs is also briefly discussed., Comment: 14 pages, 8 figures, includes the SI as an attachment Changes to v1: Added some more theory, adjusted Fig.1 for better readability

X-ray astronomy lacks high resolution spectra of interstellar dust analogues and molecules, severely hampering interstellar medium studies based on upcoming X-ray missions. Various theoretical approaches may be used to address this problem, but they must first be shown to reproduce reliable spectra compared to the experiment. In this work, we calculate the sulfur K edge X-ray absorption spectra of H2S, SO2, and OCS, whose spectra are already known from X-ray experiments and predict the X-ray spectrum of CS, which as far as we are aware has not been measured, thereby hampering its detection by X-ray telescopes. We chose these four molecules as the astrochemistry of sulfur is an unsolved problem and as the four molecules are already known to exist in space. We consider three types of methods for modeling the X-ray spectra: more accurate calculations with the algebraic-diagrammatic construction (ADC) and the CC2, CCSD, and CC3 coupled cluster (CC) approaches as well as more affordable ones with transition potential density functional theory (TP-DFT). A comparison of our computational results to previously reported experimental spectra shows that the core-valence separation (CVS) approaches CVS-ADC(2)-x and CVS-CC3 generally yield a good qualitative level of agreement with the experiment, suggesting that they can be used for interpreting measured spectra, while the TP-DFT method is not reliable for these molecules. However, quantitative agreement with the experiment is still outside the reach of the computational methods studied in this work.

Psi4 is a free and open-source ab initio electronic structure program providing Hartree–Fock, density functional theory, many-body perturbation theory, configuration interaction, density cumulant theory, symmetry-adapted perturbation theory, and coupled-cluster theory. Most of the methods are quite efficient thanks to density fitting and multi-core parallelism. The program is a hybrid of C++ and Python, and calculations may be run with very simple text files or using the Python API, facilitating post-processing and complex workflows; method developers also have access to most of Psi4’s core functionality via Python. Job specification may be passed using The Molecular Sciences Software Institute (MolSSI) QCSchema data format, facilitating interoperability. A rewrite of our top-level computation driver, and concomitant adoption of the MolSSI QCArchive Infrastructure project, make the latest version of Psi4 well suited to distributed computation of large numbers of independent tasks. The project has fostered the development of independent software components that may be reused in other quantum chemistry programs.

Knowledge of the repulsive behavior of potential energy curves $V(R)$ at $R\to0$ is necessary for understanding and modeling irradiation processes of practical interest. $V(R)$ is in principle straightforward to obtain from electronic structure calculations; however, commonly-used numerical approaches for electronic structure calculations break down in the strongly repulsive region due to the closeness of the nuclei. In the present work, we show by comparison to fully numerical reference values that a recently developed procedure [S. Lehtola, J. Chem. Phys. 151, 241102 (2019)] can be employed to enable accurate linear combination of atomic orbitals calculations of $V(R)$ even at small $R$ by a study of the seven nuclear reactions He2 Be, HeNe Mg, Ne2 Ca, HeAr Ca, MgAr Zn, Ar2 Kr, and NeCa Zn., 8 pages

We present fully numerical electronic structure calculations on diatomic molecules exposed to an external magnetic field at the unrestricted Hartree-Fock limit, using a modified version of a recently developed finite element program, HelFEM. We have performed benchmark calculations on a few low-lying states of H2, HeH+, LiH, BeH+, BH, and CH+ as a function of the strength of an external magnetic field parallel to the molecular axis. The employed magnetic fields are in the range of $B=[0,10]~B_0$ atomic units, where $B_0 \approx 2.35 \times 10^5$ T. We have compared the results of the fully numerical calculations to ones obtained with the LONDON code using a large uncontracted gauge-including Cartesian Gaussian (GICG) basis set with exponents adopted from the Dunning aug-cc-pVTZ basis set. By comparison to the fully numerical results, we find that the basis set truncation error in the gauge-including Gaussian basis set is of the order of 1 kcal/mol at zero field, that the truncation error grows rapidly when the strength of the magnetic field increases, and that the largest basis set truncation error at $B=10~B_0$ exceeds 1000 kcal/mol. Studies in larger Gaussian basis sets suggest that reliable results can be obtained in GICG basis sets at fields stronger than $B=B_0$, provided that a sufficient coverage of higher-angular-momentum functions is included in the basis set., Comment: 21 pages and 2 figures, and a further 18 pages and 16 figures of supporting information

We demonstrate that basis sets suitable for electronic structure calculations can be obtained from simple accuracy considerations for the hydrogenic one-electron ions $Y^{(Y-1)+}$ for $Y\in[1,Z]$, necessitating no self-consistent field calculations at all. It is shown that even-tempered basis sets with parameters from the commonly-used universal Gaussian basis set (UGBS) [E. V. R. de Castro and F. E. Jorge, J. Chem. Phys. 108, 5225 (1998)] reproduce non-relativistic spin-restricted spherical Hartree-Fock total energies from fully numerical calculations to better accuracy than UGBS, which is shown to exhibit huge errors for some elements, e.g. 0.19 $E_{h}$ for Th$^+$ and 0.13 $E_{h}$ for Lu as it has been parametrized for a single atomic configuration. Having shown the feasibility of the one-electron approach, partially energy-optimized basis sets are formed for all atoms in the periodic table, $1\leq Z\leq118$, by optimizing the even-tempered parameters for $Z^{(Z-1)+}$. As the hydrogenic Gaussian basis sets suggested in the present work are built strictly from first principles, also polarization shells can be obtained in the same fashion in contrast to previous approaches. The accuracy of the polarized basis sets is demonstrated by calculations on a small set of molecules by comparison to fully numerical reference values, which show that chemical accuracy can be reached even for challenging cases like SF$_6$. The present approach is straightforward to extend to relativistic calculations, and could facilitate studies beyond the established periodic table., 14 pages

The superposition of atomic potentials (SAP) approach has recently been shown to be a simple and efficient way to initialize electronic structure calculations [S. Lehtola, J. Chem. Theory Comput. 15, 1593 (2019)]. Here, we study the differences between effective potentials from fully numerical density functional and optimized effective potential calculations for fixed configurations. We find that the differences are small, overall, and choose exchange-only potentials at the local density approximation level of theory computed on top of Hartree-Fock densities as a good compromise. The differences between potentials arising from different atomic configurations are also found to be small at this level of theory. Furthermore, we discuss the efficient Gaussian-basis implementation of SAP via error function fits to fully numerical atomic radial potentials. The guess obtained from the fitted potentials can be easily implemented in any Gaussian-basis quantum chemistry code in terms of two-electron integrals. Fits covering the whole periodic table from H to Og are reported for non-relativistic as well as fully relativistic four-component calculations that have been carried out with fully numerical approaches., Comment: 12 pages, 8 figures

A uniform derivation of the self-consistent field equations in a finite basis set is presented. Both restricted and unrestricted Hartree&ndash, Fock (HF) theory as well as various density functional approximations are considered. The unitary invariance of the HF and density functional models is discussed, paving the way for the use of localized molecular orbitals. The self-consistent field equations are derived in a non-orthogonal basis set, and their solution is discussed also in the presence of linear dependencies in the basis. It is argued why iterative diagonalization of the Kohn&ndash, Sham&ndash, Fock matrix leads to the minimization of the total energy. Alternative methods for the solution of the self-consistent field equations via direct minimization as well as stability analysis are briefly discussed. Explicit expressions are given for the contributions to the Kohn&ndash, Fock matrix up to meta-GGA functionals. Range-separated hybrids and non-local correlation functionals are summarily reviewed.

The complete active space self-consistent field (CASSCF) method is the principal approach employed for studying strongly correlated systems. However, exact CASSCF can only be performed on small active spaces of ~20 electrons in ~20 orbitals due to exponential growth in the computational cost. We show that employing the Adaptive Sampling Configuration Interaction (ASCI) method as an approximate Full CI solver in the active space allows CASSCF-like calculations within chemical accuracy (, 44 pages, 5 figures

A recently developed finite element approach for fully numerical atomic structure calculations [S. Lehtola, Int. J. Quantum Chem. 119, e25945 (2019)] is extended to the description of atoms with spherically symmetric densities via fractionally occupied orbitals. Specialized versions of Hartree-Fock as well as local density and generalized gradient approximation density functionals are developed, allowing extremely rapid calculations at the basis set limit on the ground and low-lying excited states even for heavy atoms. The implementation of range-separation based on the Yukawa or complementary error function (erfc) kernels is also described, allowing complete basis set benchmarks of modern range-separated hybrid functionals with either integer or fractional occupation numbers. Finally, computation of atomic effective potentials at the local density or generalized gradient approximation levels for the superposition of atomic potentials (SAP) approach [S. Lehtola, J. Chem. Theory Comput. 15, 1593 (2019)] that has been shown to be a simple and efficient way to initialize electronic structure calculations is described. The present numerical approach is shown to afford beyond microhartree accuracy with a small number of numerical basis functions, and to reproduce literature results for the ground states of atoms and their cations for $1 \leq Z \leq 86 $. Our results indicate that the literature values deviate by up to 10 {\mu}Eh from the complete basis set limit. The numerical scheme for the erfc kernel is shown to work by comparison to results from large Gaussian basis set calculations from the literature. Spin-restricted ground states are reported for Hartree-Fock and Hartree-Fock-Slater calculations with fractional occupations for $1 \leq Z \leq 118$., Comment: 20 pages

We present the implementation of a variational finite element solver in the HelFEM program for benchmark calculations on diatomic systems. A basis set of the form $\chi_{nlm}(\mu,\nu,\phi)=B_{n}(\mu)Y_{l}^{m}(\nu,\phi)$ is used, where $(\mu,\nu,\phi)$ are transformed prolate spheroidal coordinates, $B_{n}(\mu)$ are finite element shape functions, and $Y_{l}^{m}$ are spherical harmonics. The basis set allows for an arbitrary level of accuracy in calculations on diatomic molecules, which can be performed at present with either nonrelativistic Hartree--Fock (HF) or density functional (DF) theory. Hundreds of DFs at the local spin-density approximation (LDA), generalized gradient approximation (GGA) and the meta-GGA level can be used through an interface with the Libxc library; meta-GGA and hybrid DFs aren't available in other fully numerical diatomic program packages. Finite electric fields are also supported in HelFEM, enabling access to electric properties. We introduce a powerful tool for adaptively choosing the basis set by using the core Hamiltonian as a proxy for its completeness. HelFEM and the novel basis set procedure are demonstrated by reproducing the restricted open-shell HF limit energies of 68 diatomic molecules from the $1^{\text{st}}$ to the $4^{\text{th}}$ period with excellent agreement with literature values, despite requiring \emph{orders of magnitude} fewer parameters for the wave function. Then, the electric properties of the BH and N2 molecules under finite field are studied, again yielding excellent agreement with previous HF limit values for energies, dipole moments, and dipole polarizabilities, again with much more compact wave functions than what were needed in the literature references. Finally, HF, LDA, GGA, and meta-GGA calculations of the atomization energy of N2 are performed, demonstrating the superb accuracy of the present approach., Comment: 33 pages, 2 figures. Minor changes to text

The need for accurate calculations on atoms and diatomic molecules is motivated by the opportunities and challenges of such studies. The most commonly-used approach for all-electron electronic structure calculations in general - the linear combination of atomic orbitals (LCAO) method - is discussed in combination with Gaussian, Slater a.k.a. exponential, and numerical radial functions. Even though LCAO calculations have major benefits, their shortcomings motivate the need for fully numerical approaches based on, e.g. finite differences, finite elements, or the discrete variable representation, which are also briefly introduced. Applications of fully numerical approaches for general molecules are briefly reviewed, and their challenges are discussed. It is pointed out that the high level of symmetry present in atoms and diatomic molecules can be exploited to fashion more efficient fully numerical approaches for these special cases, after which it is possible to routinely perform all-electron Hartree-Fock and density functional calculations directly at the basis set limit on such systems. Applications of fully numerical approaches to calculations on atoms as well as diatomic molecules are reviewed. Finally, a summary and outlook is given., 77 pages, 5 figures, 697 references. Added discussion on numerical methods

The spurious interaction of an electron with itself-self-interaction error-is one of the biggest problems in modern density-functional theory. Some of its most glaring effects, such as qualitatively incorrect charge separation upon dissociation, can be removed by an approximate self-interaction correction due to Perdew and Zunger (PZ) (Perdew, J.; Zunger, A. Phys. Rev. B 1981, 23, 5048). However, the correction introduces an explicit dependence on the occupied orbital densities, which makes proper minimization of the functional difficult. Previous work (Vydrov et al., J. Chem. Phys. 2006, 124, 094108) has suggested that the application of the PZ correction results in worse atomization energies than those obtained with the uncorrected parent functional. But, it has only recently been found that the correct minimization of the PZ energy functional requires complex-valued orbitals, which have not been used in previous work on atomization energies. Here, we study the effect of the proper use of complex-valued orbitals on the atomization energies of molecules in the W4-11 data set (Karton, A.; Daon, S.; Martin, J. M. Chem. Phys. Lett. 2001, 510, 165). We find that the correction has a tendency to weaken the binding of molecules. The correction using complex-valued orbitals is invariably found to yield better atomization energies than the correction with real-valued orbitals. The correction applied to the PBEsol exchange-correlation functional results in a functional that is more accurate than the (uncorrected) PBE functional.

Although many programs have been published for fully numerical Hartree--Fock (HF) or density functional (DF) calculations on atoms, we are not aware of any that support hybrid DFs, which are popular within the quantum chemistry community due to their better accuracy for many applications, or that can be used to calculate electric properties. Here, we present a variational atomic finite element solver in the HelFEM program suite that overcomes these limitations. A basis set of the type $\chi_{nlm}(r,\theta,\phi)=r^{-1}B_{n}(r)Y_{l}^{m}(\hat{\boldsymbol{r}})$ is used, where $B_{n}(r)$ are finite element shape functions and $Y_{l}^{m}$ are spherical harmonics, which allows for an arbitrary level of accuracy. HelFEM supports nonrelativistic HF and DF including hybrid functionals, which are not available in other commonly available program packages. Hundreds of functionals at the local spin density approximation (LDA), generalized gradient approximation (GGA), as well as the meta-GGA levels of theory are included through an interface with the Libxc library. Electric response properties are achievable via finite field calculations. We introduce an alternative grid that yields faster convergence to the complete basis set than commonly used alternatives. We also show that high-order Lagrange interpolating polynomials yield the best convergence, and that excellent agreement with literature HF limit values for electric properties, such as static dipole polarizabilities, can be achieved with the present approach. Dipole moments and dipole polarizabilities at finite field are reported with the PBE, PBEh, TPSS, and TPSSh functionals. Finally, we show that a recently published Gaussian basis set is able to reproduce absolute HF and DF energies of neutral atoms, cations, as well as anions within a few dozen microhartrees., Comment: 44 pages, 11 figures. Small reorganization of intro

Approximate full configuration interaction (FCI) calculations have recently become tractable for systems of unforeseen size thanks to stochastic and adaptive approximations to the exponentially scaling FCI problem. The result of an FCI calculation is a weighted set of electronic configurations, which can also be expressed in terms of excitations from a reference configuration. The excitation amplitudes contain information on the complexity of the electronic wave function, but this information is contaminated by contributions from disconnected excitations, i.e. those excitations that are just products of independent lower-level excitations. The unwanted contributions can be removed via a cluster decomposition procedure, making it possible to examine the importance of connected excitations in complicated multireference molecules which are outside the reach of conventional algorithms. We present an implementation of the cluster decomposition analysis and apply it to both true FCI wave functions, as well as wave functions generated from the adaptive sampling CI (ASCI) algorithm. The cluster decomposition is useful for interpreting calculations in chemical studies, as a diagnostic for the convergence of various excitation manifolds, as well as as a guidepost for polynomially scaling electronic structure models. Applications are presented for (i) the double dissociation of water, (ii) the carbon dimer, (iii) the {\pi} space of polyacenes, as well as (iv) the chromium dimer. While the cluster amplitudes exhibit rapid decay with increasing rank for the first three systems, even connected octuple excitations still appear important in Cr$_2$, suggesting that spin-restricted single-reference coupled-cluster approaches may not be tractable for some problems in transition metal chemistry., Comment: 15 pages, 5 figures

We describe the implementation of orbital optimization for the models in the perfect pairing hierarchy [Lehtola et al, J. Chem. Phys. 145, 134110 (2016)]. Orbital optimization, which is generally necessary to obtain reliable results, is pursued at perfect pairing (PP) and perfect quadruples (PQ) levels of theory for applications on linear polyacenes, which are believed to exhibit strong correlation in the {\pi} space. While local minima and {\sigma}-{\pi} symmetry breaking solutions were found for PP orbitals, no such problems were encountered for PQ orbitals. The PQ orbitals are used for single-point calculations at PP, PQ and perfect hextuples (PH) levels of theory, both only in the {\pi} subspace, as well as in the full {\sigma}{\pi} valence space. It is numerically demonstrated that the inclusion of single excitations is necessary also when optimized orbitals are used. PH is found to yield good agreement with previously published density matrix renormalization group (DMRG) data in the {\pi} space, capturing over 95% of the correlation energy. Full-valence calculations made possible by our novel, efficient code reveal that strong correlations are weaker when larger bases or active spaces are employed than in previous calculations. The largest full-valence PH calculations presented correspond to a (192e,192o) problem., Comment: 19 pages, 4 figures

The theory for the generation of Wannier functions within the generalized Pipek--Mezey approach [Lehtola, S.; J\'onsson, H. J. Chem. Theory Comput. 2014, 10, 642] is presented and an implementation thereof is described. Results are presented for systems with periodicity in one, two and three dimensions as well as isolated molecules. The generalized Pipek--Mezey Wannier functions (PMWF) are highly localized orbitals consistent with chemical intuition where a distinction is maintained between {\sigma}- and {\pi}-orbitals. The PMWF method is compared with the so-called maximally localized Wannier functions (MLWF) that are frequently used for the analysis of condensed matter calculations. Whereas PMWFs maximize the localization criterion of Pipek and Mezey, MLWFs maximize that of Foster and Boys and have the disadvantage of mixing {\sigma}- and {\pi}-orbitals in many cases. The PMWF orbitals turn out to be as localized as the MLWF orbitals as evidenced by cross-comparison of the values of the PMWF and MLWF objective functions for the two types of orbitals. Our implementation in the atomic simulation environment (ASE) is compatible with various representations of the wave function, including real-space grids, plane waves and linear combinations of atomic orbitals. The projector augmented wave formalism for the representation of atomic core electrons is also supported. Results of calculations with the GPAW software are described here, but our implementation can also use output from other electronic structure software such as ABINIT, NWChem and VASP., Comment: 17 pages, 5 figures

We present the generic, object-oriented C++ implementation of the completeness-optimization approach (Manninen and Vaara, J. Comput. Chem. 2006, 27, 434) in the freely available ERKALE program, and recommend the addition of basis set stability scans to the completeness-optimization procedure. The design of the algorithms is independent of the studied property, the used level of theory, as well as of the role of the optimized basis set: the procedure can be used to form auxiliary basis sets in a similar fashion. This implementation can easily be interfaced with various computer programs for the actual calculation of molecular properties for the optimization, and the calculations can be trivially parallelized. Routines for general and segmented contraction of the generated basis sets are also included. The algorithms are demonstrated for two properties of the argon atom--the total energy and the nuclear magnetic shielding constant--and they will be used in upcoming work for generation of cost-efficient basis sets for various properties.

A variational, self-consistent implementation of the Perdew-Zunger self-interaction correction (PZ-SIC), based on a unified Hamiltonian and complex optimal orbitals, is presented for finite systems and atom-centered basis sets. A simplifying approximation allowing the use of real canonical orbitals is proposed. The algorithm is based on two-step self-consistent field iterations, where the updates of the canonical orbitals and the optimal orbitals are done separately. Calculations of the energy of atoms ranging from H to Ar are presented, using various generalized gradient functionals (PBE, APBE, PBEsol) and a meta-generalized gradient functional (TPSS). While the energy of atoms is poorly described by PBEsol, which is a functional optimized to reproduce properties of solids, the PZ-SIC brings the calculations into good agreement with the best ab initio estimates. The importance of using complex optimal orbitals becomes particularly clear in calculations using the TPSS functional, where the original functional gives good results while the application of PZ-SIC with real orbitals gives highly inaccurate results. With complex optimal orbitals, PZ-SIC slightly improves the accuracy of the TPSS functional. The charge localization problem that plagues Kohn-Sham DFT functionals, including hybrid functionals, is illustrated by calculations on the CH3 + F(-) complex, where even PBEsol with PZ-SIC is found to give estimates of both energy and charge with accuracy comparable to that of coupled cluster calculations.

Semilocal approximations to the density functional for the exchange-correlation energy of a many-electron system necessarily fail for lobed one-electron densities, including not only the familiar stretched densities but also the less familiar but closely related noded ones. The Perdew-Zunger (PZ) self-interaction correction (SIC) to a semilocal approximation makes that approximation exact for all one-electron ground- or excited-state densities and accurate for stretched bonds. When the minimization of the PZ total energy is made over real localized orbitals, the orbital densities can be noded, leading to energy errors in many-electron systems. Minimization over complex localized orbitals yields nodeless orbital densities, which reduce but typically do not eliminate the SIC errors of atomization energies. Other errors of PZ SIC remain, attributable to the loss of the exact constraints and appropriate norms that the semilocal approximations satisfy, suggesting the need for a generalized SIC. These conclusions are supported by calculations for one-electron densities and for many-electron molecules. While PZ SIC raises and improves the energy barriers of standard generalized gradient approximations (GGAs) and meta-GGAs, it reduces and often worsens the atomization energies of molecules. Thus, PZ SIC raises the energy more as the nodality of the valence localized orbitals increases from atoms to molecules to transition states. PZ SIC is applied here, in particular, to the strongly constrained and appropriately normed (SCAN) meta-GGA, for which the correlation part is already self-interaction-free. This property makes SCAN a natural first candidate for a generalized SIC. Published under license by AIP Publishing.

Implentation of seminumerical stability analysis for calculations using the Perdew–Zunger self-interaction correction is described. It is shown that real-valued solutions of the Perdew–Zunger equations for gas phase atoms are unstable with respect to imaginary orbital rotations, confirming that a proper implementation of the correction requires complex-valued orbitals. The orbital density dependence of the self-interaction corrected functional is found to lead to multiple local minima in the case of the acrylic acid, H6, and benzene molecules. In the case of benzene, symmetry breaking that results in incorrect ground state geometry is found to occur, erroneously leading to alternating bond lengths in the molecule.

Novel implementations based on dense tensor storage are presented for the singlet-reference perfect quadruples (PQ) [Parkhill, Lawler, and Head-Gordon, J. Chem. Phys. 130, 084101 (2009)] and perfect hextuples (PH) [Parkhill and Head-Gordon, J. Chem. Phys. 133, 024103 (2010)] models. The methods are obtained as block decompositions of conventional coupled-cluster theory that are exact for four electrons in four orbitals (PQ) and six electrons in six orbitals (PH), but that can also be applied to much larger systems. PQ and PH have storage requirements that scale as the square, and as the cube of the number of active electrons, respectively, and exhibit quartic scaling of the computational effort for large systems. Applications of the new implementations are presented for full-valence calculations on linear polyenes (C n H n+2 ), which highlight the excellent computational scaling of the present implementations that can routinely handle active spaces of hundreds of electrons. The accuracy of the models is studied in the {\pi} space of the polyenes, in hydrogen chains (H 50 ), and in the {\pi} space of polyacene molecules. In all cases, the results compare favorably to density matrix renormalization group values. With the novel implementation of PQ, active spaces of 140 electrons in 140 orbitals can be solved in a matter of minutes on a single core workstation, and the relatively low polynomial scaling means that very large systems are also accessible using parallel computing., Comment: 13 pages, 7 figures

The Pipek-Mezey scheme for generating chemically intuitive, localized molecular orbitals is generalized to incorporate various ways of estimating the atomic charges, instead of the ill-defined Mulliken charges used in the original formulation, or Löwdin charges, which have also been used. Calculations based on Bader, Becke, Voronoi, Hirshfeld, and Stockholder partial charges, as well as intrinsic atomic orbital charges, are applied to orbital localization for a variety of molecules. While the charges obtained with these various estimates differ greatly, the resulting localized orbitals are found to be quite similar and properly separate σ and π orbitals, as well as core and valence orbitals. The calculated results are only weakly dependent on the basis set, unlike those based on Mulliken or Löwdin charges. The effect of varying the penalty exponent on the charge in the objective function was studied briefly and was found to lead to some changes in the localized orbitals when degeneracies are present. The various localization methods have been implemented in ERKALE, an open source program for electronic structure calculations.

Hydration of sulfuric acid plays a key role in new-particle formation in the atmosphere. It has been recently proposed that proton dynamics is crucial in the stabilization of these clusters. One key question is how water molecules mediate proton transfer from sulfuric acid, and hence how the deprotonation state of the acid molecule behaves as a function concentration. We address the proton transfer in aqueous sulfuric acid with O K edge and S L edge core-excitation spectra recorded using inelastic X-ray scattering and with ab initio molecular dynamics simulations in the concentration range of 0-18.0 M. Throughout this range, we quantify the acid-water interaction with atomic resolution. Our simulations show that the number of donated hydrogen bonds per Owater increases from 1.9 to 2.5 when concentration increases from 0 to 18.0 M, in agreement with a rapid disappearance of the pre-edge feature in the O K edge spectrum. The simulations also suggest that for 1.5 M sulfuric acid SO4(2-) is most abundant and that its concentration falls monotonously with increasing concentration. Moreover, the fraction of HSO4(-) peaks at ∼12 M.

Results of Perdew–Zunger self-interaction corrected (PZ-SIC) density functional theory calculations of the atomization energy of 35 molecules are compared to those of high-level quantum chemistry calculations. While the PBE functional, which is commonly used in calculations of condensed matter, is known to predict on average too high atomization energy (overbinding of the molecules), the application of PZ-SIC gives a large overcorrection and leads to significant underestimation of the atomization energy. The exchange enhancement factor that is optimal for the generalized gradient approximation within the Kohn-Sham (KS) approach may not be optimal for the self-interaction corrected functional. The PBEsol functional, where the exchange enhancement factor was optimized for solids, gives poor results for molecules in KS but turns out to work better than PBE in PZ-SIC calculations. The exchange enhancement is weaker in PBEsol and the functional is closer to the local density approximation. Furthermore, the drop inthe exchange enhancement factor for increasing reduced gradient in the PW91 functional gives more accurate results than the plateaued enhancement in the PBE functional. A step towards an optimal exchange enhancement factor for a gradient dependent functional of the PZ-SIC form is taken by constructing an exchange enhancement factor that mimics PBEsol for small values of the reduced gradient, and PW91 for large values. The average atomization energy is then in closer agreement with the high-level quantum chemistry calculations, but the variance is still large, the F2 molecule being a notable outlier.

We present a comprehensive simulation study on the solid-liquid phase transition of the ionic liquid 1,3-dimethylimidazolium chloride in terms of the changes in the atomic structure and their effect on the Compton profile. The structures were obtained by using ab initio molecular dynamics simulations. Chosen radial distribution functions of the liquid structure are presented and found generally to be in good agreement with previous ab initio molecular dynamics and neutron scattering studies. The main contributions to the predicted difference Compton profile are found to arise from intermolecular changes in the phase transition. This prediction can be used for interpreting future experiments.

Completeness-optimization is a novel method for the formation of one-electron basis sets. Contrary to conventional methods of basis set generation that optimize the basis set with respect to ground-state energy, completeness-optimization is a completely general, black-box method that can be used to form cost-effective basis sets for any wanted property at any level of theory. In our recent work [J. Lehtola, P. Manninen, M. Hakala, and K. Hämäläinen, J. Chem. Phys. 137, 104105 (2012)] we applied the completeness-optimization approach to forming primitive basis sets tuned for calculations of the electron momentum density at the Hartree-Fock (HF) level of theory. The current work extends the discussion to contracted basis sets and to the post-HF level of theory. Contractions are found to yield significant reductions in the amount of functions without compromising the accuracy. We suggest polarization-consistent and correlation-consistent basis sets for the first three rows of the periodic table, which are completeness-optimized for electron momentum density calculations.

We report on the microscopic structure of water at sub- and supercritical conditions studied using X-ray Raman spectroscopy, ab initio molecular dynamics simulations, and density functional theory. Systematic changes in the X-ray Raman spectra with increasing pressure and temperature are observed. Throughout the studied thermodynamic range, the experimental spectra can be interpreted with a structural model obtained from the molecular dynamics simulations. A spatial statistical analysis using Ripley’s K-function shows that this model is homogeneous on the nanometer length scale. According to the simulations, distortions of the hydrogen-bond network increase dramatically when temperature and pressure increase to the supercritical regime. In particular, the average number of hydrogen bonds per molecule decreases to ≈0.6 at 600 °C and p = 134 MPa.