# Optimizing excited states in quantum Monte Carlo: A reassessment of double excitations

Stuart Shepard, Anthony Scemama, Saverio Moroni, Claudia Filippi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Stuart Shepard, Anthony Scemama, Saverio Moroni, Claudia Filippi. Optimizing excited states in quantum Monte Carlo: A reassessment of double excitations. The Journal of Chemical Physics, 2025, 163, pp.024119. 10.1063/5.0272733. hal-05039397

HAL Id: hal-05039397

https://hal.science/hal-05039397v1

Submitted on 18 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Optimizing excited states in quantum Monte Carlo: A reassessment of double excitations

Stuart Shepard,¹ Anthony Scemama,² Saverio Moroni,³ and Claudia Filippi¹

¹) MESA+ Institute for Nanotechnology, University of Twente, 7500 AE Enschede,
The Netherlands

²) Laboratoire de Chimie et Physique Quantiques, Université de Toulouse, CNRS, UPS,
31062 Toulouse, France

³) CNR-IOM DEMOCRITOS, Istituto Officina dei Materiali and SISSA Scuola
Internazionale Superiore di Studi Avanzati, Via Bonomea 265, I-34136 Trieste,
Italy

(*Electronic mail: c.filippi@utwente.nl)

(*Electronic mail: moroni@democritos.it)

(*Electronic mail: scemama@irsamc.ups-tlse.fr)

(*Electronic mail: stushep8989@gmail.com)

Quantum Monte Carlo (QMC) methods have proven to be highly accurate for computing excited states, but the choice of optimization strategies for multiple states remains an active topic of investigation. In this work, we revisit the calculation of double excitation energies in nitroxyl, glyoxal, tetrazine, and cyclopentadienone, exploring different objective functionals and their impact on the accuracy and robustness of QMC. A previous study for these systems employed a penalty functional to enforce orthogonality among the states, but the chosen prefactors did not strictly ensure convergence to the target states. Here, we confirm the reliability of previous results by comparing excitation energies obtained with different functionals and analyzing their consistency. Additionally, we investigate the performance of different functionals when starting from a pre-collapsed excited state, providing insight into their ability to recover the target wave functions.
# I. INTRODUCTION

While mainly employed to study ground-state properties of molecular and condensed systems, quantum Monte Carlo (QMC) methods$^{1-3}$, specifically, variational (VMC) and fixed-node diffusion Monte Carlo (DMC) have recently been shown to be highly accurate approaches also for the computation of excited states of various characters$^{4-11}$. QMC methods provide a stochastic solution to the electronic structure problem, enabling the use of correlated wave functions with explicit dependence on inter-particle coordinates through a Jastrow correlation factor. These methods typically combine Jastrow factors with compact determinantal expansions but can also employ functional forms as complex as actual neural networks to describe the correlated nature of excited states.

When optimizing the wave function parameters for multiple states of the same symmetry, early calculations employed a state-average approach, where a common set of orbitals and Jastrow parameters were used for all states and optimized to minimize a weighted average of the energies, with the linear coefficients in the expansion ensuring orthogonality among the states$^{12}$. To overcome the potential limitations of these constrained wave functions, state-specific optimization has been used in combination with variance minimization to target individual states$^{13}$. However, using variance as the objective function was found to cause problems in some cases, leading to the loss of the target state during minimization$^{7}$. Alternatively, several approaches to multistate, state-specific optimization have resorted to minimizing objective functionals that include the energies of the relevant states along with a penalty for nonzero overlaps among them$^{9,14-18}$, performing the optimization either one state at a time or simultaneously.

Until recently$^{18}$, it was unclear how to choose the prefactors in the penalty terms to ensure optimization to the target states. In our previous work$^{9}$, we obtained accurate double excitation energies at the VMC and DMC levels for several challenging molecules by optimizing all states simultaneously using a modification of a penalty functional$^{16}$ and relatively small determinantal expansions of balanced quality among the states of interest. However, these calculations were performed with prefactor choices that do not strictly guarantee optimization to the desired states, calling into question the validity of the results and the effectiveness of the approach.

In the present work, we revisit the computation of double excitation energies of nitroxyl, glyoxal, tetrazine, and cyclopentadienone, aiming to reassess the accuracy and robustness of QMC for this class of excitations. With the use of different objective functionals in the optimization
of the multiple states, we confirm that the previous claims of accuracy fully hold, reinforcing the reliability of QMC also for excited states of double character. Furthermore, we explore the ability of different functionals to recover the correct excited states under varying orthogonality prefactor conditions. In particular, we analyze their performance in the extreme case of a pre-collapsed excited state, where the optimization process is especially challenging, gaining better insight into the limitations and strengths of different optimization strategies for excited-state QMC calculations.

