The generality of the GUGA MRCI approach in COLUMBUS for treating complex quantum chemistry

COLUMBUS is a collection of programs for high-level ab initio molecular electronic structure calculations. The programs are designed primarily for extended multi-reference (MR) calculations on electronic ground and excited states of atoms and molecules. Since its early versions,^1,2 the COLUMBUS program system^3,4 was always at the forefront of the development of proper methodology to solve chemically challenging problems, relying, of course, on the actual state-of-the-art theoretical knowledge and computer architecture. The primary focus in the 1980s was small molecules and mostly properties related to electronic energy. Besides the MR methodology based on the standard non-relativistic Hamiltonian, the treatment of relativistic effects in the form of spin-orbit (SO) configuration interaction (CI) was a unique feature.⁵ When the analytic gradient for the MRCI energy was developed and implemented in 1992,⁶ a new avenue of applications became accessible, which allowed the optimization of structures not only in the ground but also in the excited state. The next important development came in 2004 with the derivation and efficient programming of nonadiabatic couplings.^7,8 This feature opened up a new field of application toward photochemistry and photodynamics. The combination of COLUMBUS with dynamics programs such as NEWTON-X and SHARC allowed highly competitive simulations of nonadiabatic processes. The size of the molecules COLUMBUS can handle increased significantly over the years. While in the 1970s, calculations were restricted to few atoms, nowadays, the treatment of molecules with over 100 atoms is possible, depending on the reference wavefunction, symmetry, basis set, and other factors. Due to the timely response of COLUMBUS developers to the appearance of new parallel computer architectures, COLUMBUS was pioneering parallel execution with the help of the Global Array toolkit.^9–12 Today, COLUMBUS efficiently runs on mainframe computers, as well as on computer clusters, and allows applications in many fields of chemistry, materials science, and biochemistry.

COLUMBUS is dedicated to variational calculations based on multireference configuration interaction with single and double excitations (MR-CISD)¹³ and related methods, of which the MR averaged quadratic coupled-cluster approach (MR-AQCC)^14,15 is probably the most popular one because it includes size extensivity corrections. The program can also perform calculations on excited states in the form of a linear-response theory (LRT).¹⁶ A formulation optimizing the total energy (TE) in place of the correlation energy (TE-AQCC) is available as well.¹⁷ In COLUMBUS, the expansion of the wavefunction is performed in an uncontracted (uc) form in which no internal contraction (ic) of the reference wavefunction is used. For more details on this point, see Sec. II G. All wavefunction-related aspects of COLUMBUS are based on the graphical unitary group approach (GUGA), as developed by Shavitt.¹⁸ The significant advantage of this uc expansion is its flexibility, which allows the straightforward implementation of analytic energy gradients at MR-CISD and MR-AQCC levels and nonadiabatic couplings at the MR-CISD level.^6–8,19,20 An exceptional feature of COLUMBUS is the ability to perform full two-component, SO-MR-CISD calculations.^5,21

In addition to these well-established methods, several new approaches are being developed, which will appear in future releases of COLUMBUS. One of the focuses of this paper is to outline these emerging developments, which are presented in Sec. IV. Two of these techniques are the all configuration mean energy (ACME) multiconfiguration self-consistent field (MCSCF) and the graphically contracted function (GCF) methods, which allow vast configuration expansion spaces (up to 10¹⁵⁰ configuration state functions). We also describe local electron correlation schemes for MR approaches and a spin density approach within GUGA. In an extension of the nonrelativistic energy gradient approach mentioned above, analytic spin-orbit energy gradients are also being developed.

After this overview of COLUMBUS capabilities, it is worthwhile to discuss somewhat in more detail what the real focal points in terms of applications are. In a world of increasing research in materials science dedicated to the development of new compounds with interesting magneto-optical properties derived from biradicaloid character, with new demands of utilizing and understanding photodynamical processes and of dealing with complicated open-shell systems in transition metals, lanthanides and actinides, the requirements on the flexibility of programs for electronic structure theory have risen significantly. The examples mentioned, and many others, demand MR methods because of the intrinsically complicated electronic structures involved in the problems. This is the point where COLUMBUS excels. The uncontracted nature of the CI expansion provides the required high flexibility and allows precise benchmark calculations. As mentioned above, the simplicity of the variational calculations of this formulation is also the basis for the availability of analytic MR energy gradients and nonadiabatic couplings. These are features, which stand out, and are not shared by many other MR program packages. Thus, because of the availability of analytic energy gradients at the MRCI level, consistent geometry optimization at the same high-level method can be performed as the final energy calculation. This situation must be contrasted to the case where, because of the lack of analytic energy gradients for the high-level method, geometry optimizations need to be performed at a lower level.

The necessarily larger amount of computational effort can be attenuated by various selection schemes applied to the reference wavefunction and by an efficient parallelization of the MRCI step. There have been several other ways developed in the literature to reduce the computational demand of MR methods, usually achieved by introducing additional approximations or restrictions. These include low-order perturbation theory (PT), the equation-of-motion (EOM) approach, low-order PT treatment of spin-orbit, and internal contraction.^13,22,23 In the spin-orbit CI case, in particular, COLUMBUS fully treats strong correlation, weak correlation, and spin-orbit coupling; other codes make compromises and approximations to one or more aspects of those three effects. Thus, to verify the validity of these other methods and their applicability in various contexts, comparisons must be made to more accurate methods without these additional approximations. COLUMBUS has served that purpose for almost 40 years, and it will continue to do so.

Beyond this benchmark role, the available procedures in COLUMBUS are so efficient that reliable production work can be done on many interesting problems, at a precision level that is hardly achievable with other approaches. The second focus of this paper is on delivering a showcase of many examples of applications using COLUMBUS for electronic structure (Sec. II) and nonadiabatic dynamics problems (Sec. III). These examples, spanning fields of materials science, biological sciences, atmospheric chemistry, and heavy metal chemistry, should provide a practical guideline for applying the methods available in COLUMBUS.

COLUMBUS is freely available from the website https://www.univie.ac.at/columbus/. It includes executables for a simple compilation-free installation of the serial code along with the source code for the compilation of the parallel section of the CI calculation. The COLUMBUS webpage contains detailed documentation and tutorials, which introduce the user to the main application types. Moreover, the tuition material of several COLUMBUS workshops, also available from the same webpage, provides a host of information about theoretical procedures and applications.

The understanding of pancake bonding^24,25 provides a unique challenge to electronic structure theory.^26–32 It requires a theory level that can treat MR ground states, an accurate level of dynamic electron correlation, and geometry optimization at the same high computational level, all of which are available in COLUMBUS. In a typical pancake-bonded dimer, two π-conjugated radicals are bonded together in a π-stacking configuration with direct atom–atom contacts shorter than the sum of the van der Waals radii.³³ For example, C⋯C inter-radical contacts in pancake bonded dimers are often close to 3.0–3.1 Å compared to the vdW value of 3.4 Å. In addition to this characteristic geometry, pancake bonded dimers possess highly directional interactions with maximum multicenter overlaps³⁴ and low-lying singlet and triplet states. Highly conducting organic materials^35–39 often display this interaction because the strong orbital overlap between π-radicals is central to designing new organic conductors. Figure 1 illustrates a few examples of molecules forming pancake-bonded dimers and a molecular orbital (MO) diagram for the phenalenyl (PLY) dimer.

The computational challenge arises from the fact that this is not only a fundamentally multiconfigurational ground state problem, but it is burdened with the added complexity that the dispersion interactions need to be also included; the latter require a huge number of configurations to be described accurately, for which the MR-AQCC theory, combined with geometry optimization at the same level, appears to be especially appropriate. This theory has been applied to the three prototypical pancake bonded problems: the binding energy and conformational preferences in the PLY dimer,⁴⁶ PLY₂, the stability of the (TCNE⁻)₂ dimer,⁴⁷ and the problem of multiple pancake bonding⁴⁸ in a dimer of thiatriazine (TTA) and other systems. In the first two cases, we used MR-AQCC(2,2)/6-31G(d). In the latter case, we were able to use MR-AQCC(4,4)/6-311++G(2d,2p).

The strong preference for the atom-over-atom configuration, which is missing in vdW complexes, is present in the pancake bonded ones because the partial electron pairing favors this configuration. For instance, in the PLY₂ case, the torsional rotation around the C₃ axis shows that the electron pairing is completely broken at 30°, as indicated by the very large rotational barrier of 17.2 kcal/mol and by the increase in the number of effectively unpaired electrons (N_U),

in which n_i is the occupation of the ith natural orbital (NO) and M is the total number of NOs as computed using the nonlinear formula of Head-Gordon.⁴⁹ At 30° torsion, N_U becomes nearly equal to that of the triplet, indicating a broken pancake bond. Why then is the binding energy of the PLY₂ dimer only 11.5 kcal/mol? An approximate energy decomposition shows that at the short equilibrium C⋯C distance of 3.1 Å, the vdW component of the energy (including Pauli exclusion repulsion, dispersion, and a small electrostatic term) is overall repulsive at +5.7 kcal/mol, while that contribution from the SOMO-SOMO electron pairing is significantly stronger at about −17.2 kcal/mol. This approach solved the long-standing problem of explaining the strong preference in pancake bonding for the atom-over-atom overlap.

