Gator: A Python-driven program for spectroscopy simulations using correlated wave functions

The Gator program has been developed for computational spectroscopy and calculations of molecular properties using real and complex propagators at the correlated level of wave function theory. Currently, the focus lies on methods based


| INTRODUCTION
Quantum chemical calculations have become a corner stone for fundamental and applied research in molecular sciences. Simulations of spectroscopies and predictions of molecular properties are made available through various formulations of response theory and allow for the revelation of molecular structure-function relationships that are indispensable for experimental spectrum characterization and rational material design. 1,2 However, any attempt to compute molecular spectra and/or properties remains a compromise between accuracy and applicabilityparticularly so for large molecular systems embedded in complex environments. 3 While reliable quantum chemical methods exist to address the electronic ground state, electronic structure theory as a whole and efficient computer programs in particular are lagging behind the current needs in the field of computational spectroscopy. 4 A wide range of semiempirical and density functional-based methods exist, which are applicable to a wide range of large molecular systems, 5, 6 but they are not reliable and predictive without extensive benchmarking and/or comparisons with experiments. This turns out to be an increasingly difficult way forward with growing system complexity as suitable experiments are often missing and highly accurate methods are computationally unaffordable. 7 Therefore, the development of wave function-based ab initio methods and the provision of efficient computer program implementations are in high demand, not least because such methods are systematically improvable and usually possess predictive power. In contrast to semiempirical and density functional-based approximations, the physical effects contained in the ab initio Hamiltonian are a priori well defined, and from the outset, it is clear which properties and spectra can be described in a physically correct manner.
The focus of the Gator program lies on single-reference methods, currently using the algebraic diagrammatic construction (ADC) scheme for the polarization propagator up to the third order in the fluctuation potential [8][9][10][11] as implemented in the ADC connect (adcc) toolkit. 12 It had been inspired by the adcman module contained in Q-Chem 13, 14 but evolved into an independent, open-source hybrid Python/C++ module. Based on the capabilities of adcc, the Gator program enables the calculation of molecular properties using linear real and complex response functions. In addition, Møller-Plesset (MP) perturbation theory and ADC (2) have been implemented with efficient exploitation of high-performance computing (HPC) resources. This implementation is accomplished using a distributed auxiliary Fock matrix-driven tensor contraction scheme with hybrid MPI/OpenMP parallelization. Overall, the Gator program is designed to provide a flexible platform for (i) high-performance ab initio quantum chemistry, (ii) scientific extensions (methods and algorithms) with relative ease and quick turnover time, and (iii) educational training with a high degree of student interactivity. The available features of Gator are described in more detail in the corresponding subsection.
With partly different scientific objectives but an otherwise largely shared philosophy and intentions, similar program development efforts have recently been undertaken, for example: (i) the Psi4NumPy 15, 16 Python programming environment for facilitating the use of the self-consistent field (SCF) and post-Hartree-Fock (post-HF) kernels from the Psi4 program through NumPy 17 ; (ii) the PySCF program platform 18,19 for Python-based SCF and post-HF electronic structure theory calculations of finite and periodic systems; (iii) the molsturm 20 quantum chemistry framework, which provides a contraction-based and basis function-independent SCF, integrating readily with existing integral libraries and post-HF codes; and (iv) the Dalton Project 21 providing a platform written in Python for the interoperability of quantum chemical software projects. There exist several other efforts to provide massively parallel quantum chemistry programs, which do not necessarily exploit a Python-based modular program setup. For an overview, see Refs. 22, 23.