## II. METHODS

In the following, we first present the form of the trial wave function. Then, we will focus on three choices of energy functionals with an orthogonality penalty recently employed for the state-specific optimization of excited states in quantum Monte Carlo calculations⁹,¹⁶,¹⁸.

### A. Wave Functions

We employ here QMC trial wave functions of the Jastrow-Slater form,

$$ \Psi_I = J_I \sum_{i=1}^{N_{\text{CSF}}} c_{Ii} C_i[\{\phi\}_I], \quad (1) $$

where $C_i$ are spin-adapted configuration state functions (CSF), $c_{Ii}$ the expansion coefficients, $J_I$ the Jastrow factors, and $\{\phi\}_I$ the one-particle orbitals. In all calculations presented here, the ground and excited states have the same symmetry, so the CSFs, $C_i$, share the same orbital occupation patterns for all states. Note that each state, $I$, in Eq. 1 has its own optimal Jastrow factor, set of orbitals, and expansion coefficients. When each state has its own optimizable Jastrow and orbital parameters, the approach is referred to as a state-specific optimization, as opposed to a state-average optimization where all states share a common Jastrow factor and set of orbitals which are optimized by minimizing a weighted average of the energies of the individual states⁷,¹². In a state-average optimization, only the expansion coefficients are state-dependent and ensure orthogonality among the states.

The initial orbitals in the VMC optimization are obtained from a $N_{\text{state}}$-state-average [SA($N_{\text{state}}$)] CASSCF calculation. The SA-CASSCF orbitals are used as the molecular orbital basis for the CIPSI expansions, which are generated in a balanced manner for the states of interest as described in Ref.⁹. The CIPSI expansion coefficients, $c_{Ii}$, are then used as initial coefficients in the deter-
minantal components of the VMC wave functions. The subsequent state-specific optimization in VMC will lead to a different set of molecular orbitals for each state.

## B. Penalty-Based Objective Functions

Pathak et al.¹⁶ proposed the following objective function for the state-specific optimization of a given state I,

$$O_{<}[\Psi_{I}] = E_{I}[\Psi_{I}] + \sum_{J<i} \lambda_{IJ} |S_{IJ}|^{2}, \qquad (2)$$

where the normalized overlap, $S_{IJ}$, is given by,

$$S_{IJ} = \frac{\langle \Psi_I | \Psi_J \rangle}{\sqrt{\langle \Psi_I | \Psi_I \rangle \langle \Psi_J | \Psi_J \rangle}}. \qquad (3)$$

The minimization of $O_{<}$ penalizes non-zero overlaps with states lower in energy than the state under consideration. Each overlap penalty is weighted by a state-dependent prefactor, $\lambda_{IJ} > 0$. The approach proposed in Ref.¹⁶ involves optimizing states one at a time starting from the ground state. Higher energy states are then optimized using the fixed, pre-optimized, wave functions of lower energy states, or anchoring states. It was also proposed that the prefactors, $\lambda_{IJ}$, should have values on the order of $\Delta E_{IJ} = E_{I} - E_{J}$. As discussed below, with these objective functions, it is also possible to optimize all states at the same time, obtaining in principle the same optimal wave functions.

In a more recent work by the same group¹⁸, a more general optimization approach was introduced which employs a single objective function to simultaneously optimize an ensemble of states,

$$O_{e}[\{\Psi_{I}\}] = \sum_{I} w_{I} E_{I} [\Psi_{I}] + \lambda \sum_{I} \sum_{J<i} |S_{IJ}|^{2}, \qquad (4)$$

with weights, $w_{I} > 0$, and a state-independent prefactor, $\lambda > 0$, applied to the orthogonality penalty. To ensure that the target ensemble of states is the minimum of $O_{e}$, the weights (which sum to unity) must obey the condition $w_{J} > w_{I}$ for $E_{J} < E_{I}$ and the minimum value for $\lambda$ must be larger than the maximum of the critical value for each pair of states,

$$\lambda > \lambda_{c} = \max_{J<i} \left[ \Delta E_{IJ} \frac{w_{I}w_{J}}{w_{J} - w_{I}} \right]. \qquad (5)$$

The gradients of $O_{e}$ were shown to be equivalent to those of the objective functions (Eq. 2) in the limit of $w_{I}/w_{J} \to 0$ with $E_{I} > E_{J}$. In this limit, the prefactors of the penalty must simply satisfy the condition $\lambda_{IJ} > \Delta E_{IJ}$, as suggested in the first paper¹⁶.
Wheeler et al.¹⁸ proposed an iterative approach to determine the weights. We add here how these weights can be calculated exactly for a system with a single unique ΔE (i.e. a 2-state system, a 3-state system with a doubly degenerate excited state, etc. [see the Appendix]).

