Approximate analytical gradients and nonadiabatic couplings for the state-average density matrix renormalization group self-consistent field method

TL;DR

This paper introduces an approximate method for calculating analytical gradients and nonadiabatic couplings within the state-average DMRG-SCF framework, enabling accurate and efficient simulations of complex electronic states.

Contribution

It develops a new approach to compute gradients and couplings in DMRG-SCF using a single-site MPS Lagrange multiplier, avoiding sweep procedures and closely following the SA-CASSCF formalism.

Findings

Successfully applied to optimize a conical intersection in 1,2-dioxetanone.

Reproduces SA-CASSCF results with arbitrary accuracy.

Provides a computationally efficient scheme for nonadiabatic dynamics simulations.

Abstract

We present an approximate scheme for analytical gradients and nonadiabatic couplings for calculating state-average density matrix renormalization group self-consistent-field wavefunction. Our formalism follows closely the state-average complete active space self-consistent-field (SA-CASSCF) \emph{ansatz}, which employs a Lagrangian, and the corresponding Lagrange multipliers are obtained from a solution of the coupled-perturbed CASSCF (CP-CASSCF) equations. We introduce a definition of the matrix product state (MPS) Lagrange multipliers based on a single-site tensor in a mixed-canonical form of the MPS, such that a sweep procedure is avoided in the solution of the CP-CASSCF equations. We apply our implementation to the optimization of a conical intersection in 1,2-dioxetanone, where we are able to fully reproduce the SA-CASSCF result up to arbitrary accuracy.

Tables

Table 1. Table 1: Average electronic energy, gradient norm, and the maximum gradient element throughout the ( σ ​ σ ∗ 𝜎 superscript 𝜎 \sigma\sigma^{*} ),( n ​ σ ∗ 𝑛 superscript 𝜎 n\sigma^{*} ) conical intersection optimization in 1,2-dioxetanone.

Step	Electronic energy/		Avg. gradient norm/		Max. gradient element/
	Hartree		Hartree/bohr		Hartree/bohr
	CASSCF	DMRG-SCF	CASSCF	DMRG-SCF	CASSCF	DMRG-SCF
1	$-301.76392482$	$-301.76392481$	$0.098030$	$0.098090$	$-0.070811$	$-0.070874$
2	$-301.77670336$	$-301.77670325$	$0.089606$	$0.089626$	$-0.050384$	$-0.050435$
3	$-301.79392062$	$-301.79391779$	$0.159354$	$0.159318$	$0.127371$	$0.127362$
4	$-301.78862423$	$-301.78862694$	$0.073202$	$0.073216$	$-0.034099$	$-0.034124$
5	$-301.81089953$	$-301.81089903$	$0.066443$	$0.066451$	$0.037325$	$0.037307$
6	$-301.81982062$	$-301.81982236$	$0.027812$	$0.027800$	$-0.017275$	$-0.017280$
7	$-301.82254958$	$-301.82254963$	$0.017071$	$0.017075$	$0.012606$	$0.012609$
8	$-301.82398804$	$-301.82398842$	$0.011692$	$0.011693$	$-0.008136$	$-0.008114$
9	$-301.82423076$	$-301.82423099$	$0.004222$	$0.004204$	$-0.003256$	$-0.003232$
10	$-301.82426240$	$-301.82426246$	$0.000897$	$0.000897$	$-0.000552$	$-0.000549$
11	$-301.82426314$	$-301.82426325$	$0.000404$	$0.000387$	$-0.000269$	$-0.000254$

Table 2. Table 2: Mean absolute error of gradient elements (in Hartree/bohr) calculated with SA-DMRG-SCF ( m = 500 𝑚 500 m=500 ) and SA-CASSCF, both analytically and numerically. The first two entries are MAE as defined in Eq. ( 59 ) for analytical and numerical gradients, respectively, while the last two entries are the deviation between analytical and numerical SA-DMRG-SCF and SA-CASSCF gradients, respectively.

	DMRG-SCF$-$CASSCF
Analytic	$1.9\times {10}^{-3}$
Numerical	$2.2\times {10}^{-4}$
	Analytic$-$Numerical
DMRG-SCF	$2.0\times {10}^{-3}$
CASSCF	$5.6\times {10}^{-5}$

Equations

$$|\Phi ⟩=\sum _I{c}_I|{\varphi }_I⟩$$

$$\hat{H}=\sum _{\underset{\tau }{t,u}}⟨t⟩hu{a}_{t\tau }^{†}{a}_{u\tau }+\frac{1}{2}\sum _{\underset{\tau ,{\tau }^{\prime }}{t,u,v,w}}(tu|vw)a_{t\tau }^{†}{a}_{v{\tau }^{\prime }}^{†}{a}_{w{\tau }^{\prime }}{a}_{u\tau }+E_{\text{core}}$$

$$|\Phi ⟩=\Phi (\mathbf{\kappa },\mathbf{c})$$

$$E^{\Phi }=⟨\Phi ⟩\hat{H}\Phi =\sum _{pq}⟨p⟩hq{\gamma }_{pq}+\frac{1}{2}\sum _{pqrs}(pq|rs){\Gamma }_{pqrs},$$

$$\frac{d{E}^{\Phi }}{dx}$$

$$=\sum _{pq}\frac{\partial ⟨p⟩hq}{\partial x}{\gamma }_{pq}+\frac{1}{2}\sum _{pqrs}\frac{\partial (pq|rs)}{\partial x}{\Gamma }_{pqrs}.$$

$$E^{\text{SA}}=\sum _{\Psi }{\omega }_{\Psi }{E}^{\Psi }=\sum _{\Psi }{\omega }_{\Psi }⟨\Psi ⟩\hat{H}\Psi$$

$$E^{\text{SA}}=E^{\text{SA}}(\mathbf{\kappa },{\mathbf{c}}^1,{\mathbf{c}}^2,\dots ,{\mathbf{c}}^n)$$

$$\frac{\partial E^{\text{SA}}}{\partial\kappa_{pq}}=0;\quad\frac{\partial E^{\text{SA}}}{\partial c_{\Theta I}}=0;\quad\frac{\partial E^{\Theta}}{\partial c_{\Phi I}}\Bigg{|}_{\Theta\neq\Phi}\neq 0;\quad\frac{\partial E^{\Theta}}{\partial\kappa_{pq}}\neq 0;$$

$$L^{\Theta }=E^{\Theta }+\left(\sum _{pq}{\bar{\kappa }}_{pq}^{\Theta }\frac{\partial E^{\text{SA}}}{\partial {\kappa }_{pq}}-0\right)+\left(\sum _{\Psi I}{\bar{c}}_{\Psi I}^{\Theta }\frac{\partial E^{\text{SA}}}{\partial c_{\Psi I}}-0\right).$$

$$\frac{\partial L^{\Theta }}{\partial {\bar{\kappa }}_{pq}^{\Theta }}=\frac{\partial L^{\Theta }}{\partial {\bar{c}}_{\Psi I}^{\Theta }}=\frac{\partial L^{\Theta }}{\partial {\kappa }_{pq}}=\frac{\partial L^{\Theta }}{\partial c_{\Psi I}}=0.$$

$$\frac{\partial L^{\Theta }}{\partial {\kappa }_{rs}}=$$

$$\frac{\partial L^{\Theta }}{\partial c_{\Phi J}}=$$

$$\left(\begin{array}[]{cc}\mathbf{H^{\text{OO}}}&\mathbf{H^{\text{CO}}}\\ \mathbf{H^{\text{OC}}}&\mathbf{H^{\text{CC}}}\\ \end{array}\right)\left(\begin{array}[]{c}\mathbf{\bar{\bm{\kappa}}}\\ \mathbf{\bar{c}^{\Theta}}\end{array}\right)=-\left(\begin{array}[]{c}\mathbf{g^{\Theta}}\\ 0\end{array}\right)$$

$$\frac{d{E}^{\Theta }}{dx}$$

$$-\sum _{pq}{X}_{pq}^{\Theta ,\text{e}}\frac{\partial S_{pq}}{\partial x}$$

$${\mathbf{\gamma }}^{\Theta ,\text{e}}$$

$${\mathbf{\Gamma }}^{\Theta ,\text{e}}$$

$${\tilde{\gamma }}_{pq}$$

$${\bar{\gamma }}_{pq}$$

