Time-dependent linear-response variational Monte Carlo

TL;DR

This paper extends variational Monte Carlo to compute electronic excitation energies and oscillator strengths using a time-dependent linear-response approach, showing promising results over traditional methods.

Contribution

It introduces a novel linear-response VMC method that incorporates wave-function derivatives, enabling accurate excited-state calculations.

Findings

LR-VMC outperforms CIS in excitation energy accuracy

Method successfully applied to beryllium atom excitations

LR-VMC shows promise for excited-state properties of molecules

Abstract

We present the extension of variational Monte Carlo (VMC) to the calculation of electronic excitation energies and oscillator strengths using time-dependent linear-response theory. By exploiting the analogy existing between the linear method for wave-function optimisation and the generalised eigenvalue equation of linear-response theory, we formulate the equations of linear-response VMC (LR-VMC). This LR-VMC approach involves the first-and second-order derivatives of the wave function with respect to the parameters. We perform first tests of the LR-VMC method within the Tamm-Dancoff approximation using single-determinant Jastrow-Slater wave functions with different Slater basis sets on some singlet and triplet excitations of the beryllium atom. Comparison with reference experimental data and with configuration-interaction-singles (CIS) results shows that LR-VMC generally outperforms CIS for excitation energies and is thus a promising approach for calculating electronic excited-state properties of atoms and molecules.

Tables

Table 1. Table 1: 2 ​ s ​ 3 ​ s 2 𝑠 3 𝑠 2s3s ( S 1 superscript 𝑆 1 {}^{1}S ) excitation energies (in Hartree) for the beryllium atom calculated using CIS and LR-VMC/TDA including the response of the Jastrow parameters (j) and of the Jastrow and orbital parameters (j+o). The experimental value is taken from Ref. 35 .

	CIS	LR-VMC/TDA(j)	LR-VMC/TDA(j+o)	Exp.
VB1	0.378	0.2888(1)	0.2672(1)	0.249
VB2	0.228	0.2880(1)	0.2378(2)	0.249

Table 2. Table 2: 2 ​ s ​ 4 ​ s 2 𝑠 4 𝑠 2s4s ( S 1 superscript 𝑆 1 {}^{1}S ) excitation energies (in Hartree) for the beryllium atom calculated using CIS and LR-VMC/TDA including the response of the Jastrow parameters (j) and of the Jastrow and orbital parameters (j+o). The experimental value is taken from Ref. 35 .

	CIS	LR-VMC/TDA(j)	LR-VMC/TDA(j+o)	Exp.
VB1	2.639	0.6106(1)	0.5578(2)	0.297
VB2	0.470	0.6104(1)	0.321(3)	0.297

Table 3. Table 3: 2 ​ s ​ 2 ​ p 2 𝑠 2 𝑝 2s2p ( P 1 superscript 𝑃 1 {}^{1}P ) excitation energies (in Hartree) for the beryllium atom calculated using CIS and LR-VMC/TDA including the response of the orbital parameters (o). The experimental value is taken from Ref. 35 .

	CIS	LR-VMC/TDA(o)	Exp.
VB1	0.220	0.2358(1)	0.194
VB2	0.189	0.1873(2)	0.194

Table 4. Table 4: 2 ​ s ​ 2 ​ p 2 𝑠 2 𝑝 2s2p ( P 3 superscript 𝑃 3 {}^{3}P ) excitation energies (in Hartree) for the beryllium atom calculated using CIS and LR-VMC/TDA including the response of the orbital parameters (o). The experimental value is taken from Ref. 35 .

	CIS	LR-VMC/TDA(o)	Exp.
VB1	0.068	0.1064(1)	0.100
VB2	0.063	0.0929(2)	0.100