Pancake boding occurs in more complex aggregates as well, such as trimers and stacked chains, offering further challenges for the electronic structure community.

Aromatic diradicals are unusually stable as a result of their resonant aromaticity yet highly reactive because of their unpaired electrons. For instance, the open-shell para-benzyne diradical can be formed by gentle heating of enediyne [(Z)-hex-3-ene-1,5-diyne], resulting in a reactive intermediate that is only 8 kcal/mol less stable than the closed-shell reactant.⁵⁰ Aromatic molecules with proximate unpaired electrons (ortho-orientation) often display triple-bond-like features, while for many isomers, the ground state is a multiconfigurational singlet state due to through-bond coupling.⁵¹ In these cases, a single reference quantum method utilizing just one electronic configuration or Slater determinant is not able to accurately capture the physical nature of the system. Aromatic and heteroaromatic diradicals are electron-rich, and the complexity of their electronic structure results in a high density of closely spaced electronic states. This adds to the complexity of properly characterizing these systems. Often the best results are obtained by using a state-averaged approach to optimizing the multireference (MR) wavefunctions for all states nearby in energy, followed by single-state calculations at a higher, correlated level of theory.

The methods available in COLUMBUS are particularly well-suited for performing highly correlated MR calculations, including single-state and state-averaging (SA) approaches. For instance, the para isomer of benzyne is a two-configurational ground state singlet.⁵² As shown in Fig. 2(a), para-benzyne contains a very high density of close-lying electronic singlet states. Using COLUMBUS, a 32-state averaged MR-CISD/TZ calculation with 4 states from each of the 8 irreducible representations under D_2h symmetry was performed. From this, it was determined that there were 25 valence and 5 Rydberg singlet states within 10 eV of the ¹A_g ground state. The proper selection of active spaces is also essential in MR calculations. We used COLUMBUS to optimize the active space calculation for the 9,10 didehydroanthracene molecule [Fig. 2(b)].⁵³ Nine different active spaces were explored, and the MCSCF natural orbital populations were used to determine the ideal active spaces. Using the natural orbital populations, orbitals that were either unoccupied or close to doubly occupied were considered less important for inclusion in the active space than orbitals with partial occupation. Dynamical correlation between all orbitals, including those not included in the active space, was included via the MR-CISD and MR-AQCC electronic excitations. Using the optimized (8,8) active space with the MR-AQCC/TZ method revealed a ground state singlet for the anthracene diradical with 6.13 kcal/mol and 7.18 kcal/mol (0.265 eV and 0.311 eV) adiabatic and vertical singlet–triplet splitting, respectively [Fig. 2(c)]. A final example of the utility of COLUMBUS for aromatic diradicals involves the characterization of the lowest-lying singlet and triplet state energies and geometries of the three (ortho, meta, and para) didehydro isomers of pyrazine [Fig. 2(d)].⁵⁴ Single point MR-AQCC/TZ calculations with a (12,10) active space reveals that the ortho (2,3) and para (2,5) isomers are ground state singlets with adiabatic gaps of 1.78 kcal/mol and 28.22 kcal/mol (0.0771 eV and 1.224 eV), respectively, while the meta (2,6) isomer is a ground state triplet that is nearly isoenergetic with the higher-lying singlet [ΔE(S,T) = −1.40 kcal/mol (−0.061 eV)]. The singlet state of the meta isomer lies higher in energy than either the ortho or the para singlet state. A bonding analysis suggests that this is the result of unfavorable three-center-four-electron antibonding character in the 11a₁ σ orbital of the meta isomer. The relatively large adiabatic gap for the para (2,6) isomer is likely caused by through-bond coupling effects stabilizing the singlet state.

In recent years, polycyclic aromatic hydrocarbons (PAHs) with certain radical character have attracted large interest due to their exceptional optical, electronic, and magnetic properties.^55–58 Among various open-shell singlet PAHs, zethrenes characterize an attractive class of compounds, showing interesting optical properties in the near-infrared region.^59–62 Tuning the biradicaloid character and balancing it against the enhanced chemical instability is the key challenge for the development of useful materials.⁶³ Replacing sp²-carbons with sp²-coordinated nitrogen atoms while preserving the planar π-scaffold geometry has received significant attention in efforts to modify the electronic structure of PAHs.^64–68 

To monitor the achievable range in biradicaloid character, all fourteen different nitrogen-doped heptazethrene (HZ) structures were created by symmetrically replacing two C atoms by two N atoms in such a way that the original C_2h symmetry was maintained.⁶⁹ All the structures were optimized using the second-order Moller–Plesset perturbation theory.⁷⁰ State averaging (SA) over the lowest singlet (1¹A_g) and triplet (1³B_u) states was performed with four electrons in five π orbitals, denoted as SA2-CAS(4,5). The same complete active space (CAS) was used in the MR-AQCC calculations. The expansion was confined to the π-space only. It has been shown previously that freezing of the σ orbitals had a negligible effect on the singlet–triplet splitting and density of unpaired electrons.⁷¹ 

Figure 3(a) represents the density of unpaired electrons for the ground 1¹A_g state of HZ [N_U = 1.23e, Eq. (1)]. The effect of nitrogen doping on the biradicaloid character and thereby also on the chemical stability/reactivity is summarized in Fig. 3(b) in the form of a N_U color map. It signifies for each of the 14 doping positions, the corresponding total N_U values. The reference value of pristine HZ is 1.23e, where the unpaired density is mostly located at the (1,15) position [Fig. 3(a)]. Doping of the two N atoms at that position leads to a complete quenching of the radical character,⁷² indicated by the deep blue color code with a total N_U value of only 0.58e [Fig. 3(b)]. On the other hand, doping of the two N atoms at the (13,27) positions leads to a strong enhancement of radical character with a total N_U value of 2.51e.

Therefore, in summary, if one desires to stabilize the pristine HZ toward the closed-shell, one has to look at the doping positions colored in deep blue, whereas the doping positions colored in red signify the creation of large biradical character. The separation between these two regions of stability/reactivity is quite prominent. Doping of the two N atoms in the phenalenyl region significantly quenches the biradicaloid character, whereas doping in the central benzene ring of HZ enhances it. In the latter case, the importance of the radical phenalenyl system is enhanced.

Graphene nanoflakes are attracting considerable interest because their electronic properties, including their bandgaps, can be effectively tuned by varying their size and the shape of their edges.^73,74 Correlated MR computations with COLUMBUS were used to study these systems and to elucidate the formation of polyradical character with an increase in the size of the nanoflake.^75,76 These computations do not only provide accurate results, but also mechanistic insight can be obtained through the visualization of the density of unpaired electrons according to Eq. (1).

For this work, we have performed MRCI calculations on two exemplary graphene nanoflakes. In Fig. 4(a), a rectangular periacene C₇₈H₂₄ with an unperturbed zig-zag edge six phenyl rings long and an armchair edge of length five is shown. For this molecule, the number of unpaired electrons N_U was equal to 2.16e, indicating biradical character. The distribution of these unpaired electrons is shown as an isocontour plot in Fig. 4(a), highlighting that the unpaired density is concentrated around the center of the zig-zag edge. As a consequence of unpaired electrons, the singlet ground state becomes almost degenerate with the first triplet state (T₁), and a gap of only 0.055 eV is obtained via the Davidson-corrected MRCI + Q_D approach.⁷⁷ The high concentration of unpaired density at the center of the zig-zag edge suggests that a perturbation at this position might strongly affect the biradical character. To test this hypothesis, we attached an additional phenyl ring at this position for each zig-zag edge, yielding a molecule with molecular formula C₈₄H₂₆ [Fig. 4(b)]. Indeed, adding this additional phenyl ring produces an almost closed-shell structure (N_U = 0.20e) with a significantly enhanced singlet–triplet gap of 0.987 eV.

The calculations were performed using MR-CISD based on a restricted active reference (RAS) space containing eight electrons in eight orbitals. To allow for these calculations with over 100 atoms, the σ-space was frozen during the calculations, and only the π-orbitals were correlated. A split-valence basis set was used, and following Ref. 78, an underlying atomic natural orbital basis set (ANO-S-VDZ) was used.⁷⁹ 

The calculations described in Sec. II D were performed using the parallel MRCI implementation in COLUMBUS,^12,80 which is based on the Global Array toolkit^9,10 for managing the parallel environment and one-sided data access. These calculations were performed on an HPE MC990X symmetric multiprocessing (SMP) system with 2 TB RAM and 8x Intel E7-8867 v4 CPUs, each of which has 18 cores. In the following, we will use the example of the triplet computation for the system shown in Fig. 4(b), C₈₄H₂₆, to discuss the performance characteristics of the parallel calculations.