$$\displaystyle\tilde{\Gamma}_{pqrs}=\sum_{\Psi}\omega_{\Psi}\Bigg{(}\sum_{o}\Gamma^{\Psi}_{oqrs}\bar{\kappa}^{\Theta}_{op}$$

$$+{\Gamma }_{pqos}^{\Psi }{\bar{\kappa }}_{or}^{\Theta }$$

$${\bar{\Gamma }}_{pqrs}=\sum _{\Psi }2{\omega }_{\Psi }{\gamma }_{pqrs}^{\Psi {\bar{\mathbf{c}}}_{\mathbf{\Psi }}^{\mathbf{\Theta }}}.$$

$$|\Psi ⟩$$

$$=\sum _{\sigma }\sum _{a_1\dots a_{L-1}}{M}_{1a_1}^{{\sigma }_1}{M}_{a_1{a}_2}^{{\sigma }_2}\cdots M_{a_{L-1}1}^{{\sigma }_L}|\mathbf{\sigma }⟩.$$

$$|\Psi ⟩$$

$$\hat{\mathcal{W}}=\sum _{\sigma ,{\sigma }^{\prime }}\sum _{b_1\dots b_{L-1}}{W}_{1b_1}^{{\sigma }_1{\sigma }_1^{\prime }}{W}_{b_1{b}_2}^{{\sigma }_2{\sigma }_2^{\prime }}\cdots W_{b_{L-1}1}^{{\sigma }_L{\sigma }_L^{\prime }}|\mathbf{\sigma }⟩⟨{\mathbf{\sigma }}^{\prime }|.$$

$$\delta \left(⟨\Psi ⟩\hat{H}\Psi -\lambda (⟨\Psi ⟩\Psi -1)\right)=0$$

$$⟨\Psi ⟩\hat{H}\Psi$$

$$\displaystyle=\mkern-27.0mu\sum_{\underset{\sigma_{L}\sigma^{\prime}_{L}}{a_{L-1}a^{\prime}_{L-1}b_{L-1}}}\mkern-27.0muM^{\sigma_{L}\dagger}_{1a_{L-1}}W^{\sigma_{L}\sigma^{\prime}_{L}}_{b_{L-1}1}\Big{(}\ldots\sum_{\underset{\sigma_{2}\sigma^{\prime}_{2}}{a_{1}a^{\prime}_{1}b_{1}}}M^{\sigma_{2}\dagger}_{a_{2}a_{1}}W^{\sigma_{2}\sigma^{\prime}_{2}}_{b_{1}b_{2}}\Big{(}\sum_{\sigma_{1}\sigma^{\prime}_{1}}M^{\sigma_{1}\dagger}_{a_{1}1}W^{\sigma_{1}\sigma^{\prime}_{1}}_{1b_{1}}{M}^{\sigma^{\prime}_{1}}_{1a^{\prime}_{1}}\Big{)}{M}^{\sigma^{\prime}_{2}}_{a^{\prime}_{1}a^{\prime}_{2}}\cdots\Big{)}{M}^{\sigma^{\prime}_{L}}_{a^{\prime}_{L-1}1}$$

$${\mathbb{L}}_{a_0{a}_0^{\prime }}^{b_0}$$

$${\mathbb{L}}_{a_1{a}_1^{\prime }}^{b_1}$$

$$\cdots$$

$${\mathbb{L}}_{a_l{a}_l^{\prime }}^{b_l}$$

$${\mathbb{R}}_{a_{l-1}^{\prime }{a}_{l-1}}^{b_{l-1}}=\sum _{\underset{{\sigma }_l{\sigma }_l^{\prime }}{a_l{a}_l^{\prime }{b}_l}}{M}_{a_{l-1}^{\prime }{a}_l^{\prime }}^{{\sigma }_l^{\prime }}{W}_{b_{l-1}{b}_l}^{{\sigma }_l{\sigma }_l^{\prime }}{\mathbb{R}}_{a_l^{\prime }{a}_l}^{b_l}{M}_{a_l{a}_{l-1}}^{{\sigma }_l†}.$$

$${\mathcal{H}}_{IJ}$$

Full text

I Introduction

Multiconfigurational methods are tailored for molecular systems exhibiting strong electronic correlation, as found in many transition metal complexesPierloot_Calculations_TaCC_2005 ; Neese_Advanced_CCR_2007 ; Gagliardi_Transition_RiCC_2007 ; Neese_Correlated_JBIC_2011 ; Pierloot_Transition_IJQC_2011 ; Daniel_Photochemistry_CCR_2015 ; Ashley_Ironing_CCR_2017 ; Radon_Benchmarking_PCCP_2019 ; Vogiatzis_Computational_CR_2019 , bond dissociation processes and excited electronic states.Gonzalez_Progress_C_2012 ; Plasser_Electronically_TCA_2012 ; Daniel_Photochemistry_CCR_2015 ; Ghosh_Combining_CR_2018 ; Lischka_Multireference_CR_2018 Excited electronic states are of key importance for photoinduced phenomena, including light-matter interaction with DNASerrano-Andres_Are_JPPCPR_2009 ; Middleton_DNA_ARPC_2009 ; Mai_Excitation_2014 , light harvesting, photocatalysis and artificial photosynthesisGust_Solar_ACR_2009 ; Gust_Realizing_FD_2012 ; Jager_Using_CCR_2015 ; Agbe_Recent_JIEC_2019 , as well as photodynamic cancer therapyDolmans_Photodynamic_NRC_2003 ; Rose_Fiat_COCB_2008 ; Ormond_Dye_M_2013 ; Mehraban_Developments_M_2015 .

The majority of modern multiconfigurational methods are based on the complete active space self-consistent field (CASSCF) methodroos80 ; sieg81 ; shep87 ; olse11 , which, however, scales exponentially with the number of active orbitals, and as such, has been effectively limited to about 18 active orbitals.aqui15a Only recently, this limit has been raised to 20 orbitals through massive parallelization.Vogiatzis_Pushing_JCP_2017 One successful remedy to the exponential scaling problem of CASSCF has been the quantum chemical density matrix renormalization group (DMRG)White1998Ab ; mitrushenkov2001quantum ; chan2002highly ; legeza2003controlling ; chan2004algorithm ; moritz2005convergence ; moritz2005relativistic ; rissler2006measuring ; mcculloch2007density ; lege08 ; chan2009density ; mart10 ; schollwock2011density ; chan2011density ; Kurashige2009High ; Knecht2014Communication ; wout14 ; Yanai2015Density ; knec16a ; Chan2016Matrix . DMRG was originally introduced by Whitewhite1992density ; white1992real to solid state physics in 1992, but has found numerous applications in quantum chemistry since then.White_initio_JCP_1999 ; Daul_Full-CI_IJQC_2000 ; Mitrushenkov_Quantum_JCP_2001 ; Mitrushenkov_Quantum_JCP_2001 ; Mitrushenkov_Quantum_JCP_2003 ; Mitrushchenkov_importance_IJQC_2011 ; Chan_Highly_JCP_2002 ; Chan_Exact_JCP_2003 ; Chan_algorithm_JCP_2004 ; Chan_State-of-the-art_JCP_2004 ; Chan_Density-matrix_JCP_2005 ; Legeza_QC-DMRG_MP_2003 ; Moritz_Convergence_JCP_2005 ; Moritz_Relativistic_JCP_2005 ; hachmann2007radical ; Zgid_spin_JCP_2008 ; Zgid_spin_JCP_2008 ; Zgid_density_JCP_2008 ; Zgid_Obtaining_JCP_2008 ; Kurashige_High-performance_JCP_2009 . Combined with self-consistent-field orbital optimization (DMRG-SCF)zgid2008density ; ghosh2008orbital ; sun2017 ; wouters2014communication ; Ma_Second-Order_JCTC_2017 , the method is able to approximate a CASSCF wave function to arbitrary accuracy with polynomial scalingchan2002highly by reducing the size of the configurational space by means of the optimization of a so-called matrix product state (MPS) wave functionOstlund_Thermodynamic_PRL_1995 ; Rommer_Class_PRB_1997 , therefore allowing for much larger active spaces than traditional CASSCF.

Many problems in quantum chemistry tackled by DMRG and DMRG-SCFChan_Chapter_ARiCC_2009 ; mizukami2012more ; kurashige2013entangled ; Chalupsky_JAmChemSoc_Reactivity_2014 ; sharma2014low ; yana15 ; freitag2015orbital ; Olivaresamaya2015The ; Wouters_density_EPJD_2014 require only the calculation of electronic energies and properties. However, the study of photochemical phenomena often relies on optimizing excited-state structures, locating conical intersections and potential energy surface crossings, or performing ab-initio molecular dynamic simulations. These tasks require efficient calculation not only of the energy of the ground and excited electronic states, but also of energy gradients and nonadiabatic couplings.

The most efficient way to calculate gradients is through the derivation of analytical expressions, as originally introduced by PulayPulay_initio_MP_1969 ; Pulay_Direct_AoEST_1977 ; Pulay_Systematic_JACS_1979 ; Pulay_efficient_TCA_1979 . Accordingly, an analytical formulation for nonadiabatic couplings has been developed by Yarkony and co-workers.Iii1984On

Analytical gradients are easily evaluated for a fully variationally-optimized wavefunction with the help of the Hellmann-Feynman theorem.Hellmann ; Feynman_Forces_PR_1939 Single-state multiconfigurational self-consistent-field (MCSCF) wavefunctions are fully variationally optimized, and hence the state-specific CASSCF analytical gradient formulation Taylor_Analytical_JCC_1984 has appeared shortly after CASSCF has been introduced. The same holds true for single-state DMRG-SCF gradients, which have been introduced by Liu2013Multireference and Hu2015Excited and have found several applications in ground and excited state structure optimizationsLiu2013Multireference ; Hu2015Excited ; Nakatani_JChemPhys_Density_2017 and resonance Raman spectraMa2017Multiconfigurational .

However, state-specific MCSCF methods have disadvantages, especially when applied to excited states. First, a specific state may not be tracked easily. For example, the so-called root flipping, i. e. a change of the excited state order during orbital optimization may occur (although the problem can be partially alleviated if the state with the maximum overlap to the state of interest is followed during the optimization)Hu2015Excited ; Ma2017Multiconfigurational . Second, orthogonality, which greatly simplifies the calculation of excited-state properties and, in particular, of nonadiabatic couplingsIii1984On ; Yarkony_Electronic_1995 ; Schmidt_AnnuRevPhysChem_construction_2003 is not guaranteed in a state-specific optimization.

One possible remedy for the disadvantages of the state-specific optimization is the state-average MCSCF (SA-MCSCF) ansatzDocken_JChemPhys_LiH_1972 ; Werner_JChemPhys_quadratically_1981 ; Diffenderfer_JPhysChem_Use_1982 , where several states are optimized simultaneously using a single set of molecular orbitals (MOs) which yields the best average energy of the states of interest. In the state-average ansatz, the root flipping problem does not occurDocken_JChemPhys_LiH_1972 and the states are necessarily orthogonal. However, with SA-MCSCF the gradient formalism loses its simplicity as the wavefunction is no longer fully variational. Nevertheless, the Hellmann-Feynman theorem can still be applied with the help of Lagrange multipliersHelgaker_Configuration-interaction_TCA_1989 ; Jones2001Analytical that are obtained from the solution of the so-called coupled perturbed multiconfiguration self-consistent field (CP-MCSCF) equations.Dupuis_Energy_JCP_1981 ; Osamura_Generalization_JCP_1982 Jones2001Analytical were the first to describe an SA-MCSCF analytic gradient formulation based on Lagrangians. Recently, Snyder_atomic_JCP_2015Snyder_atomic_JCP_2015 ; Snyder_direct-compatible_JCP_2017 presented a GPU-based implementation of SA-MCSCF analytic gradients, following the work of Jones2001Analytical

Besides the gradients, the solution of the CP-MCSCF equations may also be used to calculate analytical nonadiabatic couplings in the SA-MCSCF formalismIii1984On ; Lengsfield_Nonadiabatic_ACP_2007 and for the second-order MCSCF orbital optimization procedure.wern80 ; wern85 ; know85 ; wern87 ; leng81 ; leng82 ; jorg81 ; yeag82 ; Olsen_Optimization_AiCP_1983

While the SA-CASSCF analytical gradient problem may be considered solved, this is not the case for SA-DMRG-SCF analytical gradients. The main challenge that remains is an adequate definition of the Lagrange parameters for the gradient calculation. For traditional MCSCF wavefunctions one usually expresses these in the configuration basis, which is impractical for MPS wavefunctions optimized by DMRG due to the exponential growth of the number of configurations with the number of active orbitals. The straightforward approach of employing all matrices in an MPS as parameters would result in a highly redundant parameterization, not least due to the gauge freedom of an MPS.Wouters_Thouless_PRB_2013 ; Nakatani_Linear_JCP_2014 . Nonredundant parametrizations of an MPS existDorando_Analytic_JCP_2009 ; Haegeman_Time-Dependent_PRL_2011 ; Wouters_Thouless_PRB_2013 ; Nakatani_Linear_JCP_2014 ; Vanderstraeten_Tangent-space_SPLN_2019 , the majority being based on the tangent-space concept. However for large active spaces (which are often found in typical applications of SA-DMRG-SCF in quantum chemistry) the number of MPS parameters, even if they are non-redundant, would still become very large, possibly rendering the gradient calculation prohibitively expensive. Therefore, it is helpful to introduce an approximate parametrization, which does not necessarily cover the full MPS parameter space, but nevertheless allows for fast and efficient estimation of gradients. In our previous work on a second-order DMRG-SCF optimization schemeMa_Second-Order_JCTC_2017 we derived the expression of the state-average energy Hessian based on an ansatz for the MPS variational parameters. In this work, we generalize this ansatz and employ it in a fast and efficient implementation to provide an efficient approximation for SA-DMRG-SCF analytical gradients and nonadiabatic couplings.

II Theory

As the SA-DMRG-SCF analytical gradient theory is closely related to the SA-CASSCF analytical gradient theory, we will begin with a brief recap of the latter. Further information can be found in Ref. 105 and especially the excellent paper by Snyder_direct-compatible_JCP_2017 . In particular, we adopt the notation of the latter work. The theory for the nonadiabatic couplings follows very closely the gradient theory – the differences between the two will be presented in Sec. II.6.

II.1 Single-state CASSCF gradients

In single-state (or state-specific) CASSCFroos1980complete , a wavefunction is defined as a linear combination of configurations $\ket{\phi_{I}}$:

[TABLE]

where the expansion coefficients $c_{I}$ are called the configuration interaction (CI) coefficients and the configurations $\phi_{I}$ are chosen in such a way that they represent all possible excitations within a predefined active orbital space and do not include any excitations outside this space. The CI coefficients are obtained by diagonalizing the matrix of the non-relativistic electronic Hamiltonian

[TABLE]

in the basis of the configurations $\phi_{I}$. Here, $\braket{t}{h}{u}$ and $(tu|vw)$ are the one- and two-electron integrals in a given orthonormal molecular orbital basis, respectively, and $a^{\dagger}_{t\tau}$ and $a_{t\tau}$ are creation and annihilation operators, respectively, for an orbital $t$ and spin $\tau$. $E_{\text{core}}$ is the sum of the energy contribution from the inactive (doubly occupied) orbitals and the nuclear repulsion energy.

The wavefunction $\ket{\Phi}$ is a function of a set of variational parameters, namely the orbital rotation parameters $\mathbf{\bm{\kappa}}$ and CI coefficients ${\bm{c}}$

[TABLE]

and is variationally optimized with respect to these parameters to yield a minimum single-state CASSCF energy

[TABLE]

with $\gamma_{pq}$ and $\Gamma_{pqrs}$ being the one- and two-particle reduced density matrices (RDMs), respectively. Since the optimized wavefunction $\ket{\Phi}$ is then fully variational with respect to all of its variational parameters, the energy gradient with respect to some perturbation $x$ (e. g. nuclear displacement) may be calculated according to the the Hellmann-Feynman theoremszalay2011multiconfiguration :

[TABLE]

II.2 SA-CASSCF analytical gradients

In SA-CASSCF, the wavefunction is determined by a variational optimization of the state-average energy of several states $\Psi$

[TABLE]

where $\omega_{\Psi}$ is the weight of state $\Psi$, and the sum of all weights equals to 1. In this paper we will only consider a situation where all weights are equal, although Jones2001Analytical as well as Snyder_direct-compatible_JCP_2017 have also considered a non-equal-weights situation. Similarly to the single-state case, the state-average energy depends on a set of variational parameters

[TABLE]

where $\mathbf{\bm{\kappa}}$ are the orbital parameters and ${\bm{c}}^{\Psi}$ are the CI coefficients for state $\Psi$.

The optimization procedure ensures that the state-average energy, but not the individual state-specific energies, is variational with respect to the orbital parameters and the CI parameters.

[TABLE]

Hence, the gradient cannot be calculated according to the Hellmann-Feynman theorem. Following the method of the Lagrange multipliers, we may construct a Lagrangian $L^{\Theta}$ for state $\Theta$ starting from the energy in Eq. (4)

[TABLE]

where $\bar{\kappa}_{pq}^{\Theta}$ and $\bar{c}_{\Psi I}^{\Theta}$ are the Lagrange multipliers corresponding to the orbital and CI parts, respectively. $L^{\Theta}$ is variational with respect to all parameters and, therefore, the Hellmann-Feynman theorem applies. The Lagrange multipliers can be obtained by exploiting the fact that the Lagrangian must be variational with respect to all of its parameters:

[TABLE]

The first two equalities in Eq. (10) are trivially fulfilled and are equivalent to the definition of the constraints. Hence, the Lagrange multipliers are obtained from the other two equalities, yielding the coupled perturbed MCSCF (CP-MCSCF) equations:

[TABLE]

or in matrix form

[TABLE]

As can be seen from Eqs. (11) and (12), $\mathbf{H}$ is the state-average Hessian matrix with its orbital-orbital ($\mathbf{H}^{\text{OO}}$), orbital-CI ($\mathbf{H}^{\text{OC}}=(\mathbf{H}^{\text{CO}})^{T}$) and CI-CI ($\mathbf{H}^{\text{CC}}$) components and $\mathbf{g^{\Theta}}$ is the state-specific orbital gradient for state $\Theta$.

After solving Eq. (13), we may construct the Lagrangian of Eq. (9), and calculate the gradient as followsSnyder_atomic_JCP_2015 ; Snyder_direct-compatible_JCP_2017 :

[TABLE]

The first two terms in Eq. (14) originate from the Hellmann-Feynman theorem. They seem identical to Eq. (5), but the one- and two-particle reduced density matrices have been replaced by their effective counterparts, $\bm{\gamma}^{\Theta,\text{e}}$ and $\bm{\Gamma}^{\Theta,\text{e}}$, that are given bySnyder_direct-compatible_JCP_2017 :

[TABLE]

where $\bm{\tilde{\gamma}}$ and $\bm{\bar{\gamma}}$, as well as $\bm{\tilde{\Gamma}}$ and $\bm{\bar{\Gamma}}$ are the orbital and the CI contributions to the RDM, respectively:

[TABLE]

[TABLE]

Here $\gamma^{\Psi\bm{\bar{c}^{\Theta}_{\Psi}}}_{pq}$ and $\Gamma^{\Psi\bm{\bar{c}^{\Theta}_{\Psi}}}_{pqrs}$ are matrix elements of the one- and two-particle transition density matrices, respectively, between state $\Psi$ and a state with the Lagrange parameters $\bm{\bar{c}^{\Theta}_{\Psi}}$ as the CI coefficients. The last term in Eq. (14) is dubbed the “connection” term by Jones2001Analytical , and is evaluated from the effective CI Lagrangian and the derivative of the MO overlap matrix.Snyder_atomic_JCP_2015 ; Snyder_direct-compatible_JCP_2017

In practice, the CP-MCSCF equations (13) are not solved directly, since the cost of evaluating and storing the full Hessian is too large. Instead, iterative solvers such as the preconditioned conjugate gradient (PCG)Jones2001Analytical or the direct inversion of the iterative subspace (DIIS)Snyder_direct-compatible_JCP_2017 are employed, which only require the computation of the product between the Hessian and a trial vector, i. e. $\mathbf{H^{\text{OO}}}\tilde{\bm{\kappa}}^{\Theta},\mathbf{H^{\text{CO}}}\tilde{\bm{\kappa}}^{\Theta},\mathbf{H^{\text{CO}}}\tilde{\bm{c}}^{\Theta}$ and $\mathbf{H^{\text{CC}}}\tilde{\bm{c}}^{\Theta}$, with $\tilde{\bm{\kappa}}^{\Theta}$ and $\tilde{\bm{c}}^{\Theta}$ as trial vectors, which become equal to the corresponding Lagrange multipliers upon convergence. Although these quantities must be recalculated each iteration (by contrast, the Hessian needs to be calculated only once), one must store only a trial vector of size $n_{\bm{\kappa}}+n_{\bm{c}}$, where $n_{\bm{\kappa}}$ is the number of orbital rotation parameters and $n_{\bm{c}}$ the number of CI coefficients, instead of a full Hessian matrix. In addition, as will become evident in Section II.5, evaluating certain matrix-vector products offers further computational advantages compared to evaluating the corresponding Hessian matrix elements.

II.3 DMRG and DMRG-SCF

Most commonly, DMRG for quantum chemistry is formulated on the basis of an MPS wavefunction parametrizationOstlund_Thermodynamic_PRL_1995 ; Rommer_Class_PRB_1997 ,

[TABLE]

Hence, the CI coefficients in Eq. (1) are encoded as a product of three-dimensional tensors $\mathbf{M}^{\sigma_{l}}$, $\ket{\bm{\sigma}}=\ket{\sigma_{1},\ldots,\sigma_{L}}$ represents the occupation number vector in analogy to the configurations in Eq. (1), and $L$ is the number of active orbitals. We may occasionally call the quantities $M^{\sigma_{l}}_{a_{l-1}a_{l}}$ the elements of an MPS tensor at site (orbital) $l$. For the optimization, it is important to express the MPS in a mixed-canonical formschollwock2011density at an arbitrary site $l$:

[TABLE]

where the tensors with elements ${A}^{\sigma_{i}}_{a_{l-1}a_{i}}$ are left-normalized and tensors with elements ${B}^{\sigma_{i}}_{a_{l-1}a_{i}}$ are right-normalized, respectively Schollwock_density-matrix_PTRSA_2011 .

In contrast with the CASSCF/CI wavefunction in Eq. (1), for which all the CI coefficients are determined in one step, DMRG optimizes MPS wavefunctions iteratively with one (or two adjacent, see below) $\mathbf{M}^{\sigma_{l}}$ MPS tensors at a time, while ensuring that maximum dimensions of these matrices do not exceed a certain value $m$, which is denoted as the bond dimension or number of renormalized block states. Due to this systematic dimension reduction, the exponential scaling of Eq. (1) is reduced to polynomial scaling. In analogy to Eq. (21), also operators may be represented as matrix product operators (MPOs):

[TABLE]

(We shall not go into detail here on how to obtain the MPO representation for various operators: this has been described in detail in our previous workkell15a .)

In that case, the optimization of a MPS wavefunction is formulated as the variational optimization of the entries of a single MPS tensor (or two adjacent MPS tensors) to minimize the expectation value of the energy $E^{\Psi}=\braket{\Psi}{\hat{H}}{\Psi}$, under the constraint that the wavefunction is normalized:

[TABLE]

Inserting the MPO expression for the Hamiltonian $\hat{H}$ into the expression for the expectation value of the energy, we obtain

[TABLE]

where we have regrouped the multiplication by indices $\sigma_{i},\sigma^{\prime}_{i}$ in the last step. We may now define left boundaries recursively asschollwock2011density

[TABLE]

and analogously right boundaries

[TABLE]

Assuming that the left and right boundaries have been constructed from left- and right-normalized tensors, respectively, inserting the definition of boundaries into Eq. (24) yields an eigenvalue equation $\mathcal{H}v=\lambda v$, where $\mathcal{H}$ is the local Hamiltonian matrix at site $l$ with matrix elements

[TABLE]

where we combined one set of indices $(a_{l-1}\sigma_{l}a_{l})$ into one composite index to form the matrix $\mathcal{H}$ and the vector $v$. After the lowest (or several lowest, see below) eigenvalue(s) $\lambda$ have been obtained, the corresponding eigenvector(s) are reshaped again into the MPS tensor $M^{\sigma_{l}}_{a_{l-1}a_{l}}$. A subsequent normalization and truncation procedure such as singular value decomposition (SVD) ensures that the maximum dimension of ${M}^{\sigma_{l}}_{a_{l-1}a_{l}}$ does not exceed $m$. A basis transformation finalizes the local optimization at site $l$, generating a new mixed-canonical form of the MPS at site $l+1$. The process (sweep) is repeated until the final site $L$ is reached and then its direction is reversed. In passing, we note that commonly two adjacent MPS tensors are optimised simultaneously (two-site DMRG)schollwock2011density , which is, however, not an important aspect for the purpose of this work. For further details on the DMRG optimization procedure we refer the reader to Refs. 43; 127; 128.

After optimization of the MPS wavefunction, we can easily obtain one- and two-particle RDMs as expectation values of operators within the MPS-MPO framework of DMRG.schollwock2011density ; Schollwock_density-matrix_PTRSA_2011 This allows us to formulate a (SA-)DMRG-SCF procedure in analogy with the (SA-)CASSCF one, with the CI coefficient optimization step in CASSCF replaced by a DMRG procedure in DMRG-SCFzgid2008density ; Ma_Second-Order_JCTC_2017 . Importantly, Eqs. (4) and (6) remain the same.

II.4 Definition of variational parameters for analytical gradients for SA-DMRG-SCF

The derivation of the state-average gradient for an MPS wavefunction by means of Eqs. (10)-(13) is far from trivial. Unlike the expansion of Eq. (1), the CI coefficients are not explicitly available for a MPS wavefunction in Eq. (21). Hence, one must define an analogous set of MPS parameters which are equivalent to the CI parameters ${\bm{c}}^{1},{\bm{c}}^{2},\ldots,{\bm{c}}^{n}$ in Eq. (7).

As mentioned in the Introduction, in this work we would like to avoid considering the set with all non-redundant entries of the MPS as free variational parameters. Instead, we introduce a much simpler, albeit an approximative parametrization of an MPS in the mixed canonical form where only the entries of the tensor associated to the canonization site are included in the variational space. We will justify the choice of our parametrization below.

Let us define two auxiliary MPSs in orbital subspaces spanned by sites 1 to $l-1$ and $l+1$ to $L$:

[TABLE]

Inserting these equations into the mixed-canonical representation of our MPS at site $l$ (Eq. (22)) yields

[TABLE]

where $\ket{\sigma_{l}}$ is the local basis state for orbital $l$, corresponding to its four possible occupations: unoccupied, spin-up, spin-down, and doubly occupied. Grouping the indices $(a_{l-1}\sigma_{l}a_{l})$ as in Eqs. (30) and (31), we arrive at an expression for $\ket{\Psi}$ equivalent to Eq. (1):

[TABLE]

(From now on indices $I$ and $J$ in equations containing MPS will refer to the grouped indices $(a_{l-1}\sigma_{l}a_{l})$, in analogy to CI configurations, which are also labelled $I$ in Eq. (1).

With a wavefunction written in terms of a CI expansion, we may define MPS Lagrange multipliers (or variational parameters) $\bar{\bm{v}}^{\Psi}$ in a basis spanned by $v_{I}$ in analogy to the CI Lagrange multipliers $\bar{\bm{c}}^{\Psi}$ and construct a Lagrangian for a DMRG-SCF wavefunction similarly to Eq. (9) employing MPS Lagrange multipliers.

Such a definition of the variational parameters places a constraint on some properties of the MPS. The configuration basis $\ket{\phi_{i}}$ in Eq. (1) is, by definition, the same for all states in the state-average description. This is not necessarily true for a set of MPS, if they are optimized individually because $\ket{a_{l-1}}$ and $\ket{a_{l}}$ (Eq. 34) will be different for each state. Moreover, in this case not even the dimensions of the individual ${M}^{\sigma_{l}}_{a_{l-1}a_{l}}$ may be the same for all states. Consequently, for our definition of MPS parameters to be valid, the basis $\ket{\phi_{i}}$ in Eq. (35) must also be the same for all states in the state-average description. This can be ensured if all MPS are optimized and brought into a mixed-canonical form simultaneously, where the truncation and the basis transformation step during the sweeps are the same for all states. In practice, this means employing common left- and right boundaries for all states by following the procedure described in Ref. 129. We have implemented such a procedure in our new DMRG MPS-MPO solver: the details of the implementation are, however, beyond the scope of this work and will be described in a future publication.

It remains to be specified at which site $l$ we write the MPS in Eq. (22) to define our MPS parameters, as the choice is, in principle, arbitrary. However, the site must be chosen a priori before solving the CP-MCSCF equations. We will call this site the linear response site for convenience. At first glance, there should be no difference, because for an optimized MPS wavefunction all mixed canonical forms at all possible linear response sites $l$ are equivalent. However, the parameters $v_{I}$ and the bases $\ket{\phi_{I}}$ for different linear response sites are not equivalent. In fact, here the approximative character of our MPS parametrization is revealed. We consider the variation of MPS parameters $v_{I}$ at a single site, but not the variation of the basis $\ket{\phi_{I}}$ that depends on the MPS tensors on the sites different from the canonization one. With sites towards the middle of the lattice, the number of parameters grows, and therefore the variational parameter space covered by the MPS parameters $v_{I}$ is enlarged. As the MPS approaches the full CI wavefunction, the variation of the basis $\ket{\Phi_{I}}$ approaches zero and the error arising from neglecting the variation of the basis (which we call basis-variation omission error (BVOE)) vanishes, although, in general, for approximate full CI wavefunctions it is nonzero. We will assess the accuracy of the gradient calculation for different linear response sites on an example in Section III.2, and quantify BVOE in Section III.5.

In passing we note that in Ref. 75 we have employed a similar definition of MPS parameters, where, in contrast to this work, i) two-site tensors at the first two sites have been employed in the definition instead of one-site MPS tensors at an arbitrary site and ii) the requirement of the same local basis $\ket{\Phi_{i}}$ for several states was not met. We consider the parametrization employed in this work as an improvement to the one from Ref. 75, as a response site in the middle of a lattice delivers larger variational flexibility of the wavefunction compared to the sites at the edge.

