Monika Dash

Jonas Feldt

Saverio Moroni
email: moroni@democritos.it

Anthony Scemama
email: scemama@irsamc.ups-tlse.fr

Claudia Filippi
email: c.filippi@utwente.nl

Excited states with selected CI-QMC: chemically accurate excitation energies and geometries

We employ quantum Monte Carlo to obtain chemically accurate vertical and adiabatic excitation energies, and equilibrium excited-state structures for the small, yet challenging, formaldehyde and thioformaldehyde molecules. A key ingredient is a robust protocol to obtain balanced ground-and excited-state Jastrow-Slater wave functions at a given geometry, and to maintain such a balanced description as we relax the structure in the excited state. We use determinantal components generated via a selected configuration interaction scheme which targets the same second-order perturbation energy correction for all states of interest at different geometries, and we fully optimize all variational parameters in the resultant Jastrow-Slater wave functions. Importantly, the excitation energies as well as the structural parameters in the ground and excited states are converged with very compact wave functions comprising few thousand determinants in a minimally augmented double-ζ basis set. These results are obtained already at the variational Monte Carlo level, the more accurate diffusion Monte Carlo method yielding only a small improvement in the adiabatic excitation energies. We find that matching Jastrow-Slater wave functions with similar variances can yield excitations compatible with our best estimates; however, the variance-matching procedure requires somewhat larger determinantal expansions to achieve the same accuracy, and it is less straightforward to adapt during structural optimization in the excited state.

I. INTRODUCTION

Quantum Monte Carlo (QMC) methods are first-principle approaches that approximately solve the Schrödinger equation in a stochastic manner. The two most commonly used variants, namely, variational (VMC) and diffusion Monte Carlo (DMC) typically employ so-called Jastrow-Slater wave functions where a determinant expansion is multiplied by a Jastrow factor which explicitly depends on the inter-particle distances and accounts for a significant portion of electronic correlation. Because of the presence of the Jastrow factor, much shorter expansions are often needed to describe the Slater component and obtain accurate VMC results, which can then be further improved with the use of DMC. Thanks to the favorable scaling with system size, these methods have been routinely employed to compute the electronic properties, particularly total energies, of relatively large molecules and solids [START_REF] Mitas | Electronic structure by quantum Monte Carlo: atoms, molecules and solids[END_REF][START_REF] Willow | Goldstone quantum Monte Carlo for correlation energies and quasiparticle energy bands of one-dimensional solids[END_REF][START_REF] Hunt | Quantum Monte Carlo calculations of energy gaps from first principles[END_REF][START_REF] Zen | Fast and accurate quantum Monte Carlo for molecular crystals[END_REF][START_REF] Chen | Quantum Monte Carlo simulations of a giant {Ni21 Gd20} cage with a S= 91 spin ground state[END_REF].

Recent methodological advances have reduced the cost per Monte Carlo step of computing energy derivatives to the one of the energy itself also for large multi-determinant wave functions [START_REF] Filippi | Simple formalism for efficient derivatives and multi-determinant expansions in quantum Monte Carlo[END_REF][START_REF] Assaraf | Optimizing the Energy with Quantum Monte Carlo: A Lower Numerical Scaling for Jastrow-Slater Expansions[END_REF]. This has enabled us not only to simultaneously optimize geometry and wave function with as many as 200,000 determinants but also to explore the dependence of the results on different lengths and types of Slater expansions, foraying into QMC wave function choices beyond conventional small-active-space definitions. Capitalizing on these developments, our recent thorough investigation [START_REF] Dash | Perturbatively Selected Configuration-Interaction Wave Functions for Efficient Geometry Optimization in Quantum Monte Carlo[END_REF] of constructing the Slater component by an automated determinant selection through a selected configuration interaction (CIPSI) approach [START_REF] Huron | Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions[END_REF][START_REF] Gouyet | SPIN-CIPSI calculations on H2CO[END_REF][START_REF] Trinquier | Silabutadienes and their silylene isomers. An ab initio study[END_REF][START_REF] Castex | Comparison between experimentally and theoretically determined potential curves of the Ar2* lowest states[END_REF][START_REF] Pelissier | Bonding between transition metal atoms. Ab initio effective potential calculations of Cu2[END_REF][START_REF] Nebot-Gil | CI calculations concerning the styrene cis-trans photoisomerization[END_REF][START_REF] Evangelisti | Convergence of an improved CIPSI algorithm[END_REF][START_REF] Cimiraglia | Second order perturbation correction to CI energies by use of diagrammatic techniques: An improvement to the CIPSI algorithm[END_REF][START_REF] Cimiraglia | Ab initio study of the photodissociation of nitrosoalkanes and nitrosamines[END_REF][START_REF] Cimiraglia | Recent advances in multireference second order perturbation CI: The CIPSI method revisited[END_REF][START_REF] Illas | Selected versus complete configuration interaction expansions[END_REF][START_REF] Harrison | Approximating full configuration interaction with selected configuration interaction and perturbation theory[END_REF][START_REF] Povill | Treating large intermediate spaces in the CIPSI method through a direct selected CI algorithm[END_REF] has lead to accurate predictions of the ground-state energies and structural parameters of butadiene with relatively short Slater expansions. Relevant determinants could be systematically introduced into the wave function that would not have been chosen otherwise based on a manual, intuitive selection.

These developments also open very interesting prospects for the application of QMC to geometry relaxation in the excited state, where most electronic structure methods either lack the required accuracy or are computationally quite expensive due to their scaling with system size. To date, there are very few studies to assess the ability of QMC to predict excited-state geometries [START_REF] Valsson | Photoisomerization of model retinal chromophores: insight from quantum monte carlo and multiconfigurational perturbation theory[END_REF][START_REF] Guareschi | Ground-and excited-state geometry optimization of small organic molecules with quantum Monte Carlo[END_REF][START_REF] Guareschi | Solvent Effects on Excited-State Structures: A Quantum Monte Carlo and Density Functional Study[END_REF][START_REF] Zulfikri | Multiple-Resonance Local Wave Functions for Accurate Excited States in Quantum Monte Carlo[END_REF], while most of the relatively limited literature on excited-state QMC calculations is primarily concerned with vertical excitation energies [START_REF] Schautz | Excitations in photoactive molecules from quantum Monte Carlo[END_REF][START_REF] Schautz | Optimized Jastrow-Slater wave functions for ground and excited states: Application to the lowest states of ethene[END_REF][START_REF] Drummond | Electron emission from diamondoids: a diffusion quantum Monte Carlo study[END_REF][START_REF] Tapavicza | Mixed time-dependent density-functional theory/classical trajectory surface hopping study of oxirane photochemistry[END_REF][START_REF] Zimmerman | Excited states of methylene from quantum Monte Carlo[END_REF][START_REF] Filippi | Absorption spectrum of the green fluorescent protein chromophore: a difficult case for ab initio methods?[END_REF][START_REF] Dubeckỳ | Ground and excited electronic states of azobenzene: A quantum Monte Carlo study[END_REF][START_REF] Send | Electronic excitations of simple cyanine dyes: reconciling density functional and wave function methods[END_REF][START_REF] Filippi | Bathochromic shift in green fluorescent protein: a puzzle for QM/MM approaches[END_REF][START_REF] Valsson | Excitation energies of retinal chromophores: critical role of the structural model[END_REF][START_REF] Valsson | Rhodopsin absorption from first principles: Bypass-ing common pitfalls[END_REF][START_REF] Scemama | Deterministic construction of nodal surfaces within quantum Monte Carlo: the case of FeS[END_REF][START_REF] Scemama | Excitation energies from diffusion Monte Carlo using selected configuration interaction nodes[END_REF][START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF][START_REF] Blunt | Excited-State Diffusion Monte Carlo Calculations: A Simple and Efficient Two-Determinant Ansatz[END_REF]. Importantly, all these studies are characterized by the use of very different wave functions ranging from the simple ansatz of a CI singles wave function to complete active space (CAS) expansions, sometimes truncated or only partially optimized in the presence of the Jastrow factor due to the limitations previously faced in sampling and optimizing large numbers of determinants. In this work, we want to overcome this empiricism in the application of QMC to excited states and a) identify the most efficient protocol to obtain a balanced and robust description of the ground and excited states in VMC at a given geometry; b) extend this protocol to the optimization of the geometry as well as the computation of energy differences between different potential energy surfaces at distinct geometries; c) demonstrate the competence of VMC in determining accurate vertical excitation energies, optimal excitedstate structures, and adiabatic excitations.

