Full wave function optimization with quantum Monte Carlo and its effect on the dissociation energy of FeS

TL;DR

This study employs advanced quantum Monte Carlo methods with orbital and guide function optimization to accurately predict the dissociation energy of FeS, revealing the importance of orbital optimization and Jastrow factors in transition metal compounds.

Contribution

First implementation of quantum Monte Carlo orbital optimization for transition metal molecules, improving agreement with experimental dissociation energies.

Findings

Orbital optimization yields the correct $^5 riangle$ ground state.

Jastrow factor size significantly affects DMC accuracy.

Optimized orbitals are essential for accurate energy predictions.

Abstract

Diffusion quantum Monte Carlo calculations with partial and full optimization of the guide function are carried out for the dissociation of the FeS molecule. For the first time, quantum Monte Carlo orbital optimization for transition metal compounds is performed. It is demonstrated that energy optimization of the orbitals of a complete active space wave function in the presence of a Jastrow correlation function is required to obtain agreement with the experimental dissociation energy. Furthermore, it is shown that orbital optimization leads to a $^5\Delta$ ground state, in agreement with experiments, but in disagreement with other high-level ab initio wave function calculations which all predict a $^5\Sigma^+$ ground state. The role of the Jastrow factor in DMC calculations with pseudo potentials is investigated. The results suggest that a large Jastrow factor may improve the DMC accuracy substantially at small additional cost.

Tables

Table 1. Table 1: FeS single determinant DMC energies in E h subscript 𝐸 ℎ E_{h} with HF, B3LYP, and VMC-optimized (opt orbs) orbitals and both Jastrow factors (Jas). All Jastrow parameters are optimized. Parentheses indicate one standard deviation.

state	Jas	HF	B3LYP	opt orbs
${}_{}{}^5\Delta$	J1	-134.0181(4)	-134.0426(4)	-134.0426(4)
${}_{}{}^5\Sigma _{}^{+}$	J1	-134.0163(4)	-134.0406(4)	-134.0386(4)
${}_{}{}^5\Delta$	J2	-134.0234(4)	-134.0486(4)	-134.0489(4)
${}_{}{}^5\Sigma _{}^{+}$	J2	-134.0223(4)	-134.0460(4)	-134.0456(3)

Table 2. Table 2: FeS single determinant VMC energies in E h subscript 𝐸 ℎ E_{h} with HF, B3LYP, and VMC-optimized (opt orbs) orbitals and both Jastrow factors (Jas). All Jastrow parameters are optimized.

state	Jas	HF	B3LYP	opt orbs
${}_{}{}^5\Delta$	J1	-133.8975(4)	-133.9164(4)	-133.9270(4)
${}_{}{}^5\Sigma _{}^{+}$	J1	-133.8961(4)	-133.9107(4)	-133.9225(4)
${}_{}{}^5\Delta$	J2	-133.9451(4)	-133.9708(4)	-133.9754(3)
${}_{}{}^5\Sigma _{}^{+}$	J2	-133.9440(4)	-133.9664(4)	-133.9681(3)

Table 3. Table 3: DMC energies in E h subscript 𝐸 ℎ E_{h} with CAS guide functions at different levels of parameter optimization as indicated. The Jastrow factor is indicated in parentheses.

	Jas(J1)	Jas(J1)+CI	all(J1)	all(J2)
${}_{}{}^5\Delta$	-134.0245(4)	-134.0348(4)	-134.0489(3)	-134.0579(4)
${}_{}{}^5\Sigma _{}^{+}$	-134.0326(4)	-134.0406(3)	-134.0485(4)	-134.0571(4)

Table 4. Table 4: Atomic DMC energies in E h subscript 𝐸 ℎ E_{h} with CAS guide function at different levels optimization and different Jastrow factors.

		Jas	Jas+MO
S	J1	-10.1298(1)	-10.1298(1)
S	J2	-10.1309(1)	-10.1314(1)
Fe	J1	-123.8073(4)	-123.8079(4)
Fe	J2	-123.8125(4)	-123.8126(4)