| Algorithmic considerations
The performance-determining tasks in post-HF spectroscopy calculations are the atomic-to-molecular orbital integral transformations followed by the tensor operations defined in the iterative solution of the underlying response equations. Conventional implementations suited for shared-memory computers typically store the full tensor of transformed antisymmetrized two-electron integrals in memory. While this approach works well for calculations of molecular systems up to about 700 one-particle basis functions on standard computers, it quickly becomes unfeasible for larger systems due to a memory bottleneck. As a complement to a conventional integral transformation, we have also implemented an MPI/OpenMP parallel auxiliary Fock matrix-driven transformation routine. The basic idea of this distributed integral transformation is to express quantities on the molecular orbital (MO) basis (such as integrals and tensor contractions) in terms of auxiliary Fock-like matrices. This has previously been exploited in an atomic orbitalbased formulation of the complete active space SCF method on graphical processing units. 24 For electron-repulsion integrals in the MO basis, we use the simple relation where Here, i, j, k, l and α, β, γ, δ refer to MO and atomic orbital (AO) indices, respectively, and N is the number of basis functions. The basic advantage is that it allows these quantities to be constructed incrementally, and by using a split MPI communicator, the work of constructing this set of auxiliary Fock matrices can be divided into the available cluster nodes. Any MO integral resulting from Equation (1) remains afterward on the node of its evaluation, being ready, for example, for the distributed calculation of the MP2 energy or ADC matrix-vector products.
For the ADC part, Gator relies on contraction-based iterative algorithms for solving response and eigenvalue equations as has been described previously for adcc. 12 In such schemes, the ADC matrix M is never explicitly constructed, but instead, the secular and response equations are solved by calculating matrix-vector products with sets of trial vectors {b i }: Formally, the ADC(2) vectors consist of a single-excitation part and a much larger double-excitation part, which can be folded into the space of single excitations. 25 We implemented this algorithm within the HPC-QC module of Gator as vectors in the single excitation space can be communicated between nodes at a negligible cost. The expressions for the ADC(2) matrix-vector product in Equation (3) have thus also been derived and implemented with the use of the auxiliary Fock matrix-driven MO transformation.
For smaller systems, the computational cost of the auxiliary Fock matrix-driven approach is not beneficial, and hence, both variants of integral transformation and calculations of σ-vectors are made available in the Gator program by VeloxChem. 26 With AO integral screening based on a combination of the Cauchy-Schwarz inequality and the density, the calculation of auxiliary Fock matrices scales somewhat better than N 3 . As there are a total of N 2 pair indices (i, j), this results in a total scaling of less than N 5 for the auxiliary Fock matrix-driven integral transformation in Equation (1).

| Modular structure of Gator
The Gator program consists of three major modules as illustrated in Figure 1: (i) HPC-QC for MP2 and ADC(2) calculations in HPC environments, exploiting the Fock matrix-driven integral transformation outlined above; (ii) Respondo for the evaluation of real and complex linear response functions; and (iii) ADC connect (adcc) 12 for excited-state calculations based on the ADC scheme up to the third order. The Hartree-Fock reference state and the required one-and twoelectron integrals are provided by the VeloxChem 26 program.
The modules of Gator are developed in a two-layered fashion. The upper layer is written in Python and includes the management of the hardware resources, input/output handling, result processing, and high-level parts of some algorithms. Numerical array operations are carried out using NumPy 17 and SciPy. 27 Parts of the lower program layer are written in C++, providing functionality for the evaluation of tensor contractions that, in adcc, are performed by libtensor. 28 The two layers communicate by means of the pybind11 library. 29

