Time-dependent linear-response variational Monte Carlo
We present the extension of variational Monte Carlo (VMC) to the calculation of electronic excitation energies and oscillator strengths using
time-dependent linear-response theory. By exploiting the analogy existing between the linear method for wave-function optimisation and the generalised eigenvalue equation of linear-response theory, we formulate the equations of linear-response VMC (LR-VMC). This LR-VMC approach involves the first- and second-order derivatives of the wave function with respect to the parameters. We perform first tests of the LR-VMC method within the Tamm-Dancoff approximation using single-determinant Jastrow-Slater wave functions with different Slater basis sets on some singlet and triplet excitations of the beryllium atom. Comparison with reference experimental data and with configuration-interaction-singles (CIS) results shows that LR-VMC generally outperforms CIS for excitation energies and is thus a promising approach for calculating electronic excited-state properties of atoms and molecules.
Keywords: excitation energies - linear method - Tamm-Dancoff approximation - oscillator strengths - beryllium
Quantum Monte Carlo (QMC) methods bk:hammond ; fou+01rmp ; tou16 are a powerful and reliable alternative to wave-function methods and density-functional theory (DFT) for quantum chemistry calculations, thanks to their favorable scaling with system size and to their suitability for high-performance computing infrastructures, such as petascale architectures. Variational Monte Carlo (VMC) bres98 combines Monte Carlo integration for computing the expectation value of the electronic Hamiltonian and the variational principle for the ground state. VMC scales as (where is the number of electrons), similar to DFT scaling. The main drawback of any QMC approach is the very large prefactor in the scaling, preventing the systematic use of QMC in quantum chemistry calculations of medium- and large-size systems. This drawback is alleviated by performing massive parallel calculations on supercomputers Coccia:2012kz ; Coccia14 .
A fundamental role is played by the trial wave function, often written as a product of a determinantal part and a bosonic Jastrow factor dru+04prb which depends on interparticle distances (with electron-nucleus, electron-electron, higher many-body terms,…). For example, one can use for the determinantal part a linear combination of configuration state functions (CSF, i.e. spatial- and spin-symmetry adapted linear combinations of Slater determinants of one-electron molecular orbitals) pet12 , or the antisymmetrised geminal power (AGP) ansatz (a single determinant of geminal pairing functions cas+03jcp ; cas+04jcp ; zen14 ). Furthermore, the optimisation of the wave function is crucial for an accurate description of both static and dynamic electron correlation. The linear method tou07 ; umr+07prl ; tou08 allows one to efficiently perform such an optimisation for all the parameters of the wave function, using only the first-order derivatives of the wave function with respect to the parameters.
The calculation of excited-state properties of molecules (from prototypical models to complex organic dyes and biochromophores) still represents an open challenge for theoreticians. The two commonly used approaches are time-dependent density-functional theory, which is not computationally demanding but sometimes lacks accuracy, and wave-function methods, which are more accurate but very computationally demanding. QMC methods were originally formulated for ground state problems and their extension to excited states is not straightforward. Relatively few applications of QMC for electronic excitations are present in literature, see e.g. the singlet and triplet energies for the benchmark CH diradical zim09 , the low-lying singlet excited states of biochromophores fili10 ; filippi2011bathochromic ; val12 ; val13 , the transition in acrolein tou12 ; flo14 , and the recent extension of the AGP ansatz for calculating excited-state energies zen+15jctc ; dup15 ; zha16a ; zha16b ; neu16 .
The basic idea of the present work stems from the formal analogy existing between the linear method for wave function optimisation and time-dependent linear-response theory bk:mc . Indeed, the generalised eigenvalue equations of linear-response theory in the Tamm-Dancoff approximation (TDA) and of the linear method at the ground-state minimum coincide. Starting from this observation, we derive and implement the linear-response equations in VMC (LR-VMC). This represents an extension of the well-established time-dependent linear-response Hartree-Fock or multiconfiguration self-consistent field methods, taking into account both static and dynamic electron correlations.
The paper is organised as follows. In Section II, VMC and linear-response theory are briefly reviewed, and the LR-VMC method is presented and discussed in detail. Results of LR-VMC calculations in the TDA of some singlet and triplet excitations of the beryllium atom are reported and discussed in Section III. Conclusions and perspectives for future work are given in Section IV.
We first briefly review the form of the wave function that we use and the linear optimisation method. We then derive the time-dependent linear-response equations and show how to implement them in VMC.
ii.1 Wave-function parametrisation
where is a Jastrow factor operator depending on a set of parameters , is the orbital rotation operator depending on a set of orbital rotation parameters , and are CSFs with associated coefficients . The CSFs are linear combinations of Slater determinants of orbitals , which are expanded in a basis of Slater functions
The Slater functions are centered on the nuclei and their spatial representation is
each function being characterized by a set of quantum numbers and an exponent , are real spherical harmonics, and a normalisation factor. The full set of parameters to consider is thus where stands for the set of exponents.
ii.2 Linear optimisation method
The linear optimisation methodtou07 ; umr+07prl ; tou08 allows one to find the optimal parameters p using an iterative procedure. At each iteration, we consider the intermediate-normalised wave function
where is the wave function for the parameters at the current iteration (taken as normalised to unity, i.e. ), and we expand it to linear order in the parameter variations ,
where are the first-order derivatives of the wave function
where are the first-order derivatives of the original wave function . Using the intermediate-normalised wave function has the advantage that the derivatives in Eq. (6) are orthogonal to , i.e. . We then determine the parameter variations by minimising the corresponding energy
we update the parameters as , and iterate until convergence.
The minimisation in Eq. (7) leads to the following generalized eigenvalue equation to be solved at each iteration
where is the current energy, and are the left and right energy gradients (identical except on a finite Monte Carlo sample), and is the Hamiltonian matrix in the basis of the first-order wave function derivatives, and is the overlap matrix in this basis. Note that in Eq. (8), 0 and stand for the zero column vector and the zero row vector, respectively.
ii.3 Linear-response theory
Starting from the previously optimised wave function, we introduce now a time-dependent perturbation (e.g, interaction with an electric field) in the Hamiltonian
where is a coupling constant. The approximate ground-state wave function evolves in time through its parameters , which become now generally complex. As before, it is convenient to introduce the intermediate-normalised wave function
where is the wave function for the initial parameters , again taken as normalised to unity (i.e., ). At each time, the time-dependent parameters can be found from the Dirac-Frenkel variational principle (see, e.g., Ref.27)
To apply Eq. (11) in linear order in , we start by expanding the wave function around to second order in the parameter variations
where are the first-order derivatives of already introduced in Eq. (6), and are the second-order derivatives of the wave function
where are the second-order derivatives of the original wave function . Again, the advantage of using the intermediate-normalised wave function is that the second-order derivatives are orthogonal to , i.e. . Plugging Eq. (12) into Eq. (11) and keeping only first-order terms in , in the limit of a vanishing perturbation (), we find
with the matrices where is the ground-state energy, , and . If we look for free-oscillation solutions of the form
where corresponds to an excitation (or de-excitation) energy, we arrive at the linear-response equation in the form of a non-Hermitian generalized eigenvalue equationbk:mc
The Tamm-Dancoff approximation (TDA) corresponds to neglecting the contributions from B, leading to
At the ground-state minimum, i.e. when the energy gradient is zero, the generalised eigenvalue equation of the linear method in Eq. (8) is equivalent to the TDA equation (17) which directly gives excitation energies .
Finally, the oscillator strength for the transition from the ground state to the excited state (with excitation energy ) can be easily extracted from the response vector
where is the vector containing the transition dipole moments for the component (, , or ) between the ground-state wave function and the wave-function derivative
and is the electronic dipole operator.
ii.4 Realisation in VMC
We now give the expressions for performing linear-response calculations in VMC, referred to as LR-VMC, i.e. for calculating the expressions in Section II.3 in a VMC run. For convenience, we also recall the expressions necessary for the linear optimisation method.
The current ground-state energy is calculated as
where is the local energy and stands for an average on a finite Monte Carlo sample of points distributed according to , with designating the electron coordinates. The left and right energy gradients are evaluated as
where is the first-order derivative of the local energy
Note that, in the limit of an infinite sample, due to the hermiticity of the Hamiltonian, and therefore the left and right gradients become identical.
The elements of the overlap matrix S are calculated as
and the elements of the matrix H are evaluated as
The elements of the matrix A are then given by
and the elements of the matrix B are
Finally, the expression of the transition dipole moment needed for calculating oscillator strengths is
where is the -component of electronic dipole moment.
In the linear optimisation method, using the non-symmetric estimator of the matrix H in Eq. (25) instead of a symmetrised one has the advantage of leading to the strong zero-variance principle of Nightingale and Melik-Alaverdian zero : in the limit where the current wave function and its first-order derivatives form a complete basis for expanding an exact eigenstate of the Hamiltonian, the parameter variations and the associated energy are found from Eq. (8) with zero variance provided that the Monte Carlo sample size is larger than the number of parameters (see discussion in Ref. 14). Unfortunately, this strong zero-variance principle does not apply when solving the linear-response equation (16). However, in the limit where is an exact eigenstate of the Hamiltonian, the left energy gradient in Eq. (LABEL:gLi) vanishes with zero variance, and thus the TDA linear equation (17) becomes equivalent to Eq. (8) for calculating excited-state energies even on a finite Monte Carlo sample. Therefore, in this case, the strong zero-variance principle applies to the calculation of the response vectors and excitation energies .
ii.5 Computational details
The calculations shown here were performed using the QMC program CHAMP CHAMP , starting from Hartree-Fock calculations done with GAMESS SchBalBoaElbGorJenKosMatNguSuWinDupMon-JCC-93 . Two Slater basis sets of different sizes were used, namely the VB1 and VB2 basis set from Ref. 31. The VB1 basis set has five and one Slater functions (), whereas the VB2 basis set has six , two , and one Slater functions (). We use a flexible Jastrow factor consisting of the exponential of the sum of electron-nucleus, electron-electron and electron-electron-nucleus terms, written as systematic polynomial and Padé expansions Umr-UNP-XX ; FilUmr-JCP-96 ; GucSanUmrJai-PRB-05 , with 4 electron-nucleus parameters, 5 electron-electron parameters and 15 electron-electron-nucleus parameters. For each VMC calculation, 10 blocks were employed with 10 steps each. One block was used for equilibration of the VMC distribution.
The beryllium atom was used as a first test of the LR-VMC approach, since accurate experimental reference values for the excitation energies are available from Ref. 35. An accurate description of the Be ground state requires a multiconfigurational wave function for accounting for the near-degeneracy between the 2 and 2 orbitals. However, for these preliminary tests, we present only results of calculations using a Jastrow-Slater single-determinant wave function for the ground state using TDA linear-response theory. This choice is motivated by the fact that a direct comparison between the LR-VMC/TDA method and configuration-interaction-singles (CIS) calculations represents a simple but essential first step for validating our approach. We expect LR-VMC/TDA to outperform CIS because the Jastrow factor in LR-VMC should account for a substantial part of electronic correlation, and we find this to be the case for most of the excitations studied. The results are presented both as errors with respect to the experimental values in Figure 1 and as detailed excitation energies in the subsequent tables.
In Table 1 results for the singlet 23 () state are reported. The effect of the Slater basis set adopted is dramatic at the CIS level, as a reasonable agreement with the reference experimental value of 0.249 Hartree is found only when the VB2 basis set is used. LR-VMC/TDA values are labelled as follow: (j) designates the response of the Jastrow parameters only, while (j+o) is the response of both the Jastrow and orbital parameters. The response of the Jastrow factor substantially improves upon the CIS VB1 estimate, going from 0.378 to 0.2888(1) Hartree. The excitation energy improves further when the response of the orbital parameters are included in the LR-VMC/TDA calculation, yielding an error of around 0.02 Hartree with respect to the experimental value. Increasing the size of the Slater basis set, i.e. moving from VB1 to VB2, we obtain a fair agreement with the experimental data when both the Jastrow and the orbital parameters are included in the response (0.2378(2) Hartree).
The singlet 24 () excitation is higher in energy, and CIS fails to recover the experimental result of 0.297 Hartree, for both basis sets, as shown in Table 2. As already mentioned for the 23 excitation, the response of the Jastrow factor plays an important role for the VB1 basis set, reducing the error in the excitation energy by around 2 Hartree. Including the orbital parameters in the response lowers the excitation energy further to 0.5578(2), but this is still a large overestimate of the experimental value. With the VB2 basis, the LR-VMC/TDA(j+o) calculation outperforms CIS, but a substantial error (0.02 Hartree) still remains for this high-lying excitation. The failure of VB1 and, to a lesser extent, of VB2 is likely due to the poor description of the 4 orbital.
The extension of our proposed approach to excitations is straightforward, with a relaxation of the spatial symmetry constraints in the orbital rotation parameters. Note that the Jastrow factor employed in this work only depends on interparticle distances, i.e. it has spherical symmetry, and therefore excited states with symmetry cannot be represented with the wave-function derivatives with respect to the Jastrow parameters. For this reason, only results concerning the response of the orbitals (o) are reported for the excitations. In Table 3, results for the singlet 22 () state are given, which is the lowest energy excitation in the beryllium atom. The CIS calculations with the VB1 and VB2 basis sets show a fair agreement with the reference value of 0.194 Hartree, the CIS calculation using the VB2 basis set being only 5 mHartree below it. The LR-VMC/TDA(o) estimate is also close to the experimental reference when the VB2 basis set is employed (0.1873(2) Hartree), while for the VB1 basis set LR-VMC/TDA(o) greatly overestimates the excitation energy.
Similarly, our implementation of linear response allows us to easily compute triplet excitations by considering triplet orbital rotation parameters. The CIS calculation underestimates the correct excitation energy by more than 30 mHartree, while the LR-VMC/TDA(o) excitation energies are very close to the reference values of 0.100 Hartree. The basis set effects are small is this case.
Finally, we computed the oscillator strength (Table 5) corresponding to the singlet () excitation, which is non zero according to selection rules. The LR-VMC(o) oscillator strengths seem more sensitive to the basis set compared to the CIS oscillator strengths. Moreover, the inclusion of the Jastrow factor does not improve the oscillator strength.
Iv Conclusions and Perspectives
In this work we have presented a formulation of time-dependent linear-response theory in the VMC framework using a Jastrow-Slater wave function. Compared to state-specific or state-average excited-state QMC methods, the advantage of this LR-VMC approach is that, after optimizing only one ground-state wave function, one can easily calculate several excitation energies of different spatial or spin symmetry. Compared to similar linear-response quantum chemistry methods, the presence of the Jastrow factor in LR-VMC allows one to explicitly treat a part of dynamical correlation. A disadvantage of the method is that the excitation energies are much more sensitive than the ground-state energy to the quality of the optimized ground-state wave function. This is true in other linear-response quantum-chemistry methods as well, but is a bigger drawback in a method that employs stochastic optimization.
Using a Jastrow-Slater single-determinant wave function and the TDA, the LR-VMC method was shown to be more accurate that CIS for most of the excitation energies of the beryllium atom that were studied. The LR-VMC approach thus seems a promising method for calculating electronic excitation energies. In the near future, a systematic study on a set of molecules will be an essential step to further validate the proposed methodology, together with calculations using the full response equation beyond the TDA. Also, we will explore using multideterminant wave functions, larger basis sets, and including the wave-function derivatives with respect to the exponents of the Slater functions.
EC thanks University of L’Aquila for financial support and the Laboratoire de Chimie Théorique for computational resources. MO and CJU were supported in part by NSF grant ACI-1534965.
- (1) B. L. Hammond, W. A. Lester, Jr., and P. J. Reynolds, Monte Carlo Methods in Ab-Initio Quantum Chemistry. World Scientific, 1994.
- (2) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, “Quantum monte carlo simulations of solids,” Rev. Mod. Phys., vol. 73, pp. 33–83, 2001.
- (3) J. Toulouse, R. Assaraf, and C. J. Umrigar, “Introduction to the Variational and Diffusion Monte Carlo methods,” Adv. Quantum Chem., vol. 73, pp. 285–314, 2016.
- (4) D. Bressanini and P. J. Reynolds, “Between classical and quantum monte carlo methods: ”variational” qmc,” Advances in Chemical Physics, Monte Carlo Methods in Chemical Physics, vol. 105, pp. 5345–5350, 1998.
- (5) E. Coccia and L. Guidoni, “Quantum monte carlo study of the retinal minimal model CHNH,” J. Comput. Chem., vol. 33, pp. 2332–2339, 2012.
- (6) E. Coccia, D. Varsano, and L. Guidoni, “Ab Initio Geometry and Bright Excitation of Carotenoids: Quantum Monte Carlo and Many Body Green’s Function Theory Calculations on Peridinin,” J. Chem. Theory Comput., vol. 10, pp. 501–506, Jan. 2014.
- (7) N. D. Drummond, M. D. Towler, and R. J. Needs, “Jastrow correlation factor for atoms, molecules, and solids,” Phys. Rev. B, vol. 70, p. 235119, 2004.
- (8) F. R. Petruzielo, J. Toulouse, and C. J. Umrigar, “Approaching chemical accuracy with quantum Monte Carlo,” J. Chem. Phys., vol. 136, p. 124116, 2012.
- (9) M. Casula and S. Sorella, “Geminal wave function with jastrow correlation: A first application to atoms,” J. Chem. Phys., vol. 119, p. 6500, 2003.
- (10) M. Casula, C. Attaccalite, and S. Sorella, “Correlated geminal wave function for molecules: An efficient resonating valence bond approach,” J. Chem. Phys., vol. 121, p. 7110, 2004.
- (11) A. Zen, E. Coccia, Y. Luo, S. Sorella, and L. Guidoni, “Static and Dynamical Correlation in Diradical Molecules by Quantum Monte Carlo Using the Jastrow Antisymmetrized Geminal Power Ansatz,” J. Chem. Theory Comput., vol. 10, pp. 1048–1061, 2014.
- (12) J. Toulouse and C. J. Umrigar, “Optimization of quantum Monte Carlo wave functions by energy minimization,” J. Chem. Phys., vol. 126, p. 084102, 2007.
- (13) C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, “Alleviation of the fermion-sign problem by optimization of many-body wave functions,” Phys. Rev. Lett., vol. 98, p. 110201, 2007.
- (14) J. Toulouse and C. J. Umrigar, “Full optimization of Jastrow-Slater wave functions with application to the first-row atoms and homonuclear diatomic molecules,” J. Chem. Phys., vol. 128, p. 174101, 2008.
- (15) P. M. Zimmerman, J. Toulouse, Z. Zhang, C. B. Musgrave, and C. J. Umrigar, “Excited states of methylene from quantum Monte Carlo,” J. Chem. Phys., vol. 131, p. 124103, 2009.
- (16) O. Valsson and C. Filippi, “Photoisomerization of model retinal chromophores: Insight from quantum Monte Carlo and multiconfigurational perturbation theory,” J. Chem Theory Comput., vol. 6, p. 1275, 2010.
- (17) C. Filippi, F. Buda, L. Guidoni, and A. Sinicropi, “Bathochromic shift in green fluorescent protein: a puzzle for QM/MM approaches,” J. Chem. Theory Comput., vol. 8, no. 1, pp. 112–124, 2011.
- (18) O. Valsson, C. Angeli, and C. Filippi, “Excitation energies of retinal chromophores: Critical role of the structural model,” Phys. Chem. Chem. Phys., vol. 14, p. 11015, 2012.
- (19) O. Valsson, P. Campomanes, I. Tavernelli, U. Rothlisberger, and C. Filippi, “Rhodopsin absorption from first principles: Bypassing common pitfalls,” J. Chem. Theory Comput., vol. 9, p. 2441, 2013.
- (20) J. Toulouse, M. Caffarel, P. Reinhardt, P. E. Hoggan, and C. J. Umrigar, “Quantum Monte Carlo calculations of electronic excitation energies: The case of the singlet (CO) transition in acrolein,” in ”Advances in the Theory of Quantum Systems in Chemistry and Physics”, Progress in Theoretical Chemistry and Physics, vol. 22, pp. 345–353, 2012.
- (21) F. M. Floris, C. Filippi, and C. Amovilli, “Electronic excitations in a dielectric continuum solvent with quantum Monte Carlo: Acrolein in water,” J. Chem. Phys., vol. 140, p. 034109, 2014.
- (22) A. Zen, E. Coccia, S. Gozem, M. Olivucci, and L. Guidoni, “Quantum monte carlo treatment of the charge transfer and diradical electronic character in a retinal chromophore minimal model,” J. Chem. Theory Comput., vol. 11, pp. 992–1005, 2015.
- (23) N. Dupuy, S. Bouaouli, F. Mauri, S. Sorella, and M. Casula, “Vertical and adiabatic excitations in anthracene from quantum Monte Carlo: Constrained energy minimisation for structural and electronic excited-state properties in the JAGP ansatz,” J. Chem. Phys., vol. 142, p. 214109, 2015.
- (24) L. Zhao and E. Neuscamman, “Equation of Motion Theory for Excited States in Variational Monte Carlo and the Jastrow Antisymmetric Geminal Power in Hilbert Space,” J. Chem. Theory. Comput., vol. 12, pp. 3719–3726, 2016.
- (25) L. Zhao and E. Neuscamman, “An efficient Variational Principle for the Direct Optimisation of Excited States,” J. Chem. Theory. Comput., vol. 12, pp. 3436–3440, 2016.
- (26) E. Neuscamman, “Variation after response in quantum Monte Carlo,” J. Chem. Phys., vol. 145, p. 081103, 2016.
- (27) R. McWeeny, Methods of Molecular Quantum Mechanics. Academic Press, 1992.
- (28) M. P. Nightingale and V. Melik-Alaverdian, “Optimisation of Ground- and Excited-State Wave functions and van der Waals Clusters,” Phys. Rev. Lett., vol. 87, p. 043401, 2001.
- (29) “CHAMP, a quantum Monte Carlo program written by C. J. Umrigar, C. Filippi and J. Toulouse, see http://www.physics,cornell.edu/cyrus/champ.html.”
- (30) M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, “General atomic and molecular electronic structure system,” J. Comput. Chem., vol. 14, p. 1347, 1993.
- (31) I. Ema, J. M. G. de la Vega, G. Ramirez, R. Lopez, J. F. Rico, H. Meissner, and J. Paldus, “Polarized basis sets of Slater-type orbitals: H to Ne atoms,” J. Comput. Chem., vol. 24, pp. 859–868, 2003.
- (32) C. J. Umrigar. unpublished.
- (33) C. Filippi and C. J. Umrigar, “Multiconfiguration wave functions for quantum Monte Carlo calculations of first-row diatomic molecules,” J. Chem. Phys., vol. 105, p. 213, 1996.
- (34) A. D. Güçlü, G. S. Jeon, C. J. Umrigar, and J. K. Jain, “Quantum Monte Carlo study of composite fermions in quantum dots: The effect of Landau-level mixing,” Phys. Rev. B, vol. 72, p. 205327, 2005.
- (35) A. Kramida and W. C. Martin, “A Compilation of Energy Levels and Wavelengths for the Spectrum of Neutral Beryllium (Be I),” J. Phys. Chem. Ref. Data, vol. 26, p. 1185, 1997.
- (36) R. Schnabel and M. Kock, “-value measurement of the Be I resonance line using a nonlinear time-resolved laser-induced-fluorescence technique,” Phys. Rev. A, vol. 61, p. 062506, 2000.