Table 5. Table 5: Spectroscopic constants of FeS for fully optimized CAS wave function with Jastrow factor J1. Equilibrium distance is given in Å, zero point energy in eV, and the harmonic frequency and anharmonicity in cm -1 .

state	$r_e$	ZPE	${\omega }_e$	${\omega }_e{x}_e$
${}_{}{}^5\Sigma _{}^{+}$	2.00(1)	0.0320(4)	518(7)	2.67(7)
${}_{}{}^5\Delta$	2.031(7)	0.0308(7)	499(11)	2.53(11)
${}_{}{}^5\Delta$	2.017a	0.0321(3)b	518(5)b	1.7(2)b

Table 6. Table 6: DMC dissociation energy D 0 subscript 𝐷 0 D_{0} in eV for FeS ( Δ 5 superscript Δ 5 {}^{5}\Delta ) with different guide functions and different optimized parameters (opt params). Spin-orbit, zero point energy, and core valence correlation corrections are included.

nodes	Jas	opt params	$D_0$
UB3LYP	J1	Jas	2.988(16)
1 Det	J1	Jas+MO	2.913(16)
1 Det	J2	Jas+MO	2.914(15)
CASSCF	J1	Jas	2.459(16)
CAS	J1	Jas+CI	2.720(16)
CAS	J1	Jas+MO+CI	3.085(15)
CAS	J2	Jas+MO+CI	3.159(15)
exp			3.31(15)a
a Drowart et al.Drowart et al. (1967)

Equations

$$\begin{array}{r}E=\frac{⟨{\Psi }_t(\text{R},\text{p})⟩H{\Psi }_t(\text{R},\text{p})}{⟨{\Psi }_t(\text{R},\text{p})⟩{\Psi }_t(\text{R},\text{p})}.\end{array}$$

$$\begin{array}{r}E=\int d\text{R}p(\text{R})E_L(\text{R})=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{i=1}^N{E}_L({\text{R}}_i)\end{array}$$

$$\begin{array}{r}-\frac{\partial }{\partial \tau }f(\text{R}.\tau )=-\frac{1}{2}{\nabla }^{2}f(\text{R},\tau )\\ +\nabla \cdot (\text{b}(\text{R})f(\text{R},\tau ))+[E_L(\text{R})-E_t]f(\text{R},\tau )\end{array}$$

$$f({\text{R}}^{\prime },{\tau }^{\prime })=\int d\text{R}K({\text{R}}^{\prime },{\tau }^{\prime };\text{R},\tau )f(\text{R},\tau )$$

$$K({\text{R}}^{\prime },{\tau }^{\prime };\text{R},\tau )\approx K_d({\text{R}}^{\prime },{\tau }^{\prime };\text{R},\tau )K_b({\text{R}}^{\prime },{\tau }^{\prime };\text{R},\tau )$$

$$K_d({\text{R}}^{\prime },{\tau }^{\prime };\text{R},\tau )=(2\pi \tau {)}^{-3N/2}{e}^{-\frac{({\text{R}}^{\prime }-\text{R}-\Delta \tau \text{b}(\text{R}){)}^2}{2\tau }}$$

$$K_b({\text{R}}^{\prime },{\tau }^{\prime };\text{R},\tau )=e^{-\tau (E_L({\text{R}}^{\prime })+E_L(\text{R})-2E_t)/2}$$

$$\begin{array}{r}|\Psi ⟩=e^U|\Phi ⟩\end{array}$$

$$\begin{array}{r}|\Phi ⟩=\sum _{i=1}^{N_{csf}}{c}_i|{\text{F}}_i⟩\end{array}$$

$$\begin{array}{r}|{\text{F}}_i⟩=\sum _k{d}_{i,k}|{\text{D}}_k^{\uparrow }⟩|{\text{D}}_k^{\downarrow }⟩.\end{array}$$

$$\begin{array}{r}|{\psi }_k⟩=\sum _{m=1}^{N_{bas}}{\alpha }_{k,m}|{\chi }_m⟩.\end{array}$$

