The Rayleigh energy function
There often exists some confusing regarding the terms variational energy, variational principle and variational parameters, so we will devote the first part of this section to clearing up these confusions.
The Rayleigh quotient: \[\begin{equation} \require{physics} E(\boldsymbol{\eta}, \vb{p}) = \frac{ \ev{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\Psi(\vb{p})} }{ \braket{\Psi(\vb{p})} } \thinspace , \end{equation}\] is very often used as an energy function for a wave function model \(\ket{\Psi(\vb{p})}\). The Ritz variational principle uses the Rayleigh quotient in conjunction with a linear expansion of the wave function model inside a Fock subspace: \[\begin{equation} \ket{\Psi(\vb{p})} = \sum_{\vb{k}} c_{\vb{k}} \ket{\vb{k}} \thinspace , \end{equation}\] where \(\ket{\vb{k}}\) represents a Fock subspace vector, often an ONV. We can show that solving the stationary conditions for the Ritz method amounts to diagonalizing the matrix representation of the Hamiltonian inside the Fock subspace. As a consequence, wave functions expressed in the same Fock subspace always have an energy higher than or equal to the lowest eigenvector of the Hamiltonian. Every wave function that is determined using the Ritz method automatically has variationally determined parameters. On the other hand, wave functions that are variationally determined, need not come from the application of the Ritz method in conjunction with the Rayleigh quotient as an energy function. We will see that there exist certain methods, such as coupled-cluster-, or in general projected Schrödinger equation-, type wave functions whose energy functions are not Rayleigh quotients. Even though we can apply some modifications such that these PSE-type parameters are variationally determined, the use of a non-Rayleigh quotient energy function allows for energy values to fall below the corresponding Rayleigh-Ritz energy value in the same Fock subspace.
Let us now proceed by deriving formulas for the first- and second-order perturbational derivatives of the energy, for methods that apply the Rayleigh quotient as an energy function. Without loss of generality, we can assume that the wave function model \(\ket{\Psi(\vb{p})}\) is normalized. For the first-order perturbational derivative, due to the Hellmann-Feynman theorem, we only require the perturbational partial derivative of the energy. For the Rayleigh energy function, this is: \[\begin{equation} \eval{ \pdv{ E (\boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0)) }{\eta_m} }_{\boldsymbol{\eta}_0} \\ = \ev**{ \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) }{ \Psi(\vb{p}^\star(\boldsymbol{\eta}_0)) } \thinspace , \end{equation}\] such that the main issue in this case is isolated into calculating the partial derivative of the Hamiltonian, which is a substantial simplification. The formulas derived in the section on perturbation-dependent Fock spaces come in handy and can be readily used to calculate first-order properties.
In order to obtain second-order molecular properties, we’ll have to solve the linear response equations. The parameter response force can in this case be calculated as: \[\begin{align} \qty[ \vb{F}_{\vb{p}} ]_{im} &= \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}) }{p_i}{\eta_m} }_{ \boldsymbol{\eta}_0, \vb{p}^\star(\boldsymbol{\eta}_0) } \\ &= \matrixel{ \eval{ \pdv{ \Psi(\vb{p}) }{p_i} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} }{ \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) }{ \Psi(\vb{p}^\star(\boldsymbol{\eta}_0)) } \\ & \hspace{12pt} + \matrixel{ \Psi(\vb{p}^\star(\boldsymbol{\eta}_0)) }{ \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) }{ \eval{ \pdv{ \Psi(\vb{p}) }{p_i} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} } \end{align}\] and the parameter response force constant as: \[\begin{align} \qty[ \vb{k}_{\vb{p}} ]_{ij} &= \eval{ \pdv{ E(\boldsymbol{\eta}_0, \vb{p}) }{p_i}{p_j} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} \\ &= \matrixel{ \eval{ \pdv{ \Psi(\vb{p}) }{p_i}{p_j} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} }{ \hat{\mathcal{H}}(\boldsymbol{\eta}_0) }{ \Psi(\vb{p}^\star(\boldsymbol{\eta}_0)) } + \matrixel{ \Psi(\vb{p}^\star(\boldsymbol{\eta}_0)) }{ \hat{\mathcal{H}}(\boldsymbol{\eta}_0) }{ \eval{ \pdv{ \Psi(\vb{p}) }{p_i}{p_j} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} } \notag \\ & \hspace{12pt} + (1 + P_{ij}) \matrixel{ \eval{ \pdv{ \Psi(\vb{p}) }{p_i} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} }{ \hat{\mathcal{H}}(\boldsymbol{\eta}_0) }{ \eval{ \pdv{ \Psi(\vb{p}) }{p_j} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} } \thinspace . \end{align}\] In order, then, to calculate the actual values for second-order molecular properties, we also need the explicit second-order perturbational derivative of the energy: \[\begin{equation} \eval{ \pdv{ E (\boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0)) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} = \ev**{ \qty( \eval{ \pdv{ \hat{\mathcal{H}}(\boldsymbol{\eta}) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} ) }{ \Psi(\vb{p}^\star(\boldsymbol{\eta}_0)) } \thinspace . \end{equation}\]