Table 5. Table 5: Oscillator strength f 𝑓 f corresponding to the 2 ​ s ​ 2 ​ p 2 𝑠 2 𝑝 2s2p ( P 1 superscript 𝑃 1 {}^{1}P ) excitation of the beryllium atom computed using CIS and LR-VMC/TDA including the response of the orbital parameters (o). The experimental value is taken from Ref. 36 .

	CIS	LR-VMC/TDA(o)	Exp.
VB1	0.648	0.435(1)	1.34(3)
VB2	0.669	0.57(2)	1.34(3)

Equations

$$|\Psi (\text{p})⟩=\hat{J}(\mathbf{\alpha })e^{\hat{\kappa }(\mathbf{\kappa })}\sum _{I=1}^{N_{\text{CSF}}}{c}_I|C_I⟩,$$

$$|{\varphi }_i⟩=\sum _{\mu =1}^{N_{\text{basis}}}{\lambda }_{i\mu }|{\chi }_{\mu }⟩.$$

$$⟨\text{r}|{\chi }_{\mu }⟩=N_n(\zeta )r^{n-1}{\text{e}}^{-\zeta r}{Y}_{\ell ,m}(\theta ,\varphi ),$$

$$|\overline{\Psi }(\text{p})⟩=\frac{|\Psi (\text{p})⟩}{⟨{\Psi }_0|\Psi (\text{p})⟩}$$

$$|{\overline{\Psi }}_{\text{lin}}(\text{p})⟩=|{\Psi }_0⟩+\sum _i\Delta p_i|{\overline{\Psi }}_i⟩,$$

$$|{\overline{\Psi }}_i⟩={\left(\frac{\partial |\overline{\Psi }(\text{p})⟩}{\partial p_i}\right)}_{\text{p}={\text{p}}^0}=|{\Psi }_i⟩-⟨{\Psi }_0|{\Psi }_i⟩|{\Psi }_0⟩,$$

$$E_{\text{lin}}=\underset{\text{p}}{\min }\frac{⟨{\overline{\Psi }}_{\text{lin}}(\text{p})|\hat{H}|{\overline{\Psi }}_{\text{lin}}(\text{p})⟩}{⟨{\overline{\Psi }}_{\text{lin}}(\text{p})|{\overline{\Psi }}_{\text{lin}}(\text{p})⟩},$$

$$\left(\begin{array}{cc}{E}_0& {\text{g}}_{\text{R}}^{\text{T}}/2\\ {\text{g}}_{\text{L}}/2& \text{H}\end{array}\right)\left(\begin{array}{c}1\\ \mathbf{\Delta }\text{p}\end{array}\right)=E_{\text{lin}}\left(\begin{array}{cc}1& {\text{0}}^{\text{T}}\\ \text{0}& \text{S}\end{array}\right)\left(\begin{array}{c}1\\ \mathbf{\Delta }\text{p}\end{array}\right),$$

$$\hat{H}(t)=\hat{H}+\gamma \hat{V}(t),$$

$$|\overline{\Psi }(\text{p}(t))⟩=\frac{|\Psi (\text{p}(t))⟩}{⟨{\Psi }_0|\Psi (\text{p}(t))⟩},$$

$$\frac{\partial }{\partial p_i^{\ast }}\frac{⟨\overline{\Psi }(\text{p}(t))|\hat{H}(t)-i\frac{\partial }{\partial t}|\overline{\Psi }(\text{p}(t))⟩}{⟨\overline{\Psi }(\text{p}(t))|\overline{\Psi }(\text{p}(t))⟩}=0.$$

$$|\overline{\Psi }(\text{p}(t))⟩$$

$$|{\overline{\Psi }}_{ij}⟩$$

$$\text{A}\mathbf{\Delta }\text{p}(t)+\text{B}\mathbf{\Delta }\text{p}(t{)}^{\ast }=i\text{S}\frac{\partial \mathbf{\Delta }\text{p}(t)}{\partial t},$$

$$\Delta \text{p}(t)=\text{X}{e}^{-i{\omega }_{n}t}+{\text{Y}}^{\ast }{e}^{i{\omega }_{n}t},$$