To this aim, we focus on the low-lying singlet n → π * excitation of formaldehyde and thioformaldehyde, and compute these excited-state properties employing compact wave functions containing between a few thousand to about 45,000 determinants obtained through different CIPSI selection schemes. These molecules are small but theoretically challenging: a recent VMC study reported an error as large as 0.2 eV on the vertical excitation energy of thioformaldehyde [START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF], and DMC calculations without a Jastrow factor required wave functions with as many as hundred thousand determinants to achieve high accuracy for formaldehyde [START_REF] Scemama | Excitation energies from diffusion Monte Carlo using selected configuration interaction nodes[END_REF]. Here, we obtain VMC vertical excitation energies within chemical accuracy (∼ 0.04 eV) of the extrapolated full CI (FCI) estimates for both molecules with expansions as small as a few thousand determinants in combination with a minimally augmented double-ζ basis set. These excellent results are obtained using CIPSI expansions constructed to yield a comparable second-order perturbation (PT2) correction in the ground and excited states, and fully reoptimized in the presence of the Jastrow factor. Matching the energy variance of the Jastrow-Slater wave functions of the two states [START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF][START_REF] Robinson | Excitation variance matching with limited configuration interaction expansions in variational Monte Carlo[END_REF] to estimate the excitation energy appears to be a more delicate procedure which, in the case of formaldehyde, recovers VMC excitation energies within 0.05 eV of our best estimates for determinantal expansions comprising at least 7000 determinants.

When optimizing the structure, we follow two different selection routes to maintain a balanced treatment of the wave function while changing the geometry: we construct the determinantal component targeting a roughly constant value of either the perturbation correction or the variance of the CIPSI expansions. We find that both schemes are viable to obtain robust VMC geometries, also in the more demanding case of formaldehyde where different correlated methods give a range of variations in the prediction of the CO bond and the out-ofplane angle in the excited state as large as 80 mÅ and 20 • , respectively [START_REF] Budzák | Accurate Excited-State Geometries: A CASPT2 and Coupled-Cluster Reference Database for Small Molecules[END_REF]. With just a minimal basis, we obtain optimal VMC structures converged with fewer than a thousand determinants and in excellent agreement with the coupled cluster estimates, namely, with deviations smaller than a couple of degrees in the angles and 3 and 10 mÅ in the ground-and excited-state bond lengths, respectively. Finally, we compute the difference between the variational minima of the VMC ground-and excited-state potential energy surfaces to evaluate the adiabatic excitation energies. For both molecules, irrespective of the determinant selection mode, we obtain VMC and DMC estimates within better than 0.05 eV of the corresponding coupled cluster values.

The paper is organized as follows. We describe the CIPSI selection scheme employed to obtain a balanced description of multiple states in Section II and present the computational details in Section III. The VMC and DMC vertical excitation energies of formaldehyde and thioformaldehyde are given in Section IV and the results of the VMC structural relaxation in the ground and excited states in Section V. We conclude in Section VI by summarizing the most important outcomes of our investigation and the future prospects of the applicability of our approach.

II. METHODS

The wave functions used in the QMC calculations are of the Jastrow-Slater form, namely, the product of a determinantal component and a positive Jastrow correlation function,

Ψ = J Ndet k c k D k , (1)
where N det is the total number of determinants and the Jastrow factor J depends here on the electron-nucleus and electron-electron distances, explicitly ensuring that the electronic cusp conditions are satisfied.

To construct the determinantal part of the wave function, we employ the CIPSI selected CI algorithm that iteratively allows us to identify energetically important determinants from the FCI space. When one is interested in multiple electronic states, it is important to obtain a description of the CI subspace which leads to a uniform and balanced treatment of all states of interest. To ensure a consistent quality of the wave functions, the selection of the determinants for the multiple states is done in a single run even if the states belong to different symmetry classes. In practice, starting from an initial reference subspace S typically given by the CI singles wave functions of all states, at every iteration, we expand the space by selecting among all singly-and doubly-excited determinants those which contribute the most to a state-average PT2 energy contribution.

If we represent a newly selected determinant with |α , we then compute its second-order energy contribution using Epstein-Nesbet perturbation theory [START_REF] Epstein | The stark effect from the point of view of Schrödinger's quantum theory[END_REF][START_REF] Nesbet | Configuration interaction in orbital theories[END_REF] as

δE (2) n = | α| Ĥ|Ψ n | 2 Ψ n | Ĥ|Ψ n -α| Ĥ |α , (2)
and the corresponding CI coefficients as

c α,n = Ψ n | Ĥ|α Ψ n | Ĥ|Ψ n -α| Ĥ|α , (3)
where n denotes a state and Ψ n is the current CI wave function in the space S. As suggested by Angeli et al. [START_REF] Angeli | Multireference perturbation CI II. Selection of the zero-order space[END_REF], the selection criterion for a determinant |α in the external space is based on the energy contribution

e α = Nstates n w n δE (2) n , (4)
where w n is the inverse of the maximum c 2 k,n with the index k running over the determinants in the current subspace S. The determinant |α is added to S if its energy contribution e α is higher than a threshold, which is automatically adjusted in such a way that the number of determinants in S is doubled at each iteration. In the applications presented in this work, the ground state is a closed-shell and the excited state an openshell singlet, both having single-reference character with the dominant configuration state function (CSF) comprising a single determinant in the ground state and two determinants in the excited state. Since the expansion is performed in the basis of the determinants and not of the CSFs, reweighting the energy contributions by w n recovers the balance for the different states. As a consequence, this criterion ensures that the multiple states resulting from the expansion will have approximately equivalent net PT2 corrections at every CIPSI iteration.

The many-body Hamiltonian, Ĥ, is diagonalized in the subspace S at every iteration to determine the current CI coefficients. One keeps expanding until the wave function has a required number of determinants or until another chosen criterion is met as for instance a desired total PT2 correction or variance of the CI wave functions as further discussed below in the Result sections [START_REF] Dash | Perturbatively Selected Configuration-Interaction Wave Functions for Efficient Geometry Optimization in Quantum Monte Carlo[END_REF][START_REF] Giner | Quantum Monte Carlo with reoptimised perturbatively selected configuration-interaction wave functions[END_REF][START_REF] Caffarel | Recent Progress in Quantum Monte Carlo[END_REF][START_REF] Giner | Fixed-node diffusion Monte Carlo potential energy curve of the fluorine molecule F2 using selected configuration interaction trial wavefunctions[END_REF][START_REF] Giner | Using perturbatively selected configuration interaction in quantum Monte Carlo calculations[END_REF]. Additionally, we impose that the selected wave functions are spin-adapted, namely, eigenfunctions of the Ŝ2 operator: all determinants corresponding to the spatial occupation patterns currently present in S are added to the reference space before the Hamiltonian is diagonalized [START_REF] Applencourt | Spin adaptation with determinant-based selected configuration interaction[END_REF][START_REF] Garniron | Quantum Package 2.0: An Open-Source Determinant-Driven Suite of Programs[END_REF]. Although wave functions for all the states are obtained through this common selection process, determinants of one symmetry character do not contribute, of course, to a state of a different character. Therefore, we separate these determinants based on space symmetry prior to proceeding with the VMC optimization.

In what we name the "expansion" scheme, we simply use subsequent wave functions of increasing length generated according to the CIPSI algorithm we just described. While the "expansion" scheme is more easily transferable to larger systems where one could experience difficulty in obtaining a large expansion to then truncate, we will also present some tests with wave functions generated from a "truncation" scheme. After having introduced a large number of determinants in the reference space and obtained the wave functions for both states by diagonalizing the Hamiltonian, we project the wave functions in a common subspace of determinants as follows: for each group of determinants, P, of the internal space sharing the same spatial occupation pattern, we compute the quantity

d(P) = k∈P Nstates n c 2 kn , (5)
All the individual spatial occupation patterns of the internal space are sorted in decreasing order of d(P) values, and the truncated determinant space is built by taking the union of the first patterns of the list. This guarantees that the truncated wave functions are spin-adapted and that the final determinants kept after truncation are the most important for all states of interest.

III. COMPUTATIONAL DETAILS

All QMC computations are carried out with the program package CHAMP [START_REF]CHAMP is a quantum Monte Carlo program package[END_REF]. We employ scalar-relativistic energyconsistent Hartree-Fock pseudopotentials and the correlationconsistent Gaussian basis sets specifically constructed for these pseudopotentials [START_REF] Burkatzki | Energy-consistent pseudopotentials for quantum Monte Carlo calculations[END_REF][START_REF]For the hydrogen atom, we use a more accurate BFD pseudopotential and basis set[END_REF]. For the majority of our calculations, we use a minimally augmented double-ζ basis set denoted here as maug-cc-pVDZ and constructed by augmenting the cc-pVDZ basis with s and p diffuse functions on the heavy atoms. Basis-set convergence tests are performed with the fully augmented aug-cc-pVTZ pseudopotential basis. In both cases, the diffuse functions are taken from the corresponding all-electron Dunning's correlation-consistent basis sets [START_REF] Kendall | Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions[END_REF]. The Jastrow factor includes two-body electron-electron and electron-nucleus correlation terms [START_REF]of two fifth-order polynomials of the electron-nuclear and the electronelectron distances, respectively, and rescale the inter-particle distances as R = (1 -exp(-κr))[END_REF].