In our previous work⁹, we employed an objective function similar to Eq. 2, namely, involving the energy of a single state but, as opposed to using anchoring states, we optimized all states at the same time to lower the computational cost and benefit from the use of correlated sampling in the optimization. In addition, we modified Eq. 2 to penalize the states for having non-zero overlap with all other states, regardless of energy ordering, which translates to one difference in the objective function,

$$O_{\neq}[\Psi_I] = E_I[\Psi_I] + \sum_{J \neq I} \lambda_{IJ} |S_{IJ}|^2, \qquad (6)$$

namely, the ≠ in the sum of the overlap penalty, The prefactors, λIJ, were chosen to be symmetric with respect to the interchange of I and J. Consequently, to exemplify this in the case of the ground state, while the objective function in Eq. 2 is simply,

$$O_{<}[\Psi_1] = E_1[\Psi_1], \qquad (7)$$

since there are no states lower in energy, Eq 6 becomes

$$O_{\neq}[\Psi_1] = E_1[\Psi_1] + \sum_{J \neq 1} \lambda_{1J} |S_{1J}|^2. \qquad (8)$$

This modification was motivated by the fact that, in circumstances where states are close in energy, their order may change throughout the optimization while, by introducing ‘≠’ in Eq. 2, the ordering of the states in energy does not need to be known in advance. Using this penalty function, we optimized the variational parameters of Jastrow-Slater wave functions (Eq. 1) built with CIPSI expansions and obtained accurate excitations energies of double character for a set of medium-sized molecules which are challenging for high-level coupled cluster methods.

However, despite our successful application of Eq. 6, its utilization can lead to potential issues in the optimization¹⁸. For a symmetric λIJ, minimizing the objective functions O< [ΨI] over the states of interest is in fact equivalent to minimizing an ensemble function, O< [ΨI], where all weights wI are equal. As a result, the states are constrained to be orthogonal while the sum of the energies of the states is minimized. Since the sum of the energies of a set of eigenstates is equal to that of the states resulting from a rotation of the same eigenstates, the ensemble functional with equal weights has a degenerate minimum. In principle, this can result in an arbitrary rotation of
the states yielding an arbitrary excitation energy smaller than that given by the eigenvalues. The presence of an issue with using equal weights can also be inferred from Eq. 5: for equal weights, $\lambda_c \to \infty$ so that, for a finite $\lambda$, the objective function will in principle not minimize to the target set of states.

The use of the objective function $O_\neq$ does of course not pose a problem if $\lambda_{IJ}$ is allowed to be non-symmetric. In this case, the equivalence with $O_e$ becomes apparent by inspecting the gradients with respect to the parameters. If we consider only two states, we can rewrite the gradients of $O_e$ with respect to a parameter, $p_1$, of $\Psi_1$ and a parameter, $p_2$ of $\Psi_2$, as

$$ \begin{aligned} \frac{\partial O_e}{\partial p_1} &= w_1 \frac{\partial}{\partial p_1} \left( E_1[\Psi_1] + \frac{\lambda}{w_1} |S_{12}|^2 \right) \\ \frac{\partial O_e}{\partial p_2} &= w_2 \frac{\partial}{\partial p_2} \left( E_2[\Psi_2] + \frac{\lambda}{w_2} |S_{12}|^2 \right), \end{aligned} \quad (9) $$

and compare them with the gradients of $O_\neq$,

$$ \begin{aligned} \frac{\partial O_\neq[\Psi_1]}{\partial p_1} &= \frac{\partial}{\partial p_1} (E_1[\Psi_1] + \lambda_{12} |S_{12}|^2), \\ \frac{\partial O_\neq[\Psi_2]}{\partial p_2} &= \frac{\partial}{\partial p_2} (E_2[\Psi_2] + \lambda_{21} |S_{21}|^2). \end{aligned} \quad (10) $$

Since $S_{12} = S_{21}$, the gradients of $O_e$ and $O_\neq$ are equivalent, up to a multiplicative factor, if $\lambda_{12} = \frac{\lambda}{w_1}$ and $\lambda_{21} = \frac{\lambda}{w_2}$. Therefore, in practice, $\lambda_{IJ}$ in $O_\neq$ should not be symmetric if optimization to the target states is desired, undermining the original motivation of introducing $O_\neq$, namely, not to have to decide a priori the ordering of the states.

In the following, we will revisit the results of Ref.⁹ employing either $O_<$ or $O_e$. In both cases, we will optimize all states simultaneously as opposed to using anchoring states in combination with $O_<$¹⁶.

### III. COMPUTATIONAL DETAILS