$$\left(\begin{array}{cc}\text{A}& \text{B}\\ {\text{B}}^{\ast }& {\text{A}}^{\ast }\end{array}\right)\left(\begin{array}{c}{\text{X}}_n\\ {\text{Y}}_n\end{array}\right)={\omega }_n\left(\begin{array}{cc}\text{S}& \text{0}\\ \text{0}& -{\text{S}}^{\ast }\end{array}\right)\left(\begin{array}{c}{\text{X}}_n\\ {\text{Y}}_n\end{array}\right).$$

$$\text{A}{\text{X}}_n={\omega }_n\text{S}{\text{X}}_n.$$

$$f_n=\frac{2}{3}{\omega }_n\sum _{\alpha =x,y,z}{\left[({\text{X}}_n+{\text{Y}}_n{)}^{\text{T}}{\mathbf{\mu }}^{\alpha }\right]}^2,$$

$${\mu }_i^{\alpha }$$

$$E_0=⟨E_{\text{L}}(\text{R})⟩,$$

$$g_{\text{L},i}$$

$$g_{\text{R},j}$$

$$E_{\text{L},j}(\text{R})=\frac{H{\Psi }_j(\text{R})}{{\Psi }_0(\text{R})}-\frac{{\Psi }_j(\text{R})}{{\Psi }_0(\text{R})}{E}_{\text{L}}(\text{R}).$$

$$S_{ij}$$

$$H_{ij}$$

$$A_{ij}$$

$$B_{ij}$$

$${\mu }_i^{\alpha }$$

Full text

Keywords: excitation energies - linear method - Tamm-Dancoff approximation - oscillator strengths - beryllium

I Introduction

Quantum Monte Carlo (QMC) methods bk:hammond ; fou+01rmp ; tou16 are a powerful and reliable alternative to wave-function methods and density-functional theory (DFT) for quantum chemistry calculations, thanks to their favorable scaling with system size and to their suitability for high-performance computing infrastructures, such as petascale architectures. Variational Monte Carlo (VMC) bres98 combines Monte Carlo integration for computing the expectation value of the electronic Hamiltonian $\hat{H}$ and the variational principle for the ground state. VMC scales as $N^{3-4}$ (where $N$ is the number of electrons), similar to DFT scaling. The main drawback of any QMC approach is the very large prefactor in the scaling, preventing the systematic use of QMC in quantum chemistry calculations of medium- and large-size systems. This drawback is alleviated by performing massive parallel calculations on supercomputers Coccia:2012kz ; Coccia14 .

A fundamental role is played by the trial wave function, often written as a product of a determinantal part and a bosonic Jastrow factor dru+04prb which depends on interparticle distances (with electron-nucleus, electron-electron, higher many-body terms,…). For example, one can use for the determinantal part a linear combination of configuration state functions (CSF, i.e. spatial- and spin-symmetry adapted linear combinations of Slater determinants of one-electron molecular orbitals) pet12 , or the antisymmetrised geminal power (AGP) ansatz (a single determinant of geminal pairing functions cas+03jcp ; cas+04jcp ; zen14 ). Furthermore, the optimisation of the wave function is crucial for an accurate description of both static and dynamic electron correlation. The linear method tou07 ; umr+07prl ; tou08 allows one to efficiently perform such an optimisation for all the parameters of the wave function, using only the first-order derivatives of the wave function with respect to the parameters.

The calculation of excited-state properties of molecules (from prototypical models to complex organic dyes and biochromophores) still represents an open challenge for theoreticians. The two commonly used approaches are time-dependent density-functional theory, which is not computationally demanding but sometimes lacks accuracy, and wave-function methods, which are more accurate but very computationally demanding. QMC methods were originally formulated for ground state problems and their extension to excited states is not straightforward. Relatively few applications of QMC for electronic excitations are present in literature, see e.g. the singlet and triplet energies for the benchmark CH2 diradical zim09 , the low-lying singlet excited states of biochromophores fili10 ; filippi2011bathochromic ; val12 ; val13 , the $n\rightarrow\pi^{*}$ transition in acrolein tou12 ; flo14 , and the recent extension of the AGP ansatz for calculating excited-state energies zen+15jctc ; dup15 ; zha16a ; zha16b ; neu16 .