Full text

I Introduction

Understanding the chemistry of transition metal compounds is the key to the comprehension of phenomena like photosynthesis and the design of highly efficient catalysts. The iron sulfide molecule FeS is the smallest member of iron and sulfur containing compounds which include important enzymes like nitrogenase. Despite its small size it is a very challenging system for ab initio wave function methods.

Variational (VMC) and diffusion (DMC) quantum Monte Carlo methods Wagner (2014); Austin et al. (2012); Lüchow (2011) are the most dominantly used stochastic methods to study real materials in physics and chemistry when there is a necessity for accuracy beyond mean field theory. DMC can provide results at least as reliable as the quantum chemistry gold standard CCSD(T) while being highly parallelizable on current hardware Grossman (2002); Nemec et al. (2010); Kannengießer et al. (2014); Klahm and Lüchow (2014); Dubecký et al. (2016).

VMC and DMC methods have successfully been applied to transition metal compounds. Most calculations to date have employed a single determinant guide function with pseudopotentials, for instance the early calculation of low-lying states and the dissociation energy of TiC by one of the present authors Sokolova and Lüchow (2000). Wagner and Mitas  Wagner and Mitas (2007, 2003) investigated a series of 3d transition metal oxides using VMC and DMC while dissociation energies were computed for transition metal carbonyl complexes by Diedrich et al. Diedrich et al. (2005). Further work includes the atomization and ionization energies for various vanadium oxide compounds Bande and Lüchow (2008). Doblhoff-Dier et al. Doblhoff-Dier et al. (2016) studied the performance of DMC for the dissociation energy of a selected set of transition metal compounds. Antisymmetrized geminal product guide functions with a Jastrow factor were employed by Casula et al. Casula et al. (2009) to describe the correct energy ordering for the iron dimer with DMC whereas Caffarel and coworkers studied the all-electron DMC spin density distribution for the $\text{CuCl}_{2}$ molecule Caffarel et al. (2014).

One of the present authors determined recently ionization and dissociation energies for 3d transition metal sulfides Petz and Lüchow (2011). Despite the overall very good performance of DMC with single determinant guide functions, a large deviation from experiment of 0.54 and 0.40 eV was observed for the dissociation energies of FeS and CrS, respectively. In this paper, FeS is therefore reexamined. We investigate in particular the effect of a multideterminant guide function and simultaneous Jastrow correlation function and orbital optimization on the quality of the DMC energies.

II Methods

Since there are comprehensive reviews on quantum Monte Carlo methods (QMC)Hammond et al. (1994); Lüchow (2011); Austin et al. (2012); Wagner (2014), VMC, DMC, and the form of the wave function employed in this work are discussed in this section only briefly.

II.1 Variational quantum Monte Carlo

In the variational quantum Monte Carlo method an upper bound to the ground state energy is calculated using a parametrized trial wave function $\Psi_{t}(\textbf{R},\textbf{p})$ with electron coordinates R and parameters p

[TABLE]

$E$ can be estimated efficiently using Monte Carlo quadrature with importance sampling,

[TABLE]

where $E_{L}(\textbf{R})=H\Psi_{t}(\textbf{R})/\Psi_{t}(\textbf{R})$ is called local energy, and $p(\textbf{R})=\Psi^{*}_{t}(\textbf{R})\Psi_{t}(\textbf{R})/\braket{\Psi_{t}(\textbf{R})}{\Psi_{t}(\textbf{R})}$ is the probability distribution. The parameters p are optimized with respect to the energy. The quality of the VMC results depends strongly on the quality and flexibility of the trial wave function.

II.2 Diffusion quantum Monte Carlo

In DMC, the imaginary time Schrödinger equation with importance sampling transformation is employed Reynolds et al. (1982)

[TABLE]

where $f(\textbf{R},\tau)=\Psi(\textbf{R},\tau)\Psi_{g}(\textbf{R})$, $\Psi_{g}(\textbf{R})$ is the guide function, $\textbf{b}(\textbf{R})=\nabla\Psi_{g}/\Psi_{g}$ is called drift or quantum force, and $E_{t}$ is the trial or reference energy. The equation is solved stochastically in its integral form