The computation included 1.3 × 10⁹ configurations, i.e., the MRCI problem corresponds to determining eigenvalues of a matrix of dimension (1.3 × 10⁹) × (1.3 × 10⁹). In the Davidson algorithm used,⁸¹ it is not necessary to store the whole CI matrix, which is far beyond the available storage. Instead, the algorithm uses a direct CI approach⁸⁸ in which it is only necessary to store a small number of CI and product vectors of the Hamiltonian matrix. If each CI vector element is stored as a double-precision number (8-byte), then storing one CI vector requires about 10 GB of memory for 1.3 × 10⁹ configurations. The present calculation uses the default value for the maximum subspace dimension of 6. Thus, in total, 12 vectors (six CI vectors and six product vectors) have to be stored, totaling 120 GB. Any other data (including the molecular orbital integrals) only require significantly less memory, meaning that 120 GB is a good guess for the global memory needed. In addition to this overall global storage space, it is also necessary to reserve local working memory for each process.

Computations used between 12 and 72 processor cores and some crucial performance data are shown in Table I. Starting with 12 cores, we find that one iteration took about 24 min, while the complete iterative MRCI convergence procedure required 284 min. For every iteration, a total of 507 GB had to be transferred, which was primarily determined by the transfer of the CI vectors and product vectors. The computer system used allowed for a transfer rate of 1.68 GB/s, meaning that the communication time, distributed over all cores, required only a fraction of the total iteration. The number of cores was varied up until 72. Table I shows that the time consistently decreases as the number of cores is increased. Using 72 cores, we find that the time for one CI iteration is sped up by a factor of 5.53, which is close to the ideal value of 6. The speedup for the overall CI program is somewhat worse, reaching only 4.80, highlighting that the preparation steps are not as well parallelized as the actual iterations. In any case, the time for the overall CI program is already below 1 h for 72 cores and further speedup is usually not required.

Table I also shows that the data volume per iteration is largely independent of the number of cores, which is generally the case if enough memory is provided to the program. The last column shows the approximate amount of core memory needed to obtain the performance shown here. Using only 12 cores, a significant amount of 15 GB per core is needed, while this value is reduced to 3 GB for 72 cores. If needed, somewhat lower memory values could be accommodated but only at the cost of significantly increased communication.

The above discussion illustrates the hardware requirements for carrying out parallel MRCI computations on systems with a CI dimension of about 1 × 10⁹. Crucially, about 200 GB of memory in total are needed for such a calculation, and the data access must occur at the rate of about 1 GB/s, requiring either an SMP machine or high-performance interconnects between nodes. Conversely, the discussion shows that the number of cores is often only a secondary concern as reasonable runtimes are already obtained with a modest number of cores.

CF₂Cl₂ is one of the most abundant chlorofluorocarbons, and its lifetime of ∼112 years⁸² makes it a subject of great concern to the ozone layer depletion. Therefore, the study of its photochemistry is crucial to the environment. In this sub-section, we present potential energy curves for 25 singlet states along one C–Cl coordinate, at the MCSCF level, as well as full geometry optimization of an ion-pair structure at the MR-CISD level. Such a structure can explain the photochemical release of chloride.⁸³ 

First, a relaxed scan along one C–Cl coordinate has been performed at the CAS (12,8) level, where two σ bonds and four Cl lone pairs define the 12 electrons. The eight active orbitals comprise the two σ_C–Cl/σ_C–Cl^* pairs and the four Cl lone pairs. The C_s symmetry has been used, and nine valence states (5A′ + 4A″, including the ground and eight nσ^* states) have been averaged at the CASSCF/aug-cc-pVDZ^84–86 level, with equal weights.

Using the optimized geometries, single-point calculations have been performed at the MCSCF level, for the nσ^* and Rydberg [3s(C) and 3p(C)] states. The configuration state functions (CSFs) have been generated through the same CAS (12,8) scheme as above along with four auxiliary (AUX) orbitals, formed by the 3s(C) and 3p(C) Rydberg orbitals. Only single CAS → AUX excitations are allowed, and 25 states have been averaged at the MCSCF level [13A′ + 12A″, including 8 nσ^*, 4 n3s(C), and 12 n3p(C) states], with equal weights. To properly account for Rydberg states, we used the mixed d-aug-cc-pVDZ(C)/aug-cc-pVDZ(F,Cl) basis set.^84,85,87

Full geometry optimization of the ion-pair state has been performed at the MR-CISD/aug-cc-pVDZ level using only four valence orbitals in the active space, that is, the σ_C–Cl/σ_C–Cl^* pair of the longer bond and two lone pairs of the dissociating Cl. A single-pointMR-CISD/aug-cc-pVDZ calculation including the eight valence orbitals in the CAS space, denoted MR-CISD(12,8), has been performed. In both cases, 3A′ + 1A″ states have been averaged at the MCSCF level.

As can be seen from Fig. 5, a cascade of nonadiabatic transitions can deactivate higher states, until the ion-pair structure (in 3¹A′) is reached, thus explaining the Cl⁻ yields in the lowest energy band of the spectrum from Ref. 83. As in other systems, the ion-pair state curve is usually well separated from the others at relatively large C–Cl distances.^88–90 

[CF₂Cl]^+δCl^−δ is bonded by 3.00 eV at the MR-CISD + Q(12,8)/aug-cc-pVDZ level, 1.00 eV lower than [CF₃]^+δCl^−δ,⁹⁰ which can be due to the longer C–Cl distance (by 0.193 Å) of the former. Mulliken charges (δ) and dipole moment are 0.60e and 7.51 D, respectively, at the MR-CISD(12,8) level.

Accurate potential energy surfaces are required for many problems in chemistry, in particular, for spectroscopy. High accuracy also means high computational expense; therefore, such calculations, even with efficient implementations as in COLUMBUS, can be utilized only for small molecules.⁹¹ Here, we show, through an application on ozone dissociation, how some unique features of COLUMBUS allow accurate treatment of complicated situations.

Accurate quantum chemical results require the consideration of not only dynamic but often also static correlation.¹³ Despite significant previous efforts,⁹² MR versions of coupled cluster (CC) methods are not yet routinely available. Therefore, MRCI methods seem to be one of the best choices, in particular, at the MR-CISD and MR-AQCC levels.¹³ 

One important shortcoming of MR-CISD is the linearly increasing size of the parameter space with the number of reference functions. If a CAS^93–95 reference is used, the cost increases exponentially with the number of active orbitals involved. To avoid this explosion in computational cost, often internally contracted (ic) functions are used^96–99 (see Ref. 13 for a detailed review). This procedure is very efficient and has been applied successfully in many cases. However, ic methods have the disadvantage that the treatment of static and dynamic correlation appears in separated steps of the calculation. Therefore, if strong interaction with another state is present and there are (avoided) crossings between different potential energy surfaces, the separation of the static and dynamic correlation is problematic because the crossing happens at different geometries for the reference and the subsequent MR-CISD calculations.¹⁰⁰ 

Such a problem is largely solved if, instead of internal contraction, an uc MRCI wavefunction is used with all reference functions individually included in the ansatz. The advantage of the uc vs ic calculations to avoid unphysical reef-like structures on dissociation curves has been shown in Refs. 101 and 102 and will be demonstrated with an example below. We note that flexible selection of reference functions in COLUMBUS allows a reduction of the cost compared to the full CAS reference calculations; it is possible to use one wavefunction in the MCSCF step to obtain orbitals and then use a different set of reference configurations in the correlated MR-CISD calculation³ for both single points and energy gradients.⁶ 

Contrary to CC methods, MRCI methods are not size-extensive,¹⁰³ i.e., the energy does not scale properly with the system size,^77,104 an error which certainly needs to be corrected in high accuracy calculations. There are several possibilities, both a priori and a posteriori, for this correction, which are summarized in Refs. 13 and 103. The most often used a posteriori correction is the Davidson (MR-CISD-Q_D) approach,⁷⁷ but we usually suggest its Pople version (MR-CISD-Q_P).¹⁰⁵ A priori corrections are more advantageous since not only the energy but also the wavefunction is corrected, and also analytic gradients are available. The two most popular versions are MR-ACPF¹⁰⁶ and MR-AQCC;^14,15 the latter seems to be more stable.¹⁰³ 