Calculations on nitroxyl, glyoxal, tetrazine, and cyclopentadienone are identical to those performed in Shepard et al.⁹ in regards to geometry, basis set, and initial trial wave function. From the referenced work, we focus on the smallest CIPSI expansions in the SA-CASSCF orbital basis which were described in the corresponding SI. To reiterate briefly for the present work, all calculations use scalar-relativistic energy-consistent Hartree Fock (HF) pseudopotentials¹⁹ with the corresponding aug-cc-pVDZ Gaussian basis sets with diffuse function exponents taken from
the corresponding all-electron Dunning's correlation-consistent basis sets²⁰. The Jastrow, orbital, and CSF parameters of each state are optimized using the stochastic reconfiguration (SR) method²¹,²². All VMC calculations are performed with the software package CHAMP²³, all CIPSI calculations with Quantum Package²⁴,²⁵, and all HF and SA($N_{\text{state}}$)-CASSCF calculations with GAMESS(US)²⁶. For more details on the present calculations, refer to Ref.⁹ and, for the CIPSI/QMC approach, refer to Refs.⁵⁻⁸.

# IV. RESULTS

## A. Checking Reliability of Previous Results

We revisit the VMC calculations of Shepard et al.⁹ on the double excitation energies of four prototypical molecules, leveraging reliable variational optimization approaches and reassessing the results. As already discussed, in our previous work, the objective functional $O_{\neq}$ was minimized using a symmetric $\lambda_{IJ}$, which can in principle be problematic and lead to an incorrect energy separation among the states. Here, we wish to show that the results from Ref.⁹ are accurately representing the actual excitation energies by comparing them with those obtained with an appropriate choice of objective function, $O_{<}$, and constraints $\lambda_{IJ}$.

For all molecules, we will repeat the optimizations of the smallest CIPSI expansions built on CASSCF orbitals (see SI of Ref.⁹). Even though the use of natural orbitals or larger expansions on CASSCF orbitals improves the excitation energies in VMC (see main text of Ref.⁹), this choice allows us to perform tests on consistently constructed expansions which are limited in size also for glyoxal and tetrazine. We compare the previous results to those obtained using $O_{<}$ as the objective functional in Table I. Since this functional imposes no constraints on the ground-state optimization, it unequivocally reflects the optimized ground-state wave function, which, according to the variational principle, is the closest to the eigenstate for the given functional form. This objective functional, in turn, ensures that the minimum-energy state orthogonal to the ground state (i.e., the first excited state) corresponds to the next higher eigenstate, provided that $\lambda_{21}$ satisfies the condition $\lambda_{21} > (E_2 - E_1)$.

The reported energies and excitation energies in Ref.⁹ all use values of $\lambda_{21}$ that are safely above this minimum limit. To maintain consistency, we adopt the same values in the present calculations, except of course for $\lambda_{IJ} = 0$ when $I < J$, as well as the same optimization procedure (i.e. number
of steps and SR damping factors). The resulting energies and excitation energies are presented in Table I for the smallest CIPSI expansions on on CASSCF orbitals of all four molecules investigated in Shepard *et al.*⁹.

<table><caption>TABLE I. VMC total energies (E, Ha) and excitation energies (ΔE, eV) obtained in Ref.⁹ ($O_\ddagger$ with symmetric $\lambda_{IJ}$) and in the present work ($O_<$). In all 2-state optimizations, all non-zero entries of $\lambda_{IJ}$ are 1.0 Ha. The same optimization parameters as in Ref.⁹ are used.</caption><thead><tr><th></th><th>E (1<sup>1</sup>A')</th><th>E (2<sup>1</sup>A')</th><th>ΔE</th></tr></thead><tbody><tr><td>Nitroxyl</td><td></td><td></td><td></td></tr><tr><td>O<sub>‡</sub></td><td>-26.4891(2)</td><td>-26.3271(2)</td><td>4.41(1)</td></tr><tr><td>O<sub>&lt;</sub></td><td>-26.4896(2)</td><td>-26.3264(2)</td><td>4.44(1)</td></tr><tr><td>Glyoxal</td><td></td><td></td><td></td></tr><tr><td>O<sub>‡</sub></td><td>-44.5729(2)</td><td>-44.3633(2)</td><td>5.70(1)</td></tr><tr><td>O<sub>&lt;</sub></td><td>-44.5729(2)</td><td>-44.3630(2)</td><td>5.71(1)</td></tr><tr><td>Tetrazine</td><td></td><td></td><td></td></tr><tr><td>O<sub>‡</sub></td><td>-52.1826(2)</td><td>-52.0008(2)</td><td>4.95(1)</td></tr><tr><td>O<sub>&lt;</sub></td><td>-52.1830(2)</td><td>-52.0008(2)</td><td>4.96(1)</td></tr><tr><td>Cyclopentadienone</td><td>E (1<sup>1</sup>A<sub>1</sub>) E (2<sup>1</sup>A<sub>1</sub>) E (3<sup>1</sup>A<sub>1</sub>) ΔE<sub>12</sub> ΔE<sub>13</sub> ΔE<sub>23</sub></td><td></td><td></td></tr><tr><td><sup>a</sup>O<sub>‡</sub></td><td>-46.7532(2) -46.5281(2) -46.4930(2) 6.13(1) 7.08(1) 0.96(1)</td><td></td><td></td></tr><tr><td><sup>b</sup>O<sub>&lt;</sub></td><td>-46.7534(2) -46.5278(2) -46.4937(2) 6.14(1) 7.07(1) 0.93(1)</td><td></td><td></td></tr></tbody></table>