[TABLE]

with the propagator $K$. In the short time limit $\Delta\tau=\tau^{\prime}-\tau\to 0$ the propagator can be approximated by the product of a drift diffusion propagator $K_{d}$ and a weighting/branching term $K_{b}$

[TABLE]

where

[TABLE]

and

[TABLE]

In practice, a sample is in turn propagated with $K_{d}$ and weighted/branched with $K_{b}$ while enforcing the nodes of $\Psi_{g}$ onto $\Psi$ (fixed node approximation). The ground state energy is calculated with a mixed estimator Reynolds et al. (1982). After removal of the time step error by extrapolation $\Delta\tau\to 0$, the accuracy of the DMC energy depends only on the accuracy of the nodes of $\Psi_{g}$ which is a $3n-1$ dimensional hypersurface for $n$ electrons.

II.3 Wave functions and pseudo potentials

The same Slater Jastrow form is used for trial and guide function

[TABLE]

where $e^{U}$ is the symmetric Jastrow correlation factor describing the dynamic electron correlation, and $\ket{\Phi}$ is an antisymmetric Slater determinant or a linear combination thereof. $U(\boldsymbol{\beta})$ is expanded in many body terms depending on parameters $\boldsymbol{\beta}$. In this work, a 14-term Jastrow factor with electron-electron, electron-nucleus, and electron-electron-nucleus terms in the Boys/Handy form as suggested for QMC by Schmidt and Moskowitz  Schmidt and Moskowitz (1990) and a more accurate form with cusp-less electron-electron-nucleus terms which is identified as sm666 (69 terms) in reference  Lüchow et al. (2015) are employed. In the following, the two forms are denoted by J1 and J2, respectively.

The antisymmetric part $\ket{\Phi}$ is a Slater determinant or a linear combination of configuration state functions (CSF) $\ket{\textbf{F}_{i}}$, spin adapted linear combination of Slater determinants,

[TABLE]

[TABLE]

Spin up and spin down determinants $\ket{\textbf{D}^{\uparrow}_{k}}$ and $\ket{\textbf{D}^{\downarrow}_{k}}$ are built from orthogonal single particle molecular orbitals (MO) expanded in an atomic basis

[TABLE]

In this work, the triple zeta basis set by Burkatzki, Filippi, and Dolg has been employed  Burkatzki et al. (2008, 2007). The corresponding small core relativistic pseudo potentials have been used to improve the efficiency of QMC and to include scalar relativistic effects. As usual in QMC, the nonlocal pseudo potentials are localized on a spherical grid using the trial wave function Mitáš et al. (1991).

Jastrow parameters $\boldsymbol{\beta}$, CI coefficients c, and MO coefficients $\boldsymbol{\alpha}$ are variational parameters of the wave function. All parameters have been optimized with respect to the VMC energy. The MO coefficients are not optimized directly but rather using appropriate orbital rotations to keep the orthogonality of the orbitals Toulouse and Umrigar (2007). To the best of our knowledge, this is the first time orbital optimization is applied to transition metal compounds.

The perturbative and the linear method are used for wave function optimization, both of which make use of efficient statistical estimators for derivatives with respect to the parametersToulouse and Umrigar (2007). Special care has to be taken for the pseudo potential contribution to the local energy because the localized pseudo potential depends on all wave function parameters. Analytical derivatives with respect to all parameters have been implemented into the QMC code amolqc Manten and Lüchow (2001), developed in our group.

Analytical derivatives of the localized pseudo potential are particularly expensive for the orbital rotations. Therefore, the perturbative method was employed for the MO optimization. In this method, only the energy denominator $\Delta\epsilon_{i}$ requires derivatives of the pseudo potential (see Eq. (43) and (44) in Ref. Toulouse and Umrigar (2007)). Attempts to replace $\Delta\epsilon_{i}$ by guess values were not successful. Instead, we calculated $\Delta\epsilon_{i}$ with low accuracy using a small subset of the sample and fixed the energy denominators for subsequent iterations.