We optimize all wave function parameters (Jastrow, orbital, and linear coefficients) by energy minimization in VMC using the stochastic reconfiguration (SR) method [START_REF] Sorella | Weak binding between two aromatic rings: Feeling the van der Waals attraction by quantum Monte Carlo methods[END_REF] in a conjugate gradient implementation [START_REF] Neuscamman | Optimizing large parameter sets in variational quantum Monte Carlo[END_REF]. We optimize the ground and excited states separately since the two states have different symmetries both at the ground-and the excited-state optimal structures. We relax the geometry in Z-matrix coordinates and simply follow the direction of steepest descent, appropriately rescaling the interatomic forces and using an approximate constant diagonal hessian. After convergence, we perform additional optimization steps and average the last 40 structures to estimate the structural parameters presented below. To remove occasional spikes in the forces, we use an improved estimator of the forces obtained by sampling the square of a modified wave function close to the nodes [START_REF] Attaccalite | Stable liquid hydrogen at high pressure by a novel ab initio molecular-dynamics calculation[END_REF]. In the DMC calculations, we treat the pseudopotentials beyond the locality approximation using the T-move algorithm [START_REF] Casula | Beyond the locality approximation in the standard diffusion Monte Carlo method[END_REF] and employ an imaginary time step of 0.02 a.u. As shown in the Supporting Information (SI), this time step yields DMC excitation energies converged to better than 0.01 eV for the smallest wave function in formaldehyde and is therefore appropriate for all wave functions of higher quality considered in this work.

The CIPSI calculations are performed with Quantum Package [START_REF] Scemama | Quantum Package v1.0[END_REF] using orbitals obtained from small complete active space self-consistent field (CASSCF) calculations in the program GAMESS(US) [START_REF] Schmidt | General atomic and molecular electronic structure system[END_REF][START_REF] Gordon | Theory and applications of computational chemistry[END_REF]. As explained above, the CIPSI expansions are constructed to be eigenstates of Ŝ2 and the selected determinants are subsequently mapped into the basis of configuration state functions, thereby effectively reducing the number of linear optimization parameters for QMC.

IV. VERTICAL EXCITATION ENERGIES

We begin our investigation by computing the lowest singlet vertical excitations energies of n → π * character of formaldehyde (CH 2 O) and thioformaldehyde (CH 2 S). In the ground state, both molecules posses C 2v symmetry and the relevant ground (S0) and excited (S1) states belong to the A 1 and the A 2 irreducible representation, respectively. The QMC vertical excitation energies are computed on the ground-state structures optimized at the CC3/aug-cc-pVTZ level without the frozen-core approximation [START_REF] Scemama | Excitation energies from diffusion Monte Carlo using selected configuration interaction nodes[END_REF][START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF]. Importantly, as we show below, the CC3 geometries are identical to the optimal ground-state structures obtained at the VMC level.

A. Formaldehyde

The lowest n → π * excited state of formaldehyde has been the subject of a recent DMC investigation [START_REF] Scemama | Excitation energies from diffusion Monte Carlo using selected configuration interaction nodes[END_REF], where Jastrow-free CIPSI wave functions with as many as 300,000 determinants were employed to recover a vertical excitation energy within 0.08(3) eV of the best theoretical estimate of 3.97 eV [START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF]. In fact, calculations involving about 9000 and 75,000 determinants yielded excitation energies about 0.2 and 0.1 eV higher, respectively, than the reference value. Given the size of the molecule and the single-reference character of the excitation, this finding is rather surprising and warrants an investigation where we generate optimal wave functions in the presence of the Jastrow factor prior to the DMC step.

Here, we revisit the VMC and DMC computation of this vertical excitation energy with compact CIPSI expansions in the Jastrow-Slater wave functions, containing between 1000 and 40,000 determinants. These determinant components are the result of subsequent expansions at the CIPSI level constructed to achieve a balanced description of the two states of interest at each step by selecting determinants that yield comparable PT2 contributions for both states as discussed above. The determinants corresponding to A 1 and A 2 symmetry are then isolated and the Jastrow-Slater wave functions are separately fully optimized by energy minimization in VMC. The resulting VMC and DMC total and vertical excitation energies are listed in Table I and plotted in Fig. 1.

The VMC vertical excitation energies computed with the maug-cc-pVDZ basis are already within 0.05 eV of the CC3/aug-cc-pVTZ value with about 1000 and 2400 determinants in the ground and the excited state, respectively. Further increasing the size of the wave functions by a factor of 10, we observe a rapid convergence to within ∼ 0.02 eV of the CC3 value. The DMC calculations with the VMC optimized wave functions gain about 17-30 mHartrees for the range of expansions studied, so the gain is uniform for both states. Consequently, the gap is quite stable also at the DMC level and consistently compatible within statistical error (to better than 0.02 eV) with the corresponding VMC values. A similar behavior is observed when employing the fully-augmented triple-ζ basis set with the VMC excitation energy being rather stable as a function of the expansion size and in agreement within statistical error with the corresponding DMC value. From these calculations, we can therefore make three important observations: a) the weighted CIPSI algorithm used here yields an automated, balanced determinant selection of the two states; b) in combination with VMC optimization, we obtain reliable estimates for the vertical excitation energy with rather small and compact Jastrow-Slater wave functions; c) to estimate the excitation energy, there is no need to perform a DMC calculation as the main correlation effects on the energy difference have already been captured at the VMC level.

The extrapolated FCI estimate of the excitation energy computed with the same pseudopotential maug-cc-pVDZ basis set is 3.99 eV (see Fig. S1 in the Supporting Information) in perfect agreement with the all-electron value of 3.99 eV obtained with the corresponding aug-cc-pVDZ basis set. In the all-electron calculations, the use of the larger aug-cc-pVTZ only reduces the FCI value to 3.98 eV, that can be further corrected for basis set and frozen-core errors to yield the best theoretical estimate of 3.97 eV [START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF]. Our results of 3.99 and 4.02 eV obtained with the maug-cc-pVDZ and aug-cc-pVTZ, respectively, are in excellent agreement with the refer- ence value.

So far, we have obtained a balanced description of both states by selecting determinants in the CIPSI expansion that yield comparable net PT2 contribution to both states at every step. In Fig. 2, we plot the resulting excitations energies together with the estimates obtained as difference of the fits of the ground-and excited-state VMC energies against the number of total determinants as E fit EX (N det ) -E fit GS (N det ) (see Section SI.III). Although we optimize our wave functions separately for the two states in VMC, we use the total number of determinants as common index since the determinants were obtained from the same CIPSI selection [65]. Furthermore, since we now have at hand multiple wave functions for the ground and excited states, we can alternatively follow another physically appropriate approach to compute the excitation energy by matching the variances of the states of interest [START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF][START_REF] Robinson | Excitation variance matching with limited configuration interaction expansions in variational Monte Carlo[END_REF]. To do so, we fit the VMC energies separately for the ground and the excited states against the corresponding variances (see Fig. S3) and then subtract the two fits as

E fit EX (σ 2 ) -E fit GS (σ 2
). The resulting variance-matched excitation energies are plotted in Fig. 2, where smaller variances correspond to larger expansions.

For our fully optimized Jastrow-CIPSI wave functions, both schemes to estimate the excitation energy yield reasonable values, which become of course more compatible for the larger expansions. The excitation energies computed by fitting the energies against the number of determinants display a clear and fast convergence, and are closer to the extrapolated FCI estimate over the whole range of expansions explored here. For the variance-matching scheme, the fit at high variances deviate from the FCI value by 0.08-0.12 eV and, only beyond about 7000 determinants and a variance of about 0.31 a.u, the fit starts giving reasonable estimates within 0.05 eV of our best estimate. Therefore, to make a reliable variancematched prediction, one needs to employ larger wave functions containing several thousand determinants. We note that we are not in a regime [66] where a reliable estimate of the exact energy can be obtained by linear extrapolation to zero variance.