| Implementation and efficiency of HPC-QC
In addition to the standard implementations of MP2 and ADC(2) available within the adcc module, we have implemented new and highly parallel versions of MP2 and ADC(2) that exploit the auxiliary Fock matrix-driven MO transformation within the HPC-QC module. VeloxChem is responsible for the construction of the auxiliary Fock matrices, including the underlying evaluation of two-electron integrals on the AO basis using a hybrid MPI/OpenMP parallelized module for efficient execution on HPC clusters. The module accepts, depending on memory limits, an arbitrary number of density matrices in a batch and handles those in parallel. With a single evaluation of the set of two-electron integrals, the corresponding complete batch of auxiliary Fock matrices is constructed. Such a multiauxiliary Fock matrix construction leads to additional costs for integral distribution. As a larger batch size comes at the cost of reduced OpenMP parallel efficiency in the integral evaluation, it may become necessary to balance the two forms of parallelism in the integral transformation depending on the size of the system and the adopted basis set.
In the same way, HPC-QC is hybrid parallelized. With regard to MO integral storage, a split MPI communicator is used across cluster nodes, and any given integral naturally resides in the memory of the node executing Equation (1) to enable usage of the aggregated memory of the available cluster nodes, as well as to minimize the need for communication. For example, the storage of an integral block hoojjvvi can be distributed over the occupied-occupied (oo) or the virtual-virtual (vv) pair indices. The distributed storage of MO integrals allows for a parallel implementation of ADC (2), where the tensor operations take place on the individual compute nodes that include a list of distributed pair indices and the corresponding auxiliary Fock matrices on the MO basis. The construction of matrix-vector products (σ-vectors) in ADC(2) is evaluated by looping over the distributed pair indices and contracting the corresponding auxiliary Fock matrices on the MO basis with the trial vectors of single excitations. Storage of the large vectors of double excitations is avoided by packing them into the space of single excitations. 25 As such, parallel ADC(2) calculations can be carried out efficiently on cluster nodes with only moderate amounts of memory.
In Figure 2, we show results for guanine oligomers. We investigated the scaling of the MO integral calculation compared to the construction of σ-vectors with respect to the number of contracted basis functions (gray curve in Figure 2). The MO integral calculation involves all those necessary at the level of ADC(2) theory. For guanine oligomers with the def2-SVP basis set, the MO integral calculation scales as N 4.9 in line with the discussion above. The evaluation of MO integrals therefore uses the majority of the computational time (core-hours) in the parallel implementation of ADC (2) as the construction of σ-vectors shows a lower scaling of N 3.8 .

| Capabilities
Calculations using the ADC(2), 8 ADC(2)-x, 9 and ADC(3) 10 methods can be performed using the adcc module for shared-memory architectures. For spectrum calculations in the soft X-ray region, the core valence separation (CVS) is F I G U R E 1 Gator program capabilities and module overview available. [30][31][32] All excited state and transition properties implemented in adcc are also available in the Gator program. Absorption cross sections and C 6 dispersion coefficients can be obtained through the complex polarization propagator (CPP) approach, [33][34][35][36][37][38] implemented in the Respondo library. Calculations with explicit solvents are currently available using the polarizable embedding (PE) model [39][40][41] through the CPPE library 42 as interfaced with VeloxChem.
MP2 energies and ADC(2) excitation energies can additionally be obtained by employing the HPC-QC module using the distributed aggregated memory on HPC clusters.

| HOW TO USE GATOR
To illustrate the modular, low-barrier nature of our Python-based environment, we outline the calculation of the X-ray absorption spectrum of water in Figure 3. This calculation was performed with the Gator program in a Jupyter notebook, 43 which is an open, web-based interactive environment of increasing popularity. By dividing the code into multiple blocks, Jupyter enables the calculation of individual code sections, keeping results in local memory. This simplifies rewriting and recomputing code snippets, such that the user does not need to recalculate all quantities after a modification is performed. This approach allows for an interactive path to programming, analysis, and education. In Gator, we take advantage of these concepts by combining ab initio computational spectroscopy with intuitive analysis and plotting routines through NumPy 17 and Matplotlib. 44 Molecular systems can be investigated in an iterative manner by varying parameters such as basis sets, convergence criteria, numerical solvers, and convolution/broadening functions.
Gator also offers traditional job submission using input files. An example input file is shown in Figure 4.