Due to the complete active space (CAS) nature of the Slater part of the wave function, many spin-up and spin-down determinants in the CSF expansion are identical and calculated only once. The 630 Slater determinants of the CAS wave function for the FeS ${}^{5}\Sigma^{+}$ are built from only nine different spin-up and 126 spin-down determinants. This leads to substantial speed up compared to a naive calculation of all Slater determinants Scemama et al. (2016). The multideterminant DMC calculations in this work require only about four times the CPU time for corresponding single determinant calculation.

III Results

FeS has been studied with many experimental and theoretical techniques. The ground state was assigned $X\ ^{5}\Delta$ by Zhang et el. Zhang et al. (1996) using mass spectrometry. This assignment was firmly established by Takano et al. using microwave spectroscopy Takano et al. (2004). The dissociation energy is not very accurately known. Drowart et al. Drowart et al. (1967) determined $D_{0}=3.31(15)$ eV following an earlier measurement of an upper bound of 3.34 eV by Marquart and Berkowitz Marquart and Berkowitz (1963). An experimental bond length of 2.017 Å is obtained for the $X\ ^{5}\Delta$ state from the rotational constant measured by Takano and coworkers Takano et al. (2004). DeVore et al. performed IR spectroscopy on FeS in different matrices and obtained a harmonic frequency of 540 cm*-1* while Zhai et al.  reported 520(30) cm*-1* employing photo electron spectroscopyDeVore and Franzen (1975); Zhai et al. (2003). Using dispersed fluorescence spectroscopy, Wang et al.  obtained more accurately 518 $\pm$ 5 cm*-1*Wang et al. (2011).

This molecule has also been studied extensively using density functional theory (DFT) and wave function methods. Surprisingly for such a small molecule, there is disagreement about the ground state. Many DFT calculations correctly predict a ${}^{5}\Delta$ ground state Hübner et al. (1998); Bridgeman and Rothery (2000); Schultz et al. (2005); Wu et al. (2007); Clima and Hendrickx (2007); Liang et al. (2009); Li et al. (2013), but some functionals including B3LYP and VSXC predict a ${}^{5}\Sigma^{+}$ ground state Hübner et al. (1998); Wu et al. (2007); Schultz et al. (2005) with several basis sets. High-level multireference calculations resulted all in a ${}^{5}\Sigma^{+}$ ground state. Sauer and coworkers used CASSCF/ACPF and found the ${}^{5}\Sigma^{+}$ state to lie 0.14 eV below the ${}^{5}\Delta$ state Hübner et al. (1998). Clima and Hendrickx employed CASPT2 and also found a ${}^{5}\Sigma^{+}$ ground state Clima and Hendrickx (2007). An earlier CASSCF/ICACPF calculation by Bauschlicher and Maitre considered only the ${}^{5}\Delta$ state Bauschlicher and Maitre (1995).

In this work, DMC with single determinant and CAS multideterminant guide functions is used to determine the ground state symmetry, the dissociation energy, and the spectroscopic constants. The electronic configuration for the ${}^{5}\Delta$ state is $\sigma^{2}\pi^{4}\sigma^{2}\delta^{3}\sigma^{1}\pi^{2}$ and $\sigma^{2}\pi^{4}\sigma^{2}\delta^{2}\sigma^{2}\pi^{2}$ for ${}^{5}\Sigma^{+}$. In the CASSCF calculations, 12 electrons were distributed in the active space formed by the 3d and 4s orbitals of the iron atom and the 3p orbitals of the sulfur atom. Dissociation energies are calculated by subtracting the atomic energies in their ground states ${}^{3}P$ and ${}^{5}D$ for sulfur and iron, respectively. State averaged CASSCF calculations were carried out for the atoms to preserve symmetry. Initial wave functions for QMC were created using Hartree-Fock (HF), B3LYP, and CASSCF with the GAMESS code Schmidt et al. (1993). The QMC code amolqcManten and Lüchow (2001) developed in our group is used for all QMC calculations.