Finally, while the "expansion" scheme represents the most practical route to generate the determinantal component in a Jastrow-CIPSI wave function for larger systems, given the small size of formaldehyde, we can also investigate the effect of using the "truncation" scheme where a large CIPSI expansion is truncated either to yield wave functions with a similar norm for all states (Table S1, truncation from 40,000 total determinants) or to include the most relevant patterns (Table S2, truncation from 3 million total determinants) as discussed in Section II. We find that both truncated wave functions yield energetically equivalent estimates of the vertical excitation energies at the VMC and DMC levels as those obtained with comparable expansion sizes with the "expansion" scheme of Table I. We note that the "truncation" scheme should be applied to a relatively large expansion since, during the CIPSI iterations, the largest coefficients vary much at the beginning but tend to stabilize as the size of the expansion grows. If the size of the initial expansion is such that the determinants kept after truncation have converged coefficients relative to 2. VMC vertical excitation energies (circles) of formaldehyde obtained with the "expansion" scheme. The excitation energy is also estimated as the difference of the fits of the energies against the total number of determinants (bottom x-axis) and against the variance of the two states (top x-axis). We also show the extrapolated FCI and CC3 values obtained with the aug-cc-pVTZ basis set [START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF]. the largest coefficient, then truncating the wave function will be independent of the size of the starting expansion and, therefore, will be equivalent to truncating the FCI wave function. In general, fully reoptimizing the trial wave functions appears to yield accurate and robust estimates of the excitation energies of formaldehyde for both the "expansion" and the "truncation" scheme. This small molecule has been the subject of many recent investigations [START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF][START_REF] Loos | Theoretical 0-0 Energies with Chemical Accuracy[END_REF][START_REF] Wiberg | Comparison of CIS-and EOM-CCSD-calculated adiabatic excited-state structures. Changes in charge density on going to adiabatic excited states[END_REF][START_REF] De Mare | Ab initio investigation of the effect of geometry optimization and configuration interaction on the relative energies and structures of the ground and lowest excited states of thioformaldehyde[END_REF] and is a variant of formaldehyde, wherein the O atom is replaced by a S atom. Interestingly, the very recent study by Flores et. al [START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF] employed Jastrow-CIPSI wave functions to compute the excitation energy of the n → π * transition but obtained a VMC excitation energy of 2.07(2) eV, that is, about 0.2 eV lower than the reference extrapolated stochastic-heat-bath CI value of 2.31(1) eV on their geometry. They optimized all wave function parameters, albeit in variance minimization. In view of our success with the computation of the vertical excitation energy of formaldehyde at the VMC level, it is somewhat puzzling that the simple substitution of oxygen with sulfur would worsen so much the performance of the method. Therefore, we revisit here the same molecule using fully optimized Jastrow-CIPSI wave functions and both the CIPSI selection procedure of multiple states with similar PT2 energy correction and the variancematching scheme.

We summarize our QMC results obtained from the CIPSI "expansion" scheme and the maug-cc-pVDZ basis set in Table II and also present them in Fig. 3. As in the case of formaldehyde, we find that we are able to obtain a stable estimate of the vertical excitation energy both at the VMC and DMC level with relatively little effort. The excitation energy is essentially the same within statistical error at both the VMC and the DMC level when increasing the total number of determinants from 800 to 44,000. The use of larger expansions yields of course a gain in the total energies, which is however less than 1 and 0.5 mHartree in VMC and DMC, respectively, when doubling the expansion size between the last two entries of Table II. Computation of the DMC energies on top of the final VMC optimized wave functions uniformly lowers the energies of both states by about 18-20 mHartrees, thereby unaltering the energy separation between them. Therefore, the CIPSI selection scheme combined with wave function optimization ensures a balanced description of the two states and yields a VMC excitation energy of 2.23 eV in excellent agreement with the CC3 value of 2.23 eV and the extrapolated FCI estimate of 2.22 eV obtained with a larger aug-cc-pVTZ basis [START_REF] Loos | Theoretical 0-0 Energies with Chemical Accuracy[END_REF].

In Fig. 4, we plot the excitation energies obtained for the different CIPSI expansions together with the two estimates one obtains by computing the difference of the fits of the ground-and excited-state energies against the total number of determinants and against the variances of the two states as done for formaldehyde. We find that both estimates are compatible well within 0.05 eV over the whole range of expansion sizes/variances. We therefore do not reproduce the large error of 0.2 eV reported in Ref. [START_REF] Pineda Flores | Excited State Specific Multi-Slater Jastrow Wave Functions[END_REF], which cannot therefore be due to the use of the variance-matching recipe with reoptimized Jastrow-CIPSI wave functions.

V. OPTIMAL EXCITED-STATE STRUCTURES AND ADIABATIC EXCITATION ENERGIES

Upon relaxation of formaldehyde in the first singlet excited state, the oxygen atom moves out of the plane. The optimal excited-state structure therefore possesses C s symmetry with the symmetry plane being perpendicular to the initial molecular plane and passing through the CO bond. The excited-state optimization therefore follows the A state as the symmetry of the molecule is lowered from C 2v to C s . In the first singlet excited state, thioformaldehyde remains instead planar and preserves C 2v symmetry so that the excited state maintains A 2 character.

A. Formaldehyde

We first optimize the structure of formaldehyde in the ground state for some of the Jastrow-Slater CIPSI wave functions that we have used to compute the vertical excitation energies of Table I. Similar to our previous findings for butadiene [START_REF] Dash | Perturbatively Selected Configuration-Interaction Wave Functions for Efficient Geometry Optimization in Quantum Monte Carlo[END_REF], we obtain an accurate geometry in the ground state with relatively short expansions. As shown in Table III, the optimal bond lengths and bond angles computed with the compact aug-cc-pVDZ basis are within 1 mÅ and 1 • of the corresponding CC3/aug-cc-pVTZ values already with the smallest wave function considered. Energetically, we gain about 1 mHartree upon structure optimization with the smallest wave function of 580 determinants and the average VMC energies of the two larger cases are identical to the initial en- FIG. 3. Convergence of the VMC and DMC energies with the number of CIPSI determinants of the ground and excited states (top) and of the excitation energy (bottom) of the n → π * excitation of thioformaldehyde. The maug-cc-pVDZ basis is used. We also show the extrapolated FCI and CC3 values obtained with the aug-cc-pVTZ basis set [START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF].

ergies on top of the starting CC3 geometry.

To relax the geometry in the excited state, we can simply start from the ground-state structure, with the oxygen slightly displaced out of plane, and use the same excited-state expansion as in the calculation of the vertical excitation. Even though such a wave function misses all determinants of B 2 TABLE III. Optimal VMC ground-state bond lengths (Å) and bond angles (deg) of formaldehyde obtained using Jastrow-CIPSI wave functions with the maug-cc-pVDZ basis. The VMC energy (a.u.) is averaged over the last 40 iteration of geometry optimization. We also list the CC3/aug-cc-pVTZ and the experimental values. 4. VMC vertical excitation energies (circles) of thioformaldehyde obtained with the "expansion" scheme. The excitation energy is also estimated as the difference of the fits of the energies against the total number of determinants (bottom x-axis) and against the variance of the two states (top, x-axis). We also show the extrapolated FCI and CC3 values obtained with the aug-cc-pVTZ basis set [START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF]. character which will acquire non-zero weight as the molecule moves out of plane, this quick procedure gives us already a good estimate of the excited-state geometry with a CO bond of about 1.334 Å and an out-of-plane angle of 31.5 • (see Table S4). These structural parameters compare very favorably with the corresponding values obtained in CC3 and a full valence CASPT2 calculation [START_REF] Budzák | Accurate Excited-State Geometries: A CASPT2 and Coupled-Cluster Reference Database for Small Molecules[END_REF]. We stress that the excited-state geometry of formaldehyde has been found to depend significantly on the level of theory, with highly-correlated methods spanning differences of as much as 0.08 Å in the CO bond length and 20 degrees in the out-of-plane angle subtended by this bond on the HCH plane [START_REF] Budzák | Accurate Excited-State Geometries: A CASPT2 and Coupled-Cluster Reference Database for Small Molecules[END_REF].

No

Starting from such an out-of-plane excited-state geometry, we generate a series of CIPSI expansions of symmetry A and continue the structural relaxation in the excited state to further investigate the dependence of the resultant geometry on the complexity of the wave function. Since we also want to compute the adiabatic excitation energy, we have to evaluate differences of the ground/excited-state energies on the ground/excited-state optimal geometries. To obtain a balanced description, we therefore generate the CIPSI wave functions at different geometries either targeting a given value of the PT2 correction or of the CI variance for all four states. For convenience, we shall hereafter refer to them as the iso-PT2 and the iso-variance procedure, respectively. The results of these calculations are shown in Table IV.