We show the importance of both the uncontracted ansatz and the size-extensivity correction on the accurate potential energy surface of ozone. Here, the potential along the minimum energy path (MEP) leading to the dissociation of one O atom is important for accurate prediction of highly excited vibration levels as well as to describe the scattering of an O atom by an O₂ molecule.¹⁰⁷ The latter process shows an unusual isotope effect,^108,109 which is difficult to explain theoretically. Theoretical methods often predict a “reef-like” structure with a small barrier and a van-der-Waals minimum along the MEP (see, e.g., Refs. 107 and 110 for reviews). In the case of ozone, Holka et al.¹⁰⁷ showed the significant effect of size-extensivity correction on the size of the barrier, while Dawes et al.¹¹⁰ demonstrated that including several internally contracted reference functions in a multistate ic-MRCI calculation causes the reef to disappear. Figure 6 shows how the effects mentioned above, i.e., uncontraction of the reference space (left panel) and the inclusion of size-extensivity corrections in the form of the Davidson (Q_D) and Pople (Q_P) corrections, as well as by MR-AQCC (right panel), lower the barrier and lead essentially at the uc-MR-AQCC level to its disappearance. As discussed by Tyuterev et al.^111,112 and Powell et al.,^102,113 only the barrierless potential is capable of reproducing the experimental findings both for the vibrational levels and the scattering.

Spin-orbit MRCI calculations in COLUMBUS are based on spin-averaged molecular orbitals (i.e., the polarization of spinors is recovered at the CI level) and an effective one-electron spin-orbit (SO) operator scheme.⁵ Scalar relativistic effects enter through modified one-electron integrals. By a suitable definition of the spin functions, the Hamiltonian is real; the generally complex odd-electron case is embedded within an artificial, real N + 1 electron case, doubling the size of the Hamiltonian matrix. While spin-orbit relativistic effective core potentials (RECPs), such as the RECPs of Cologne and Stuttgart^115–117 and Ermler et al.,^118,119 are natively supported, eXact-2-Component (X2C)¹²⁰ and arbitrary order Douglas–Kroll–Hess (DKH)¹²¹ for scalar relativistic contributions along with the Atomic Mean-Field approximation (AMFI)¹²² for the spin-orbit interaction are accessible via the COLUMBUS/MOLCAS interface.¹²³ For a review on spin-orbit coupling, cf. Ref. 124.

COLUMBUS supports variational, uncontracted spin-orbit MRCI calculations treating electron correlation and spin-orbit coupling on the same footing (one-step method).⁵ Since this procedure expands the wavefunction in terms of CSFs with multiple spin multiplicities including all (2S + 1) components of each, the CSF space of SO-MRCI is several times larger than that of a conventional MRCI calculation, which requires only a single component of a single spin multiplicity. The flexible MRCI paradigm and an effective parallelization scheme adapted to current computer architecture enable the application of COLUMBUS to heavy element science, e.g., the spectroscopy of lanthanides and actinides.

The chemistry of the early actinide (An) elements features high oxidation state species such as the linear actinyl ions, OAnO⁺ (V) and OAnO²⁺ (VI), which exist for uranium through americium.¹²⁵ The actinyl ions are well studied experimentally, and significant contributions by theory and computation have been made.^126–129 Actinide containing systems are prime candidates for computational approaches due to their radioactivity and short lifetimes, especially for the latter members. The methods must treat relativistic effects and nondynamical correlation. Herein, the SO coupling in the actinyls was explored using the one-step, variational, uncontracted, RECP-based two-component formalism, wherein SO and electron correlation are computed simultaneously and treated equally.

The archetypal uranyl (VI) ion is a closed shell ion. The uranium 6s and 6p orbitals in combination with the oxygen 2s orbital form the inner valence MOs containing 12 electrons, as shown in Fig. 7. Another 12 electrons fill the bonding MOs formed from the uranium 5f and 6d orbitals and the oxygen 2p orbital. As the actinyl series progresses, the additional electrons occupy the nonbonding 1δ_u and 1ϕ_u orbitals, stemming from the 5f manifold, producing a multitude of low-lying electronic states and culminating for AmO₂⁺ in a high spin ground state of ⁵Σ⁺_0+g from the 1δ²_u 1ϕ²_u configuration. Beginning with americium, the chemistry of the actinides becomes more lanthanide-like and the III oxidation state dominates.

Quasi-actinyls, further potential members of the actinyl series, are under investigation. An interesting question is whether the known weak-field coupling in the 1δ_u, 1ϕ_u subspace and significant antibonding character of the 3π_u orbitals in the actinyls will continue into the quasi-actinyls to yield low spin states. Compact correlation-consistent double zeta plus polarization basis sets developed for use with RECPs and SO operators were employed.¹³⁰ Large reference spaces consisting of the fully occupied 1π_u, 2π_u, 3σ_u MOs and the partially occupied 1δ_u, 1ϕ_u, 3π_u nonbonding MOs were used in SO-MR-CISD calculations. Table II lists the ground occupations and states. CmO₂²⁺ is isoelectronic with AmO₂⁺ as are the successive pairs shown in Table II. High spin results are obtained, suggesting that the 3π_u orbital is predominantly 5f and confirming the lanthanoid nature of these actinide elements. The non-octet ground state of CfO₂⁺ is anomalous and may be due to an inadequate reference space.

One of the first tasks in the one-step, variational, uncontracted, two-component formalism is to develop Gaussian basis sets for use with the RECPs, a key element of which is the COLUMBUS version of the atomic self-consistent-field program by the basis set expansion method, known as ATMSCF.¹³¹ The ATMSCF program is a modernized and enhanced version of the Chicago atomic self-consistent-field (Hartree–Fock) program of 1963.¹³² Energy-expression coefficients now treat the ground states of all atoms to the extent that Russell–Saunders (LS) coupling applies. Excited states with large angular-momentum orbitals can be handled. Relativistic effects can be included to the extent possible with RECPs.

A common problem in basis set exponent optimization is exponent collapse, where two exponents approach each other very closely and their corresponding coefficients become very large in magnitude with opposite signs. The current code manages this problem by expressing the natural logarithms of all the exponents for each l-value as a series of Legendre polynomials and then constraining the number of independent coefficients.¹³³ 

Employing this software, basis sets have been developed for various elements (e.g., Refs. 134, 130, and 135). An example application of basis sets developed using ATMSCF is the electronic structure and spectra of actinyl ions.¹²⁹ 

Nonadiabatic nuclear dynamics requires the electronic structure data (ESD), energies, energy gradients, derivative couplings, and, for more sophisticated calculations, dipole and transition dipole moments and the spin-orbit interaction. There are two ways to present these data: on-the-fly¹³⁶ (as discussed in Sec. III A) or as coupled diabatic representations fit to functional forms, most recently neural network forms.^137,138 In the on-the-fly approach, the wavefunctions are determined when and where needed so that the above-noted quantities can always be calculated. The alternative approach, the fit surface approach, uses quasi-diabatic functional forms fitted to reliably reproduce adiabatic ESD. The strengths of these approaches are complementary. In the on-the-fly approach, ESD is always available since the electronic wavefunctions are determined at each nuclear geometry sampled. For this approach, high accuracy electronic structure techniques cannot be used as their single point evaluation is too costly. Fitted surface techniques, on the other hand, precalculate and represent as functional forms the energy, energy gradients, and derivative couplings with the locus of points at which the ESD is determined and fit, guided by the regions sampled by surface hopping trajectories.¹³⁹ However, fit surface techniques do not usually include interactions with an electric field or the spin-orbit interaction, which can be a significant deficiency. Either way, the electronic structure method employed in the calculations must be able to describe large sections of the configurational space of the nuclei, including multireference regions. Such a requirement makes the MR method in COLUMBUS ideal to deal with these problems. The combination of both approaches, fitting and on-the-fly dynamics, is also possible, using ESD of high accuracy, as, e.g., obtained with COLUMBUS, to fit diabatic model potentials,¹⁴⁰ also based on the vibronic coupling models¹⁶⁵ or machine learning potentials,¹⁴¹ and to run subsequent on-the-fly trajectories (see also Secs. III B 4 and III C).

The electronic structure methods included in COLUMBUS can be used for running nonadiabatic dynamics simulations with mixed quantum-classical (NA-MQC) approaches¹⁴² and simulating spectra (absorption and emission) based on the nuclear ensemble approach.¹⁴³ The NEWTON-X^136,144,145 and SHARC^146,147 programs specialize in these types of dynamics and spectrum simulations and have dedicated interfaces to COLUMBUS. NEWTON-X, in particular, was first developed to work with COLUMBUS and only later gained interfaces to other quantum-chemical programs.

Both NEWTON-X and SHARC perform dynamics with variants of the fewest switches surface hopping.¹⁴⁸ In this method, potential energy surfaces, their gradients, and the couplings between the states do not need to be known a priori. These quantities are calculated on-the-fly at the nuclear position defined by classical trajectories. Thus, at every time step of integration of the Newton equations, the dynamics program automatically invokes COLUMBUS to calculate the required electronic quantities. After COLUMBUS execution, the dynamics program automatically extracts the electronic information and uses that data to continue the propagation. Alternatively, the wavefunctions produced by COLUMBUS can be used to compute time-dependent wavefunction overlaps,^149,150 which are used either in Hammes-Schiffer–Tully time-derivative couplings¹⁵¹ or in local diabatization procedures.¹⁵² 