III.1 Effect of orbital optimization in single determinant QMC calculations

VMC and DMC energies for single determinant wave functions were calculated for HF, B3LYP Kohn-Sham orbitals, and energy-optimized orbitals for both the ${}^{5}\Delta$ and the ${}^{5}\Sigma^{+}$ state. The bond lengths were fixed at the minima obtained from the potential curves calculated at the multideterminant DMC level (see below), 2.031 Å for ${}^{5}\Delta$ and 2.000 Å for ${}^{5}\Sigma^{+}$. The Jastrow parameters were optimized with respect to the energy in all calculations. If both Jastrow parameters and orbitals were optimized the parameters were optimized simultaneously. The DMC results, extrapolated to zero time step, are shown in Table 1.

The ${}^{5}\Delta$ energy is in all variants below the ${}^{5}\Sigma^{+}$ energy. The large difference between the DMC calculations with HF and B3LYP orbitals is very unusual. A comparison of HF and B3LYP orbitals shows a significant difference only in the quadruply occupied $\pi$ orbital which has substantially more bonding character in B3LYP compared to HF, see Figure 1. The same effect has been reported by Wagner and Mitas for some transition metal oxides Wagner and Mitas (2003, 2007). The B3LYP orbitals appear to be optimal because orbital optimization in addition to Jastrow parameter optimization does not lower the DMC energy, the ${}^{5}\Sigma^{+}$ energy is even slightly above the DMC energy with B3LYP orbitals.

When comparing the DMC results with HF or B3LYP guide functions and both Jastrow factors one observes a reduction of the DMC energy by approximately 6 m$E_{h}$ when using the larger J2. Since the nodes of the guide function are unchanged this difference is due to the different localization of the pseudo potential. These results highlight the importance of using an accurate Jastrow factor for the localization of the pseudo potential. Note that J1 is not a small Jastrow factor as it already contains three-body electron-electron-nucleus terms. When the orbitals are optimized in the presence of the Jastrow factor the more accurate J2 improves the nodes of the guide functions further in addition to reducing the pseudo potential localization error.

In Table 2, the corresponding VMC energies are shown. The last column containing the results for energy-optimized orbitals (and simultaneous optimization of the Jastrow factor) demonstrates the effectiveness of the optimization. Substantial systematic reduction of the VMC energy is obtained when using optimal orbitals instead of B3LYP orbitals or when using the larger Jastrow factor.

In Figure  2, the quadruply occupied $\pi$ orbital is shown after energy optimization with both Jastrow factors. Although the VMC energy is substantially lower compared to VMC with B3LYP orbitals, the $\pi$ orbitals display much less bonding character and for the more accurate J2 even less than for J1. We observe that the bonding character of the orbitals does not correlate with the stability of FeS. This fact may hint to the importance of many-electron effects in transition metal bonding.

These calculations show a preference for the ${}^{5}\Delta$ state but do not result in a dissociation energy close to the experimental value because both states show significant multireference character and require a multideterminant guide function for accurate determination of the dissociation energy.

III.2 Multideterminant QMC calculations

The most common way to include nondynamical correlation is a CAS-type wave function. The CAS space described above results in 630 determinants for ${}^{5}\Sigma^{+}$ and 500 for ${}^{5}\Delta$. For efficient evaluation, the wave function is expressed in terms of determinants while CSFs are employed for the parameter optimization. The Slater-Jastrow wave function is now capable of covering nondynamical correlation with the determinant part and simultaneously dynamical correlation with the Jastrow function. In Table 3, DMC results are shown for both states with some or all parameters optimized with respect to the energy. If not optimized, the CI and MO parameters are taken from the ab initio CASSCF wave function. The last column in Table 3 contains DMC energies with the fully optimized guide function employing the J2 Jastrow factor. All entries are a result of time step extrapolation. The DMC results are remarkable, in particular when comparing with the single determimant results in Table 1. For both states, all DMC energies are below the single determinant results with HF orbitals as expected, but surprisingly the single determinant results with B3LYP orbitals are substantially below the CAS-based DMC results even with optimized CI coefficients. Only after optimization of the orbitals the energy is substantially improved compared with the DMC/B3LYP result.

