The Chronus Quantum (ChronusQ) Software Package

TL;DR

ChronusQ is an open-source, parallel C++ software package designed for advanced ab initio electronic structure calculations, emphasizing time dependence, spin, relativistic effects, and magnetic field inclusion.

Contribution

It introduces a modern, parallel C++ framework with comprehensive features for electronic structure theory, including relativistic and magnetic field treatments, advancing computational chemistry tools.

Findings

Supports time-dependent and spin-dependent electronic structure calculations

Includes relativistic effects and magnetic field perturbations

Utilizes MPI and OpenMP for efficient parallel performance

Abstract

The Chronus Quantum (ChronusQ) software package is an open source (under the GNU General Public License v2) software infrastructure which targets the solution of challenging problems that arise in ab initio electronic structure theory. Special emphasis is placed on the consistent treatment of time dependence and spin in the electronic wave function, as well as the inclusion of relativistic effects in said treatments. In addition, ChronusQ provides support for the inclusion of uniform finite magnetic fields as external perturbations through the use of gauge-including atomic orbitals (GIAO). ChronusQ is a parallel electronic structure code written in modern C++ which utilizes both message passing (MPI) and shared memory (OpenMP) parallelism. In addition to the examination of the current state of code base itself, a discussion regarding ongoing developments and developer contributions will also be provided.

Equations

$$\text{RS / US:}{F}^{\sigma }(P^{\alpha },P^{\beta })C^{\sigma }=S{C}^{\sigma }{E}^{\sigma },\sigma \in \{\alpha ,\beta \}.$$

$$\text{GS:}{F}^{2C}(P^{2C})C^{2C}=S{C}^{2C}{E}^{2C},$$

$$X^{2C}=\left[\begin{array}{cc}{X}^{\alpha \alpha }& X^{\alpha \beta }\\ X^{\beta \alpha }& X^{\beta \beta }\end{array}\right],X\in \{F,P\}.$$

$$X^{2C}=\sum _{I=0}^3{X}^I\otimes {\sigma }_I,$$

$$2X^0$$

$$2X^1$$

$$-2i{X}^2$$

$$2X^3$$