In the calculation of spectra based on nuclear ensembles, the procedure is similar. Having a collection of nuclear geometries, for instance, generated from a Wigner distribution¹⁵³ for the quantum harmonic oscillator or even extracted from a previously run dynamics, the third-party program invokes COLUMBUS to calculate transition energies and transition dipole moments for each of these geometries. The spectrum is then generated by postprocessing of these results.¹⁴³ 

Many features make COLUMBUS attractive for these types of simulations: first, the availability of analytical potential energy gradients and analytical nonadiabatic (and spin-orbit) couplings at the MCSCF and MRCI levels. Having such analytical quantities turns out to be essential for the quality of results and computational efficiency; this happens because dynamics simulations based on numerical gradients tend to be unstable due to the switching of the dominant diabatic character in the adiabatic states near state crossings. Moreover, the computational cost of numerical gradients becomes prohibitively large, even for medium-sized molecules.

A second feature making COLUMBUS a well-suited platform for on-the-fly dynamics is its extreme flexibility to build active and reference spaces. Dynamics based on MCSCF are known for being unstable due to the frequent switching of orbitals in these spaces.¹⁵⁴ For example, a CAS(2,2) initial set of π and π^* orbitals may switch into a π and σ^* set as soon as a molecular bond is stretched. Such orbital switching causes discontinuities in the total energy, ruining the trajectory. Dynamics simulations based on post-MSCF methods inherit the same problems from their MCSCF step. The brute-force solution for such a problem, to increase the active and reference spaces until they contain all orbitals needed for dynamics, is usually computationally prohibitive. With COLUMBUS, however, it is possible to include these additional orbitals in separated subspaces, which may be completely isolated or interact with each other via single excitations. This approach keeps the number of generated configurations manageable. In Sec. III E, we discuss an example of such specialized spaces in dynamics simulations.

COLUMBUS has been used for dynamics simulations of many molecules. Some illustrative case studies using the COLUMBUS/NEWTON-X interface are the dynamics of adenine in the gas phase at the MR-CIS level¹⁵⁵ and the dynamics of azomethane in solution at the generalized valence bond perfect-pairing (GVB-PP) MCSCF level with molecular mechanics to describe the solvent using a quantum mechanical/molecular mechanical (QM/MM) approach.¹⁵⁶ 

The COLUMBUS/SHARC interface allows the treatment of arbitrary spin-multiplicities and their interactions via spin-orbit and nonadiabatic couplings on equal footings. This approach was used to study the photodissociation dynamics of SO₂ at the MR-CIS level that reproduced the results of simulations using exact quantum dynamics.¹⁵⁷ As expected from selection rules based on symmetry, only one out of three triplet states is populated, concurring with population flow between the ¹B₁ and ¹A₂ singlet states.

Recently, the COLUMBUS/NEWTON-X interface has been explored to simulate the dynamics of transient anions.¹⁵⁸ This new approach has been developed to describe the attachment of a low-energy electron to a molecule and the subsequent dynamics, including the possibility of electron autodetachment. The use of multiconfigurational methods in such a problem is mandatory due to the dissociative character of the temporary anions.

The direct GUGA based configuration interaction (CI) in COLUMBUS provides an important advantage for fitting/diabatizing ESD, its ability to determine the ab initio (first) derivative coupling

(and for that matter energy gradients) using highly efficient analytic gradient techniques.^7,8 As a consequence, the key vector requirement for diabatization is readily formulated. Here, $ΨJa,(ab)$ is the ab initio (ab), adiabatic (a) wavefunction from COLUMBUS. R, (q) are the nuclear (electronic) coordinates. The requirement is that the energy-difference scaled ab initio derivative coupling

equals, in a least squares sense, the energy difference scaled derivative coupling

obtained from the fit diabatic electronic Schrödinger equation

The algorithm we employ^139,159 to simultaneously fit and diabatize high quality ab initio ESD is outlined below. It is based on the following form for H^d, the (fit) diabatic potential energy matrix,

Here, B is an N^state × N^state (N^state is small) symmetric matrix with 1 in the (u, v) and (v, u) elements, g(R) are the fitting functions, and P is a projector in a subgroup of a complete nuclear permutation inversion group. The linear parameters V_l provide a least-squares fit of the ab initio determined adiabatic energies, energy gradients, and energy-difference scaled derivative couplings to those obtained from Eqs. (5) and (6). The algorithm obtains the array of points Rⁿ at which to perform the least-squares fit from quasi-classical surface hopping trajectories based on Tully’s fewest switches method.¹⁶⁰ 

Figure 8 illustrates the potential of this method to diabatize ab initio ESD and describe the vicinity of conical intersections and reports earlier^159,161 results from a full 33-dimensional, 4 diabatic state representation of the S₀, S₁, and S₂ states of phenol, C₆H₅OH. The use of analytic gradient techniques and the inclusion of the derivative coupling make it possible to provide the diabatic representation in its full dimensionality even for systems this large. Depicted in Fig. 8(a) are the adiabatic potential energy curves obtained ab initio and from H^d. The diabatic potentials are also given. Figure 8(b) depicts the derivative coupling along a linear synchronous transit path passing through two conical intersections. The path changes at r(OH) ∼ 1.25 Å, producing the observed discontinuity.

The fit coupled diabatic state approach is particularly well suited for quantum nuclear dynamics and for nonadiabatic processes where high accuracy is needed. Laser control of the electronic state produced in a nonadiabatic process has long been a holy grail of physical chemistry. The photodissociation of ammonia, NH₃(⁠$X̃$⁠, v_i) + hϖ → NH₃(⁠$Ã$⁠) → NH₂(⁠$X̃$⁠, $Ã$⁠), was thought to provide such an example¹⁴⁰ where the excitation of v₃ preferentially produced NH₂(⁠$Ã$⁠) by manipulating the near conical intersection dynamics. Using an accurate H^d,¹⁶² obtained by carrying out the required electronic structure calculations with COLUMBUS, and full six dimensional quantum dynamics, we convincingly demonstrated that this was not the case.¹⁶³ 

An option to avoid the computational cost of on-the-fly computations (Sec. III A) and the algorithmic challenges of a diabatic potential fit (Sec. III B) amounts to the construction of a vibronic coupling model, using only minimal electronic structure data. A vibronic coupling model¹⁶⁴ operates by creating a Taylor series expansion of the diabatic Hamiltonian matrix. The diagonal terms in the linear vibronic coupling (LVC) matrix derive from the gradients of the different excited states, while the off-diagonal terms are related to the nonadiabatic couplings. Spin-orbit terms can be naturally added as constants in a diabatic picture.^164–167 COLUMBUS provides all three types of terms analytically and, thus, allows us to parameterize an entire LVC model using only one single-point computation.

In Ref. 167, it was examined whether LVC models constructed in this way can indeed be used to obtain a rough description of the dynamics to be expected. Linear spin-vibronic coupling models were constructed for five molecules: SO₂, adenine, 2-aminopurine, 2-thiocytosine, and 5-azacytosine; these were subsequently used for surface-hopping dynamics. The correct qualitative time-resolvedpicture was obtained for four of the molecules, while only in the case of 5-azacytosine was an inconsistency observed as the decay to the closed-shell ground state could not be reproduced. These results suggest that the combination of vibronic coupling models with surface hopping dynamics could prove to be a handy tool for obtaining a rough estimate of the photodynamical behavior of a molecule. The approach can naturally be extended to a quadratic vibronic coupling model using the strategy described in Ref. 168.

The availability of analytic gradients and nonadiabatic couplings at the MRCI level in COLUMBUS, in addition to the interface with NEWTON-X, provides a unique opportunity to do nonadiabatic dynamics of complicated systems including dynamical correlation for excited states at a multireference level. An additional advantage of COLUMBUS is the possibility of varying the level of electronic structure theory [by comparing CASSCF and MRCI] to determine how the results depend on the methodological approach. Uracil is an attractive test system for such comparison since static electronic structure calculations have shown a great dependence of the barrier on the level of theory used to describe the excited state.^169,170