After the discussion of the single determinant energies this result may not be fully surprising as the strong lowering of the DMC energy with B3LYP orbitals indicates a significant influence of the dynamical electron correlation on the nodal structure. This influence is missing when CASSCF orbitals are employed. Only the orbital optimization in the presence of a Jastrow correlation function is capable of catching this influence. Optimization of the orbital parameters simultaneously with J1 Jastrow and CI parameters leads to an improvement of the energy by 14 and 8 m$E_{h}$ for the ${}^{5}\Delta$ and ${}^{5}\Sigma^{+}$ state, respectively, compared to optimization of J1 Jastrow and CI parameters only. Note that the DMC energies with the CASSCF wave function show a ${}^{5}\Sigma^{+}$ ground state. Only after orbital optimization a ${}^{5}\Delta$ ground state is identified although the energy difference is not statistically significant (as long as spin-orbit corrections are neglected, see below).

This discussion allows us to understand why DFT calculates mainly a ${}^{5}\Delta$ ground state whereas high-level wave function methods yield a ${}^{5}\Sigma^{+}$ state. DMC with CASSCF orbitals covers the dynamical correlation energy in a similar way (although more complete) as CASPT2 or ACPF. While the DMC accuracy is limited by the CASSCF nodes, CASPT2 and ACPF are limited by the excitation level out of the CASSCF reference and the basis set. When CAS orbitals optimized in the presence of a Jastrow factor (i.e. dynamic correlation) are used we observe a stabilization of ${}^{5}\Delta$ with respect to ${}^{5}\Sigma^{+}$. We expect the same result for CASPT2 if the optimization of the CAS orbitals for the full CASPT2 energy is done.

III.3 Calculation of the dissociation energy

The electronic dissociation energy is obtained by subtracting the atomic DMC energies shown in Table  4. There is only one CSF for each atom with the active space given above. The time step is extrapolated to zero. In case of the atoms the changes of the DMC energy for different levels of optimization is not significant and is of the order of the statistical error. The comparison of the two Jastrow factors demonstrates once more the improved pseudopotential localization with the larger Jastrow factor.

The experimental dissociation energy of $D_{0}=3.31(15)$ eV measured by Drowart et al. refers to 0 K Drowart et al. (1967). For comparison with the experimental value the zero point energies (ZPE) and the spin-orbit contribution are required. Scalar relativistic effects are taken into account by the pseudo potentials. Finally, the core-valence correlation energy is estimated.

Potential energy curves are calculated for both states at the DMC level with fully optimized multideterminant guide functions with the J1 Jastrow factor. For efficiency reasons, a fixed time step resulting in a 99% acceptance ratio was employed. From a fit of the curves to Morse functions we obtain the spectroscopic constants and the zero point energy shown in Table 5.

The equilibrium bond length of 2.031(7) Å is slightly longer than the experimental value 2.017 Å obtained by microwave spectroscopy Takano et al. (2004). The harmonic frequencies and anharmonicities agree with experimental values $\omega_{e}$ = 518(5) cm*-1* and $\omega_{e}x_{e}$ = 1.7(2) cm*-1* measured by Wang et al.  Wang et al. (2011). The zero point energy is slightly larger for the ${}^{5}\Sigma^{+}$ state.

While spin-orbit effects cannot be included directly in the fixed node DMC framework, Melton at al. Melton et al. (2016) were able to calculate spin-orbit interactions in DMC in fixed phase approximation. In this work, spin-orbit corrections are taken from the literature. For the atoms, -0.024 eV for sulfur and -0.050 eV for iron are derived as first-order spin-orbit corrections from experimental splittings Kramida et al. (2015). The ${}^{5}\Sigma^{+}$ state shows no first-order spin-orbit splitting while for the ${}^{5}\Delta$ state we obtain from Schultz et al. Schultz et al. (2005) a spin-orbit correction of -0.022 eV. Thus, the spin-orbit correction reduces the dissociation energy of the ${}^{5}\Sigma^{+}$ state by 0.074 eV and the dissociation energy of the ${}^{5}\Delta$ state by 0.052 eV.