| NUMERICAL EXAMPLES
Within the electric dipole approximation, the absorption cross section of a randomly oriented molecular system is obtained from the isotropic average of the complex polarizability 33,35,37 : F I G U R E 3 Gator notebook mode: The notebook input for the calculation of the X-ray absorption spectrum of water (upper left) and program output (lower left), focusing only on information for the requested five excited states. The resulting spectrum plot is shown on the right F I G U R E 4 Gator program input file for performing the core valence separation-complex polarization propagator-algebraic diagrammatic construction (CVS-CPP-ADC(2)) spectrum calculation shown in Figure 3 allowing for the calculation of an absorption spectrum at arbitrary energy intervals by solving the associated linear response equations over a range of frequencies.
To illustrate the capabilities of the Gator program, we computed the UV/vis and X-ray absorption cross sections, the static polarizability, and the polarizability at imaginary frequencies, yielding C 6 dispersion coefficients 38 for furan and thiophene as shown in Figure 5. A damping term γ = 0.24 eV was used for the calculation of absorption cross sections. For comparison, the excitation energies and intensities obtained from explicit calculations of 25 eigenstates are convoluted using a Lorentzian function with the same broadening parameter. All property calculations were performed using the 6-311++G** 45 basis set. The structures of the molecules were optimized at the B3LYP 46 /cc-pVTZ 47 level of theory using the Q-Chem 5.1 program. 14 With imaginary frequency arguments, we obtain a smooth monotonous polarizability as illustrated in the left panels of Figure 5, resulting in an isotropic static polarizability ( α 0 ) and a C 6 coefficient of 47.07 and 906.3 a.u. for furan and 62.34 and 1438 a.u. for thiophene. Experimental static polarizabilities are reported to be 48.59 and 65.18 a.u. for furan and thiophene, respectively-see Ref. 51 and references therein.
For the calculation of valence excitation spectra, Gator offers the possibility to either solve explicitly for eigenstates or to scan over relevant frequency intervals using damped response theory. Using a Lorentzian broadening for the transitions in the first approach will result in a spectrum that is practically identical to that obtained in the second approach, as illustrated in the middle and right panels of Figure 5. The linear absorption cross section obtained from (CVS-)CPP-ADC and the broadened eigenvalue results agree perfectly for lower energies. However, for higher energies, differences between the two curves occur because the resolved 25 eigenstates are not sufficient to obtain the full spectrum window.
The carbon K-edge spectra of furan and thiophene obtained with CVS-ADC(2) are given together with the experimental gas-phase spectra as insets in the right panels of Figure 5. Here, we include the first use of CVS-CPP-ADC, which has been developed and implemented in Gator. The use of the CVS approximation in addition to CPP here serves multiple purposes: (i) reduced matrix dimensions and a resulting decrease in computational effort; (ii) removal of spurious valence-continuum features that may contaminate the resolved spectrum window; and (iii) reducing difficulties in converging pure CPP-ADC calculations in the energy region of core excitations. Convergence issues have also been observed for CPP calculations of core excitations in the context of coupled cluster theory, and we note that the use of the CVS scheme introduces only a minor scalar shift in absolute energies. The overall discrepancy with respect to experiment amounts to a blue-shift of 0.5 eV. 52 To demonstrate the capabilities of Gator in computing the absorption spectra of solvated molecules at the ADC (2) level of theory, the spectrum of noradrenaline has been computed in the gas phase and in the aqueous solution using the standard eigenvalue solver, as well as the CPP approach. The vacuum structure of noradrenaline was optimized at the CAM-B3LYP/def2-svp 53, 54 level of theory as implemented in ORCA 4, 55 which served as input for the subsequent calculation of the seven energetically lowest singlet excited states at the ADC(2)/cc-pVDZ level of theory. Then, CPP-ADC(2) calculations were carried out on the frequency range obtained via the regular ADC(2) approach with a spacing of 0.1 eV, yielding the one-photon absorption cross section for each frequency; see Figure 6.
To model the UV/vis spectrum of noradrenaline in solution, it was placed in a box of size 48 × 44 × 43 Å 3 of water molecules using packmol, 56 and molecular dynamics (MD) simulations were carried out with NAMD 2.12 57 using the CHARMM 36 force field 58 and parameters for noradrenaline obtained with CGenFF. 59 The system was equilibrated for 5 ns at 310 K and a pressure of 1 atm (NPT ensemble). A subsequent hybrid QM/MM minimization was performed with NAMD, 60 where noradrenaline represented the QM region, and the MM region was assigned to the solvent. The QM computations were run at the CAM-B3LYP/def2-SVP level of theory with ORCA 4 using an electrostatic embedding scheme. The final snapshot of the quantum mechanics/molecular mechanics (QM/MM) minimization was extracted and parameterized using PyFramE 61 with a standard potential model. In the spectrum calculation, the solvent effects are modeled self-consistently with the PE model applied for the Hartree-Fock reference state. [39][40][41] The absorption spectra of noradrenaline in vacuum and solution are shown in Figure 6, where we note that energies in the latter calculation are not perturbatively corrected. 41