The basic idea of the present work stems from the formal analogy existing between the linear method for wave function optimisation and time-dependent linear-response theory bk:mc . Indeed, the generalised eigenvalue equations of linear-response theory in the Tamm-Dancoff approximation (TDA) and of the linear method at the ground-state minimum coincide. Starting from this observation, we derive and implement the linear-response equations in VMC (LR-VMC). This represents an extension of the well-established time-dependent linear-response Hartree-Fock or multiconfiguration self-consistent field methods, taking into account both static and dynamic electron correlations.

The paper is organised as follows. In Section II, VMC and linear-response theory are briefly reviewed, and the LR-VMC method is presented and discussed in detail. Results of LR-VMC calculations in the TDA of some singlet and triplet excitations of the beryllium atom are reported and discussed in Section III. Conclusions and perspectives for future work are given in Section IV.

II Theory

We first briefly review the form of the wave function that we use and the linear optimisation method. We then derive the time-dependent linear-response equations and show how to implement them in VMC.

II.1 Wave-function parametrisation

We consider Jastrow-Slater-type wave functions parametrised as tou07 ; tou08

[TABLE]

where $\hat{J}({\bm{\alpha}})$ is a Jastrow factor operator depending on a set of parameters $\bm{\alpha}$, $e^{\hat{\kappa}(\bm{\kappa})}$ is the orbital rotation operator depending on a set of orbital rotation parameters $\bm{\kappa}$, and $\ket{C_{I}}$ are CSFs with associated coefficients ${\textbf{c}}=\{c_{I}\}$. The CSFs are linear combinations of Slater determinants of orbitals $\ket{\phi_{i}}$, which are expanded in a basis of Slater functions $\{\ket{\chi_{\mu}}\}$

[TABLE]

The Slater functions are centered on the nuclei and their spatial representation is

[TABLE]

each function being characterized by a set of quantum numbers $n,\ell,m$ and an exponent $\zeta$, $Y_{\ell,m}(\theta,\phi)$ are real spherical harmonics, and $N(\zeta)$ a normalisation factor. The full set of parameters to consider is thus ${\textbf{p}}=\left\{\bm{\alpha},{\textbf{c}},\bm{\kappa},\bm{\zeta}\right\}$ where $\bm{\zeta}$ stands for the set of exponents.

II.2 Linear optimisation method

The linear optimisation methodtou07 ; umr+07prl ; tou08 allows one to find the optimal parameters p using an iterative procedure. At each iteration, we consider the intermediate-normalised wave function

[TABLE]

where $\ket{\Psi_{0}}=\ket{\Psi({\textbf{p}}^{0})}$ is the wave function for the parameters ${\textbf{p}}^{0}$ at the current iteration (taken as normalised to unity, i.e. $\langle\Psi_{0}|\Psi_{0}\rangle=1$), and we expand it to linear order in the parameter variations $\bm{\Delta}{\textbf{p}}={\textbf{p}}-{\textbf{p}}^{0}$,

[TABLE]

where $\ket{\overline{\Psi}_{i}}$ are the first-order derivatives of the wave function $\ket{\overline{\Psi}({\textbf{p}})}$

[TABLE]