Surface hopping dynamics on all nucleobases using COLUMBUS were performed several years ago¹⁷¹ providing significant progress in our understanding of how the bases dissipate energy. More recently, nonadiabatic dynamics were compared to pump–probe experiments, which provide the ultimate test of the accuracy of the results.¹⁷² Trajectory surface hopping dynamics on uracil were performed using CASSCF and MR-CIS.¹⁷² The aim was to see how the two different electronic structure methods affect the excited state dynamics and compare to pump–probe strong-field and weak field ionization experiments by Horton and co-workers.¹⁷² An absorption spectrum was first calculated using a Wigner distribution of the ground state, and geometries with excitation energies within the experimental pump pulse window were selected for excited-state dynamics. The combination of COLUMBUS/NEWTON-X provides all tools for initiating, carrying out, and analyzing the dynamics. Figure 9 (taken from Ref. 172) shows the theoretical and experimental results superimposed. The probe is a vacuum ultraviolet (VUV) pulse (156 nm, 7.95 eV), which can ionize the molecule when the population is on the excited states S₁ and S₂ of uracil, but the signal (ion yield) disappears when the population decays to the ground state, so it provides an experimental verification of the time constants for radiationless decay to the ground state. For comparison with the experiment, the calculation sums the S₁ and S₂ populations as a function of time. Both the CASSCF and MR-CIS results agree well with experiment, which shows that significant population remains on the excited states even after 1 ps. Interestingly, however, the detailed dynamics predicted at the CASSCF and MR-CIS levels differ with each other. CASSCF predicts that significant population will remain on S₂, while MR-CIS predicts that the population will remain on S₁, but since the signal comes from both states, it is hard to distinguish them experimentally in this case.

Other recent relevant surface hopping calculations using NEWTON-X/COLUMBUS were performed on radical cations of uracil, cyclohexadiene, and hexatriene.¹⁷³ These studies were very successful in extracting insight into the role of the derivative coupling in the decay rates since it was shown that, in addition to the magnitude of the derivative coupling, the direction is also crucial. Again, the flexibility of COLUMBUS was instrumental in performing these calculations since specific RASSCF wavefunctions had to be used in order to obtain the best description of the excited states.

COLUMBUS is extremely flexible in building advanced configurational and reference spaces composed of uncorrelated or inter-correlated sub-spaces. In this section, we exemplify such features through a case study, the nonadiabatic dynamics of ethylene.¹⁷⁴ 

It is a well-established experimental fact that, after UV excitation into the V state (ππ^*), ethylene returns to the ground state via internal conversion within a few tens of femtoseconds.^175–177 Although dynamics simulations can quantitatively predict this ultrafast deactivation,^{174,178–184} there is a quantitative disagreement between experiments and diverse theoretical models, which has been hard to explain fully.

Distinct reasons lead to this disagreement. Some of the reasons relate to experimental effects such as the detection window setting in time-resolved spectroscopy. Some of the reasons are computational in nature mainly associated with the electronic structure level employed. Briefly, there are three main challenges in the computational simulations of ethylene. First, the description of the V state, which requires a re-optimization of the σ orbitals after the inclusion of dynamic electron correlation in the π system,^185,186 which is not possible within the standard MCSCF/MRCI procedure. As a consequence, a conventional CASSCF calculation of the V state overestimates its potential energy by more than 1 eV, and extended CI expansions are necessary to compensate.¹⁸⁷ Second, the description of the dynamics of ethylene is further complicated by the many Rydberg states lying below the V state, which naturally requires large orbital spaces and diffuse orbital basis sets. Finally, the significant amount of photoenergy released into the vibrational modes quickly leads to large CH stretching and even CH dissociation, requiring the inclusion of σ^* orbitals into the configurational space. Moreover, as discussed in Sec. III A, on-the-fly dynamics simulations based on CASSCF and post-CASSCF methods suffer from instabilities caused by orbital rotations, changing the character of the orbitals in the configurational and reference spaces, causing discontinuities in the total energy.

Thus, to deal with all these issues, on-the-fly nonadiabatic dynamics simulations of ethylene require large valence-Rydberg orbital spaces, with all valence electrons, and diffuse basis sets. Such simulations based on standard configurational and reference spaces would be too costly. Nevertheless, a chemically motivated ensemble of disjointed sub-spaces can affordably help alleviate these problems.

In Ref. 174, the nonadiabatic dynamics of ethylene was simulated with surface hopping at the MR-CISD level, including all valence electrons, representing valence and Rydberg states, and using the jun-cc-pVDZ basis set.¹⁸⁸ The configurational space included the σ_CC, four σ_CH, π, π^*, four σ_CC^*, and four Rydberg molecular orbitals.

The configurational space for the MCSCF calculations was built as follows:

In this symbolic expansion, the perfect pairing (PP) for CX (PP_CX) within the framework of GVB^189,190 is a subspace containing a σ_CX and a σ_CX^* orbital, allowing only doubly occupied configurations. CAS(2,2)_π is a subspace containing the π and π^* orbitals allowing, as usual, singly and doubly occupied configurations. Single excitations from CAS(2,2)_π into the four Rydberg orbitals (Ryd) are allowed. Single and double excitations between CAS(2,2)_π and PP_CH are also allowed.

The reference space for the MRCI calculations was built as

where each CAS(4,4)_πσ(CH) subspace contains the π and π^*, and the σ_CH and σ_CH^* orbitals for each hydrogen atom.

Although there are still open questions concerning the ultrafast deactivation of ethylene, dynamics simulations based on these involved spaces appreciably improved the computational results, bringing the excited-state lifetime closer to the experimental findings.

As of the time of the submission of this paper, several exciting program and method developments are in progress, which are described in Subsections IV A–IV D. These examples are included in this report to demonstrate upcoming new features in COLUMBUS. Only the local correlation (LC) treatment of Sec. IV D is routinely available in the current distribution version of COLUMBUS. Developmental versions of ACME MCSCF and GCF (Sec. IV A) can be made available to interested users on request.

The conventional MCSCF method is limited to about n ≲ 16 active electrons and active orbitals because the Hamiltonian diagonalization and reduced density matrix (RDM) computation effort increases dramatically with the increase in active orbital dimension n (e.g., as $nNe$ for full-CI type expansions with N_e electrons). A new orbital optimization approach¹⁹¹ is being developed that eliminates this restriction. Within each iteration of a conventional MCSCF approach, the symmetric eigenvalue equation, HV = VE, is solved in the configuration state function (CSF) basis of dimension N_CSF for the state of interest E_k or in a state-averaged (SA) calculation for the weighted sum of several N_av states of interest, $Ē=∑k=1NavwkEk$⁠. Consider first the special case: N_av = N_CSF and w_k = 1/N_CSF ∀ k. These conditions, combined with the trace identities, Tr(E) = Tr(V^THV) = Tr(HVV^T) = Tr(H), result in

This SA energy does not depend on the wavefunction expansion coefficients V_jk, and it can be computed using only the diagonal H_kk matrix elements,

These H_kk elements depend only on the small subset of Hamiltonian integrals h_pp, g_pppp, g_ppqq, and g_pqpq and require only O(n²) effort each to compute. This special case of state averaging is called the All Configuration Mean Energy (ACME) conditions with equal weights. In principle, this would allow $Ē$ to be computed with O(N_CSFn²) effort. The COLUMBUS MCSCF code is based on GUGA, in which the CSF expansion space is represented as walks within a Shavitt graph. A recursive procedure based on this graphical representation allows $Ē$ to be computed instead with only O(ωn²) effort, where ω is a factor that ranges from $ONe0$ up to $ONe2$⁠, depending on the complexity of the Shavitt graph. Since ω ≪ N_CSF, the effort for this recursive procedure is independent of N_CSF. The sparse ACME RDM, consisting of unique nonzero elements $D¯pp$⁠, $d¯pppp$⁠, $d¯ppqq$⁠, and $d¯pqqp$⁠, can also be computed recursively and requires a comparable amount of effort.

Given the RDM elements, the orbital optimization gradient and hessian can be constructed, and efficient first- or second-order convergent methods may be employed to optimize the orbitals. The ACME RDM sparsity may be exploited to reduce the computational effort for this optimization compared to the typical state-specific optimization. Due to the favorable effort scaling of the ACME algorithms and to the elimination of the H eigenvalue equation, essentially an unlimited number of active orbitals can be accommodated. Timings are given in Ref. 191 for up to 256 active orbitals and up to N_CSF ≈ 10¹⁵⁰. The most significant remaining computational effort in each iteration is, depending on the MCSCF implementation, either the four-index AO to MO transformation of the integrals or the MO to AO transformation of the ACME RDM. There is no need to artificially limit the number of active orbitals or to artificially restrict valence orbitals to be doubly occupied as is typically required in traditional MCSCF implementations. The price for this flexibility is that the state-specific MCSCF wavefunctions and energies are not available during the orbital optimization procedure.

Given the ability to optimize orbitals with the ACME conditions, subsequent state-specific high-level electronic structure calculations can be performed using these orbitals. Analytic geometry gradients can be computed for the high-level methods that use these orbitals. This analytic gradient procedure is based on the efficient and flexible successive orbital transformation formulation that has been previously developed for MRCI wavefunctions.²⁰ As in the MCSCF optimization step itself, the individual MCSCF state-specific wavefunctions and energies are not required or referenced; only the inexpensive ACME RDMs are required in this procedure. The successive orbital transformation formulation can also be applied to the efficient computation of nonadiabatic coupling between individual states computed with high-level electronic structure methods.^7,8 Future efforts will focus on the elimination of the equal-weight condition while maintaining the same computational effort.