II.5 Implementation of SA-DMRG-SCF analytical gradients

The definition of MPS parameters in Subsection II.4 and the identification of a local optimization of an MPS tensor ${M}^{\sigma_{l}}_{a_{l-1}a_{l}}$ as CI problem allow us to formulate the SA-DMRG-SCF analytical gradient theory very closely following the SA-CASSCF theory, i. e. by applying Eqs. (10)-(13) to SA-DMRG-SCF wavefunctions (which have been optimized as described in Subsection II.4). This has another advantage, namely that we may derive our implementation from an existing SA-CASSCF analytical gradient one. In this work, we have based our implementation on the SA-CASSCF analytical gradient implementation in OpenMOLCAS.Bernhardsson_MolPhys_direct_1999 ; Jones2001Analytical ; OpenMOLCAS It employs the PCG method to solve the CP-MCSCF equations (13) iteratively, evaluating the Hessian-trial vector products in each iteration.

Below we shall outline the quantities required for the analytical gradient calculation, whose evaluation is specific to the MPS wavefunction.

Transition density matrices.

The evaluation of transition density matrix elements may be performed with the help of Eq. (25)kell15a . For an MPS

[TABLE]

an expectation value $\braket{\Psi}{\hat{\mathcal{W}}}{\bar{\Psi}}$ may be calculated as