Overall, the two sets of results are in good agreement, confirming the reliability of the QMC calculations of Shepard *et al.*⁹ on the double excitation energies of these molecules. Notably, the excitation energy in the nitroxyl calculation is in disagreement for the reported precision, the excitation obtained with the $O_\ddagger$ functional being slightly lower (we note that the third state of cyclopentadienone also differs, but with a higher excitation energy, which is not attributable to the issue we are analyzing). Nitroxyl is in fact the only system for which we observe a growing overlap when the $O_\ddagger$ objective function is employed as also shown in Figure 1. The optimization of the nitroxyl wave functions were however particularly long, allowing us to observe the progressive mixing between the states and the corresponding reduction in the excitation energy, an effect which

<footnote><sup>a</sup> $\lambda_{12,21} = 1.5$ Ha, $\lambda_{13,31} = 2.0$ Ha, $\lambda_{23,32} = 1.0$ Ha.</footnote>

<footnote><sup>b</sup> $\lambda_{21} = 1.5$ Ha, $\lambda_{31} = 2.0$ Ha, $\lambda_{32} = 1.0$ Ha, with the rest zero.</footnote>
could have remained hidden in the shorter optimizations for the larger molecules. Increasing the weighting factor, $\lambda$, in the penalty term mitigates the issue but does not entirely eliminate it: as $\lambda$ is increased from 1 to 2 to 4 a.u, state mixing appears at progressively larger step counts. On the other hand, the use of the $O_<$ function in combination with the smallest value of $\lambda$ leads to a plateau in the overlap and prevents the reduction of the excitation energy, confirming the robustness of the approach.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Variational optimization of the ground and excited states of nitroxyl using the objective functions (red) $O_<$ and (purple) $O_\neq$. We plot (top) the overlap between the two states and (bottom) the excitation energy as a function of the optimization step and for increasing values of $\lambda$ factor (a.u.) in the penalty.</figcaption></figure>

Therefore, in practice, when the $O_\neq$ with a symmetric $\lambda_{IJ}$ is employed as in Ref. $^9$, the states might remain close to the target ones in optimization runs of typical length if a sufficiently large weight in the penalty is adopted and if the starting trial wave functions are of good quality as in our case. However, in general, this practice will lead to state admixture and an underestimated excitation energy, and one should instead use appropriate constraints in combination with the $O_<$ or $O_e$ functionals. That said, the excitation energies obtained with $O_<$ in Table I confirm that the results in Ref. $^9$ reflect the predictive capability of the QMC/CIPSI approach for the computation of double excitations which was the primary focus of that work.
## B. Testing Objective Functions with Pre-Collapsed States

To challenge each approach and further test the prescribed constraints on $\lambda_{IJ}$, we pre-collapse the small CIPSI wave function for the excited state of nitroxyl by optimizing it with no overlap penalty, and then continue the optimization imposing different orthogonality constraints. The lengths and parameters of the SR optimization are identical in all runs, and only the objective functional and penalty prefactor are changed. We represent the prefactors as multiples of a critical value, $\lambda = a\lambda_c$, where we assume $\lambda_c$ to take the extrapolated full CI value of 0.159 a.u. (4.32 eV) for the excitation energy of nitroxyl.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Excitation energy (eV) of nitroxyl during optimization runs started from a pre-collapsed excited-state wave function. The functionals $O_<$ (left), $O_e$ (or $O_\mp$ with non-symmetric $\lambda_{IJ}$) (center), and $O_\mp$ with symmetric $\lambda_{IJ}$ (right) are employed corresponding to different choices of $(\lambda_{12}, \lambda_{21})$, all scaled with $a\lambda_c$.</figcaption></figure>