We perform four iso-PT2 optimization tests. Starting from a VMC geometry of Table III, we generate four ground-state wave functions targeting different values of PT2 correction and further optimize the structure with these wave functions. We follow a similar procedure for the excited state starting from an out-of-plane VMC excited-state geometry of Table S4. For the ground state, the resulting bond lengths and an- gles remain within 1 mÅ and 1 • , respectively, of the CC3 values for all four wave functions. In the excited state, the CO bond undergoes a marginal lengthening and the out-of-pane angle slightly increases by 2-3 • becoming closer to the CC3 angle. The adiabatic excitation energies are computed as differences between iso-PT2 ground and excited-state energies on the corresponding equilibrium structures, and are highly stable irrespective of the size of the wave function: the VMC values lie between 3.55-3.56 eV and the DMC ones are somewhat higher and compatible within statistical error with the CC3 estimate. The iso-variance tests yield very similar results as the iso-PT2 ones and the procedure is therefore also a viable route to compute excited-state structures and adiabatic energies.

B. Thioformaldehyde

Also for thioformaldehyde, we investigate the VMC convergence of the ground-state optimization with three groundstate wave functions used in the computation of the vertical excitation energies. As shown in Table VI, we are able to obtain converged bond lengths and bond angles within less than 3 mÅ and 1 • of the CC3 values already with the smallest set of 702 determinants and a maug-cc-pVDZ basis. Energetically, we gain much less than 1 mHartee in all three cases upon optimization starting from the CC3 structure.

In the first singlet excited state, unlike its oxygen counterpart, thioformaldehyde remains planar and the excited state does not change character. We therefore start by simply relaxing the structure with one of the excited-state wave functions employed in the estimation of the vertical excitation energy. To proceed in the structural optimization, we then generate new CIPSI expansions following the iso-PT2 scheme and matching the PT2 corrections used in the corresponding ground-state calculations. As shown in Table VI, despite testing a large range of expansions, there is no significant variation in the length of the CS bond, which is 10 mÅ longer than the CC3 value. The adiabatic excitation energy is again rather stable across all cases and within 0.01-0.04 eV of the CC3 estimate.

VI. DISCUSSION

We have demonstrated the excellent performance of QMC in the accurate computation of the vertical excitation energies of two small but theoretically challenging systems, formaldehyde and thioformaldehyde. Using fully optimized Jastrow-CIPSI wave functions where the determinantal components are constructed to yield a similar PT2 correction for both states, we are able to obtain VMC excitation energies compatible within less than 0.02 eV of the extrapolated FCI values using compact expansions of a few thousand determinants Next, we have investigated the ability of QMC to obtain accurate ground-and excited-state structures. In the ground state, we easily obtain geometries in excellent agreement with those produced by other high-level correlated approaches such as CC3 and full valence CASPT2, using very compact Jastrow-CIPSI wave functions. The maximum deviation is a meager 3 mÅ in the CS bond of thioformaldehyde. During the relaxation of the excited-state structure, we regenerate the determinantal component in our Jastrow-CIPSI wave functions following two different determinant selection schemes: i) keeping a roughly constant PT2 energy correction during the optimization and ii) targeting a fixed value of the CI energy variance. Like the vertical excitation energies, the excitedstate structural parameters show relatively low sensitivity to the size of the wave function and both selection schemes yield similar structures with no clear distinction in the convergence properties. The largest deviation with respect to the CC3 bond lengths is of about 10 mÅ when using the double-ζ basis. Geometrical relaxation of the excited state additionally allows us to estimate the adiabatic excitation energies, which we find to be compatible within less than 0.05 eV of the CC3 values for all wave function sizes.

To summarize, the use of a CIPSI selection scheme targeting similar PT2 corrections for the states of interest in combination with the full optimization of the Jastrow-Slater wave function enables us to obtain extremely stable and accurate estimates of the vertical excitation energies already at the VMC level, namely, without the need of DMC. This can be achieved with very compact wave functions, demonstrating the accuracy of the "expansion" scheme, which is viable also for larger systems. Furthermore, our iso-PT2 protocol to regenerate the wave function during relaxation along a potential energy surface leads to consistent and high-quality structures for both ground and excited state and to accurate adiabatic excitation energies. Our robust estimates of the structural parameters and the vertical and adiabatic excitation energies with moderate Slater expansions open important prospects for the use of VMC structural optimization with Jastrow-CIPSI wave functions as an efficient and reliable approach to characterize excited-state potential energy surfaces.

ACKNOWLEDGEMENT

This work is part of the Industrial Partnership Programme (IPP) "Computational sciences for energy research" of the Netherlands Organisation for Scientific Research (NWO-I, formerly FOM). This research programme is co-financed by Shell Global Solutions International B.V. This work was carried out on the Dutch national supercomputer Cartesius with the support of SURF Cooperative, and using HPC resources from CALMIP (Toulouse) under allocation 2018-0510. J.F. acknowledges the Deutsche Forschungsgemeinschaft (DFG) for financial support (Grant FE 1898/1-1). The authors thank Dr. Pierre-Francois Loos for useful discussion and data on the coupled cluster geometries. The authors declare no competing financial interest. SUPPINFO CIPSI energies and FCI/maug-cc-pVDZ extrapolation for formaldehyde and thioformaldehyde; performance of "truncation" scheme in VMC; fits of ground-and excited-state VMC energies as a function of the total number of determinants and

FIG. 1 .

 1 FIG.1. Convergence of the VMC and DMC energies of the ground and excited states (top) and of the excitation energy (bottom) of the n → π * excitation of formaldehyde with the combined number of determinants in the CIPSI expansions. The maug-cc-pVDZ basis is used. We also show the extrapolated FCI and CC3 values obtained with the aug-cc-pVTZ basis set[START_REF] Loos | A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks[END_REF].

TABLE I .

 I VMC and DMC ground-and excited-state energies, and vertical excitation energies (eV) of formaldehyde obtained with fully optimized Jastrow-CIPSI wave functions, where a series of increasing determinantal expansions are obtained with the "expansion" scheme.

	Basis	No. det	No. param		VMC			DMC
		S0	S1 S0	S1	E(S0)	E(S1)	∆E	E(S0)	E(S1)	∆E
	maug-cc-pVDZ 343 436 848 614 -22.88852(32) -22.74015(33) 4.037(13) -22.91702(26) -22.76903(26) 4.027(10)
		580 1124 946 961 -22.89762(32) -22.74929(31) 4.036(12) -22.92182(25) -22.77393(24) 4.024(09)
		994 2360 1104 1357 -22.90177(30) -22.75459(27) 4.005(11) -22.92362(23) -22.77635(23) 4.008(09)
		1703 4182 1375 1937 -22.90512(25) -22.75753(25) 4.016(10) -22.92530(22) -22.77772(22) 4.016(08)
		2747 7110 1762 2805 -22.90717(25) -22.76046(24) 3.992(09) -22.92675(22) -22.77947(22) 4.008(08)
		3050 8320 1874 3141 -22.90719(20) -22.76051(20) 3.991(08) -22.92642(18) -22.77916(17) 4.007(07)
		5932 16871 2915 5550 -22.90852(20) -22.76233(20) 3.978(08) -22.92687(18) -22.77997(17) 3.997(07)
		10854 29786 4681 9177 -22.90961(20) -22.76303(20) 3.989(08) -22.92728(14) -22.78072(16) 3.988(06)
	aug-cc-pVTZ	675 912 3237 2256 -22.90151(19) -22.75361(19) 4.024(07) -22.92260(16) -22.77499(16) 4.017(06)
		1488 3214 3592 3474 -22.91364(17) -22.76638(17) 4.007(06) -22.92778(14) -22.78068(13) 4.015(05)
		3058 7672 4210 5050 -22.92021(15) -22.77242(15) 4.021(06) -22.93112(12) -22.78359(12) 4.014(05)
		5849 15338 5222 7401 -22.92050(16) -22.77333(14) 4.005(06) -22.93136(11) -22.78360(10) 4.021(04)
		12987 31710 7731 12040 -22.92188(15) -22.77466(14) 4.006(06) -22.93145(11) -22.78393(13) 4.014(05)

TABLE II .

 II VMC and DMC ground-and excited-state energies, and vertical excitation energies (eV) of thioformaldehyde obtained with fully optimized Jastrow-CIPSI wave functions obtained with the "expansion" scheme. The maug-cc-pVDZ basis set is used.

	No. det	No. param		VMC			DMC
	S0	S1	S0 S1	E(S0)	E(S1)	∆E	E(S0)	E(S1)	∆E
	353 475 854 607 -17.04424(26) -16.96212(25) 2.234(10) -17.06968(25) -16.98671(24) 2.258(9)
	702 1488 995 1044 -17.04998(25) -16.96814(24) 2.227(09) -17.07246(23) -16.99023(23) 2.237(9)
	1165 2702 1168 1461 -17.05239(24) -16.97019(24) 2.237(09) -17.07284(22) -16.99086(22) 2.231(9)
	1834 4692 1419 2034 -17.05388(23) -16.97100(23) 2.255(09) -17.07394(22) -16.99171(22) 2.237(9)
	2500 6316 1662 2514 -17.05494(23) -16.97260(23) 2.241(09) -17.07434(21) -16.99211(22) 2.238(8)
	3432 8338 1997 3112 -17.05622(23) -16.97313(23) 2.261(09) -17.07482(22) -16.99219(22) 2.248(8)
	5712 14562 2810 4848 -17.05698(22) -16.97497(22) 2.232(08) -17.07542(22) -16.99314(26) 2.239(9)
	14218 30142 5650 9162 -17.05780(15) -16.97590(19) 2.236(07) -17.07566(18) -16.99345(26) 2.237(8)
		B. Thioformaldehyde