[TABLE]

To exploit this for the one- and two-particle reduced transition density matrix elements $\gamma^{\Psi\bar{\Psi}}_{pq}$ and $\Gamma^{\Psi\bar{\Psi}}_{pqrs}$, we replace $\hat{\mathcal{W}}$ in Eq. (37) with the MPO form of the operators

[TABLE]

respectively, and evaluate Eq. (37). (For details on obtaining the MPO form of various operators, see Ref. 128).

If the bra and ket in Eq. (37) differ only by a single tensor, then, instead of recalculating the full Eq. (37) for every expectation value one must calculate all contractions for all sites $\neq l$ only once for all expectation values and recalculate only the contractions involving tensors at site $l$. During the solution of the CP-MCSCF equations, only the tensor at site $l$ is updated. This can be seen in Eqs. (41) and (72) in the work of Snyder_direct-compatible_JCP_2017 , who update the CI coefficients but retain the basis intact; by analogy, we update only the vector $v_{I}$ in Eq. (35) and therefore only the MPS tensor at site $l$. Furthermore, in the framework of SA-DMRG-SCF, $\ket{\bar{\Psi}}$ always differs from $\ket{\Psi}$ only by a single tensor, as multiple states are optimized with a common boundary. However, after a complete sweep the differing site is always at the beginning of the lattice. Hence, to bring the MPS into a usable form, we must perform a simultaneous canonization of all states. Details on this procedure are beyond the scope of this work and will be described in a future publication.