For the $O_<$ functional, the orthogonality constraint is applied only to the excited state and the
optimization of the ground state is thus decoupled from the excited-state one, the only influence
being a lower sampling efficiency due to the use of a guiding function⁹. As shown in Fig. 2 (left),
the states separate in energy for $\lambda > \lambda_c$ and, at convergence, yield the same excitation energy. In
the case of $\lambda = \lambda_c$, the excited state does not recover and slowly collapses toward the ground state,
and more rapidly for half the value of $\lambda_c$. There is also an acceleration in the optimization for
larger and larger $\lambda$ although the relative speed up reduces after a certain value. For example, a
factor of $a = 6.3$ does not provide a faster convergence of the excited state than $a = 4$ and one must
also take into account that a larger $\lambda$ will in principle introduce more noise into the optimization.
A factor of $a = 6.3$ is roughly equivalent to 1.0 Ha, which was used in Section IV A to validate the
calculations in Shepard *et al*. This verifies that $\lambda = 1.0$ Ha is more than large enough to ensure
optimization to the target states in the tests we performed above.

When using the $O_e$ functional (or equivalently $O_\ddagger$ with non-symmetric $\lambda_{IJ}$), we must set two
prefactors given by $\lambda_{12} = \lambda/w_1$ and $\lambda_{21} = \lambda/w_2$ (Eqs. 9 and 10) where, again, $\lambda = a\lambda_c$ and
the weights $w_1$ and $w_2$ are set as in the Appendix. Looking at Fig. 2 (middle), we find that
the excitation energy follows the expected behavior with increasing scaling of the penalty. For
$a \ge 1.0$, the states separate in energy while, for $a = 0.5$, the excitation energy remains severely
underestimated. The factors of $a = 4.0$ and $a = 10$ show a similar convergence of the excitation
energy, yielding a faster separation between the states than $a = 1$.

Finally, in Figure 2 (right), we observe that, when the prefactors are symmetric ($\lambda_{12} = \lambda_{21}$),
values of $a < 4$ do not recover the target states but lead to a collapse of the excited state. For
$a = 4.0$, the overlap suggests that the states are orthogonal (not shown), yet the excitation energy
is far lower than expected for the optimal wave functions.

V. CONCLUSIONS

In our previous work⁹, QMC calculations were performed for the double excitation energies of
prototypical challenging molecules, where the wave functions were optimized in a state-specific
fashion penalizing the states for having non-zero overlap to all other states. The resulting func-
tional, named here $O_\ddagger$, was chosen to have symmetric prefactors, $\lambda_{IJ}$, for the orthogonality con-
straints between two states so that the ordering of the states does not need to be known *a priori*.
The use of symmetric prefactors, however, guarantees that the states optimize to the target ones
only in the limit of $\lambda_{IJ}$ going to infinity¹⁸.
Here, we revisit our calculations imposing orthogonality only to the lower states using what we
call the $O_<$ functional, and optimizing all states simultaneously. We confirm the accuracy of our
previous results and the predictive capability of the approach outlined therein$^{5-8}$ for computing
excitation energies in QMC. For nitroxyl, however, the excitation energy obtained with the $O_\ddagger$
functional$^9$ is somewhat lower. We attribute this discrepancy to the particularly long optimization
performed with $O_\ddagger$ for this molecule, which leads to a growing overlap for the chosen $\lambda$ prefactor
in the orthogonality constraint. While increasing $\lambda$ delays this growth, appropriate functionals$^{16,18}$
such as $O_<$ or the ensemble $O_e$ variant should be employed to prevent state mixing.

Finally, we test the use of the functionals $O_\ddagger$, $O_<$, and $O_e$ in a state-specific optimization for nitroxyl, starting from an excited-state wave function that was intentionally collapsed toward the ground state. While $O_\ddagger$ fails to recover the target states, both $O_<$ and $O_e$ quickly respond to the presence of the orthogonality constraint for sufficiently large prefactors. The availability of reliable functionals$^{16,18}$, along with effective optimization approaches, enhances the robustness of QMC for excited states, especially when combined with the use of smart multi-determinant wave functions$^{5-9}$.

## ACKNOWLEDGMENTS

We thank Lucas Wagner for useful discussions. This work was in part supported by the Euro-
pean Centre of Excellence in Exascale Computing TREX — Targeting Real Chemical Accuracy
at the Exascale. This project has received funding from the European Union’s Horizon 2020 —
Research and Innovation program — under grant agreement no. 952165. This work used the
Dutch national e-infrastructure with the support of the SURF Cooperative using grant no. NWO-
2022.032.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding authors
upon reasonable request.
# Appendix: Weights for a two-state system

As shown in Ref. ¹⁸, when using the ensemble functional, $O_e$, the pairwise $\lambda_{IJ}$ should be larger than a critical value,