| GATOR FOR TRAINING AND EDUCATION
The Gator program is a convenient platform for undergraduate training and education in computational and theoretical chemistry with students ranging from the introductory bachelor to the advanced master level. It offers insight into the underlying general theory of quantum chemistry and the adopted electronic structure theory methods with direct and concrete access to the intermediate quantities (tensor elements and wave functions) that are addressed in the workflow of quantum chemical calculations. This is accomplished by interfacing the Gator program with the Jupyter notebook as described above. This allows the student to steer, control, and analyze the course of a desired quantum chemical calculation in a modular manner. Students who are initially not familiar with programming are naturally motivated to overcome barriers involved with entering the field of theory and program development. With basic knowledge of Python programming, it is thus possible for students to design and implement key algorithms such as, for example, the SCF F I G U R E 6 UV/vis spectra of noradrenaline in vacuum and solution. Continuous lines show the convoluted stick spectra, obtained from a Lorentzian function with γ = 0.124 eV. Dots depict the one-photon absorption cross section σ(ω) computed with complex polarization propagator (CPP) optimization of electron densities or the response theory calculations of excitation energies. Following this philosophy, it is possible to incorporate such practical exercises into existing courses, which has been successfully demonstrated already by the Psi4 16 community for example. The deepened algorithmic understanding gained by the students will be advantageous when entering the research level in this field.

| OUTLOOK
The Gator platform will undergo development following several major lines: (i) The HPC-QC module will be extended to include implementations for the calculation of transition and excited-state properties, as well as the inclusion of external one-particle potentials exploiting the intermediate state representation. 62 (ii) The complex linear response function will be implemented in the HPC-QC module. (iii) The capabilities of the Respondo library will be further extended by including higher-order response functions to allow for simulations of advanced linear and nonlinear spectroscopies. (iv) Quantum chemical methods other than the ADC approaches will be included with the primary aim at approximate coupled cluster theory of second order (CC2) 63 and also equation-of-motion and linear-response coupled cluster methods with singles and doubles. [64][65][66] (v) Modular libraries for environment models other than the PE model are under development. (vi) The engine of the HPC-QC module, that is, the efficient construction of the large number of auxiliary Fock matrices on the AO basis, will be adapted for execution on heterogeneous cluster nodes with varying combinations of central processing units and graphics processing units. These extensions will be realized as developments of the already existing modules (i)-(iii), while (iv) will require the development of a new module of its own.

| SOFTWARE AVAILABILITY
Gator is freely available under the GPLv3 license. The source code can be found on GitHub (https://github.com/gatorprogram/gator). Input files and notebook examples can be found on the GitHub website. Furthermore, the package can be conveniently installed using the conda package manager (https://anaconda.com) by running the command conda install gator -c gator, which will also automatically install all other dependencies.

RELATEDWIREs ARTICLES The ORCA program system
Molpro: A general-purpose quantum chemistry program package VeloxChem: A Python-driven density-functional theory program for spectroscopy simulations in high-performance computing environments