Sigma vectors.

Calculation of Hessian-trial vector products requires calculation of sigma vectors. For a CI wavefunction Eq. (1), a sigma vector is defined as the left-hand side of the CI eigenvalue equation, i. e.

[TABLE]

where $H_{ij}$ is the matrix element of the Hamiltonian. (Note that in contrast to other literature, e. g. Refs. 105, 109, 132 we denote the sigma vector with a capital $\Sigma$ to avoid confusion with MPS indices $\sigma$). An expression for the sigma vector for MPS wavefunctions is easily obtained from Eqs. (30) and (31)

[TABLE]

Also in DMRG this equation constitutes the left-hand side of the local eigenvalue equation for the optimization of an MPS tensor, and, as such, belongs to the core part of our DMRG MPS-MPO implementation.kell15a

Diagonal of the Local Hamiltonian.

Also diagonal elements of the local Hamiltonian $\mathcal{H}^{\sigma_{l}\sigma^{\prime}_{l}}_{a_{l-1}a_{l}a^{\prime}_{l-1}a^{\prime}_{l}}$ from Eq. (30) are required in our gradient implementation to calculate the preconditioner for the PCG algorithm.

With the above quantities, we may setup the CP-MCSCF equations and evaluate the gradients. We will closely follow the work of Snyder_direct-compatible_JCP_2017 on SA-CASSCF gradients and present here only DMRG-specific steps in order to be as brief as possible.