$$ \lambda_{IJ} > \lambda_c^{IJ} = (E_I - E_J) \frac{w_I w_J}{w_J - w_I}, \qquad (A.1) $$

where $w_I < w_J$ if $E_I > E_J$. In practice, given a set of weights $\{w_L\}$, a single critical $\lambda_c$ can be used corresponding to the maximum of the critical values. A possible choice of weights is the one that gives the same critical lambda for all pairs of states, yielding one free parameter, $\lambda_c$, and an iterative approach to determine the weights.

For a 2-state system, if we set $\lambda_c^{21} = \lambda_c = (E_2 - E_1)$, we can rewrite Eq. A.1 as

$$ 1 = \frac{w_1 w_2}{w_1 - w_2}, \qquad (A.2) $$

which, when used with the condition on the sum of the weights, $w_1 + w_2 = 1$, leads to the quadratic equation,

$$ w_1^2 + w_1 - 1 = 0, \qquad (A.3) $$

which has as solutions $w_1 = (\sqrt{5} - 1)/2$ and $w_2 = (3 - \sqrt{5})/2$.

This can be readily generalized to the case of a 3-state system where the energies of the two excited states are equal, etc. For a 3-state system of this type, using that $w_2 = w_3$ and that $w_1 + w_2 + w_3 = 1$, one obtains $w_1 = (\sqrt{2} - 1)$ and $w_2 = w_3 = (1 - 1/\sqrt{2})$.

# REFERENCES

¹W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, “Quantum monte carlo simulations of solids,” Rev. Mod. Phys. **73**, 33–83 (2001).

²A. Lüchow, “Quantum monte carlo methods,” Wiley Interdiscip. Rev.: Comput. Mol. Sci. **1**, 388–402 (2011), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.40.

³B. M. Austin, D. Y. Zubarev, and W. A. Lester, “Quantum monte carlo and related approaches,” Chem. Rev. **112**, 263–288 (2012), pMID: 22196085, https://doi.org/10.1021/cr2001564.

⁴J. Feldt and C. Filippi, *Excited-state calculations with quantum Monte Carlo. Quantum chemistry and dynamics of excited states: methods and applications* (John Wiley & Sons, Ltd, 2020) pp. 247–275.
⁵M. Dash, J. Feldt, S. Moroni, A. Scemama, and C. Filippi, "Excited states with selected configuration interaction-quantum monte carlo: Chemically accurate excitation energies and geometries," J. Chem. Theory Comput. **15**, 4896-4906 (2019), pMID: 31348645, https://doi.org/10.1021/acs.jctc.9b00476.

⁶M. Dash, S. Moroni, C. Filippi, and A. Scemama, "Tailoring cipsi expansions for qmc calculations of electronic excitations: The case study of thiophene," J. Chem. Theory Comput. **17**, 3426-3434 (2021).

⁷A. Cuzzocrea, A. Scemama, W. J. Briels, S. Moroni, and C. Filippi, "Variational principles in quantum monte carlo: The troubled story of variance minimization," J. Chem. Theory Comput. **16**, 4203-4212 (2020), pMID: 32419451, https://doi.org/10.1021/acs.jctc.0c00147.

⁸A. Cuzzocrea, S. Moroni, A. Scemama, and C. Filippi, "Reference excitation energies of increasingly large molecules: A qmc study of cyanine dyes," Journal of Chemical Theory and Computation **18**, 1089-1095 (2022), pMID: 35080893, https://doi.org/10.1021/acs.jctc.1c01162.

⁹S. Shepard, R. L. Panadés-Barrueta, S. Moroni, A. Scemama, and C. Filippi, "Double excitation energies from quantum monte carlo using state-specific energy optimization," Journal of Chemical Theory and Computation **18**, 6722-6731 (2022), pMID: 36314602, https://doi.org/10.1021/acs.jctc.2c00769.

¹⁰D. Pfau, S. Axelrod, H. Sutterud, I. von Glehn, and J. S. Spencer, "Accurate computation of quantum excited states with neural networks," Science **385**, eadn0137 (2024), https://www.science.org/doi/pdf/10.1126/science.adn0137.

¹¹P. B. Szabó, Z. Schätzle, M. T. Entwistle, and F. Noé, "An improved penalty-based excited-state variational monte carlo approach with deep-learning ansatzes," Journal of Chemical Theory and Computation **20**, 7922-7935 (2024), pMID: 39213603, https://doi.org/10.1021/acs.jctc.4c00678.

¹²C. Filippi, M. Zaccheddu, and F. Buda, "Absorption spectrum of the green fluorescent protein chromophore: A difficult case for ab initio methods?" J. Chem. Theory Comput. **5**, 2074-2087 (2009), pMID: 26613149, https://doi.org/10.1021/ct900227j.