where $\ket{{\Psi}_{i}}=\left(\partial\ket{\Psi({\textbf{p}})}/\partial p_{i}\right)_{{\textbf{p}}={\textbf{p}}^{0}}$ are the first-order derivatives of the original wave function $\ket{\Psi({\textbf{p}})}$. Using the intermediate-normalised wave function has the advantage that the derivatives in Eq. (6) are orthogonal to $\ket{\Psi_{0}}$, i.e. $\langle\Psi_{0}|\overline{\Psi_{i}}\rangle=0$. We then determine the parameter variations $\bm{\Delta}{\textbf{p}}$ by minimising the corresponding energy

[TABLE]

we update the parameters as ${\textbf{p}}^{0}\to{\textbf{p}}^{0}+\bm{\Delta}{\textbf{p}}$, and iterate until convergence.

The minimisation in Eq. (7) leads to the following generalized eigenvalue equation to be solved at each iteration

[TABLE]

where $E_{0}=\bra{\Psi_{0}}\hat{H}\ket{\Psi_{0}}$ is the current energy, $g_{\text{L},i}=2\bra{\overline{\Psi}_{i}}\hat{H}\ket{\Psi_{0}}$ and $g_{\text{R},j}=2\bra{{\Psi}_{0}}\hat{H}\ket{\overline{\Psi}_{j}}$ are the left and right energy gradients (identical except on a finite Monte Carlo sample), and $H_{ij}=\bra{\overline{\Psi}_{i}}\hat{H}\ket{\overline{\Psi}_{j}}$ is the Hamiltonian matrix in the basis of the first-order wave function derivatives, and $S_{ij}=\langle\overline{\Psi}_{i}|\overline{\Psi}_{j}\rangle$ is the overlap matrix in this basis. Note that in Eq. (8), 0 and ${\textbf{0}}^{\text{T}}$ stand for the zero column vector and the zero row vector, respectively.

II.3 Linear-response theory

Starting from the previously optimised wave function, we introduce now a time-dependent perturbation (e.g, interaction with an electric field) in the Hamiltonian

[TABLE]

where $\gamma$ is a coupling constant. The approximate ground-state wave function $\ket{\Psi({\textbf{p}}(t))}$ evolves in time through its parameters ${\textbf{p}}(t)$, which become now generally complex. As before, it is convenient to introduce the intermediate-normalised wave function

[TABLE]

where $\ket{\Psi_{0}}=\ket{\Psi({\textbf{p}}^{0})}$ is the wave function for the initial parameters ${\textbf{p}}^{0}$, again taken as normalised to unity (i.e., $\langle\Psi_{0}|\Psi_{0}\rangle=1$). At each time, the time-dependent parameters ${\textbf{p}}(t)$ can be found from the Dirac-Frenkel variational principle (see, e.g., Ref.27)

[TABLE]

To apply Eq. (11) in linear order in $\gamma$, we start by expanding the wave function $\ket{\overline{\Psi}({\textbf{p}}(t))}$ around ${\textbf{p}}^{0}$ to second order in the parameter variations $\Delta{\textbf{p}}(t)={\textbf{p}}(t)-{\textbf{p}}^{0}$

[TABLE]

where $\ket{\overline{\Psi}_{i}}$ are the first-order derivatives of $\ket{\overline{\Psi}({\textbf{p}})}$ already introduced in Eq. (6), and $\ket{\overline{\Psi}_{ij}}$ are the second-order derivatives of the wave function $\ket{\overline{\Psi}({\textbf{p}})}$

[TABLE]

where $\ket{{\Psi}_{ij}}=\left(\partial^{2}\ket{\Psi({\textbf{p}})}/\partial p_{i}\partial p_{j}\right)_{{\textbf{p}}={\textbf{p}}^{0}}$ are the second-order derivatives of the original wave function $\ket{\Psi({\textbf{p}})}$. Again, the advantage of using the intermediate-normalised wave function is that the second-order derivatives are orthogonal to $\ket{\Psi_{0}}$, i.e. $\langle\Psi_{0}|\overline{\Psi}_{ij}\rangle=0$. Plugging Eq. (12) into Eq. (11) and keeping only first-order terms in $\bm{\Delta}{\textbf{p}}(t)$, in the limit of a vanishing perturbation ($\gamma\rightarrow 0$), we find