Due to the pseudopotentials, the valence correlation energy is calculated with QMC. Rather than attempting expensive all-electron QMC calculations, the core-valence correlation contribution is estimated with multireference perturbation theory known as MC-QDPT2 or MR-MP2  Nakano (1993a, b) and implemented in GAMESS Schmidt et al. (1993). While the same active space as in QMC is employed the core-valence basis set TK+NOSeC-V-QZP with all diffuse functions is used Koga et al. (1999); Noro et al. (2000). The dissociation energy is calculated with and without correlating the core electrons to estimate the core-valence correlation contribution to the dissociation energy. This contribution is surprisingly large. For the ${}^{5}\Delta$ state, the dissociation energy increases by the core-valence correlation contribution by about 0.143 eV and by about 0.146 eV for the ${}^{5}\Sigma^{+}$ state.

With the DMC energies from fully optimized wave functions using the J2 Jastrow factor and all of the above corrections the energy gap between ${}^{5}\Delta$ and ${}^{5}\Sigma^{+}$ increases to 0.042(15) eV, confirming the ${}^{5}\Delta$ state as ground state in agreement with experiment.

Adding all corrections to the electronic dissociation energy for the ${}^{5}\Delta$ state we obtain +0.060 eV. Table  6 shows the DMC dissociation energy $D_{0}$ for the different guide wave functions with this correction added. Moderate agreement with the experimental value is obtained when unrestricted B3LYP orbitals are employed for the atoms and the molecule. The agreement deteriorates with both Jastrow factors slightly when using restricted determinants with optimized orbitals (where the symmetry of the orbitals is preserved during the optimization). Substantially stronger deviation from the experimental value is obtained not only with an unoptimized CASSCF wave function where only the Jastrow part is optimized, but also when optimizing the CI coefficients of the CAS wave function. Finally, improvement over the unrestricted single determinant energy is calculated with full optimization of the CAS wave function, i.e. optimization of orbitals, CI coefficients, and Jastrow parameters simultaneously. With the large J2 Jastrow function, excellent agreement with the experimental value is obtained. The significant effect of the Jastrow factor on the DMC dissociation energy was surprising to us. As discussed above it is due to the pseudo potential localization and the effect of the Jastrow factor on the orbital optimization. The additional CPU cost for the substantially larger J2 Jastrow factor is, on average, only 10%.

III.4 Conclusion

Single and CAS-based multideterminant DMC calculations have been carried out with partial and full optimization of Jastrow, CI, and orbital parameters. The ${}^{5}\Delta$ state was identified as ground state, in agreement with experiment but in disagreement with other wave function methods. ${}^{5}\Delta$ is obtained as ground state only after optimization of the CAS orbitals simultaneously with the Jastrow and CI parameters. The CASSCF reference seems to be biased toward the ${}^{5}\Sigma+$ state which can explain why CASPT2 yields a ${}^{5}\Sigma^{+}$ ground state. The optimization of the CAS orbitals in the presence of the Jastrow correlation function allows thus to construct compact accurate wave functions yielding excellent energies even for challenging molecules such as FeS. In addition to the dissociation energy, a potential curve was calculated, yielding harmonic frequencies and anharmonicities in good agreement with experiment. Furthermore, the influence of the Jastrow factor on the DMC results has been investigated. While most DMC calculations on transition metal compounds employed small Jastrow factors, the above results indicate that a large Jastrow factor allows for higher DMC accuracy through improved pseudo potential localization and improved optimized orbitals. The additional cost for a larger Jastrow factor is small in pseudo potential calculations.

Acknowledgements.

The authors gratefully acknowledge computing resources granted by Jülich-Aachen research alliance (JARA) and RWTH Aachen University under projects jara0119 and rwth0128.