Second-order response derivatives for the variational Lagrangian

In order to calculate second-order molecular properties using a variational Lagrangian, we take the perturbational derivative of the partial derivative of the Lagrangian and invoke the chain rule: \[\begin{equation} \require{physics} \begin{split} & \eval{ \frac{ \dd{^2}{E(\boldsymbol{\eta})} }{\dd{\eta_m} \dd{\eta_n}} }_{\boldsymbol{\eta}_0} \\ & \hspace{12pt} = \eval{ \pdv{ \mathscr{L} ( \boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0), \boldsymbol{\lambda}^\star(\boldsymbol{\eta}_0) ) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} + \sum_i^x \qty( \eval{ \pdv{ \mathscr{L} ( \boldsymbol{\eta}, \vb{p}, \boldsymbol{\lambda}^\star(\boldsymbol{\eta}_0) ) }{p_i}{\eta_n} }_{ \boldsymbol{\eta}_0, \vb{p}^\star(\boldsymbol{\eta}_0) } ) \qty( \eval{ \pdv{ p_i^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \\ & \hspace{24pt} + \sum_a^y \qty( \eval{ \pdv{ \mathscr{L} ( \boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}), \boldsymbol{\lambda} ) }{\lambda_a}{\eta_n} }_{ \boldsymbol{\eta}_0, \boldsymbol{\lambda}^\star(\boldsymbol{\eta}_0) } ) \qty( \eval{ \pdv{ \lambda_a^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \thinspace , \end{split} \label{eq:molecular_Hessians} \end{equation}\] which means that for second-order properties, we again require the first-order response of the wave function parameters: \[\begin{equation} \qty[ \vb{x} ]_{im} = \eval{ \pdv{ p_i^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} \thinspace . \end{equation}\] Furthermore, we also require the first-order response of the Lagrange parameters: \[\begin{equation} \qty[ \vb{y} ]_{am} = \eval{ \pdv{ \lambda_a^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} \thinspace . \end{equation}\]

In order to find the first-order wave function and multiplier responses, we have to solve the linear response equations. Once found, we are able to calculate second-order properties by using: \[\begin{equation} \eval{ \frac{ \dd{^2}{\mathcal{E}(\boldsymbol{\eta})} }{\dd{\eta_m} \dd{\eta_n}} }_{\boldsymbol{\eta}_0} = \eval{ \pdv{ \mathscr{L} ( \boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0), \boldsymbol{\lambda}^\star(\boldsymbol{\eta}_0) ) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} + \qty[ \vb{x}^{\text{T}} \vb{G}_{\boldsymbol{\lambda}} ]_{mn} + \qty[ \vb{y}^{\text{T}} \vb{F}_{\vb{p}} ]_{mn} \thinspace , \end{equation}\] in which the explicit second-order perturbational derivative is given by: \[\begin{equation} \eval{ \pdv{ \mathscr{L} ( \boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0), \boldsymbol{\lambda}^\star(\boldsymbol{\eta}_0) ) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} = \eval{ \pdv{ E( \boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0) ) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} + \sum_a^y \lambda_a^\star(\boldsymbol{\eta}_0) \qty( \eval{ \pdv{ f_a ( \boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0) ) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} ) \thinspace . \end{equation}\]