[TABLE]

with the matrices $A_{ij}=\bra{\overline{\Psi}_{i}}\hat{H}-E_{0}\ket{\overline{\Psi}_{j}}=H_{ij}-E_{0}S_{ij}$ where $E_{0}$ is the ground-state energy, $B_{ij}=\bra{\overline{\Psi}_{ij}}\hat{H}\ket{\Psi_{0}}$, and $S_{ij}=\langle\overline{\Psi}_{i}|\overline{\Psi}_{j}\rangle$. If we look for free-oscillation solutions of the form

[TABLE]

where $\omega_{n}$ corresponds to an excitation (or de-excitation) energy, we arrive at the linear-response equation in the form of a non-Hermitian generalized eigenvalue equationbk:mc

[TABLE]

The Tamm-Dancoff approximation (TDA) corresponds to neglecting the contributions from B, leading to

[TABLE]

At the ground-state minimum, i.e. when the energy gradient is zero, the generalised eigenvalue equation of the linear method in Eq. (8) is equivalent to the TDA equation (17) which directly gives excitation energies $\omega_{n}=E_{\text{lin}}-E_{0}$.

Finally, the oscillator strength $f_{n}$ for the transition from the ground state to the excited state $n$ (with excitation energy $\omega_{n}$) can be easily extracted from the response vector $({\textbf{X}}_{n},{\textbf{Y}}_{n})$

[TABLE]

where $\bm{\mu}^{\alpha}$ is the vector containing the transition dipole moments for the component $\alpha$ ($x$, $y$, or $z$) between the ground-state wave function $\ket{\Psi_{0}}$ and the wave-function derivative $\ket{\overline{\Psi}_{i}}$

[TABLE]

and $\hat{\mu}^{\alpha}$ is the electronic dipole operator.

II.4 Realisation in VMC

We now give the expressions for performing linear-response calculations in VMC, referred to as LR-VMC, i.e. for calculating the expressions in Section II.3 in a VMC run. For convenience, we also recall the expressions necessary for the linear optimisation method.

The current ground-state energy is calculated as

[TABLE]

where $E_{\text{L}}({\textbf{R}})=[H\Psi_{0}({\textbf{R}})]/\Psi_{0}({\textbf{R}})$ is the local energy and $\left\langle...\right\rangle$ stands for an average on a finite Monte Carlo sample of points ${\textbf{R}}_{k}$ distributed according to $\Psi_{0}({\textbf{R}})^{2}$, with ${\textbf{R}}=({\textbf{r}}_{1},{\textbf{r}}_{2},...,{\textbf{r}}_{N})$ designating the electron coordinates. The left and right energy gradients are evaluated as

[TABLE]

and

[TABLE]

where $E_{\text{L},j}({\textbf{R}})$ is the first-order derivative of the local energy

[TABLE]

Note that, in the limit of an infinite sample, $\left\langle E_{\text{L},j}({\textbf{R}})\right\rangle=0$ due to the hermiticity of the Hamiltonian, and therefore the left and right gradients become identical.

The elements of the overlap matrix S are calculated as

[TABLE]

and the elements of the matrix H are evaluated as

[TABLE]

The elements of the matrix A are then given by

[TABLE]

and the elements of the matrix B are

[TABLE]

Finally, the expression of the transition dipole moment needed for calculating oscillator strengths is

[TABLE]

where $\mu^{\alpha}({\textbf{R}})=-\sum_{k=1}^{N}r_{k,\alpha}$ is the $\alpha$-component of electronic dipole moment.