¹³L. Otis and E. Neuscamman, "A promising intersection of excited-state-specific methods from quantum chemistry and quantum monte carlo," WIREs Computational Molecular Science **13**, e1659 (2023), https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1659.

¹⁴N. S. Blunt, S. D. Smart, G. H. Booth, and A. Alavi, "An excited-state approach within
full configuration interaction quantum monte carlo,” J. Chem. Phys. **143**, 134117 (2015), https://aip.scitation.org/doi/pdf/10.1063/1.4932595.

<footnote n="15">¹⁵K. Choo, G. Carleo, N. Regnault, and T. Neupert, “Symmetries and many-body excitations with neural-network quantum states,” Phys. Rev. Lett. **121**, 167204 (2018).</footnote>

<footnote n="16">¹⁶S. Pathak, B. Busemeyer, J. N. B. Rodrigues, and L. K. Wagner, “Excited states in variational monte carlo using a penalty method,” J. Chem. Phys. **154**, 034101 (2021), https://doi.org/10.1063/5.0030949.</footnote>

<footnote n="17">¹⁷M. T. Entwistle, Z. Schätzle, P. A. Erdman, J. Hermann, and F. Noé, “Electronic excited states in deep variational monte carlo,” Nature Communications **14**, 274 (2023).</footnote>

<footnote n="18">¹⁸W. A. Wheeler, K. G. Kleiner, and L. K. Wagner, “Ensemble variational monte carlo for optimization of correlated excited state wave functions,” Electronic Structure **6**, 025001 (2024).</footnote>

<footnote n="19">¹⁹M. Burkatzki, C. Filippi, and M. Dolg, “Energy-consistent pseudopotentials for quantum monte carlo calculations,” J. Chem. Phys. **126**, 234105 (2007), https://doi.org/10.1063/1.2741534.</footnote>

<footnote n="20">²⁰R. A. Kendall, T. H. Dunning Jr, and R. J. Harrison, “Electron affinities of the first-row atoms revisited. systematic basis sets and wave functions,” J. Chem. Phys. **96**, 6796–6806 (1992).</footnote>

<footnote n="21">²¹S. Sorella, “Generalized lanczos algorithm for variational quantum monte carlo,” Phys. Rev. B **64**, 024512 (2001).</footnote>

<footnote n="22">²²M. Casula, C. Attaccalite, and S. Sorella, “Correlated geminal wave function for molecules: An efficient resonating valence bond approach,” J. Chem. Phys. **121**, 7110–7126 (2004), https://doi.org/10.1063/1.1794632.</footnote>

<footnote n="23">²³Shinde, R.; Landinez Borda, E. J.; Shepard, S.; Slootman, E.; Cuzzocrea, A.; Azizi, V.; Lopez-Tarifa, P.; Renaud, N.; Umrigar, C. J.; Moroni, S.; Filippi, C. *Cornell-Holland Ab-Initio Materials Package (CHAMP-EU)*, 2024; https://github.com/filippi-claudia/champ (accessed 2024-7-31).</footnote>

<footnote n="24">²⁴Y. Garniron, T. Applencourt, K. Gasperich, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. As-saraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J.-P. Malrieu, M. Véril, M. Caffarel, P.-F. Loos, E. Giner, and A. Scemama, “Quantum package 2.0: An open-source determinant-driven suite of programs,” J. Chem. Theory Comput. **15**, 3591–3609 (2019), pMID: 31082265, https://doi.org/10.1021/acs.jctc.9b00176.</footnote>

<footnote n="25">²⁵A. Scemama, M. Caffarel, A. Benali, D. Jacquemin, and P.-F. Loos, “Influence of pseudopotentials on excitation energies from selected configuration interaction and diffusion monte carlo,” Res. Chem. **1**, 100002 (2019).</footnote>
<sup>26</sup>G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. E. Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, T. Harville, S. Irle, J. Ivanic, K. Kowalski, S. S. Leang, H. Li, W. Li, J. J. Lutz, I. Magoulas, J. Mato, V. Mironov, H. Nakata, B. Q. Pham, P. Piecuch, D. Poole, S. R. Pruitt, A. P. Rendell, L. B. Roskop, K. Ruedenberg, T. Sattasathuchana, M. W. Schmidt, J. Shen, L. Slipchenko, M. Sosonkina, V. Sundriyal, A. Tiwari, J. L. Galvez Vallejo, B. Westheimer, M. Wloch, P. Xu, F. Zahariev, and M. S. Gordon, “Recent developments in the general atomic and molecular electronic structure system,” J. Chem. Phys. **152**, 154102 (2020).