Equivalently to SA-CASSCF theory, we must set up and solve the CP-DMRG-SCF equations

[TABLE]

where the only difference to Eq. (13) is that the Lagrange CI parameters $\mathbf{\bar{c}^{\Theta}}$ have been replaced with the MPS parameters $\bar{\bm{v}}^{\Theta}$. As we employ the PCG solverJones2001Analytical for the CP-DMRG-SCF equations, we must be able to compute the Hessian-trial vector products $\mathbf{H^{\text{OO}}}\tilde{\bm{\kappa}}^{\Theta},\mathbf{H^{\text{CO}}}\tilde{\bm{v}}^{\Theta},\mathbf{H^{\text{CO}}}\tilde{\bm{\kappa}}^{\Theta}$ and $\mathbf{H^{\text{CC}}}\tilde{\bm{v}}^{\Theta}$ with $\tilde{\bm{\kappa}}^{\Theta}$ and $\tilde{\bm{v}}^{\Theta}$ as trial vectors (which become equal to $\mathbf{\bar{\bm{\kappa}}}$ and $\bar{\bm{v}}^{\Theta}$ at convergence). $\mathbf{H^{\text{OO}}}\tilde{\bm{\kappa}}^{\Theta}$ and $\mathbf{g^{\Theta}}$ are evaluated exactly as for SA-CASSCF (see, e. g., Eqs. (13)-(21) in Snyder_direct-compatible_JCP_2017 ). To express the remaining Hessian-trial vector products, we introduce the generalized orbital gradient matrix $\mathbf{T}(\Psi,\bar{\Psi})$. This matrix has been introduced first in Snyder_direct-compatible_JCP_2017 (see Eqs. (35)-(37) therein). However, in our implementation we followed Refs. 130; 105 and the standard textbook by Helgaker, Jørgensen and Olsenhelgaker2014molecular . The equations for the calculation of the matrix elements $T_{pq}(\Psi,\bar{\Psi})$ as employed in our implementation are provided below:

[TABLE]

where the generalized Fock matrix $F_{pq}(\Psi,\bar{\Psi})$ is calculated as

[TABLE]

if the first index is virtual. ${}^{\text{I}}F_{pq}$ and ${}^{\text{A}}F_{pq}(\Psi,\bar{\Psi})$ are the inactive and active Fock matrices, respectively:

[TABLE]

Eqs. (44)-(48) are generalizations of Eqs. (10.8.27)-(10.8.32) from Ref. 132 to transition density matrices between states $\Psi$ and $\bar{\Psi}$. In particular, the state-specific orbital gradient is obtained with $g_{pq}^{\Theta}=T_{pq}(\Theta,\Theta)$.

Similarly to Snyder_direct-compatible_JCP_2017 , the $\mathbf{H^{\text{CO}}}\tilde{\bm{v}}^{\Theta}$ product is evaluated as

[TABLE]

where $T_{pq}(\Psi,\tilde{\bm{v}}_{\Psi}^{\Theta})$ indicates a matrix element of the $\mathbf{T}$ matrix calculated for a state $\Psi$, where the MPS tensor at the linear response site $l$ has been replaced by the corresponding block of the trial vector of MPS parameters $\tilde{\bm{v}}_{\Psi}^{\Theta}$, and an unmodified state $\Psi$.

The $\mathbf{H^{\text{CO}}}\tilde{\bm{\kappa}}^{\Theta}$ product is evaluated from

[TABLE]

with

[TABLE]

The modified sigma vector $\tilde{\Sigma}^{\Psi}_{I}$ is evaluated from Eq. (41), where, however, a modified Hamiltonian

[TABLE]

with transformed integralshelgaker2014molecular

[TABLE]

[TABLE]

is employed instead of the Hamiltonian in the form of Eq. (2).

We should also note that calculating the $\mathbf{H^{\text{CO}}}\tilde{\bm{\kappa}}^{\Theta}$ product yields another computational advantage over the calculation of the full Hessian block $\mathbf{H^{\text{CO}}}$ (in addition to the obvious advantage of not storing the full Hessian). The calculation of $\mathbf{H^{\text{CO}}}$ would require one- and two-particle RDM derivatives with respect to MPS parameterswern80 ; Ma_Second-Order_JCTC_2017 , which we have shown in previous workMa_Second-Order_JCTC_2017 , has the cost of an evaluation of a transition density matrix for each single MPS parameter. As the number of MPS parameters may become extremely large, albeit limited by $m$, the cost of evaluating RDM derivatives would easily become a bottleneck of the calculation. However, calculating $\mathbf{H^{\text{CO}}}\tilde{\bm{\kappa}}^{\Theta}$ by means of Eq. (50) avoids the calculation of RDM derivatives altogether.

The $\mathbf{H^{\text{CC}}}\tilde{\bm{v}}^{\Theta}$ product, similarly to the CI case (cf. Eq. (39) in Ref. 109 or Eq. (30) in Ref. 105), is evaluated from

[TABLE]

Here we must again calculate a sigma vector $\bar{\Sigma}^{\Psi}_{I}$ with the help of Eq. (41), where, similarly to Eq. (49), the MPS tensor $\{M^{\sigma_{l}}_{a_{l-1}a_{l}}\}$ is replaced with the block of the trial vector of MPS parameters $\tilde{\bm{v}}^{\Theta}_{\Psi}$, reshaped according to the MPS tensor structure.

After solving the CP-DMRG-SCF equations we obtain a set of the orbital and MPS parameters $\mathbf{\bar{\bm{\kappa}}}$ and $\bar{\bm{v}}^{\Theta}$, which allow us to construct the effective one- and two-particle RDMs from Eqs. (15)-(20). The transition density matrices $\gamma^{\Psi\bm{\bar{v}^{\Theta}_{\Psi}}}_{pq}$ and $\Gamma^{\Psi\bm{\bar{v}^{\Theta}_{\Psi}}}_{pqrs}$, analogs of their CI counterparts $\gamma^{\Psi\bm{\bar{c}^{\Theta}_{\Psi}}}_{pq}$ and $\Gamma^{\Psi\bm{\bar{c}^{\Theta}_{\Psi}}}_{pqrs}$ from Eqs. (18) and (20), are obtained as in Eq. (49), namely by replacing the MPS tensor $\{M^{\sigma_{l}}_{a_{l-1}a_{l}}\}$ in the MPS $\Psi$ at the linear response site $l$ with the vector of MPS parameters $\bar{\bm{v}}^{\Theta}_{\Psi}$, reshaped accordingly, and evaluating the transition density matrices between the unmodified $\Psi$ and the modified state. Finally, the effective density matrices allow us to construct the Lagrangian in Eq. (9) and evaluate the gradient according to Eq. (14). This is completely identical to the SA-CASSCF procedure and has been described extensively in Refs. 105 and 109.