A novel expansion basis for electronic wavefunctions (see Ref. 192 and the references therein) is being developed within the COLUMBUS Program System. In this approach, the wavefunction is written as a linear combination of graphically contracted functions (GCF), and each GCF, in turn, is formally equivalent to a linear combination of CSFs that comprise an underlying full-CI linear expansion space of dimension N_CSF. The CSF coefficients that define the GCFs are nonlinear functions of arc factors, which themselves depend on a smaller number of essential variables, N_φ = N_CSF; these essential variables are optimized in order to minimize either a ground or excited state E_k or a state-averaged energy $Ē$⁠.

The initial implementation of the GCF method relied on the nonlinear basis dimension N_GCF to extend the wavefunction flexibility and to converge molecular properties toward the full-CI limit. An alternative, and potentially more efficient, approach to enhance the wavefunction flexibility has been implemented that consists of allowing multiple partially contracted wavefunctions to be associated with each Shavitt graph node within each GCF. The initial approach is now called the single-facet GCF (SFGCF) method, and this new approach is called the multifacet GCF (MFGCF) method. An MFGCF basis function is a matrix product state (MPS), and the ground- and excited-state wavefunctions are linear combinations of these MPSs. Several properties and algorithms previously developed for the SFGCF method have been implemented within the MFGCF formulation: state-averaging, Hamiltonian matrix construction, RDM construction, Slater determinant overlaps, graph density computation and display, and spin-density matrix computation.

MFGCF expansions with facet counts in the range f_MAX ≈ 4–10 have been shown to approach the full-CI potential energy surfaces (PESs) to within chemical accuracy (1 kcal/mol or better). The GCF method scales much better with orbital basis function dimension and with the number of electrons than traditional high-level electronic structure methods; e.g., an energy expectation value may be computed with O(⁠$NGCF2$ωn⁴) computational effort. No intrinsic restrictions are imposed on the orbital occupations, and there are no artificial excitation-level or occupation restrictions with respect to a reference function or reference space; in this sense, the method is more correctly characterized as a multiconfigurational method rather than a multireference method. Because the wavefunction is a linear combination of GCF basis functions rather than a single expansion term, the method may be used for both ground and excited electronic states, the increased wavefunction flexibility leads to higher accuracy, and this expansion form facilitates the computation of transition moments, nonadiabatic coupling, and other properties that at present can only be computed reliably with multireference approaches.

The ongoing developments with the GCF method include increasing the wavefunction flexibility and robustness by allowing for symmetry-dependent arc factors, exploring both fixed-facet and variable-facet optimization approaches, incorporating multiheaded Shavitt graphs to describe simultaneously molecular states with various charges and spin multiplicities, allowing state-averaging over multiple point group symmetry irreducible representations, and incorporating spin-orbit coupling. One goal is to combine accurate, high-level, state-specific, GCF wavefunctions with the molecular orbitals from the ACME MCSCF method to compute properties of large molecular systems.

The spin-density matrix with elements $Dqp(1,0;M)=ψ;S,Mapα†aqα−apβ†aqβψ;S,M$ is usually computed in a Slater-determinant approach as the difference of the α and β blocks of the 1-RDM. Since COLUMBUS employs the GUGA formalism,¹⁹³ this spin-orbital information is normally not available. A straightforward procedure for the spin density computation would be to transform the wavefunction to a determinantal basis and then compute the spin-density matrix in that representation. However, this may be inefficient and even unfeasible, for large CSF expansions, and an alternative approach is indicated.

It was shown¹⁹⁴ that the maximum, M = S, spin-density matrix can be calculated as

where $Dqp(0)$ is a 1-RDM element and $dqkkp=ψêpkkqψ$ are the (unsymmetrized) 2-RDM elements.^192,194 These two quantities are independent of M and are available within the GUGA formalism as part of the standard computation of the 1-RDM and 2-RDM. Once the spin-density matrix is obtained for the M = S case, the Wigner–Eckart relation^195,196$Dqp(1,0;M)=MSDqp(1,0;S)$ can be used to obtain the spin-density matrix for any of the 2S + 1 members of the degenerate multiplet (M = −S, −S + 1, …, S).

The second term on the right-hand side of Eq. (11) is a sum over a subset of the elements of a four-dimensional array, and a loop over these matrix elements should be as optimized as possible. This would require O(n³) effort. Alternatively, it is convenient to compute the spin-density matrix alongside the 1-RDM and 2-RDM, in which case the additional effort is minimal.

The spin-density matrix calculation is being implemented within the MCSCF and MRCI codes in COLUMBUS and will be available in the next stable version release. It is our hope that this implementation can be used to shed some light on problems such as how the unpaired electrons of graphene nanoflakes (Sec. II D) may give rise to spintronics applications⁷³ of these systems.

The application of the Born–Oppenheimer approximation to separate the nuclear and electronic Hamiltonians is insufficient to model many systems appropriately. We have already discussed several examples in Sec. III. Another example of the limitations of the Born–Oppenheimer approximation occurs in collisionally induced changes to the states split by the spin-orbit interaction. Such phenomena are of importance in calculating the dynamics of many systems, including Diode-Pumped Alkali Lasers (DPALs), in which alkali metals interact with noble gases to produce laser radiation.^197–201 

The coupling between electronic and nuclear states is quantified by the nonadiabatic derivative coupling terms (DCTs), which have the form

and the resultant off-diagonal energy coupling surface.^202,203 Analytic gradients, of interest in geometry optimization problems, are calculated in a similar fashion,

and thus, the methods to calculate one analytically parallels the other. The analytic calculation of DCTs and gradients, which can be performed at a single geometry R, is faster than calculation via finite difference methods, which require at least 3N_atom + 1 geometries for the evaluation (where N_atom is the number of nuclei).

COLUMBUS implements analytic DCT calculations for non-relativistic MRCI wavefunctions via the GUGA approach¹⁸ using the formalism of Lischka et al.⁷ and Dallos et al.,⁸ which is based upon the methodology of calculating analytic gradients of MRCI wavefunctions by Shepard.²⁰ By broadening those methods to include the treatment of relativistic MRCI wavefunctions of Yabushita et al.,⁵ which includes the explicit consideration of the single electron spin-orbit operator and the implicit consideration of other relativistic effects in the Relativistic Effective Core Potential (RECP),²⁰⁴ we are now able to calculate analytic DCTs and gradients for relativistic MRCI wavefunctions.

As an example, consider the radial derivative coupling in the LiHe system. This coupling arises due to the change in distance between the nuclei between the $A2Π1/2$ and $B2Σ1/2$ potential energy surfaces (PESs). Ordinarily, these surfaces are each doubly-degenerate for the m_j = ±$12$ states. The Yabushita implementation, a convenient way to utilize the real arithmetic infrastructure already present in COLUMBUS, requires the addition of an artificial, non-interacting electron for odd electron systems.⁵ This constraint turns the double degeneracy of each set of surfaces into a quadruple degeneracy. As with any degeneracy, the four electronic wavefunctions that COLUMBUS converges upon for either surface will be an arbitrary mixture of the four pure canonical states which do not mix the z-component of total angular momentum, m_j, or the z-component of the ghost electron’s spin, m_{s ghost},

(with an equivalent relationship for the $Π1/2i$states). Here, the unprimed wavefunctions are the arbitrary mixtures of m_j and m_{s ghost} (with 1 ≤ i ≤ 4), and the primed wavefunctions are the pure canonical states, each j representing a specific pairing of m_j and m_{s ghost}. Using the expansion in Eq. (14), we have shown that²⁰⁵ 

that is, the norm of the canonical DCT is the root sum square of the mixed DCTs involving any single electronic wavefunction $Σ1/2i$and all electronic wavefunctions $Π1/2j$ (the argument is equally valid for the DCTs of a single wavefunction $Π1/2i$ with all $Σ1/2j$⁠).

From these DCTs, it is possible to calculate the coupling angle between the $A2Π1/2$and $B2Σ1/2$ PESs.⁷ That angle is used to calculate an off-diagonal coupling surface and the correct mixing of the adiabatic surfaces to produce the nonadiabatic surfaces. Figure 10 shows these MRCI surfaces for LiHe as calculated with a Stuttgart basis.^70,206

This methodology has been implemented in an experimental version of COLUMBUS at the Department of Defense Shared Resource Center (DSRC). Due to the current program structure, COLUMBUS is not natively capable of producing the RECP gradients including the spin-orbit core potential gradient integrals within the same program. For this implementation, the program system NWChem²⁰⁷ was leveraged to produce these terms in preparation for the COLUMBUS calculations.

When these calculations were performed with KHe, the results compared favorably with approximation methods,²⁰⁵ indicating that integration into a future production version of COLUMBUS is feasible.