In the linear optimisation method, using the non-symmetric estimator of the matrix H in Eq. (25) instead of a symmetrised one has the advantage of leading to the strong zero-variance principle of Nightingale and Melik-Alaverdian zero : in the limit where the current wave function $\ket{\Psi_{0}}$ and its first-order derivatives $\ket{\overline{\Psi}_{i}}$ form a complete basis for expanding an exact eigenstate of the Hamiltonian, the parameter variations $\bm{\Delta}{\textbf{p}}$ and the associated energy $E_{\text{lin}}$ are found from Eq. (8) with zero variance provided that the Monte Carlo sample size is larger than the number of parameters (see discussion in Ref. 14). Unfortunately, this strong zero-variance principle does not apply when solving the linear-response equation (16). However, in the limit where $\ket{\Psi_{0}}$ is an exact eigenstate of the Hamiltonian, the left energy gradient $g_{\text{L},i}$ in Eq. (LABEL:gLi) vanishes with zero variance, and thus the TDA linear equation (17) becomes equivalent to Eq. (8) for calculating excited-state energies even on a finite Monte Carlo sample. Therefore, in this case, the strong zero-variance principle applies to the calculation of the response vectors ${\textbf{X}}_{n}$ and excitation energies $\omega_{n}$.

II.5 Computational details

The calculations shown here were performed using the QMC program CHAMP CHAMP , starting from Hartree-Fock calculations done with GAMESS SchBalBoaElbGorJenKosMatNguSuWinDupMon-JCC-93 . Two Slater basis sets of different sizes were used, namely the VB1 and VB2 basis set from Ref. 31. The VB1 basis set has five $s$ and one $p$ Slater functions ($[5s,1p]$), whereas the VB2 basis set has six $s$, two $p$, and one $d$ Slater functions ($[6s,2p,1d]$). We use a flexible Jastrow factor consisting of the exponential of the sum of electron-nucleus, electron-electron and electron-electron-nucleus terms, written as systematic polynomial and Padé expansions Umr-UNP-XX ; FilUmr-JCP-96 ; GucSanUmrJai-PRB-05 , with 4 electron-nucleus parameters, 5 electron-electron parameters and 15 electron-electron-nucleus parameters. For each VMC calculation, 104 blocks were employed with 104 steps each. One block was used for equilibration of the VMC distribution.

III Results

The beryllium atom was used as a first test of the LR-VMC approach, since accurate experimental reference values for the excitation energies are available from Ref. 35. An accurate description of the Be ground state requires a multiconfigurational wave function for accounting for the near-degeneracy between the 2$s$ and 2$p$ orbitals. However, for these preliminary tests, we present only results of calculations using a Jastrow-Slater single-determinant wave function for the ground state using TDA linear-response theory. This choice is motivated by the fact that a direct comparison between the LR-VMC/TDA method and configuration-interaction-singles (CIS) calculations represents a simple but essential first step for validating our approach. We expect LR-VMC/TDA to outperform CIS because the Jastrow factor in LR-VMC should account for a substantial part of electronic correlation, and we find this to be the case for most of the excitations studied. The results are presented both as errors with respect to the experimental values in Figure 1 and as detailed excitation energies in the subsequent tables.

In Table 1 results for the singlet 2$s$3$s$ (${}^{1}S$) state are reported. The effect of the Slater basis set adopted is dramatic at the CIS level, as a reasonable agreement with the reference experimental value of 0.249 Hartree is found only when the VB2 basis set is used. LR-VMC/TDA values are labelled as follow: (j) designates the response of the Jastrow parameters only, while (j+o) is the response of both the Jastrow and orbital parameters. The response of the Jastrow factor substantially improves upon the CIS VB1 estimate, going from 0.378 to 0.2888(1) Hartree. The excitation energy improves further when the response of the orbital parameters are included in the LR-VMC/TDA calculation, yielding an error of around 0.02 Hartree with respect to the experimental value. Increasing the size of the Slater basis set, i.e. moving from VB1 to VB2, we obtain a fair agreement with the experimental data when both the Jastrow and the orbital parameters are included in the response (0.2378(2) Hartree).