II.6 Nonadiabatic couplings

As in the SA-CASSCF case Fdez.Galvan_Analytical_JCTC_2016 ; Snyder_direct-compatible_JCP_2017 , the calculation of the nonadiabatic couplings for two states $\Theta$ and $\Lambda$, $\Braket{\Theta}{{\partial}/{\partial x}}{\Lambda}$, is very similar to the calculation of the gradients, with the following differences:

• •

The orbital gradient $\mathbf{g^{\Theta}}$ in the CP-MCSCF equations (13) or (42) is replaced by the generalized orbital gradient $\mathbf{T}(\Theta,\Lambda)$ with the average being calculated over two different states ($\Theta$ and $\Lambda$).

• •

In the expressions for the effective one- and two-particle RDMs (Eqs. (15) and (16)) the state-specific density matrices $\bm{\gamma^{\Theta}}$ and $\bm{\Gamma^{\Theta}}$ are replaced by the symmetrized transition density matrices $\bm{\gamma^{\Theta\Lambda}}$ and $\bm{\Gamma^{\Theta\Lambda}}$, respectively.

• •

For the calculation of the nonadiabatic couplings from the Lagrangian, Eq. (14) is scaled by $(E_{\Theta}-E_{\Lambda})^{-1}$ and an additional term enters, leading to

[TABLE]

where ${}^{\text{AS}}f_{pq}$ is the contribution that arises from the antisymmetric transition density matrix $\gamma^{\Theta\Lambda,\text{ns}}_{pq}$ and the antisymmetric derivative overlap $\left((\frac{\partial\phi_{p}}{\partial x}|q)-(p|\frac{\partial\phi_{q}}{\partial x})\right)$Snyder_direct-compatible_JCP_2017 . In SA-CASSCF nonadiabatic coupling theory, this contribution is termed “configuration state function (CSF) contribution”Fdez.Galvan_Analytical_JCTC_2016 , however here we prefer not to use this term as we do not operate in CSF basis in DMRG.

III Numerical examples

III.1 General computational details

In the numerical examples below all SA-CASSCF calculations have been performed with the OpenMOLCASOpenMOLCAS software package, and the SA-DMRG-SCF calculations with the development version of the QCMaquiskell15a ; kell16 DMRG program and the QCMaquis-OpenMOLCAS interface. In all examples, we have chosen the number of sweeps and the maximum number of renormalized block states $m$ in DMRG-SCF to approximate the CASSCF energy to an accuracy of $10^{-8}$ Hartree. Although DMRG-SCF makes energy and gradient calculations for large active spaces accessible, in this work we restrict ourselves to a problem that can still be treated with the standard CASSCF for the sake of comparison.

For our first example, cyclobutadiene (Section III.2), a full-valence active space, consisting of 12 electrons in 12 orbitals, with Dunning’s cc-pVDZdunn89 basis set was chosen. The optimized orbitals in both SA-CASSCF and SA-DMRG-SCF calculations turned out to be identical, as expected, and the maximum number of renormalized block states $m$ in DMRG-SCF was set to 2000.

For our second example, the 1,2-dioxetanone (Sections III.3-III.5), an active space of 16 electrons in 13 orbitals, along with the ANO-RCC-VDZP basis setRoos_New_JPCA_2005 was employed as suggested in Ref. 137. The maximum number of renormalized block states $m$ was raised to 5000, again to reproduce the CASSCF energy to an accuracy of $10^{-8}$ Hartree.

III.2 Cyclobutadiene: Dependence of the gradient on the choice of the linear response site

In this subsection we examine how the choice of the linear response site determines the error in gradients and nonadiabatic couplings, with cyclobutadiene as an example. The automerization reaction of cyclobutadiene is the most prominent example of a process primarily driven by heavy-atom tunnellingCarpenter_Heavy-atom_JACS_1983 ; Carsky_Heavy-atom_TCA_1992 ; Schoonmaker_Quantum_JCP_2018 . As such, it has been subject to ab-initio dynamic simulationsEckert-Maksic_Automerization_JCP_2006 , which necessitate the usage of nuclear gradients. A single molecular structure along the automerization reaction path from Ref. 141 has been chosen for our example below.

Figure 1 shows the mean absolute error of gradients of the S1 state and the nonadiabatic couplings between the S0 and S1 states for various linear response sites, along with the corresponding number of MPS parameters. The mean absolute error is defined as

[TABLE]

where $g_{i}^{\text{DMRG-SCF}}$ and $g_{i}^{\text{CASSCF}}$ are DMRG-SCF and CASSCF gradients or nonadiabatic couplings, respectively, and $n$ is their total number. Both errors clearly show their minimum in the middle of the lattice, where the number of MPS variational parameters of the wavefunction is the largest. Note that here, i. e., for sites 6 and 7, the number of MPS Lagrange multipliers is equal to the number of CI Lagrange multipliers in the CASSCF reference calculations, and therefore one would expect that in case of a perfect match of CASSCF and DMRG-SCF wavefunctions the gradients should also be identical. Still, we observe a small error of approximately $10^{-7}$ Hartree/bohr, which is attributed to the small numerical differences between CASSCF and DMRG-SCF wavefunctions, caused by nonzero convergence thresholds for the orbital optimization, DMRG sweeps, and the PCG method in CP-DMRG-SCF equations. By tightening the convergence thresholds the error decreases. Consequently, for the subsequent examples, we chose a linear response site in the middle of the lattice for the best accuracy.

With linear response sites towards the edge of the lattice, the gradient error increases by up to three orders of magnitude. This error arises from our single-site MPS parametrization and the BVOE because, as we already remarked above, tensors at the edges of the lattice have fewer elements leading to fewer MPS parameters. However, by varying the linear response site we may trade the gradient accuracy against the computational cost, as CP-DMRG-SCF equations with fewer Lagrange multipliers converge faster. In the middle of the lattice in our example the number of MPS parameters is equal to the number of CSFs in a CASSCF calculation, hence spanning the same variational space as in the CASSCF wavefunction: the BVOE approaches zero in this case. For SA-DMRG-SCF wavefunctions that approximate SA-CASSCF, as we will see in Section III.5, BVOE is nonzero even for linear response sites in the middle of the lattice, but the error is still small enough for practical purposes.

Avoiding BVOE completely to minimize the DMRG-SCF gradient error would require a parametrization that comprises MPS parameters from all sites: either a tangent-space based approachWouters_Thouless_PRB_2013 ; Nakatani_Linear_JCP_2014 or an extension of our current scheme to a sweep-like procedure, where CP-DMRG-SCF equations are solved for every site and the nonvariational contributions to the RDMs are obtained at the end of the sweep. Both procedures are subject of our future work, but we would expect both of them to have a very high computational cost either due to the sweep procedure or employing a very large number of MPS parameters. Yet, the current algorithm based on a single linear response site, as we see here, already achieves a very good accuracy at a rather moderate computational cost by solving the CP-DMRG-SCF equations only once.

Surprisingly, the nonadiabatic couplings show errors up to four orders of magnitude larger than the gradients. However, the error in nonadiabatic couplings can be attributed to their larger sensitivity to the wavefunction quality in general: in a DMRG-SCF calculation with an extremely tight energy convergence threshold for the sweep of $10^{-12}$ Hartree the average deviation of the couplings from their SA-CASSCF counterparts could be reduced to $7\times 10^{-6}$ Hartree/bohr. The sensitivity of nonadiabatic couplings to the quality of the reference wavefunction is not limited to DMRG-SCF but can be observed between several reference CASSCF calculations: a SA-CASSCF calculation with a lower Davidson diagonalization threshold of $10^{-9}$ Hartree yields nonadiabatic couplings which differ on average by $1.3\times 10^{-7}$ Hartree/bohr from the original CASSCF calculation. Note that these deviations in the nonadiabatic couplings have no practical effect on a conical intersection optimization, as we will see below.

III.3 1,2-Dioxetanone: Optimization of a conical intersection