Due to their high computational cost, MR methods are rarely used to predict the electronic structures and reaction energetics of large molecules. Multireference local correlation (LC) treatment is an appealing approach to reduce efficiently the computational cost as molecular size increases. The strategy developed for COLUMBUS consisted of localizing only the reference occupied orbital space,²⁰⁸ using the concepts of the weak pair (WP) approximation of Sæbø and Pulay^209,210 and a geometrical analysis of Carter and co-workers^211–213 developed for their TIGERCI code (see Refs. 214–216 for the code description and characteristic application). This approach reduced the number of configurations significantly while keeping the active space unchanged, thereby simplifying comparison with standard calculations. More details beyond the summary of the selection scheme presented here can be found in Ref. 208.

The present approach restricts the general formalism to the doubly occupied orbitals, which are localized in a first step according to the Pipek–Mezey procedure.²¹⁷ Then, a Mulliken population analysis is performed for each localized orbital to determine the atoms contributing the most to the orbital. A charge-weighted average position r_c and a maximum distance r_max are used to draw a sphere with radius αr_max (α is an adjustable parameter, default value of 1.0) centered at r_c. An example is shown in Fig. 11. Localized orbitals whose assigned spheres overlap are referred to as strong pairs. Weak pairs are those for which the assigned spheres do not overlap. Following the work of Carter and co-workers, all double excitations, which arise from simultaneous single excitations from weak orbital pairs, are neglected. This scheme leads to a straightforward program implementation with a conceptual simplicity in terms of well-defined localized orbitals.

As an example, this scheme was applied to Diels–Alder-type, strain-promoted, oxidation-controlled, cycloalkyne-1,2-quinone cycloaddition (SPOCQ) reactions,²¹⁸ which is an interesting option in the quest for faster metal-free click cycloaddition reactions. MR calculations using the MR-AQCC, MRCI, and MRCI+Q (including size-extensivity corrections computed with the Davidson–Silver method^77,219) were performed for both the reactant complex and transition state (TS) of the 1,2-benzoquinone (QUIN) plus bicyclo[6.1.0]non-4-yne (BCN) system²¹⁸ (see Fig. 12). Four active spaces were used at the MR level, in which eight active electrons were distributed in (a) eight active orbitals [CAS(8,8)]; (b) seven active orbitals [CAS(8,7)]; (c) eight active orbitals with restrictions of single excitations from a restricted active space (RAS) and at most one electron in an auxiliary (AUX) space, denoted [RAS(2)/CAS(4,4)/AUX(2)-1ex] or AS(2/4/2,8)-1ex in short; and (d) the same space as in (c) only using double excitations, denoted AS(2/4/2,8)-2ex. In most cases, the 6-31G^* basis set was used, but for the best-estimate calculation, the 6-311G(2d) basis set replaced the 6-31G^* basis for the atoms directly involved in the reaction (C6 to C15, O16, and O17), denoted as 6-311G(2d)-red (see Fig. 12 for atom numbering). The core orbitals were always frozen.

Table III lists different choices investigated for freezing localized orbitals in the correlation treatment, in combination with the WP approximation. In most of the cases, the σ orbitals most distant from the reaction center (seven C–C bonds within the range of C1 to C7 and eight C–H bonds within the range of H18 to H25 of BCN) (denoted as σ-freeze) were frozen. Both the orbital freezing (Table III, calc. 1 and calc. 2) and weak pairs (see Table III, calc. 3 and calc. 4) had a relatively small effect on the reaction’s activation energy (AE). The number of CSFs for MR-AQCC(8,7) is 18.161 × 10⁹ (Table III, calc. 1); the WP approximation reduces the CSFs nearly by a factor of 4. Freezing the σ orbitals (σ-freeze) (Table III, calc. 2) reduces the number of CSFs even more, in total by almost seven times. Although for most cases listed in Table III the reduction in the number of CSFs due to the WP approximation is nearly a factor of three, without the σ-freeze, the reduction is even more. The MR-AQCC best estimate activation energy is 10.5 kcal/mol. Corrections include zero-point energy and thermal contributions at 298 K (−0.55 kcal/mol) and solvent effects in 1,2 dichloroethane (−1.8 kcal/mol), both taken from M06-2X DFT calculations. This gives an enthalpy of activation of 8.2 kcal/mol, somewhat higher than the experimental value of 4.5 kcal/mol. Different DFT calculations using different functionals yield values ranging from 4.9 kcal/mol to 10.1 kcal/mol.²¹⁸ 

The examples presented above illustrate the wide variety of “difficult” chemical problems to which ab initio multireference correlation methods should be applied. These applications include the important class of π-conjugated biradical compounds, the treatment of excited potential energy surfaces at highest levels including nonadiabatic couplings, the combination of COLUMBUS with surface hopping dynamics software (NEWTON-X and SHARC), and a fully variational spin-orbit MRCI. The straightforward implementation of local electron correlation by means of localized orbitals and the weak pair approximation leads to a significant enhancement of the accessible molecular size for MRCI and MR-AQCC calculations. The emerging new capabilities include the ACME MCSCF and GCF methods, which allow practically unlimited CSF expansion sets, analytic spin-orbit MRCI energy gradients, and standardized connections to surface hopping dynamics program packages.

Issues concerning accessible sizes of molecules, basis sets, and MRCI dimensions cannot be answered easily because too many factors play a crucial role, and not any combination thereof in terms of accessible limits will be valid. However, there are many tools available in COLUMBUS to mitigate drastic increases in computational resources. The exponential increase in the CAS dimension is well known and is undoubtedly a major factor to be considered, especially because of the uncontracted character of the wavefunction. However, many examples show that in the variational approaches used here, in comparison with perturbational ones, the size of the active space can be restricted to the truly open-shell orbitals (static electron correlation), as has also been discussed in the same spirit by Pulay.²²⁰ However, this is not the only way to find relief from the factorial increase in the CAS. Alternatives exist due to the flexibility of construction schemes in the wavefunctions. Restrictions can be imposed by specifying cumulative occupation limits for the construction of the GUGA configuration space as well as cumulative spin restrictions for each orbital. These allow straightforward construction of direct-product expansion spaces,²²¹ generalized valence bond^189,190 (GVB) type expansions, restricted active spaces²²² (RAS), and many other possibilities. Concerning CSF expansion sizes, several billion are accessible on standard parallel computer systems and basis set sizes up to 1023 are possible. The applications discussed here demonstrate the wide scope of possibilities and the generality and flexibility of GUGA as applied within COLUMBUS.

R.S. and S.R.B. were supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, Gas Phase Chemical Physics Program, through Argonne National Laboratory under Contract No. DE-AC02-06CH11357. E.A.C. is grateful for support from the U.S. Department of Energy, Office of Science, Offices of Basic Energy Sciences and Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing, via Award No. DE-AC02-05CH11231. S.M. was funded by the Department of Energy, Award No. DEFG02-08ER15983. L.T.B. was funded by the High-Energy Laser Joint Technology Office, Albuquerque, NM. D.R.Y was supported by the US Department of Energy (Grant No. DE-SC0015997). C.P. acknowledges support from the Department of Energy (Grant No. DE-SC0001093), the National Science Foundation (Grant Nos. CHE-1213271 and CHE-18800014), and the donors of the American Chemical Society Petroleum Research Fund. P.G.S. was supported by the National Research, Innovation and Development Fund (NKFIA), Grant No. 124018. H.L. and A.J.A.A. are grateful for support from the School of Pharmaceutical Science and Technology (SPST), Tianjin University, Tianjin, China, including computer time on the SPST computer cluster Arran. M.B. and F.K. thank the support of the Excellence Initiative of Aix-Marseille University (A^*MIDEX), the project Equip@Meso (Grant No. ANR-10-EQPX-29-01), and the WSPLIT project (Grant No. ANR-17-CE05-0005-01). D.N. acknowledges support from the Czech Science Foundation (Grant No. GA18-09914S). S.A.d.M., E.V., and M.M.A.d.N. were funded by the Brazilian agencies: Coordination for the Improvement of Higher Education Personnel (CAPES); Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Project Nos. 303884/2018-5 and 423112/2018-0; and Financier of Innovation and Research (FINEP). The authors also acknowledge computer time at the Supercomputer Center of the Federal University of Rio Grande do Sul (CESUP-UFRGS). R.F.K.S. acknowledges Fundacão de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under Grant No. 2019/07671-4 and CNPq under Grant Nos. 407760/2018-0 and 305788/2018-3. I.B. thanks CNPq through research Grant Nos. 304148/2018-0 and 409447/2018-8. F.B.C.M. gratefully acknowledges the financial assistance of the Brazilian agency CNPq under Project Nos. 307052/2016-8 and 404337/2016-3. F.B.C.M., A.J.A.A., and H.L. thank the FAPESP/Tianjin University SPRINT program (Project No. 2017/50157-4) for travel support. H.L., F.P., M.O., and L.G. acknowledge gratefully computer time at the computer cluster of the Vienna Scientific Cluster, Austria, under Project Nos. 70376, 70726, and 70264.

The authors declare no competing financial interest.