The singlet 2$s$4$s$ (${}^{1}S$) excitation is higher in energy, and CIS fails to recover the experimental result of 0.297 Hartree, for both basis sets, as shown in Table 2. As already mentioned for the 2$s$3$s$ excitation, the response of the Jastrow factor plays an important role for the VB1 basis set, reducing the error in the excitation energy by around 2 Hartree. Including the orbital parameters in the response lowers the excitation energy further to 0.5578(2), but this is still a large overestimate of the experimental value. With the VB2 basis, the LR-VMC/TDA(j+o) calculation outperforms CIS, but a substantial error ($>$0.02 Hartree) still remains for this high-lying excitation. The failure of VB1 and, to a lesser extent, of VB2 is likely due to the poor description of the 4$s$ orbital.

The extension of our proposed approach to $P$ excitations is straightforward, with a relaxation of the spatial symmetry constraints in the orbital rotation parameters. Note that the Jastrow factor employed in this work only depends on interparticle distances, i.e. it has spherical symmetry, and therefore excited states with $P$ symmetry cannot be represented with the wave-function derivatives with respect to the Jastrow parameters. For this reason, only results concerning the response of the orbitals (o) are reported for the $P$ excitations. In Table 3, results for the singlet 2$s$2$p$ (${}^{1}P$) state are given, which is the lowest energy excitation in the beryllium atom. The CIS calculations with the VB1 and VB2 basis sets show a fair agreement with the reference value of 0.194 Hartree, the CIS calculation using the VB2 basis set being only 5 mHartree below it. The LR-VMC/TDA(o) estimate is also close to the experimental reference when the VB2 basis set is employed (0.1873(2) Hartree), while for the VB1 basis set LR-VMC/TDA(o) greatly overestimates the excitation energy.

Similarly, our implementation of linear response allows us to easily compute triplet excitations by considering triplet orbital rotation parameters. The CIS calculation underestimates the correct excitation energy by more than 30 mHartree, while the LR-VMC/TDA(o) excitation energies are very close to the reference values of 0.100 Hartree. The basis set effects are small is this case.

Finally, we computed the oscillator strength $f$ (Table 5) corresponding to the singlet $2s2p$ (${}^{1}P$) excitation, which is non zero according to selection rules. The LR-VMC(o) oscillator strengths seem more sensitive to the basis set compared to the CIS oscillator strengths. Moreover, the inclusion of the Jastrow factor does not improve the oscillator strength.

IV Conclusions and Perspectives

In this work we have presented a formulation of time-dependent linear-response theory in the VMC framework using a Jastrow-Slater wave function. Compared to state-specific or state-average excited-state QMC methods, the advantage of this LR-VMC approach is that, after optimizing only one ground-state wave function, one can easily calculate several excitation energies of different spatial or spin symmetry. Compared to similar linear-response quantum chemistry methods, the presence of the Jastrow factor in LR-VMC allows one to explicitly treat a part of dynamical correlation. A disadvantage of the method is that the excitation energies are much more sensitive than the ground-state energy to the quality of the optimized ground-state wave function. This is true in other linear-response quantum-chemistry methods as well, but is a bigger drawback in a method that employs stochastic optimization.

Using a Jastrow-Slater single-determinant wave function and the TDA, the LR-VMC method was shown to be more accurate that CIS for most of the excitation energies of the beryllium atom that were studied. The LR-VMC approach thus seems a promising method for calculating electronic excitation energies. In the near future, a systematic study on a set of molecules will be an essential step to further validate the proposed methodology, together with calculations using the full response equation beyond the TDA. Also, we will explore using multideterminant wave functions, larger basis sets, and including the wave-function derivatives with respect to the exponents of the Slater functions.

Acknowledgments

EC thanks University of L’Aquila for financial support and the Laboratoire de Chimie Théorique for computational resources. MO and CJU were supported in part by NSF grant ACI-1534965.