1,2-dioxetanone (Fig. 2a) is a simple chemoluminescent compound, whose thermal dissociation mechanism, along with those of its substituted derivatives have been extensively studied in chemo- and bioluminescence studies.Schmidt_Kinetics_JACS_1978 ; Adam_Cyclic_JACS_1979 ; Turro_Chemiluminescent_JACS_1980 ; Liu_Theoretical_JACS_2009 ; Liu_CASSCFCASPT2_CPL_2009 ; Navizet_Chemistry_C_2011 ; Yue_Theoretical_JPCL_2015 ; Farahani_combined_RA_2017 ; Vacher_Chemi-_CR_2018 ; Yue_Two_JCTC_2019 Liu_Theoretical_JACS_2009 performed an extensive computational study on the dissociation pathway of 1,2-dioxetanone and located two conical intersections between the $\sigma\sigma^{*}$ and the $n\sigma^{*}$ states along the pathway. In the vicinity of these conical intersections both states share a significant biradical character, which mandates a multiconfigurational treatment such as CASSCF. We performed an optimization of the first ($\sigma\sigma^{*}$),($n\sigma^{*}$) conical intersections along the reaction path, named by Liu_Theoretical_JACS_2009 as “${}^{1}\sigma\sigma^{*}$-TS”, both with CASSCF and DMRG-SCF. The optimization run for both methods is presented in Table 1. At each step the energies and the gradients are nearly identical, and the discrepancies are very small – below $2.8\times 10^{-6}$ Hartree for the energies, $6\times 10^{-5}$ Hartree/bohr for gradient norms and $6.3\times 10^{-5}$ Hartree/bohr for the maximum gradient element – in all cases well below the convergence thresholds for the optimization. These small discrepancies again decrease even further by tightening the convergence thresholds. The optimized structures are also essentially identical for both methods: the maximum absolute difference in the Cartesian coordinates is $2\times 10^{-4}$ Å for both optimized structures. The most important bond lengths and angles of the optimized structure are presented in Fig. 2b. As we see, a conical intersection optimization with our SA-DMRG-SCF analytical gradient and nonadiabatic coupling ansatz is able to accurately reproduce the SA-CASSCF optimization result.

III.4 1,2-Dioxetanone: Gradient convergence with number of renormalized block states

Having established that our SA-DMRG-SCF gradient ansatz is able to reproduce SA-CASSCF gradients to an accuracy that is only dependent on convergence thresholds in cases where an $m$ value is sufficiently large to reproduce the SA-CASSCF wavefunction to arbitrary accuracy, we now study the gradient error for smaller $m$ values, i. e., when the approximation of the SA-CASSCF wavefunction by a SA-DMRG-SCF wavefunction is of a reduced quality. For this, a gradient calculation has been performed for a structure in the vicinity of the ($\sigma\sigma^{*}$),($n\sigma^{*}$) conical intersection with SA-DMRG-SCF(16,13), but with several $m$ values varying from 200 to 6000.

Fig. 3 shows the mean and maximum absolute gradient error for each $m$, along with the average energy error. We recognize that for large $m$ values ($m\geq 2000$), the gradient error is one to two orders of magnitude larger than the energy error, and that an $m$ value which is sufficient to converge the energy up to an error of $10^{-8}$ Hartree is not sufficient to converge the gradient to the same accuracy as the energy. Even a much larger value ($m=6000$) does not converge the gradient to the same accuracy, although the error decreases by another order of magnitude. We expect that tightening the convergence threshold for the energy along with larger $m$ values would improve the convergence even further. However (as we see both from this result and the conical intersection optimization in the previous section) even at $m=2000$, the gradient error is well below typical convergence thresholds for optimizations and can therefore be neglected for practical applications.

However, for small $m$ values we see an interesting development: although the maximum gradient error is approximately the same or larger than the energy error, the mean gradient error is actually smaller. The crossover of the mean energy and gradient errors occurs around $m=1000$. While we believe that this result is caused by a fortituous error cancellation, it implies that for typical DMRG-SCF calculations, which are performed for active spaces inaccessible by standard CASSCF the accuracy of the gradient error is of the same order of magnitude or smaller than the energy error. In other words, in these calculations we may extract energies and gradients of a similar accuracy from the same wavefunction without further refinements such as imposing tighter convergence thresholds.

III.5 1,2-Dioxetanone: Comparison with numerical gradients

To further assess the errors of the SA-DMRG-SCF analytical gradient for small $m$ values, in particular to quantify the BVOE, we compared the accuracy of the analytical gradient to a numerical gradient obtained for the same $m$ value. For this comparison, again a (16,13) active space and an $m$ value of 500 was chosen. Although it is beyond the scope of this work to provide many examples, we attempted to employ several testcases, and, as such, here we calculated the S1 excited state gradient for the S0 optimized structure (rather than a structure in the vicinity of the conical intersection as in Section III.4). For the numerical gradient, a convergence with the step size has been verified by repeating the numerical gradient calculations with step sizes of 0.0025 and 0.00125 bohr and ensuring that the difference in every gradient component does not exceed $10^{-7}$ Hartree/bohr. Since even a numerical gradient that is fully converged with respect to the step size may deviate from the analytical gradient, a comparison between the CASSCF analytical and numerical gradient was considered to quantify this error. Table 2 lists the mean absolute deviations between numerical and analytical gradients, both for DMRG-SCF and CASSCF. From the comparison of the differences between analytical and numerical gradients we estimate the BVOE to be about $1.9\times{}10^{-3}$ Hartree/bohr (see first entry in Table 2). It is an order of magnitude larger than the gradient error which arises due to the small $m$ value alone (equal to the mean absolute deviation between numerical DMRG-SCF and CASSCF gradients, $2.2\times{}10^{-4}$ Hartree/bohr).

It is noteworthy that the BVOE in this example is approximately one order of magnitude larger than the analogous error for $m=500$ in Fig. 3. This is not surprising as one would expect the error to decrease when the non-variational MPS contribution to the gradient decreases. In the vicinity of a conical intersection $S_{0}$ and $S_{1}$ wavefunctions may mix and thus their electronic character becomes more similar. Under these circumstances one would expect a smaller non-variational contribution to the gradient. However, at the S0 minimum the S0 and S1 wavefunctions are very different, which in turn would lead to a larger non-variational contribution, and, consequently, a larger BVOE. While this shows that our approximation is particularly well suited for the optimization of conical intersections, it is the scope of future work to assess our scheme for other photochemical applications, such as excited state structure optimizations or surface hopping dynamics.

IV Conclusions

In this paper, we presented an implementation of a Lagrangian-based ansatz for the analytical state-average DMRG-SCF gradients and nonadiabatic couplings which requires construction and solution of CP-MCSCF equations. Our ansatz generalizes the SA-CASSCF gradients theory first presented by Jones2001Analytical in 2001 to wavefunctions encoded as MPSs. We derive the Lagrange multipliers in the CP-MCSCF equations from the mixed-canonical representation of the MPS wavefunction at one particular site, here called the linear response site, which can be chosen arbitrarily. The choice of the linear response site can be made a tradeoff between accuracy and computational cost. As it is, such an MPS parametrization constitutes an approximation, because we consider the variation of the MPS at one particular site only. However, as the MPS approaches the full CI wavefunction, the error associated with this approximation approaches zero. We further argue that for the validity of our definition of MPS Lagrange multipliers the MPS for all states must be optimized simultaneously with common left and right boundaries. Finally, we showed that our SA-DMRG-SCF gradient and nonadiabatic coupling ansatz will exactly reproduce SA-CASSCF gradients and nonadiabatic couplings, respectively, if the SA-DMRG-SCF wavefunction reproduces the SA-CASSCF wavefunction exactly. If this is not the case, both gradients and, in particular, nonadiabatic couplings show errors that are larger than the errors in the energies. This is largely due to the approximate nature of our MPS parametrization, however, this error is still small enough for practical purposes, especially for conical intersection optimizations. We have demonstrated the feasibility of our method in practical applications by performing a conical intersection optimization of 1,2-dioxetanone, reproducing the result of a SA-CASSCF conical intersection optimization. Our development, in addition to applications for conical intersection optimization for systems with large active spaces, paves the way for surface hopping studies and other excited state studies with DMRG.

Acknowledgments

L.F. acknowledges the Austrian Science Fund for a Schrödinger fellowship (Project No. J 3935). This work was supported by the Schweizerischer Nationalfonds (SNF project 200021_182400), National Natural Science Foundation of China (NO. 21703260), ETH Zürich (ETH Fellowship No. FEL-49 18-1) and the Informatization Program of the Chinese Academy of Sciences (NO. XXH13506-403). We thank Prof. Roland Lindh and Dr. Christopher Stein for fruitful discussions.