TABLE IV .

 IV Optimal VMC ground-and excited-state bond lengths (Å) and bond angles (deg), and VMC and DMC adiabatic excitation energies (eV) of formaldehyde obtained using Jastrow-CIPSI wave functions with targeted PT2 correction (iso-PT2). The VMC energy (a.u.) is averaged over the last 40 iterations of geometry optimization. o.o.p denotes the out of plane angle of the CO bond.

	State δEPT2 No. det. No. par.	CO	CH	H-C-H	o.o.p	EVMC	∆E VMC adia ∆E DMC adia
	S0	0.20	515	1183 1.20663(05) 1.09851(03) 116.84(4) -0.098(25) -22.89892(5)
		0.15	1184	1496 1.20795(08) 1.09832(05) 116.96(6)	0.075(15) -22.90411(6)
		0.10	2784	2208 1.20752(19) 1.09947(05) 117.00(5)	0.073(54) -22.90723(5)
		0.07	4799	2946 1.20710(14) 1.09959(10) 117.01(8) -0.030(29) -22.90805(4)
		CC3 [42]		1.208	1.100	116.4	0.000
		Expt. [70]		1.208	1.116	116.3	0.000
	S1	0.20	1088	1363 1.33948(07) 1.08468(02) 119.83(03) 32.916(055) -22.76787(5) 3.566(2) 3.612(9)
		0.15	3106	2037 1.33971(36) 1.08630(06) 119.03(20) 34.778(051) -22.77342(5) 3.556(2) 3.614(8)
		0.10	8058	3520 1.33606(21) 1.08647(12) 119.34(17) 34.027(233) -22.77626(3) 3.564(2) 3.597(8)
		0.07 15278 5589 1.33597(15) 1.08645(08) 119.56(09) 33.294(070) -22.77746(3) 3.553(1) 3.592(8)
		CC3 [42, 67]		1.326	1.089	118.3	36.8	3.602
		Expt. [71-73]		1.321-1.323 1.092-1.103 118.1-121.5 20.5-34

TABLE V .

 V Same as Table IV using Jastrow-CIPSI wave functions with targeted CI variance (iso-variance).

	State σ 2 CI No. det. No. par.	CO	CH	H-C-H	o.o.p	EVMC	∆E VMC adia ∆E DMC adia
	S0 0.80 515	1183 1.20663(05) 1.09851(03) 116.84(4) -0.098(025) -22.89892(5)
	0.60 1459	1605 1.20761(14) 1.09886(04) 116.98(5) -0.091(023) -22.90489(6)
	0.40 3755	2566 1.20732(23) 1.09886(05) 116.99(4) -0.057(235) -22.90761(6)
	S1 0.80 782	1167 1.33888(11) 1.08562(06) 119.40(07) 34.061(104) -22.76474(6) 3.651(2) 3.652(8)
	0.60 3106	2037 1.33971(36) 1.08630(06) 119.03(20) 34.778(051) -22.77342(5) 3.577(2) 3.624(8)
	0.40 9820	4031 1.33713(31) 1.08641(06) 119.23(11) 33.522(079) -22.77682(6) 3.563(2) 3.611(8)

TABLE VI .

 VI Optimal VMC ground-and excited-state bond lengths (Å) and bond angles (deg), and VMC and DMC adiabatic excitation energies (eV) of thioformaldehyde obtained using Jastrow-CIPSI wave functions with targeted PT2 correction (iso-PT2). The VMC energy (a.u.) is averaged over the last 40 iterations of geometry optimization. o.o.p denotes the out of plane angle of the SO bond. augmented double-ζ basis. Performing DMC calculations on top of the fully optimized VMC wave functions leads to a uniform gain across both states, thereby not affecting the VMC estimate of the excitation. If we compute the excitation energy by matching instead wave functions with similar variances, we obtain a less robust procedure in the case of formaldehyde, where relatively large expansions of about 7000 determinants are required to recover results within 0.05 eV of our best estimates.

	State δEPT2 No. det. No. par.	CS	CH	H-C-H	o.o.p	EVMC	∆E VMC adia ∆E DMC adia
	S0	0.17	702	995 1.62218(14) 1.08427(04) 116.869(53) -0.042(50) -17.05044(4)
		0.10	2500	1662 1.62250(16) 1.08469(05) 116.862(32) 0.070(36) -17.05570(4)
		0.07	5712	2810 1.62200(11) 1.08464(06) 116.786(34) -0.058(46) -17.05715(6)
		0.045 14218 5650 1.62207(10) 1.08480(07) 116.756(26) -0.041(53) -17.05817(3)
		CC3 [42]		1.619	1.083	116.1	0.000
		Expt. [70, 74, 75]		1.611-1.614 1.093-1.096 116.2-116.9	0.000
	S1	0.17	858	764 1.72393(27) 1.07879(07) 121.309(163) -0.024(46) -16.97123(4) 2.155(2) 2.152(6)
		0.10	3982	1820 1.72514(04) 1.07895(04) 121.388(27) -0.106(79) -16.97809(5) 2.112(2) 2.131(6)
		0.07	8418	3110 1.72309(07) 1.07909(05) 121.149(26) 0.115(79) -16.97998(4) 2.100(2) 2.119(6)
		0.045 16866 5475 1.72182(09) 1.07921(03) 121.144(28) 0.041(85) -16.98090(3) 2.103(1) 2.112(6)
		CC3 [42, 67]		1.709	1.078	120.2	0.0	2.112
		Expt. [70, 74, 76, 77]	1.682-1.708 1.077-1.093 116.8-121.6	0.0-8.9
	and a minimally

as a function of the energy variance; DMC time-step extrapo-lation; intermediate excited-state geometry optimizations.


Nicolas Peyret
email: nicolas.peyret@supmeca.fr

Marco Rosatello
email: marco.rosatello@supmeca.fr

Gaël Chevallier
email: gael.chevallier@univ-fcomte.fr

Jean-Luc Dion
email: jean-luc.dion@supmeca.fr

A Mindlin derived Dahl friction model

Friction modeling is essential in a vast amount of engineering applications. In this paper Mindlin contact friction model, a Masing rule compliant model, is manipulated in order to find a correlation with Dahl dynamic friction model, a Masing rule not compliant model. The resulting equations are then numerically compared and the correlation demonstrated with several loading functions. The result is a modified Dahl model that captures the partial slip behavior following the contact theory. The Dahl model, classically used in control engineering as an empirical model due to its implementation facility, finds in this new formulation a physical meaning to its parameters and finally identified physical parameters (according from Mindlin model) for the Dahl model.

Introduction

Friction plays a fundamental role in a wide variety of physical systems and the research community around this topic is always very active. The creation of models to accurately describe the phenomenon of friction has always been one of the main objectives of this research field, but the underlying reasons were different. For example, in the tribology field, friction models are needed in order to evaluate the fatigue life of bodies in contact, while for other disciplines, such as control engineering, friction models are necessary in order to develop control strategies that are able to compensate for the friction behavior.

Friction models can be divided into two main categories: macro-models and constitutive models. Macro-models are based on empirical and experimental observations of the friction behavior and are divided in three main categories: quasi-static, dynamic and hysteretic models. Quasi-static models are the classical friction models, which assumes the friction force to be a function of the relative velocity of the contact surfaces. The simplest and also best-known quasi-static model is Coulomb's, from which other models were derived, in order to take into account phenomena such as viscous friction, Stribeck effect, stiction and kinetic friction.

Dynamic models are particularly used by control engineers and the first dynamic model was developed by Dahl [START_REF] Dahl | A Solid Friction Model[END_REF] , who derived it from the materials strain-rate equation. Further developments on Dahl's model brought to the creation of other models such as the Bliman-Sorine model [START_REF] Bliman | Friction modeling by hysteresis operators. application to dahl, sticktion and stribeck effects[END_REF] , the Haessig and Friedland bristle model [START_REF] Haessig | On the modeling and simulation of friction[END_REF] and the one conceived by Canudas de Wit and al. in [START_REF] De Wit | A new model for control of systems with friction[END_REF] , best known as the LuGre model [START_REF] Dion | A new multicontact tribometer for deterministic dynamic friction identification[END_REF][START_REF] Piatkowski | Dahl and lugre dynamic friction models -t h e analysis of selected properties[END_REF] . In automatic control domain, and especially in programmed handling activities, friction is an important physical phenomenon to master in order to ensure reliability and precision [START_REF] Piatkowski | Analysis of translational positionning of unit loads by directionally-oriented friction force fields[END_REF] .