$$F=\{\begin{array}{ll}{F}^0\otimes {\sigma }_0+F^1\otimes {\sigma }_1+F^2\otimes {\sigma }_2+F^3\otimes {\sigma }_3& \mathrm{GS}\\ F^0\otimes {\sigma }_0+F^3\otimes {\sigma }_3& \mathrm{US}\\ F^0\otimes {\sigma }_0& \mathrm{RS}\end{array}$$

$$⟨O^{2C}⟩=\mathrm{Tr}\left[O^{2C}{P}^{2C}\right]=2\sum _{I=0}^3\mathrm{Tr}\left[O^I{P}^I\right]$$

$$F^0=h^0+J(P^0)-K(P^0),$$

$$F^k=h^k-K(P^k),k\in \{1,2,3\},$$

$${\chi }^{\mathrm{CGTO}}(\mathbf{r}-\mathbf{R})\mapsto \exp \left(\frac{i}{2}(\mathbf{R}\times \mathbf{B})\cdot (\mathbf{r}-\mathbf{R})\right){\chi }^{\mathrm{CGTO}}(\mathbf{r}-\mathbf{R}),$$

$$i\hslash {\partial }_{t}P(t)=[F(P(t),t),P(t)].$$

$$P(t)=U(t,t_0)P(t_0)U(t,t_0{)}^{†},U(t,t_0)=\exp \left(\Omega (t,t_0)\right),$$

$$P(t_{k+1})=U_{\mathrm{MMUT}}(t_k)P(t_{k-1})U_{\mathrm{MMUT}}^{†}(t_k),U_{\mathrm{MMUT}}(t_k)=\exp \left(-\frac{2i\Delta t}{\hslash }F(P(t_k),t_k)\right).$$

$$P(t_{k+1})=U_{EM2}(t_k)P(t_k)U_{EM2}^{†}(t_k),$$

$$F_m(t_k,t_{k+1})=\frac{1}{2}(F(P(t_k),t_k)+F(P_{\mathrm{FE}}(t_{k+1}),t_{k+1})),$$

$$P_{\mathrm{FE}}(t_{k+1})=U_{\mathrm{FE}}(t_k)P(t_k)U_{\mathrm{FE}}^{†}(t_k),U_{\mathrm{FE}}(t_k)=\exp \left(-\frac{i\Delta t}{\hslash }F(P(t_k),t_k)\right).$$

$$F(P(t),t)=F_0(P(t))+\sum _{\xi }^{N_{\mathrm{pert}}}{V}_{\xi }(t),$$

$$V_{\xi }(t)=V_{\xi }\cdot f(t)(\Theta (t-t_{\mathrm{on}})-\Theta (t-t_{\mathrm{off}})),$$

$$f(t)=\kappa \cdot \{\begin{array}{ll}1& \text{Step Function}\\ t-t_{\mathrm{on}}& \text{Linear Ramp}\\ \cos (\omega (t-t_{\mathrm{on}}))& \text{Plane Wave}\\ \exp \left(-\alpha (t-t_{\mathrm{on}}{)}^2\right)& \text{Gaussian}\end{array}$$

$$⟨O_{\zeta }(t)⟩=\mathrm{Tr}\left[O_{\zeta }P(t)\right].$$

$$HT=ST\Omega ,$$

$$⟨O_{\zeta }⟩=⟨⟨O_{\zeta };V_{\xi }⟩{⟩}_{\omega }^{\gamma }\approx {\mathbf{O}}_{\zeta }^{†}{\mathbf{X}}_{\xi }(\omega ),$$

$$(H-(\omega +i\gamma )S){\mathbf{X}}_{\xi }(\omega )={\mathbf{V}}_{\xi }.$$

$$I_{\xi }(\omega )=\underset{-\infty }{\overset{\infty }{\int }}dt{e}^{-i\omega t}⟨{\Phi }_0|(\hat{1}+\hat{\Lambda }){\bar{\mu }}_{\xi }{e}^{i{\bar{H}}_{N}t}{\bar{\mu }}_{\xi }|{\Phi }_0⟩,\xi \in \{x,y,z\}.$$

$${\bar{\mu }}_{\xi }=e^{-\hat{T}}{\hat{\mu }}_{\xi }{e}^{\hat{T}}$$

$${\bar{H}}_N=e^{-\hat{T}}\hat{H}{e}^{\hat{T}}-E_{\mathrm{CC}}.$$

$$f(\omega )=\frac{2}{3}\omega \sum _{\xi }\Re \{I_{\xi }(\omega )\}.$$

Full text

The Chronus Quantum (ChronusQ) Software Package

David B. Williams–Young

 , Alessio Petrone , Shichao Sun 22footnotemark: 2 ,

Torin F. Stetina 22footnotemark: 2 , Patrick Lestrange 22footnotemark: 2 , Chad E. Hoyer 22footnotemark: 2 ,

Daniel R. Nascimento , Lauren Koulias 22footnotemark: 2 , Andrew Wildman 22footnotemark: 2 ,

Joseph Kasper 22footnotemark: 2 , Joshua J. Goings , Feizhi Ding ,

A. Eugene DePrince III 55footnotemark: 5 , Edward F. Valeev , Xiaosong Li 22footnotemark: 2 Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA, 94720Department of Chemistry, University of Washington, Seattle, WA, 98195, United States Corresponding Author: dbwy@lbl.govDipartimento di Scienze Chimiche, Università di Napoli ‘Federico II’, Complesso Universitario di M.S. Angelo, via Cintia, I-80126 Napoli, ItalyDepartment of Chemistry and Biochemistry, Florida State University, Tallahassee, FL 32306, United StatesDepartment of Chemistry, Yale University, 225 Prospect Street, New Haven, CT, 06520, United StatesDivision of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, CA, 91125, United StatesDepartment of Chemistry, Virginia Tech, Blacksburg, VA 24061, United States Corresponding Author: xsli@uw.edu

Article Type:

Software Focus

Abstract

The Chronus Quantum (ChronusQ) software package is an open source (under the GNU General Public License v2) software infrastructure which targets the solution of challenging problems that arise in ab initio electronic structure theory. Special emphasis is placed on the consistent treatment of time dependence and spin in the electronic wave function, as well as the inclusion of relativistic effects in said treatments. In addition, ChronusQ provides support for the inclusion of uniform finite magnetic fields as external perturbations through the use of gauge-including atomic orbitals (GIAO). ChronusQ is a parallel electronic structure code written in modern C++ which utilizes both message passing (MPI) and shared memory (OpenMP) parallelism. In addition to the examination of the current state of code base itself, a discussion regarding ongoing developments and developer contributions will also be provided.

INTRODUCTION

The development of the Chronus Quantum (ChronusQ) software package [1] began in the research group of Xiaosong Li at the University of Washington in the beginning of 2014. The primary purpose of ChronusQ at that time was to serve as an outlet for the dissemination of some of the more “unconventional” electronic structure method development taking place in the group. Chief among these developments was methodological experimentation in the fields of time-dependent and relativistic electronic structure theory, and specifically at their intersection. Over the years, ChronusQ has received contributions from many developers, and new and exciting functionality continues to be added to its code base to this day. ChronusQ is free software released under the GNU General Public License (v2), and is made publicly available at http://www.chronusquantum.org. In this context, free software is meant to indicate that ChronusQ is gratis, as well as the fact that it grants its recipients the freedom to use, modify and redistribute the code base under the same license.

The set of electronic structure methods implemented in the ChronusQ program spans both the novel, and the tried and tested. Contemporary software which implement similar pieces of quantum chemistry functionality include the ReSpect[2], DIRAC[3], DALTON[4], and BAGEL[5] programs. However there are a number of electronic structure methods implemented in ChronusQ which are, at the time that this article has been written and to the best of the authors’ knowledge, not currently available in any other public, open source quantum chemistry codes. These methods will be highlighted in the main text.

As the name might imply (Chronus being the latinization of Time personified in Greek mythology), the core of the functionality in ChronusQ centers around time-dependent electronic structure theory. However, it is through the lens of time-dependent phenomena that ChronusQ also provides a slew of related capabilities such the treatment of electronically excited states and molecular properties. ChronusQ provides functionality to solve time-dependent electronic structure problems both in the time and frequency domains. In conjunction with the time-dependent electronic structure methods, special emphasis is placed on the treatment of electronic spin throughout the code base, specifically in the context of broken spin–symmetry solutions. As such, all of the electronic structure methods implemented in ChronusQ fully support the use of two-component spinor wave functions. In addition to these core functionalities, ChronusQ supports a number of complementary features which are not typically found in general purpose quantum chemistry software. Relativistic effects are supported throughout the code through the use of the exact two–component (X2C) method [6, 7, 8, 9, 10, 11]. Magnetic fields may also be added as an external perturbation through the use of gauge including atomic orbitals (GIAO) [12, 13, 14, 15, 16, 17, 18]. All electronic structure methods implemented in ChronusQ are based on the linear combination of atomic orbitals (LCAO) paradigm to construct the electronic wave function using contracted Gaussian-type orbital (CGTO) basis set expansions. Both generally contracted and segmented CGTO basis sets are supported. ChronusQ is a single-reference quantum chemistry code which utilizes single Slater determinant descriptions of the ground state electronic wave function based on either Hartree-Fock (HF) or Kohn-Sham density functional theory (KS). Using the Slater determinant as a reference wave function, ChronusQ also supports various post-SCF methods such as the particle hole [19, 20, 21] and particle particle [22, 23, 24] random phase approximations. Additionally, the development of many body methods, such as coupled cluster theory [25], are currently underway.

Since its inception, the primary design goal of ChronusQ has been to enable rapid development of novel research in electronic structure theory without sacrificing performance on modern parallel computing architectures. With the exception of a few external dependencies, ChronusQ is a standalone software infrastructure completely written in modern C++. The current release of ChronusQ relies on the C++14 standard, while a number of the development streams utilize features outlined in C++17. To this end, extensive use of modern C++ paradigms (SFINAE, polymorphism, etc.) and template metaprogramming are leveraged throughout the code base to enable a level of abstraction which allows developers to write expressive and efficient code. For an examination of some of the more salient aspects of the C++ design of ChronusQ, we refer the reader to the Appendix. ChronusQ utilizes a hybrid parallelism scheme to improve scalability on large computing clusters which leverages both message passing and shared memory protocols to achieve inter- and intranode concurrency, respectively. Calculations involving thousands of CPU cores and thousands of basis functions are regularly performed with ChronusQ on some of the world’s largest supercomputers [26, 27]. ChronusQ is compatible with any MPI-2 compliant message passing implementation, such as MPICH, OpenMPI and Intel(R) MPI. To allow for the use of ChronusQ on work stations which may not possess an MPI implementation, message passing parallelism is an optional feature which may be enabled or disabled by the build system. Shared memory parallelism is achieved through a combination of the explicit use of OpenMP threads and reliance on threaded external libraries such as those which implement the BLAS [28, 29, 30, 31]  and LAPACK [32] standards for numerical linear algebra.

After several years of development, ChronusQ has become a full-fledged, standalone, ab initio quantum chemistry software package which has proven itself useful both as a development platform to rapidly prototype electronic structure methods, as well as a robust production code capable of tackling complex science problems [26, 18, 33, 34, 35]. In this article, we briefly review the design, features and current developments of the ChronusQ software package.

DESIGN AND NOVEL STRUCTURAL FEATURES

Input Specification

ChronusQ utilizes the INI format for the specification of its input. An example ChronusQ input file specification for the time propagation of a single water molecule at the cc-pVTZ / X2C-B3LYP [36, 37] level of theory in the presence of an electric dipole field may be found in Fig. 1. ChronusQ packages a representative subset of commonly used CGTO basis sets, including the Sapporo [38, 39] basis sets commonly used in relativistic calculations, with its source code. Arbitrary CGTO basis sets may be utilized by specifying the full path for the desired basis set file in the Gaussian94 file format. For a detailed examination of the anatomy and options supported by the ChronusQ input file, we refer the reader to the ChronusQ user manual which may be found at http://www.chronusquantum.org.

Software Dependencies

While ChronusQ is best described as a standalone quantum chemistry software suite, this description is primarily limited to the electronic structure capabilities of the software. ChronusQ relies on a number of external libraries to provide complementary functionality which facilitate efficient implementations of quantum mechanical methods. Where possible and practical, these software dependencies are packaged with the source code itself; however a number of these dependencies must be installed separately and discovered by the build system upon configuration. A list and description of these dependencies is outlined in the ChronusQ user manual.

To ensure portability on many different software configurations across different architectures, ChronusQ utilizes a CMake build system. The CMake build system attempts to determine the existence and linkage of the various dependencies as well as provide specifications of various build options which influence the performance and configuration of ChronusQ such as debug and optimization flags, enabling MPI bindings, etc. A comprehensive guide to the CMake build system may be found in the ChronusQ user manual.

ChronusQ relies extensively on optimized and threaded implementations of BLAS and LAPACK numerical linear algebra standards to ensure high-levels of performance on modern computing architectures. The open source OpenBLAS [40, 41] library is packaged with ChronusQ to ensure that such performance is portable to a large number of different architectures. ChronusQ is also compatible with vendor tuned linear algebra software such as the Intel(R) Math Kernel Library (MKL) and IBM Engineering and Scientific Software Library (ESSL) for specific architectures. If MPI has been enabled, ChronusQ may be optionally configured to utilize the ScaLAPACK [42] library to perform high-performance distributed memory linear algebra suitable for large computing clusters.

With the exception of those integrals which are specific to the X2C and GIAO methods implemented in ChronusQ, molecular integrals over contracted Gaussian basis functions are evaluated using the Libint library [43]. ChronusQ was one of the early adopters of the modern C++ interface of the Libint library which has enabled the development of expressive and reusable integral drivers. For the evaluation of exchange-correlation (XC) functionals (and their derivatives) in the KS implementation, ChronusQ utilizes the Libxc library [44, 45]. To ensure portability of standardization of the ChronusQ checkpoint and binary files, we utilize the HDF5 library.

General Infrastructure for the Unified Treatment of Electronic Spin

Perhaps the most common motif in the design of ChronusQ is in the general manner with which it treats electronic spin throughout the code base. It is in this context that the development of ChronusQ has benefited from the era in which it has been written: spin was considered a priority from the beginning. Traditionally, quantum chemistry codes are generally written from a “bottom-up” approach as it relates to the breaking of spin symmetries: spin restricted methods (RS: $S^{2}$ and $S_{z}$ eigenfunctions) are written before spin unrestricted (US: $S_{z}$ eigenfunctions) methods which are written before generalized spin methods (GS: broken spin-symmetry), if such GS treatments are implemented at all. If GS is to be neglected in a particular implementation, this “bottom-up” approach is logical as the Roothaan-Hall equations separate into distinct (but coupled) eigenvalue problems in the binary spin basis [46],

[TABLE]

Here $\alpha$ and $\beta$ are spin indices, $S$ is the basis overlap matrix, and $F^{\sigma}$, $C^{\sigma}$, $P^{\sigma}$, and $E^{\sigma}$ are the spin dependent Fock matrix, molecular orbital matrix, density matrix, and diagonal matrix of canonical orbital eigenenergies, respectively. Round brackets, $(\cdot)$, indicate a functional dependence. However, GS takes the form of a single eigenvalue problem of twice the dimension [47],

[TABLE]

where

[TABLE]

Generalizing a code base from RS to US requires the manipulation of one additional density/Fock matrix whereas generalization of US to GS requires the manipulation of one fewer density/Fock matrix. GS almost seems artificial in this treatment as it more closely resembles RS in the context of the manipulation of individual matrices but resembles US in the treatment of multiple spin degrees of freedom. Such a structure requires a significant amount of logic in the code base to control the proper workflow and specialization of many aspects of the code to either work with RS/US or GS. In ChronusQ, the treatment of spin was a primary consideration from the initial development, thus some time was spent in developing a unified code infrastructure that treats RS, US, and GS on an equal footing. To this end, we adopt a “top-down” approach to spin manipulation such that GS is a generalization of US which is a generalization of RS, as would be implied by the successive restriction of the wave function to be an eigenfunction of $S_{z}$ and $S^{2}$. Such a structure is possible through a change of spin basis to the Pauli spin basis [48, 27],

[TABLE]

where

[TABLE]

and $\{\sigma_{I}\}$ is the set of Pauli matrices with $\sigma_{0}$ representing the 2x2 identity matrix. $X^{0}$ is referred to as the scalar component of the spinor operator $X^{2C}$, while $X^{1}$, $X^{2}$ and $X^{3}$ are referred to as its spin components. In this basis, specialization from GS to US to RS simply implies the removal of certain spin components of density and Fock matrices. Thus, in ChronusQ, spin-dependent operators are stored internally as collections of scalar and spin components according to

[TABLE]

The use of this scheme for spin separation is primarily a mechanism to enable a readable and expressive code structure void of fringe logic to handle the manipulation of the different spin cases. Internally, the eigenvalue problems solved for the Roothaan-Hall equations are exactly those given in Eqs. 1 and 2. In the case of the molecular orbital matrix and canonical orbital eigenvalues, no such simplification can be made. However, due to the fact that ChronusQ is primarily a single-reference quantum chemistry code, the majority of the functionality may be implemented through the use of the density matrix as opposed to the orbitals directly. Thus the vast majority of the code base benefits from this spin separation scheme. The manipulation of the cases for the eigenvalue problem is facilitated through templated matrix utilities which efficiently convert between the two formats, which yields only minimal effect on the overall structure of the code base. In practice, we have found the following points to be the primary motivation for the chosen approach for the treatment of spin in ChronusQ:

Generalization of a method to support two-component spinors is systematic. If a method may be expressed in the Pauli spin basis, generalization is typically as simple as performing some set of mathematical operations on the scalar component followed by another set of (possibly different) operations on the spin components (e.g. Eq. 8) 2. 2.

The off-diagonal matrices $X^{\alpha\beta}$ and $X^{\beta\alpha}$ are generally non-Hermitian, while $X^{2}$ and $X^{3}$ are Hermitian as long as $X^{2C}$ is Hermitian. Thus when constructing moieties such as $K(X^{i})$ (see Eq. 8), the symmetry of $X^{i}$ may be utilized; reducing the number of floating point operations in direct contraction schemes. 3. 3.

The manipulation of general (possibly spin dependent) perturbations as they relate to the Hamiltonian and the evaluation of molecular properties is straight forward. All perturbations which are spin independent (such as electric multipoles, electron-nuclear interaction, etc.) only appear in the scalar Fock matrix, while those which are spin dependent (such as magnetic fields) only appear in the spin components of the Fock matrix. Further, due to the fact that the Pauli matrices are traceless, the expectation value of a general one particle spinor operator, $O^{\mathrm{2C}}$, with respect to the one particle spinor density matrix, $P^{\mathrm{2C}}$, may be expressed as [48, 27]

[TABLE]

CURRENT PROGRAM FEATURES

Ground State Methods

Reference Wave Functions

ChronusQ implements a general self–consistent field (SCF) module which provides a unified treatment of RS, US, and GS spin symmetries to obtain ground state wave functions described by either Hartree–Fock (HF) or Kohn–Sham density functional theory (KS) approximations. Within the general framework for the treatment of spin described in the previous section, the HF Fock matrix may be expressed as [48, 27]

[TABLE]

where $h$ is the core Hamiltonian, and $J$ and $K$ are the Coulomb and exchange matrices, respectively. Thus the treatment of RHF, UHF and GHF may be expressed simply by forming additional exchange terms for the required spin components of the density matrix. The formation of the Fock matrix is based on the direct SCF scheme of Almlöf, et al. [49] using a integral screening scheme resembling that of the LinK method [50].

The formation of the KS Fock matrix in ChronusQ adopts a form analogous to Eqs. 6 and 8 for the RKS, UKS and GKS methods [48, 27]. For the use of two-component spinor density, ChronusQ utilizes a set of torque free auxiliary variables which allow the adaptation of collinear XC functionals for use with non-collinear densities [27, 51]. This non-collinear DFT formalism preserves the self-consistent, global zero-torque theorem for the XC magnetic field [52] while admitting the ability to yield non-zero local torque. In practice, it is the local torque which allows for the accurate description of spin dynamics and time-evolution of the spin density [52, 53, 54, 55, 56, 57]. The XC integration and batching scheme adopted in ChronusQ closely resembles that of Stratmann, et al. [58]. The assembly of the XC potential in ChronusQ generalizes the scheme of Burow, et al. [59] for use with open-shell and spinor densities [27]. ChronusQ supports local density (LDA), generalized gradient (GGA), and hybrid approximations for the XC functional. The formation of both the HF and KS Fock matrix is efficiently parallelized using both inter- and intranode communication schemes. Ref 27 demonstrates the linear scaling (strong and weak) of the KS implementation in ChronusQ.

Relativistic Effects and External Magnetic Fields

Relativistic effects may be included variationally in the SCF procedure for both HF and KS wave functions through the exact two-component (X2C) method [6, 7, 8, 9, 10, 11]. Practically, the extension of the SCF module to implement the X2C-HF and X2C-KS methods simply amounts to the manipulation of additional spin components for the core Hamiltonian in Eq. 8. To ensure numerical stability, relativistic calculations are performed using Gaussian charge models [60] for nuclear charge distributions rather than the point charge distributions typically used in non-relativistic theory. For calculations including one-electron spin orbit contributions to the core Hamiltonian, Boettger scaling [61] is used to approximately account for the higher-order two-electron terms. ChronusQ also supports SCF optimization in the presence of finite magnetic fields through the use of gauge including atomic orbitals (GIAO) [12, 13] at the GIAO-RHF [15], GIAO-UHF, and GIAO-GHF levels of theory [17, 18]. For the formation of the GIAO Fock matrix, the magnetic field manifests both as explicit terms added to the spin component of the Fock matrix as well as an implicit dependence of the CGTO basis function through the addition of a complex phase factor,

[TABLE]

where $\chi^{\mathrm{CGTO}}$ is a CGTO centered at $\mathbf{R}$ and $\mathbf{B}$ is the external magnetic field. In addition, ChronusQ includes orbital–Zeeman, spin–Zeeman, and diamagnetic field interactions of the non-relativistic Pauli Hamiltonian in the GIAO methods [14, 15, 17, 18]. GIAO Fock matrices may be formed directly in the AO basis through a generalization of the direct SCF method of Almlöf, et al. [49] for complex integrals. To the best of the authors’ knowledge, ChronusQ contains the only publicly available implementation of two-component spinor GIAO methods [18] at this time.

Self-Consistent Field Optimization

ChronusQ implements both first and second order optimization schemes for the ground state SCF wave function. First order SCF optimization utilizes the widely adopted [62, 63, 64, 65] scheme involving a combination of damping and direct inversion in the iterative subspace (DIIS) extrapolation of the Fock and density matrices to accelerate convergence. Second order optimization may proceed by either the quadratically convergent SCF method of Bacskay [66], the stability check and perturbation method [67, 68], or some combination of the two. Typically, the best results for the optimization of GSCF wave functions in ChronusQ involves first order optimization of the wave function to a stationary saddle point followed by second order optimization to the minimum energy solution. For a full description of the SCF control options supported by ChronusQ, we refer the reader to the ChronusQ user manual.

Real-Time Electronic Dynamics

One of the hallmarks of the functionality contained in ChronusQ is the implementation of real-time propagation for quantum systems in the presence of external perturbations, as well as its general implementation of such functionality with respect to the treatment of electronic spin. Real-time propagation is important in the study of matter-field interaction due to the fact that it is non-perturbative relative to the strength of the external perturbation; both linear and non-linear effects may be captured through proper analysis [69, 70, 71, 72, 73]. This class of real-time functionality is encapsulated in the RT module of ChronusQ. For a complete description of the specification of options for the RT module in ChronusQ, we refer the reader to the ChronusQ user manual.

Fundamentally, the RT module solves the non-linear Louiville-von Neumann equation [74, 75],

[TABLE]

Formally, Eq. 10 may be solved (propagated) exactly in the time domain through the Magnus expansion of the time-domain propagator [76],

[TABLE]

where $\Omega(t,t_{0})$ is a non-terminating series expansion which must be truncated in practice. Many explicit and implicit time integration methods may be derived from the Magnus expansion of the propagator [77, 78, 72] of which ChronusQ implements a functional subset. The default integration scheme in ChronusQ is the modified midpoint unitary transformation (MMUT) method of Li, et al. [79]. MMUT is a symplectic multi-step (leap-frog), explicit integration scheme based on the Magnus expansion with error formally $\mathcal{O}(\Delta t^{2})$,

[TABLE]

The leap-frog character of MMUT exhibits an excellent trade off between accuracy and efficiency over other Magnus-based expansions for the oscillatory solutions which are characteristic of these types of simulations. Because MMUT is a leap-frog method, a non-leap-frog step must be used to seed and restart the integration every so often to maintain accuracy and stave off energy drift [79]. For this purpose, the default seed and restart step in ChronusQ is the explicit second order (trapezoidal) Magnus (EM2) propagator given by [78]

[TABLE]

where

[TABLE]

For both the MMUT and EM2 integration methods, the matrix exponential may be formed using either exact diagonalization based methods, or approximate Taylor series and Chebyshev series expansions [80, 77, 81]. The latter approximate approaches yield drastic performance improvements for large quantum systems/basis sets due to the poor scaling and parallelization of eigenvalue decomposition.

ChronusQ admits the ability to express an arbitrary number of (one-body) time-dependent external perturbations, $\{V_{\xi}(t)\}$, of the form

[TABLE]

where $F_{0}$ is only implicitly time-dependent through the time dependence of the density matrix. Generally, a particular perturbation takes the functional form

[TABLE]

where $\Theta$ is the Heaviside step function, $t_{\mathrm{on/off}}$ are reference times to turn the perturbation on and off, $V_{\xi}$ is a physical, time-independent one-body operator, and $f(t)$ is an envelope function which describes the time dependence and amplitude of the perturbation. For $V_{\xi}$, ChronusQ supports both the length and velocity gauge electric multipole operators up through the electric octupole, and the magnetic dipole operator. Additionally, ChronusQ supports the following envelope functions:

[TABLE]

where $\kappa$ is a chosen parameter to control the amplitude of the perturbation. An example field specification is given in Fig. 1.

Throughout the time-propagation, ChronusQ evaluates the time-dependent expectation values of various one body operators, $\{O_{\zeta}\}$, according to

[TABLE]

Depending on the spin characteristics of $O_{\zeta}$, Eq. 18 may be efficiently evaluated according to Eq. 7. For example, evaluation of a time-dependent component (i.e. the $z$ component) of the electric dipole amounts to only evaluating the trace of the $z$-dipole matrix with respect to the scalar density matrix, while evaluating a component of the ($z$) component of the magnetic dipole requires evaluating the trace of the $z$-magnetic dipole with respect to $P^{3}(t)$. The supported choices for $O_{\zeta}$ are the same as those for $V_{\xi}$ outlined above. The time-series for the desired time-dependent expectation values may be obtained from the ChronusQ checkpoint file after the calculation for post processing at the users’ discretion. Using these time-series, users may simulate quantities such as absorption spectra [72, 82], electric and magnetic circular dichroism spectra [83, 35], charge transfer [84, 85, 86], and second harmonic generation [71, 87].

Frequency Domain Methods

In addition to time-domain methods, ChronusQ implements a suite of frequency domain methods to solve time-dependent problems in the RESP module. For a full description of the RESP module and its input specification, we refer the reader to the ChronusQ user manual. As opposed to explicit, non-perturbative solution of Eq. 10 in the time domain, these methods solve Eq. 10 in the frequency domain under a series of approximations, namely that the applied perturbation is weak in field strength and harmonic in the time domain. Under these approximations, solution of Eq. 10 (or more specifically Eq. 18) instead becomes a large linear algebra problem to obtain the quantum propagator in the frequency domain. Frequency domain methods fall into one of two classes: residue methods and response methods. In this section, we review the salient aspects of the implementation of these methods in ChronusQ.

Residue Methods

The quantum propagator diverges in the frequency domain at frequencies which correspond to transition energies of the quantum system [19]. Residue methods aim to determine the poles and residues of the propagator at these frequencies, thus obtaining information pertaining to the excited states of the quantum system such as transition energies and densities. For the single determinant reference wave functions implemented in ChronusQ, residue methods generally take the form of a large generalized eigenvalue problem of a matrix pencil $(H,S)$,

[TABLE]

where the form of $H$ and $S$ are dependent on the representation of the reference wave function (i.e. in the case of HF/KS response, $H$ represents the orbital Hessian and $S$ is an indefinite metric [19, 20, 21]). In general, the columns of the eigenvectors $T$ represent the transition density pertaining to the transition energies of the corresponding diagonal element of $\Omega$. ChronusQ supports residue calculations using two different approximations of the quantum propagator, namely particle-hole (ph-RPA) [19, 20, 21] and particle-particle (pp-RPA) [22, 23, 24] random phase approximations. Both HF or KS reference wave functions may be used for either of these propagator approximations. Determining propagator poles and residues under the ph-RPA yields information relating to neutral (particle number conserving) excitations of the quantum system in question, while those corresponding to the pp-RPA yield information pertaining to two electron addition (double electron affinity) and removal (double ionization). Additionally, ChronusQ supports inclusion of both relativistic effects (the X2C-ph-RPA [88] and X2C-pp-RPA [89] methods) and magnetic fields (the GIAO-ph-RPA method [33, 90]) in its residue calculations. The GIAO-ph-RPA method using GIAO-RHF and GIAO-GHF reference wave functions have enabled the calculation of linear response approximations for the magnetic circular dichroism spectra of closed shell systems [33]. To the best of the authors’ knowledge, ChronusQ contains the only publicly available implementation of the GIAO-ph-RPA method at this time. We note for clarity that the GIAO-ph-RPA method, as implemented in ChronusQ, treats the linear response of molecular systems to external perturbations in the presence of explicit finite, static magnetic fields. We are careful to make this distinction in this work to avoid confusion with work on perturbative GIAO response methods such as those developed by Kjærgaard, et al [91, 92].

Generally, the dimension of $H$ is too large to directly diagonalize using methods such as those provided in (Sca)LAPACK. For these large problems, ChronusQ implements the iterative generalized preconditioned locally harmonic residual (GPLHR) method [93, 94] to perform the diagonalization. Rather than obtaining all possible eigenvalues of $(H,S)$, the GPLHR method performs a partial diagonalization for a select number of desired eigenvalues. Unlike the more typically utilized Davidson (Lanczos) method [95, 96, 97], GPLHR is able to directly (and robustly) solve for eigenpairs in the spectral interior of $(H,S)$, thus making it an ideal solver for high-energy spectroscopies such as X-Ray absorption [98]. Like the Davidson method, GPLHR only requires the specification of action (product) of $(H,S)$ onto a trial set of vectors to perform the diagonalization. ChronusQ implements direct, density based contraction schemes to perform the matrix products for both the ph-RPA [99] and pp-RPA [24] methods. Because transition densities exhibit the same spin factorization as the density matrix (Eq. 4) [88], the direct contractions closely resemble those of Eq. 8 in the AO basis, only with a non-Hermitian density operator. As such, the generalization of the residue calculations to use RS, US and GS reference wave functions closely resembles the structure of the Fock matrix formation.

Linear Response Methods

Evaluation of the residues may be avoided entirely if the moiety of interest is Eq. 18 in the frequency domain. Within the linear response (weak perturbation) approximation, frequency dependent properties may be evaluated directly through the solution of a shifted system of linear equations given by [100, 101, 102]

[TABLE]

Here, $O_{\zeta}$ and $V_{\xi}$ are defined as they were in Eqs. 18 and 16, respectively, and $(H,S)$ is the same pencil as Eq. 19. $\gamma$ is a complex damping parameter which ensures convergence of the propagator at all frequencies (including poles) [101]. The boldface $\mathbf{O}_{\zeta}$ and $\mathbf{V}_{\xi}$ are property gradients defined by the selection of the propagator approximation and operators $O_{\zeta}$ and $V_{\xi}$, respectively [19]. ChronusQ currently supports linear response calculations using the ph-RPA only. Much like the residue methods, the dimension of $H$ is often too large to factorize or invert directly using (Sca)LAPACK. As such, ChronusQ provides an implementation of the generalized minimum residual (GMRES) method [103] to iteratively solve large linear systems of the form of Eq. 20b.

In addition to providing solvers which solve Eq. 20b in the full dimension of $H$, ChronusQ implements a series of model order reduction (MOR) techniques [26, 104] to improve the efficiency of evaluating Eq. 20b for multiple frequencies in a spectral region of interest. The MOR technique constructs a rational Krylov subspace which interpolates Eq. 20a in the frequency domain. The subspace is constructed as the linear span of solutions to Eq. 20b at multiple, carefully chosen interpolation frequencies. Using Galerkin projection under the constructed subspace, the full dimensional problem is projected onto a reduced space problem which is rapidly evaluatable. The MOR method implemented in ChronusQ has been demonstrated to provide excellent approximations to the X-Ray absorption spectrum for large molecular systems with quadratic computational scaling. ChronusQ was the first software to provide such MOR functionality and remains the sole publicly available implementation at this time to the best of the authors’ knowledge.

CURRENT DEVELOPMENTS

As is the case with many research codes, there are a number of development streams which have not yet made it into the primary release of ChronusQ. In this section, we examine some of the most recent development in the ChronusQ ecosystem which will be folded into forthcoming releases.

Kohn–Sham Density Functional Theory in Magnetic Fields

Magnetic density functional theory (BDFT) strikes a balance between cost and accuracy for simulating molecules in magnetic fields [105]. Currently, the released version of ChronusQ only supports GIAO basis functions for the treatment of finite magnetic fields within the GIAO-HF level of theory. GIAO-RKS reference wave functions [106, 107, 108] are currently supported in the development version of ChronusQ, and will be present in the forthcoming release.

Quantum Embedding Methods

In recent years, quantum embedding capabilities have been developed in ChronusQ [34]. Currently, users may perform an embedded SCF calculation by specifying both an embedding potential and embedding magnetic field on a real space grid using the XSF file format [109]. Such moieties have regularly been obtained using a combination of the Abinit [110] software and the software presented in previous work by Huang, et al. on embedding [111]. CGTO basis representations of the embedding perturbations are produced by numerical integration of the interpolated real space representations using a similar scheme as the one used for the XC integration for KS methods in ChronusQ (for details, see [34]). Embedded SCF calculations are supported for RSCF, USCF, and GSCF optimizations for both HF and KS levels of theory. Embedding X2C-HF/KS wave functions is currently under development.

Two-Component Post-SCF Methods

The development version of ChronusQ includes working implementations of both ground- and excited-state post-SCF methods. Among the ground-state post-SCF methods, standard second-order Møller-Plesset perturbation theory (MP2) [112], approximate second-order coupled-cluster (CC2) [113], and coupled-cluster with single and double excitations (CCSD) [25] have been implemented. For post-SCF excited states, ChronusQ includes an implementation of the time-dependent (TD) equation-of-motion coupled-cluster (EOM-CC) [114, 115, 116] formalism of Nascimento and DePrince [117, 118], in both the TD-EOM-CC2 and TD-EOM-CCSD flavors.

Unlike other reported time-dependent coupled-cluster implementations [119, 120, 121, 122, 123], the TD-EOM-CC formalism in ChronusQ is based on the time evolution of a coupled-cluster dipole moment function from which one can obtain a spectral line shape function, $I_{\xi}(\omega)$, defined as

[TABLE]

Here, $\bar{\mu}_{\xi}$ represents the $\xi$-component of the similarity-transformed electric dipole operator

[TABLE]

and $\bar{H}_{N}$ represents the normal-ordered similarity-transformed Hamiltonian

[TABLE]

The symbols $E_{\rm CC}$, $\hat{T}$, and $\hat{\Lambda}$ represent the energy, excitation operator, and de-excitation operator associated with the ground-state CC wavefunction, respectively. The real part of the line shape function defines oscillator strengths

[TABLE]

which may be evaluated over an arbitrarily broad energy window [117, 118].

A notable feature of the ground- and excited-state implementations discussed here is that they rely on the same template metaprogramming techniques employed throughout the rest of ChronusQ (see Appendix) and thus are compatible with both real and complex RHF, UHF, and GHF (including X2C-HF) reference wavefunctions. Hence, like many other reported implementations of time-independent relativistic EOM-CC theory [124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138], the X2C-based TD-EOM-CC2 and TD-EOM-CCSD algorithms in ChronusQ incorporate spin-orbit and scalar relativistic effects from first principles. To the best of the authors’ knowledge, this implementation is the only one that incorporates spin-orbit coupling within a time-dependent EOM-CC formalism at this time. Active effort is also being undertaken toward the development of standard (frequency domain) EOM-CC methods [116] within ChronusQ.

CONCLUSIONS

ChronusQ is a standalone, high-performance quantum chemistry software package specializing in time-dependent and relativistic electronic structure methods. Special emphasis is placed on the general treatment of electronic spin throughout the code base, allowing for the implementation of electronic structure methods with various spin symmetries while requiring only minor code modifications. Further, the modern C++ code design adopted by ChronusQ allows for the expression of complex quantum chemistry algorithms in a readable, generic, and elegant manner. Thus, as a development code, ChronusQ has proven to be a near optimal platform for the rapid development of novel electronic structure methods. In addition, with extensive exploitation of high performance numerical linear algebra software, and use of both shared and message passing parallelism schemes, ChronusQ is a performant, production quality quantum chemistry code which is suitable for both workstations and supercomputers alike.

In this article, we reviewed much of the salient aspects of the functionality implemented in the ChronusQ program as well as a number of the current developments which will accompany its upcoming release. The core of the functionality in ChronusQ revolves around the solution of time-dependent electronic structure problems in either the time or frequency domain. ChronusQ supports the inclusion of relativistic effects in the electronic wave function through the X2C method throughout the code. This support in enabled in part by the general manner with which ChronusQ treats electronic spin, as outlined in the main text. Further, ChronusQ supports the inclusion of finite, uniform magnetic fields as external perturbations in all HF based electronic structure methods. Implementations of methods such as GIAO-GHF and the GIAO-ph-RPA in ChronusQ are, at the time this article has been written and to the best of the authors’ knowledge, the sole implementations of these methods in publicly available, open source quantum chemistry codes. Further, ChronusQ provides implementations of cutting edge numerical linear algebra algorithms such as the GPLHR method, and MOR approximations for problems arising in linear response theory. To the best of the authors’ knowledge, ChronusQ provides the sole implementation of MOR functionality for linear response problems available to the public.

ChronusQ is a free (as in software freedom) quantum chemistry software package which is open to the public. Contributions from the quantum chemistry community as a whole are not only allowed, but encouraged. ChronusQ may be obtained online through our website http://www.chronusquantum.org. This website also contains the ChronusQ Wiki and user manual which serve as the most up-to-date references for the functionality contained in ChronusQ.

ACKNOWLEDGMENTS

The development of Chronus Quantum ChronusQ computational software is supported by the National Science Foundation (OAC-1663636 to X.L. and A.E.D.). The method development of magneto-optical spectroscopy is supported by National Science Foundation (CHE-1856210 to X.L.), and the development of relativistic electronic structure theory is supported by the Office of Science, US Department of Energy (DE-SC0006863 to X.L.). E.V. was partially supported by the U.S. National Science Foundation (award 1550456).

AUTHOR CONTRIBUTIONS

The overall structure and modern C++ design of ChronusQ is due to David B. Williams-Young (DBWY). Parallelism in ChronusQ is also due to DBWY. The name “Chronus Quantum” is due to Joshua J. Goings (JJG). The development of memory management utilities in ChronusQ is due to DBWY, Torin F. Stetina (TFS), and Andrew Wildman (AW). The development of the ground state SCF module is due to DBWY, Patrick Lestrange (PL), TFS, and Xiaosong Li (XL). The development of DFT related methods is due to Alessio Petrone (AP), DBWY and TFS. ChronusQ primarily utilizes the Libint library to evaluate molecular integrals over contracted Gaussian basis functions. The Libint library is due to Edward F. Valeev (EFV), and its integration into ChronusQ is due to DBWY and EFV. ChronusQ also contains an in house integral code which performs molecular integral evaluations for those integrals which are not provided by Libint, such as those required by the GIAO and X2C methods. The in house integral code is due to Shichao Sun (SS) and XL. The development of GIAO capabilities in ChronusQ are due to SS and XL. The development of the X2C method in ChronusQ is due to DBWY and Joseph Kasper.

The RT module is due to DBWY, AP, AW, SS, JJG, and Feizhi Ding. The MMUT algorithm implemented in the RT module is due to XL. The RESP module and the development of the MOR capabilities in ChronusQ are due to DBWY. The extension of DFT and GIAO capabilities to the RESP module are due to AP and SS, respectively. Validation of the RESP module is in large part due to TFS and AW.

The development of embedding methods in ChronusQ is due to Chad E. Hoyer. The development of post-SCF methods such as Coupled Cluster and Møller-Plesset perturbation theory is due to DBWY, Daniel R. Nascimento, Lauren Koulias, and A. Eugene DePrince III. The tensor contractions required by these post-SCF methods are efficiently performed by the TiledArray library, which is due to EFV.

APPENDIX: MODERN C++ SOFTWARE DESIGN

The adoption of modern C++ in ChronusQ has allowed for a level of abstraction throughout the code base which has enabled the development of expressive, reusable and efficient code. ChronusQ was originally written in C++11, but more recently we have adopted a number of features in the C++14 standard (such as generic lambda expressions, variable templates, etc.) to decrease verbosity and improve expressiveness throughout the code. Current development streams of ChronusQ require features from C++17 (such as polymorphic memory resources, template deduction, etc.); however, these features will not be merged into the release stream until the standard is more widely supported in vendor C++ compilers. In this appendix we review a number of the more salient aspects on the modern C++ design of ChronusQ.

Template Metaprogramming

The use of template metaprogramming in ChronusQ is primarily motivated by the fact that we wish to develop code which supports both real and complex wave functions with maximal code reuse and minimal duplication of type-generic implementation details. To this end, the use of templates allows for the specialization of a certain number type-specific implementation details while allowing the developer to express their algorithms in a type-generic manner.

In electronic structure theory, the majority of the type-specific implementation details are encapsulated in the use of linear algebra software. Such operations include matrix multiplication, diagonalization, etc. Implementations of the BLAS and LAPACK standards are typically expressed in either the C or FORTRAN programming languages such that their calling APIs are type-specific, i.e. the multiplication of two real (double precision) matrices (DGEMM) is performed using a different function than the one that performs the multiplication of two complex (double precision) matrices (ZGEMM). However, the generic structure of these APIs (i.e. the general information one needs to pass to perform the operation such as data pointers, dimensions, etc.) are invariant between the function calls. To this end, ChronusQ implements templated wrappers around a functional subset of the BLAS and LAPACK standards to provide uniform calling API which allows for their integration into type-generic algorithmic structures. In addition, ChronusQ also provides type-generic templated wrappers for the MPI API to a similar effect. We note for completeness that similar functionality is currently being developed and standardized in the BLAS++, LAPACK++, and SLATE software libraries [139]. There are currently development efforts to integrate such software into ChronusQ.

Type Independent API and Dynamic Polymorphism

Due to the fact that ChronusQ supports many combinations of methods throughout the code base (i.e. non-relativistic, time-dependent, X2C, GIAO, etc.), the majority of the core functionality must be templated over both a type which describes the wave function and a type which describes the basis functions to maintain a type-generic development platform. However, because C++ is strongly typed, providing a type-generic API for type-agnostic operations such as performing the self consistent field (SCF) optimization would require the calling program to implement a combinatorial number of paths to perform these operations for every allowable combination of template parameters. To solve this problem in ChronusQ, we have adopted an object oriented code structure which separates type agnostic interfaces from type dependent implementations through the use of dynamic polymorphism.

To demonstrate this state of affairs, we examine the code snippet in Fig. 2 which describes the implementation of the SingleSlater class in ChronusQ. SingleSlater implements the minimum necessary functionality (i.e. formation and storage of the Fock matrix, contractions with the orbital Hessian for response calculations, etc.) for the description of the electronic wave function as a single Slater determinant. In this example, SingleSlaterBase provides a type erased API which is accessible from the calling program to perform type agnostic operations. Additionally, SingleSlaterBase provides the required API specification for functionality in which the underlying implementation must be type dependent, i.e. Fock matrix formation and storage, etc. These functions are left as pure virtual to ensure that any derived class must provide their implementations. Using the derived implementations of the type dependent operations, one may implement operations such as the SCF optimization in a type agnostic manner.