Hysteretic models come from the theory of elasticity and they are used to model friction and energy dissipation in materials, geomechanics and jointed structures. Well-known hysteretic models are Masing model [START_REF] Masing | Zur heynschen theorie der verfestigung der metalle durch verborgen elastische spannungen[END_REF] and Iwan model [START_REF] Iwan | A distributed-element model for hysteresis and its steady-state dynamic response[END_REF] , which are based on distributions of Jenkins elements.

Differently from macro-models, the constitutive models are derived from Hertz contact theory (see [START_REF] Johnson | Contact Mechanics[END_REF] ), in order to describe the stress-strain relationship in the contact region. First Cattaneo [START_REF] Cattaneo | Sul contatto di due corpi elastici: distribuzione locale degli sforzi[END_REF] and then Mindlin [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF] , developed Hertz theory to take into account a tangential alternative load. This model is classically used in the study of fretting and wear in contact physics and in tribology domain.

The aim of this paper is, firstly, to find a correlation between the Mindlin model, which is a constitutive theoretical model that satisfies Masing rule [START_REF] Masing | Eigenspannungen und verfestigung beim messing[END_REF] , and Dahl's model, which instead is a dynamic empirical friction model that doesn't satisfy Masing rule, and secondly, the paper addresses to both Dahl and Mindlin model users: the goal is to provide physical parameters for Dahl model and an easier numerical implementation contact model.

In Section 2 the theory and the equations of the concerned friction models will be described . In Section 3 , Mindlin's equations will be manipulated in order to find the correlation with Dahl's model and finally, in Section 4 , a verification of the obtained correlation will be carried out. This correlation allows to give to the Dahl model, a new formulation with an analogy between the empirical model and the contact physics.

Theory

Mindlin contact model

In [START_REF] Cattaneo | Sul contatto di due corpi elastici: distribuzione locale degli sforzi[END_REF] , Cattaneo studied the problem of the elastic contact between two bodies pressed on each other by a constant normal force N 0 , with a tangential load T monotonically growing from zero to a maximum value T * . Mindlin extended the model for a general oscillatory force of amplitude ± T * and considered the contact between a sphere and a plane (see [START_REF] Johnson | Contact Mechanics[END_REF] ).

From Hertz theory the radius a of the contact area is defined as:

a = 3 3 4 N 0 ρ E * (1)
where

1 E * = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 (2)
and the generalized bending radius ρ is defined as:

1 ρ = 1 R 1 + 1 R 2 (3) being E 1 , E 2 the Young moduli; ν 2 1 , ν 2
2 the squares of the Poisson ratios and R 1 , R 2 the bending radii of the two surfaces in contact.

In order to describe the friction behavior at the contact surface, both Cattaneo and Mindlin used the Coulomb's law.

Considering µ the friction coefficient, if the tangential load T reaches the critical value

T c = µN 0 (4)
the two surfaces starts to slide relatively to each other entering in what is usually called the sliding or macro-slip phase, during which the tangential load T will remain constant.

During the monotonic tangential loading phase, Mindlin and Cattaneo demonstrated that the contact area is divided into two parts: a stick zone and a slip zone. The stick zone is a circle of radius c , while the slip zone is the annulus of internal radius c and external radius a . In Fig. 1 b is possible to see the two regions, while in Fig. 1 a are shown the normal and tangential stress distribution at the contact interface. The relationship between the stick zone radius c and the contact radius a was found by Cattaneo to be:

c a = 1 - T µN 0 1 3 (5)
It's easy to see that, when the tangential load T increases, the stick region becomes smaller, until T = T c , when c = 0 and the macro-slip phase takes place. Therefore, when 0 < T < T c , the contact is in a phase called partial slip. If we consider A c as the contact area of radius a and A s as the slip annular area of internal radius c and external radius a , some researchers prefer to divide this phase into a micro-slip phase, where A s < < A c and a partial slip phase, where A s < A c . We now consider an oscillating loading cycle of amplitude ± T * . The hysteresis cycle has the shape of the one shown in Fig. 2 , with the tangential displacement δ on the x-axis and the tangential force T on the y-axis. The cycle can be divided in three phases: initial loading from point O to point A (phase 0), unloading and reloading from point A to point C (phase 1), unloading and reloading back from point C to point A (phase 2). According to Mindlin model in [10] , it's possible to write the tangential displacement δ as a of the applied tangential force T . During the initial loading the tangential displacement δ 0 can be written as:

δ 0 = 3 µN 0 G * 16 a 1 -1 - T µN 0 2 3 ( 6
)
where

G * = 2 -ν 1 G 1 + 2 -ν 2 G 2 (7)
in which G 1 and G 2 are the shear moduli of the materials in contact. During the unloading and reloading phases, the respective tangential displacements δ 1 and δ 2 can be written as:

δ 1 = 3 µN 0 G * 16 a 2 1 - T * -T 2 µN 0 2 3 -1 - T * µN 0 2 3 -1 Phase 1 (8) δ 2 = - 3 µN 0 G * 16 a 2 1 - T * + T 2 µN 0 2 3 -1 - T * µN 0 2 3 -1 Phase 2 (9)
It's important to point out that, if the tangential force T reaches the critical value T c = µN 0 , according to Coulomb's friction model, it will remain constant until the load changes direction.

Dahl friction model

Dahl friction model [START_REF] Dahl | A Solid Friction Model[END_REF] was introduced in 1968 to study the rolling friction in ball bearings. Dahl's main idea was to use the materials stress-strain curve of classical solid mechanics in order to model friction. Considering T to be the friction force and T c = µN its critical value according to Coulomb's law, Dahl's model general form is:

dT dx = σ 1 - T T c sign ( ˙ x ) α ( 10
)
where σ is the stiffness coefficient, which represents the slope of the curve at T = 0 , while α is a parameter that determines the shape of the curve and it's usually taken as α = 1 . Furthermore, the force | T | will never be larger than the critical value

T c .
It's possible to see that in this model the friction force is only function of the displacement x and of the sign of the velocity ˙

x .

Masing models

Masing models are simple tools for describing the hysteretic behavior of many systems. The model was firstly introduced to model the Bauschinger effect in materials (see [START_REF] Skelton | The bauschinger effect, masing model and the ramberg-osgood relation for cyclic deformation in metals[END_REF] ), which is a reduction in compression strength after a plastic deformation in tension.

Masing [START_REF] Masing | Eigenspannungen und verfestigung beim messing[END_REF] was the first to create a model for elasto-plastic material during cycle loading. His theory provides also for the so called Masing rule or hypothesis, which says that the hysteresis paths during cycle loading can be obtained from the monotonic initial loading, also called "backbone" curve, with an expansion factor of two.

By keeping the notation shown in Fig. 2 the mathematical representation of Masing rule is the following:

δ 1 (T ) = δ * -2 δ 0 T * -T 2 ( 11
)
where δ * represents the maximum tangential displacement amplitude reached in the last cycle. A practical consequence of this rule, which is also a material property during cycle loading, is that the slope right after the load inversion has the same value of the initial slope during the monotonic initial loading phase. It's important to clarify that the original Masing hypothesis is only valid for steady-state cycles or loading between fixed limits. At this point it's possible to check if the previously described Mindlin and Dahl models satisfy Masing rule. By applying Eqs. ( 11)-( 13) it's possible to obtain Eq. ( 14) , proving that Mindlin contact model is a Masing model.

If now we consider Dahl model, since Eq. ( 10) represents the instantaneous slope of the hysteresis curve in the plane x -F , it's easy to verify that:

dT dx (0 , ˙ x > 0) = σ = dT dx (T * , ˙ x < 0) = σ 1 + T * T c α (12)
which means that Dahl model does not satisfy Masing hypothesis. The Masing assumption allows to deduce the full hysteresis cycle only by knowing the first quarter of it, instead of having to compute the complete and closed cycle. For example, this is particularly useful when the dissipated energy per cycle has to be calculated. If the model satisfies Masing hypothesis, you can just compute the dissipated energy during the initial loading and multiply this value by four to obtain the dissipated energy in the closed hysteresis cycle.

Correlation between Mindlin and Dahl models

After having described the two models in the previous sections, a correlation procedure is now shown. The first step is to find the expression of the tangential force T as a function of the displacement δ from Eqs. ( 6) to (9) .

The equations of the tangential force T during the different phases of the hysteresis cycle can then be written as:

T 0 = µN 0 1 -1 - 16 a 3 µN 0 G * δ 0 3 2 (13) T 1 = T * -2 µN 0 1 - 1 2 16 a 3 µN 0 G * δ + 1 - T * µN 0 2 3 + 1 3 2 (14) T 2 = -T * + 2 µN 0 1 - 1 2 16 a 3 µN 0 G * δ + 1 - T * µN 0 2 3 + 1 3 2 (15)
At this point the derivatives of Eqs. ( 13) -( 15) with respect to the variable δ are found as:

dT 0 dδ 0 = 8 a G * 1 - 16 a 3 µN 0 G * δ 1 2
(16)

dT 1 dδ = 8 a √ 2 G * 1 + 16 a 3 µN 0 G * δ + 1 - T * µN 0 2 3 1 2
(17)

dT 2 dδ = 8 a √ 2 G * 1 - 16 a 3 µN 0 G * δ + 1 - T * µN 0 2 3 1 2 (18)
Finally Eqs. ( 6) -( 9) are replaced respectively in Eqs. ( 16) -( 18) , obtaining:

dT dδ 0 = 8 a G * 1 - T T c 1 3
Phase 0 (19)

T * + T sign ( ˙ δ) 2 T c 1 3 Phase 1 & Phase 2 ( 20
)
If we analyze these last two equations, it's possible to see that they are equivalent to the original Dahl friction model of Eq. (10) with the two parameters σ and α chosen as(See Peyret in [START_REF] Peyret | Dissipation de l'énergie mécanique dans les assemblages: effet du frottement en sollicitation dynamique[END_REF] :

σ = 8 a G * α = 1 3 (21)
Therefore, Eqs. ( 19) -( 20) represent the Mindlin-Dahl correlated model.

Results and discussion

In order to validate the model and compare the model, a tangential displacement

δ = λδ c sin (2 πt) (22)
is imposed to the model, in which δ c is the critical displacement at which the friction force becomes equal to the Coulomb friction and macro-slip occurs:

δ c = 3 µNG * 12 a ( 23
)
where λ is an dimensionless coefficient. If λ ≥ 1, macro-slip occurs. The values of the parameters set used in simulations are, for Dahl's model: σ = 4 . 839 e + 08 , α = 1 / 3 ; for Mindlin's model :

µ = 0 . 3 , N = 10 0 [N], a = 0 . 0027 [m], G * = 4 . 42 e -11
[1/Pa] Fig. 3 shows the two excitation functions applied to the models: in Fig. 3 a the imposed displacement remains inside the micro-slip regime, while in Fig. 3 b macro-slip occurs.

Concerning the numerical resolution of the models, for Mindlin's model there are no particular problems, while for the correlated Dahl model, the derivation in respect of time must be found. The latter is obtained as:

dT dt = dT dδ dδ dt = σ ˙ δ 1 - T T c α (24)
where ˙ δ = 2 πλδ c cos (2 πt)

At this point the displacement functions shown in Fig. 3 , one for micro-slip with λ = 0 . 8 and one for macro-slip with λ = 1 . 2 , are applied to three different models. The first model is, of course, Mindlin's model, represented by Eqs. ( 13) - [START_REF] Peyret | Dissipation de l'énergie mécanique dans les assemblages: effet du frottement en sollicitation dynamique[END_REF] .

The second is Dahl's correlated model integrated in time and represented by the following equations: ˙

˙ T = σ ˙ δ 1 - T T c α Phase 0 (26)
T = σ ˙ δ 1 - T * + T sign ( ˙ δ) 2 T c α Phase 1 & Phase 2 (27)
This model will be referred as Full-Dahl . Ultimately, the third model is again Dahl's correlated model, but this time ignoring Eq. ( 26) for the initial loading, thus using only Eq. ( 27) . In order to distinguish it from the Full-Dahl , this model will be referred as Simplified-Dahl .

In the model comparison process, firstly only one full hysteresis cycle as represented in Fig. 2 is considered: Fig. 4 a shows the tangential friction force as a function of the displacement for the three models previously described, while Fig. 4 b shows the percentage error of the correlated Dahl's models relative to the original Mindlin model. Full-Dahl model correctly follows Mindlin's model, with a percentage error always smaller than 2%. Concerning the Simplified-Dahl model, after the initial loading phase which obviously produces a consistent error, the latter decreases throughout the unloading and reloading phases, but it remains unacceptable.

However, if we consider the experiments in the tribology and material fatigue domain, the number of cycles considered is high. Fig. 5 shows the behavior of the models for a higher number of cycles: it's possible to notice that, after about 3 cycles, the percentage error of the Simplified-Dahl model reaches the level of the Full-Dahl one, thus it becomes acceptable, being always under 2%.

The results discussed until now are referred to the imposed displacement function shown in Fig. 3 a, which implies an hysteresis cycle entirely inside the micro-slip domain. In order to obtain the behavior during the macro-slip phase, the imposed displacement function of Fig. 3 b is applied, and the results are shown in Fig. 6 , where the tangential effort and the percentage error are calculated over one full cycle. It's possible to see that, for the loading phase, the behavior of the different models is the same, with a consistent error for the Simplified-Dahl model. After this phase, however, the behavior changes, assisting to a fast decrease in the Simplified-Dahl error, which becomes even lower than the Full-Dahl one: this aspect can be simply due to numerical errors. Thus, in this case, the error of the simplified model becomes acceptable right after the loading phase, and not after a few cycles as it was shown for the micro-slip case of Fig. 5 a.

After having obtained good results with a simple sinus excitation, which is also one of the condition imposed by Mindlin, the exponential growing sinus represented in Fig. 7 , is applied to the models. It's possible to see that at the beginning the function is within the micro-slip and after some cycles enters in the macro-slip phase. Since the amplitude δ * of the excitation is not constant, it was calculated at each step as the envelope of the loading function, represented by the green dot lines in Fig. 7 . The results of the application of this particular excitation produces the results shown in Fig. 8 : it's possible to say that there is a good matching even if the assumptions of Mindlin model are not respected(steady sine excitation of amplitude T * ). The percentage error shows the behavior already seen for the simple sinus excitation, but in this case there is a residual varying error that probably depends on the way the maximum amplitude δ * is calculated.

Conclusions

In the paper, two well-known friction models, belonging to two different categories of models, were presented. It was demonstrated that, by mathematically manipulating the Mindlin's model equations, it's possible to find an alternative formulation of the Dahl dynamic friction model. The models were then validated through the application of different excitation functions and the error in respect to Mindlin's model was evaluated.

The interesting aspects of the obtained results are several: firstly, considered that Mindlin's model observes the Masing rule, while this is not satisfied by Dahl's model, the correlation was able to find a Masing rule compliant Dahl model. Secondly, Mindlin's model is composed by three different equations, each one for a different phase of the hysteresis cycle: The Full-Dahl model is also composed by three differential equations, but it was shown to be precise enough with only two (Simplified Dahl). The latter aspect implies an easier numerical implementation and it can be used as a friction model, not only in tribology, but also in other contexts, when it's necessary to accurately represent accurately the micro-slip phase, and probably most importantly, this work allowed to connect Dahl model, which began as an empirical global behavior law, to the contact physics.

Finally, on one hand, Dahl model is an easy-to-use model with only 3 empirical parameters, without any physical meaning. On the other hand, Mindlin model is realistic and based on clear physical parameters, but these are numerous, making it still difficult to use with reasonable computation costs. Our work allows to exploit the advantages of both models, while avoiding their drawbacks. The bottom line is that we have linked the Dahl parameters to the Mindlin physical ones, which allows us to use the easy-to-compute Dahl model with parameters that are linked to the actual physics of parts and surfaces in contact.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Stress distribution at the contact interface (a) and contact area (b) for a contact between a sphere and a plane with a tangential force T < T c .

Fig. 3 .

 3 Fig. 3. Imposed displacement excitation.

Fig. 4 .

 4 Fig. 4. Comparison Mindlin, Full-Dahl and Simplified-Dahl during a full cycle of micro-slip excitation.

Fig. 5 .

 5 Fig. 5. Comparison between Mindlin, Full-Dahl and Simplified-Dahl during 5 cycles of micro-slip excitation.

Fig. 6 .

 6 Fig. 6. Comparison between Mindlin, Full-Dahl and Simplified-Dahl for one cycle in a macro-slip situation.

Fig. 7 .

 7 Fig. 7. Exponential growing sinus excitation.

Fig. 8 .

 8 Fig. 8. Comparison between Mindlin, Full-Dahl and Simplified-Dahl during 5 cycles of exponential excitation.