{
"cells": [
{
"cell_type": "markdown",
"id": "9d35a667",
"metadata": {
"tags": []
},
"source": [
"# Energy gradients\n",
"(sec:energy-gradients)=\n",
"\n",
"In order to perform a geometry optimization and find a minimum energy structure or a transition state,\n",
"one needs to calculate the derivative of the energy $E$ with respect to the **nuclear coordinates**.\n",
"This is usually referred to as the *molecular gradient* and it needs to vanish at any stationary geometry.\n",
"The second derivatives of the energy with respect to nuclear displacement (the molecular **Hessian**) can be used to characterize\n",
"the stationary structure, i.e., to confirm whether a true minimum or a transition state has been found.\n",
"But not only molecular gradients and Hessians can be calculated as derivatives of the energy,\n",
"also other properties such as permanent and induced (dipole) moments, polarizabilities, and magnetizabilites\n",
"when taking the derivative with respect to external **electromagnetic field**, or NMR and EPR parameters\n",
"when additionally **nuclear magnetic moments** are involved {cite}`jensen2006`.\n",
"\n",
"The energy derivatives themselves can either be calculated **numerically** by using finite differences\n",
"or **analytically**. The former is usually simple to implement, but suffers from difficulties\n",
"in numerical accuracy and computational efficiency.\n",
"For the latter, considerable programming effort is requiered, but it has the advantages of greater speed, precision, and convenience.\n",
"\n",
"(sec:numerical-gradients)=\n",
"## Numerical Gradients\n",
"\n",
"The simplest method to calculate the derivative of the energy $E(\\chi)$ with respect to some parameter $\\chi$ is to use finite difference approximations.\n",
"Choosing a small change $\\chi_0$ in $\\chi$, a two-point estimation is given by\n",
"```{math}\n",
"%:label: eq:divided_difference\n",
"\\frac{\\mathrm{d} E}{\\mathrm{d} \\chi} \\approx \\frac{E(\\chi + \\chi_0) - E(\\chi)}{\\chi_0} \\, ,\n",
"```\n",
"which is known as a first-order **divided difference** and its error to the true derivative is approximately proportional to $\\chi_0$.\n",
"The true derivative of $E$ at $\\chi$ is given by the limit\n",
"```{math}\n",
"%:label: eq:true_derivative\n",
"\\frac{\\mathrm{d} E}{\\mathrm{d} \\chi} = \\lim_{\\chi_0 \\to 0} \\frac{E(\\chi + \\chi_0) - E(\\chi)}{\\chi_0} \\, .\n",
"```\n",
"\n",
"(eq:symmetric-difference-quotient)=\n",
"Another two-point formula is the **symmetric difference quotient** given by\n",
"```{math}\n",
"%:label: eq:symmetric_difference_quotient\n",
"\\frac{\\mathrm{d} E}{\\mathrm{d} \\chi} \\approx \\frac{E(\\chi + \\chi_0) - E(\\chi - \\chi_0)}{2\\chi_0} \\, ,\n",
"```\n",
"where it can be shown that the first-order error cancels, so it is approximately proportional to $\\chi_0^2$,\n",
"which means that for small $\\chi_0$ this is a more accurate approximation to the true derivative\n",
"and it is commonly used in numerical derivative codes. VeloxChem uses a default value of $0.001~a_0$.\n",
"However, in both cases two calculations of the energy $E(\\chi)$ need to be performed in order to obtain\n",
"the derivative with respect to one variable $\\chi$.\n",
"\n",
"There are also higher-order methods for approximating the first derivative, such as the \"five-point method\" given by\n",
"```{math}\n",
"%:label: eq:five_point_method\n",
"\\frac{\\mathrm{d} E}{\\mathrm{d} \\chi} \\approx \\frac{E(\\chi - 2 \\chi_0) - 8 E(\\chi - \\chi_0) + 8 E(\\chi + \\chi_0) - E(\\chi + 2 \\chi_0)}{12\\chi_0} \\, ,\n",
"```\n",
"where the error is approximately proportional to $\\chi_0^4$ and it hence gives a very accurate approximation to the gradient.\n",
"However, since four individual energy calculations need to be performed for each perturbation,\n",
"this expression is usually only employed for debugging of the analytical derivative expressions.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "6ce0f6ce-f9a8-485f-8fd3-67cda0832003",
"metadata": {},
"source": [
"\n",
"(sec:analytical-gradients)=\n",
"## Analytical Gradients\n",
"%To add: Introduction, approaches to derive analytical gradients; Lagrangian approach.\n",
"\n",
"To determine the energy gradient, we first need to identify the non-variational components of the energy functional. In the case of self-consistent field (SCF)\n",
"methods, i.e., HF and Kohn--Sham DFT, these are the molecular-orbital (MO) coefficients, which exhibit an implicit dependence on the nuclear coordinates when atom-centred basis functions such as Gaussian or Slater-type atomic orbitals (AOs) are used. Considering this implicit dependence, the total energy derivative with respect to a particular nuclear coordinate $\\xi$ is obtained via the chain rule {cite}`Rehn2015,Levchenko2005`:\n",
"%\n",
"```{math}\n",
"%:label: eq:energy_functional_general\n",
"\\frac{\\mathrm{d}E}{\\mathrm{d}\\xi}=\\frac{\\partial E}{\\partial \\xi}+\\frac{\\partial E}{\\partial\\mathbf{C}}\\frac{\\mathrm{d}\\mathbf{C}}{\\mathrm{d}\\xi} \\, .\n",
"```\n",
"%\n",
"Here, $\\mathbf{C}$ is the MO coefficient matrix, which transforms from a set of AOs $\\{\\chi_\\mu\\}$ to a set of MOs $\\{\\phi_p\\}$ via \n",
"\n",
"$$\n",
"\\phi_p=\\sum_\\mu C_{\\mu p}\\chi_\\mu \\, .\n",
"$$\n",
"\n",
"The first term, $\\partial E/\\partial\\xi$, is the Hellmann--Feynman contribution which describes the explicit dependence of the energy on the nuclear coordinate $\\xi$ through the nuclear-electron and nuclear-nuclear interaction terms of the Hamiltonian {cite}`Levchenko2005_thesis, Helgaker1988_analytical`. The second term stems from the implicit dependence of the energy on $\\xi$ due to the fact that the molecular orbitals are expanded in a finite atomic-centred basis set {cite}`Helgaker1988_analytical` \n",
"%This term would vanish if a basis set independent on atomic positions, for example a plane wave basis, were used.\n",
"\n",
"It may seem at first surprising that the derivative ${\\partial E}/{\\partial\\mathbf{C}}$ has to be computed. If the MO coefficients are obtained variationally for a specified molecular geometry, how is it that this derivative is not zero? The key to this conundrum lies in the phrase \"for a specified molecular geometry\". Since the DFT energy and density are constructed using a constrained LCAO parametrization, if we perform a nuclear displacement, the \"old\" MO coefficients no longer correspond to the minimum energy and must be re-optimized. Thus the partial derivative with respect to the MO coefficients, as well as the derivative of the MO coefficients with respect to $\\xi$ are required. The explicit computation of the latter is complicated, but can be avoided by making use of a new functional, the Lagrangian, for which the partial derivative $\\partial L / \\partial \\mathbf{C}$ is zero and constrained to the HF/DFT configuration space by construction {cite}`Levchenko2005, Helgaker2014`:\n",
"%ensuring that we re still in the H Thus, the derivative\n",
"%Thus, even though the local points of the PES are variational with respect to the MO coefficients, the \"global\" energy functional is not. We must therefore, include the partial derivative with respect to the MO coefficients, as well as the derivative of the MO coefficients with respect to $\\xi$, which is quite complicated to compute. However, we can avoid explicitly computing $\\mathrm{d}\\mathbf{C}/\\mathrm{d}\\xi$ by using a trick. The idea is to create a new functional, the Lagrangian ($L$), which is by construction equivalent to the energy functional, but for which the partial derivative $\\partial L / \\partial \\mathbf{C}$ is zero \\cite{Levchenko2005, Helgaker2014}:\n",
"%\n",
"```{math}\n",
"%:label: eq:Lagrangian_general\n",
"L(\\mathbf{C},\\boldsymbol{\\Lambda})=E(\\mathbf{C})+\\boldsymbol{\\Lambda}f_c(\\mathbf{C}) \\, ,\n",
"```\n",
"%\n",
"where $\\boldsymbol{\\Lambda}$ are a set of undetermined Lagrange multipliers and $f_c(\\mathbf{C})=0$ define the constraints for the non-variational parameters $\\mathbf{C}$. These constraints ensure that we are moving only in the HF/DFT configuration space, rather than in the infinite space of all orthogonally equivalent combinations of orbital bases {cite}`Helgaker1988_analytical`.\n",
"%From Helgaker and Jorgensen:\"At each geometry we have a set of AOs from which an infinite set of orthogonally equivalent orbital bases can be constructed. As the geometry changes we must pick out exactly one of these orbital bases at each geometry X. In this way an orthogonal orbital connection is established. (A connection is called orthogonal if it preserves orthonormality between the orbitals.) We further require that the connection is continuous and differentiable. One may also wish to impose an additional requirement on the connection, namely that it is translationally and rotationally invariant. This may seem to be a trivial requirement. However, a connection is conveniently defined in terms of atomic Cartesian displacements rather than in terms of a set of nonredundant internal coordinates. This implies that each molecular geometry may be described in an infinite number of translationally and rotationally equivalent ways\"\n",
"\n",
"By using the Lagrangian, we have shifted the difficult problem of computing ${\\mathrm{d} \\mathbf{C}}/{\\mathrm{d}\\xi}$ to the much simpler problem of determining the unknown Lagrange multipliers which satisfy $\\partial L / \\partial \\mathbf{C}=0$. Equations for these are derived by imposing that the explicit form of $\\partial L / \\partial \\mathbf{C}$ is zero {cite}`Helgaker2014`. Once the Lagrange multipliers have been obtained, the total derivative of the energy functional with respect to the nuclear coordinate $\\xi$ can be computed as:\n",
"\\begin{equation*}\n",
"\\frac{\\mathrm{d}E}{\\mathrm{d}\\xi}=\\frac{\\mathrm{d}L}{\\mathrm{d}\\xi}=\\frac{\\partial L}{\\partial\\xi}\\, .\n",
"\\end{equation*}\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b6a40259-80a4-4782-b9aa-6d5032ceac76",
"metadata": {},
"source": [
"\n",
"### Ground state\n",
"\n",
"(sec:hf-gradients)=\n",
"#### HF\n",
"As an example, we will derive in the following the analytical expression for the Hartree--Fock (HF) energy gradient.\n",
"The electronic ground-state energy described at the HF level of theory can be written as\n",
"(eq:HF_energy_fct)=\n",
" ```{math}\n",
"%:label: eq:HF_energy_fct\n",
"E_\\mathrm{HF}=\\sum_{i}f_{ii}-\\frac{1}{2}\\sum_{ij}\\braket{ij||ij}\\, ,\n",
"```\n",
"where $f_{pq}$ are Fock matrix elements, $\\braket{pq||rs}$ are anti-symmetrized two-electron integrals in physicists' (\"1212\") notation {cite}`Szabo2012`. As usual, the indices $i,j,\\ldots$ denote occupied molecular orbitals,\n",
"$a,b,\\ldots$ denote unoccupied (virtual) ones, and $p,q,\\ldots$ stand for both occupied and virtual, while Greek letters $\\mu,\\nu,\\ldots$ are used for AO indices.\n",
"\n",
"The Lagrangian for this energy functional is constructed as follows:\n",
"\n",
"(eq:hf-lagrangian)=\n",
"```{math}\n",
"%:label: eq:HF_Lagrangian\n",
"L(\\mathbf{C},\\boldsymbol{\\Lambda},\\boldsymbol{\\Omega})=E_\\mathrm{HF}+\\sum_{p,q}\\lambda_{pq}\\left(f_{pq}-\\delta_{pq}\\epsilon_p\\right)+\\sum_{p,q}\\omega_{pq}\\left(S_{pq}-\\delta_{pq}\\right) \\, ,\n",
"```\n",
"with $\\boldsymbol{\\Lambda}=\\{\\lambda_{pq}\\}$ and $\\boldsymbol{\\Omega}=\\{\\omega_{pq}\\}$ being the Lagrange multipliers,\n",
"$\\mathbf{F}=\\{f_{pq}\\}$ the Fock matrix, $\\{\\epsilon_p\\}$ the HF orbital energies, and $\\mathbf{S}=\\{S_{pq}\\}$ the overlap matrix.\n",
"Here, the conditions for the Lagrange multipliers $\\{\\lambda_{pq}\\}$ and $\\{\\omega_{pq}\\}$ ensure that the Fock matrix is diagonal and, respectively, the overlap matrix is unity for the HF state.\n",
"\n",
"To calculate the total derivative of the energy with respect to $\\xi$ we must now only calculate the partial derivative of the Lagrangian with respect to the same variable:\n",
"(eq:partial_L)=\n",
"```{math}\n",
"%:label: eq:partial_L\n",
"\\frac{\\partial L}{\\partial \\xi}=\\frac{\\partial E_\\mathrm{HF}}{\\partial \\xi}+\\sum_{p,q}\\lambda_{pq}\\frac{\\partial f_{pq}}{\\partial \\xi}+\\sum_{p,q}\\omega_{pq}\\frac{\\partial S_{pq}}{\\partial \\xi}\\, .\n",
"```\n",
"We want to re-write the above derivative of the Lagrangian in terms of effective density matrices.\n",
"For this purpose, we express the energy in terms of the one- and two-particle density matrices, $\\boldsymbol{\\gamma} = \\{\\gamma_{pq}\\}$\n",
"and $\\boldsymbol{\\Gamma} = \\{\\Gamma_{pqrs}\\}$, respectively:\n",
"(eq:EHF_DM)=\n",
"```{math}\n",
"%:label: eq:EHF_DM\n",
"E_\\mathrm{HF}=\\sum_{p,q}f_{pq}\\gamma_{pq}+\\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}\\braket{pq||rs}\\, .\n",
"```\n",
"\n",
"With this definition, the [above equation](eq:partial_L) becomes:\n",
"(eq:partial_L_final)=\n",
"```{math}\n",
"%:label: eq:partial_L_final\n",
"\\frac{\\partial L}{\\partial \\xi}&=\\sum_{p,q}(\\lambda_{pq}+\\gamma_{pq}) f^{(\\xi)}_{pq}+\\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}\\braket{pq||rs}^{(\\xi)}+\\sum_{p,q}\\omega_{pq}S^{(\\xi)}_{pq}\\nonumber\\\\\n",
"&=\\sum_{p,q}(\\lambda_{pq}+\\gamma_{pq}) h^{(\\xi)}_{pq}+\\sum_{p,q}(\\lambda_{pq}+\\gamma_{pq})\\sum_{i}\\braket{pi||qi}^{(\\xi)}+\\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}\\braket{pq||rs}^{(\\xi)}+\\sum_{p,q}\\omega_{pq}S^{(\\xi)}_{pq}\\,,\n",
"```\n",
"where $h_{pq}$ represents a matrix element of the core-Hamiltonian operator. The superscript $(\\xi)$ indicates a partial derivative with respect to variable $\\xi$,\n",
"i.e., with fixed MO coefficients $\\mathbf{C}$. Explicitly, they are given by {cite}`Levchenko2005`:\n",
"(eq:hpq_Spq)=\n",
"```{math}\n",
"%:label: eq:hpq\n",
"h_{pq}^{(\\xi)}&=\\frac{\\partial h_{pq}}{\\partial \\xi}=\\sum_{\\mu,\\nu}C_{\\mu p}h^{\\xi}_{\\mu\\nu}C_{\\nu q}\\, ,\\\\\n",
"h^{\\xi}_{\\mu\\nu}&=\\bra{\\phi_\\mu}\\frac{\\partial \\hat{h}}{\\partial \\xi}\\ket{\\phi_\\nu}+\\bra{\\frac{\\partial\\phi_\\mu}{\\partial \\xi}}\\hat{h}\\ket{\\phi_\\nu}+\\bra{\\phi_\\mu}\\hat{h}\\ket{\\frac{\\partial\\phi_\\nu}{\\partial\\xi}}\\, ,\\\\\n",
"\\braket{pq||rs}^{(\\xi)}&=\\frac{\\partial \\braket{pq||rs}}{\\partial\\xi}=\\sum_{\\mu,\\nu,\\theta,\\sigma}C_{\\mu p}C_{\\nu q}\\braket{\\phi_{\\mu}\\phi_{\\nu}||\\phi_{\\theta}\\phi_{\\sigma}}^{\\xi}C_{\\theta r}C_{\\sigma s}\\, ,\\\\\n",
"\\braket{\\phi_{\\mu}\\phi_{\\nu}||\\phi_{\\theta}\\phi_{\\sigma}}^{\\xi}&=\\braket{\\frac{\\partial\\phi_{\\mu}}{\\partial\\xi}\\phi_{\\nu}||\\phi_{\\theta}\\phi_{\\sigma}}+\\braket{\\phi_{\\mu}\\frac{\\partial\\phi_{\\nu}}{\\partial\\xi}||\\phi_{\\theta}\\phi_{\\sigma}}+\\braket{\\phi_{\\mu}\\phi_{\\nu}||\\frac{\\partial\\phi_{\\theta}}{\\partial\\xi}\\phi_{\\sigma}}\\nonumber + \\braket{\\phi_{\\mu}\\phi_{\\nu}||\\phi_{\\theta}\\frac{\\partial\\phi_{\\sigma}}{\\partial \\xi}}\\,\\\\\n",
"S_{pq}^{(\\xi)} &=\\frac{\\partial S_{pq}}{\\partial\\xi}=\\sum_{\\mu,\\nu}C_{\\mu p}S_{\\mu\\nu}^\\xi C_{\\nu q}\\, ,\\\\\n",
"S_{\\mu\\nu}^\\xi &=\\braket{\\frac{\\partial\\phi_\\mu}{\\partial\\xi}|\\phi_\\nu}+\\braket{\\phi_\\mu|\\frac{\\partial\\phi_\\nu}{\\partial\\xi}} \\,.\n",
"```\n",
"\n",
"We also made use of the definition of the Fock matrix {cite}`Szabo2012`: \n",
"```{math}\n",
"%:label: eq:FockMatEl\n",
"f_{pq}=h_{pq}+\\sum_i\\braket{pi||qi}\\, .\n",
"```\n",
"[This equation](eq:partial_L_final) is the working equation to compute the partial derivative of the Lagrangian with respect to variable $\\xi$,\n",
"and implicitly the equation for the total derivative of the energy with respect to the same variable.\n",
"Two ingredients are required to calculate $\\partial L/\\partial\\xi$: (1) the derivatives of the core-Hamiltonian matrix,\n",
"of the anti-symmetrized two-electron integrals, and of the overlap matrix,\n",
"and (2) finding the density matrices $\\boldsymbol{\\gamma}$ and $\\boldsymbol{\\Gamma}$ as well as the Lagrange multipliers $\\boldsymbol{\\Lambda}$ and $\\boldsymbol{\\Omega}$.\n",
"\n",
"By comparing [this equation](eq:EHF_DM) to [that one](eq:HF_energy_fct), we can immediately identify the nonvanishing blocks of the density matrices:\n",
"```{math}\n",
"%:label: eq:HF_gamma_ij\n",
"\\gamma_{ij} &= \\delta_{ij} \\, , \\\\\n",
"\\Gamma_{ijkl} &= -2\\delta_{ik}\\delta_{jl}\\, .\n",
"```\n",
"Equations for the $\\{\\lambda_{pq}\\}$ and $\\{\\omega_{pq}\\}$ multipliers are obtained by imposing the Lagrangian to be stationary with respect to the orbital transformation matrix $\\{C_{\\mu t}\\}$:\n",
"```{math}\n",
"%:label: eq:OrbitalRspCond\n",
"\\frac{\\partial L}{\\partial C_{\\mu t}}=0\\, ,\n",
"```\n",
"or the programmable version {cite}`Levchenko2005,Rehn2019` :\n",
"```{math}\n",
"%:label: eq:OrbitalRspCond_program\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial L}{\\partial C_{\\mu t}}=0\\, .\n",
"```\n",
"To calculate the partial derivative of the Lagrangian with respect to $C_{\\mu t}$, we will need the following three expressions:\n",
"```{math}\n",
"%:label: eq:derivFock\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial f_{pq}}{\\partial C_{\\mu t}}&=\\sum_{\\mu}C_{\\mu u}\\left(\\frac{\\partial h_{pq}}{\\partial C_{\\mu t}}+\\sum_i\\frac{\\partial \\braket{pi||qi}}{\\partial C_{\\mu t}}\\right)\\nonumber\\\\\n",
"&=h_{uq}\\delta_{pt}+h_{pu}\\delta_{qt}+\\sum_i\\braket{ui||qi}\\delta_{pt}+\\sum_i\\braket{pi||ui}\\delta_{qt}\\nonumber\\\\\n",
"&\\quad + \\sum_i\\braket{pu||qi}\\delta_{it}+\\sum_i\\braket{pi||qu}\\delta_{it}\\nonumber \\\\\n",
"&=f_{uq}\\delta_{pt}+f_{pu}\\delta_{qt}+\\braket{pu||qt}\\delta_{t\\epsilon_o}+\\braket{pt||qu}\\delta_{t\\epsilon_o} \\, ,\\label{eq:derivFock} \\\\\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial \\braket{pq||rs}}{\\partial C_{\\mu t}} &= \\braket{uq||rs}\\delta_{pt}+\\braket{pu||rs}\\delta_{qt}+\\braket{pq||us}\\delta_{rt}+\\braket{pq||ru}\\delta_{st} \\, , \\\\\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial S_{pq}}{\\partial C_{\\mu t}} &= S_{uq}\\delta_{pt}+S_{pu}\\delta_{qt}\\, , \\label{eq:derivOverlap}\n",
"```\n",
"where we have used the definitions of the Fock matrix, two-electron integrals and overlap matrix in terms of the orbital transformation matrix $\\{C_{\\mu p}\\}$ -- see [here](eq:hpq_Spq). The Kronecker delta $\\delta_{t\\epsilon_\\mathrm{o}}$ is equal to one if $t$ is an occupied orbital and is zero otherwise.\n",
"\n",
"Using the Lagrangian expressed in terms of the density matrices:\n",
"```{math}\n",
"%:label: eq:L_DMs\n",
"L=\\sum_{p,q}\\gamma_{pq}f_{pq}+\\sum_{p,q}\\lambda_{pq}(f_{pq}-\\delta_{pq}\\epsilon_p)+\\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}\\braket{pq||rs}+\\sum_{p,q}\\omega_{pq}(S_{pq}-\\delta_{pq})\\,,\n",
"```\n",
"\n",
"the partial derivative of the Lagrangian with respect to $\\mathbf{C}$ can be written as:\n",
"```{math}\n",
"%:label: eq:OrbRsp1\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial L}{\\partial C_{\\mu t}}&=\\sum_{p,q}\\left(\\gamma_{pq}+\\lambda_{pq}\\right)\\left[f_{uq}\\delta_{pt}+f_{pu}\\delta_{qt}+\\braket{pu||qt}\\delta_{t\\epsilon_o}+\\braket{pt||qu}\\delta_{t\\epsilon_o}\\right]\\nonumber \\\\\n",
"&+\\frac{1}{4}\\sum_{p,q,r,s}\\Gamma_{pqrs}\\left[\\braket{uq||rs}\\delta_{pt}+\\braket{pu||rs}\\delta_{qt}+\\braket{pq||us}\\delta_{rt}+\\braket{pq||ru}\\delta_{st}\\right]\\nonumber\\\\\n",
"&+\\sum_{p,q}\\omega_{pq}\\left(S_{uq}\\delta_{pt}+S_{pu}\\delta_{qt}\\right)\\, .\n",
"```\n",
"By using the conditions $f_{pq}=\\epsilon_p\\delta_{pq}$ and $S_{pq}=\\delta_{pq}$, the above equation becomes:\n",
"```{math}\n",
"%:label: eq:OrbRsp2\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial L}{\\partial C_{\\mu t}}&=\\sum_{p,q}\\left(\\gamma_{pq}+\\lambda_{pq}\\right)\\left[\\epsilon_u\\delta_{uq}\\delta_{pt}+\\epsilon_u\\delta_{pu}\\delta_{qt}+\\braket{pu||qt}\\delta_{t\\epsilon_o}+\\braket{pt||qu}\\delta_{t\\epsilon_o}\\right]\\nonumber \\\\\n",
"&+\\frac{1}{4}\\sum_{p,q,r,s}\\Gamma_{pqrs}\\left[\\braket{uq||rs}\\delta_{pt}+\\braket{pu||rs}\\delta_{qt}+\\braket{pq||us}\\delta_{rt}+\\braket{pq||ru}\\delta_{st}\\right]\\nonumber\\\\\n",
"&+\\sum_{p,q}\\omega_{pq}\\left(\\delta_{uq}\\delta_{pt}+\\delta_{pu}\\delta_{qt}\\right) \\nonumber \\\\\n",
"&=2\\left(\\gamma_{ut}+\\lambda_{ut}\\right)\\epsilon_u+2\\sum_{p,q}\\left(\\gamma_{pq}+\\lambda_{pq}\\right)\\braket{pu||qt}\\delta_{t\\epsilon_o}+\\sum_{p,q,r}\\Gamma_{tpqr}\\braket{up||qr}+2\\omega_{ut}\\, ,\n",
"```\n",
"where we have used $\\gamma_{pq}=\\gamma_{qp}$, $\\braket{pu||qt}=\\braket{qt||pu}$, $\\Gamma_{pqrs}=\\Gamma_{qpsr}=\\Gamma_{srpq}$ (real orbitals), and we have imposed that $\\lambda_{pq}=\\lambda_{qp}$, and $\\omega_{pq}=\\omega_{qp}$ (symmetric representation). Some of the indices have been renamed.\n",
"\n",
"To obtain equations for the orbital response Lagrange multipliers, we first have to decouple $\\boldsymbol{\\Lambda}$ from $\\boldsymbol{\\Omega}$ by taking the difference\n",
"```{math}\n",
"%:label: eq:decouple\n",
"\\sum_{\\mu}C_{\\mu u}\\frac{\\partial L}{\\partial C_{\\mu t}}-\\sum_{\\mu}C_{\\mu t}\\frac{\\partial L}{\\partial C_{\\mu u}}&=2(\\gamma_{ut}+\\lambda_{ut})(\\epsilon_u-\\epsilon_t)\\nonumber\\\\&+2\\sum_{p,q}(\\gamma_{pq}+\\lambda_{pq})\\left(\\braket{pu||qt}\\delta_{t\\epsilon_o}-\\braket{pt||qu}\\delta_{u\\epsilon_o}\\right)\\nonumber\\\\\n",
"&+\\,\\,\\,\\,\\sum_{p,q,r}\\left(\\Gamma_{tpqr}\\braket{up||qr}-\\Gamma_{upqr}\\braket{tp||qr}\\right).\n",
"```\n",
"The system of equations for $\\boldsymbol{\\Lambda}$ is then obtained by choosing $u$ and $t$ from different orbital spaces in the following equation,\n",
"(eq:OrbRspEq)=\n",
"```{math}\n",
"%:label: eq:OrbRspEq\n",
"&2(\\gamma_{ut}+\\lambda_{ut})(\\epsilon_u-\\epsilon_t)+2\\sum_{p,q}(\\gamma_{pq}+\\lambda_{pq})\\left(\\braket{pu||qt}\\delta_{t\\epsilon_o}-\\braket{pt||qu}\\delta_{u\\epsilon_o}\\right)\\nonumber\\\\\n",
"&+\\sum_{p,q,r}\\left(\\Gamma_{tpqr}\\braket{up||qr}-\\Gamma_{upqr}\\braket{tp||qr}\\right)=0\\,.\n",
"```\n",
"Once $\\boldsymbol{\\Lambda}$ is determined, $\\boldsymbol{\\Omega}$ can be calculated in a similar way, using the following equation\n",
"(eq:omega)=\n",
"```{math}\n",
"%:label: eq:omega\n",
"2\\left(\\gamma_{ut}+\\lambda_{ut}\\right)\\epsilon_u+2\\sum_{p,q}\\left(\\gamma_{pq}+\\lambda_{pq}\\right)\\braket{pu||qt}\\delta_{t\\epsilon_o}+\\sum_{p,q,r}\\Gamma_{tpqr}\\braket{up||qr}+2\\omega_{ut}=0\n",
"```\n",
"If we explicitly write the equations for different blocks of $\\boldsymbol{\\Lambda}$, we find that they are all zero. This simplifies the equation for $\\boldsymbol{\\Omega}$ to:\n",
"```{math}\n",
"%:label: eq:omega_HF\n",
"2\\gamma_{ut}\\epsilon_u+2\\sum_{p,q}\\gamma_{pq}\\braket{pu||qt}\\delta_{t\\epsilon_o}+\\sum_{p,q,r}\\Gamma_{tpqr}\\braket{up||qr}+2\\omega_{ut}=0\n",
"```\n",
"The only non-zero block of $\\boldsymbol{\\Omega}$ is the occupied-occupied block:\n",
"```{math}\n",
"%:label: eq:omega_hf_oo\n",
"&u=i,\\, t=j \\nonumber\\\\\n",
"&2\\gamma_{ij}\\,\\epsilon_i+2\\sum_{p,q}\\gamma_{pq}\\braket{pi||qj}\\delta_{j\\epsilon_o}+\\sum_{p,q,r}\\Gamma_{jpqr}\\braket{ip||qr}+2\\omega_{ij} = 0 \\nonumber \\\\\n",
"&\\Leftrightarrow\n",
"2\\delta_{ij}\\epsilon_i+2\\sum_{k,l}\\delta_{kl}\\braket{ki||lj}+\\sum_{klm}(-2\\delta_{jl}\\delta_{km})\\braket{ik||lm}+2\\omega_{ij} =0\\nonumber \\\\\n",
"&\\Leftrightarrow\n",
"\\omega_{ij}=-\\epsilon_i\\delta_{ij},\n",
"```\n",
"\n",
"(eq:HF_grad_final)=\n",
"Using the expressions for the density matrices and the non-zero Lagrange multipliers, we get the following expression for the electronic HF energy derivative:\n",
"```{math}\n",
"%:label: eq:HF_grad_final\n",
"\\frac{\\mathrm{d}E_{\\text{HF}}}{\\mathrm{d}\\xi}&=\\sum_{i,j}\\delta_{ij} h^{(\\xi)}_{ij}+\\sum_{i,j}\\delta_{ij}\\sum_{k}\\braket{ik||jk}^{(\\xi)} - \\frac{1}{2}\\sum_{i,j,k,l}\\delta_{ik}\\delta_{jl}\\braket{ij||kl}^{(\\xi)}-\\sum_{i,j}\\epsilon_{i}\\delta_{ij}S^{(\\xi)}_{ij}\\nonumber \\\\\n",
"&=\\sum_{i}h^{(\\xi)}_{ii}+\\frac{1}{2}\\sum_{i,j}\\braket{ij||ij}^{(\\xi)}-\\sum_{i}\\epsilon_{i}S^{(\\xi)}_{ii}\n",
"```\n",
"Finally, the derivative of the total energy is obtained by adding the trivial contribution from the nuclear repulsion energy term {cite}`Szabo2012`,\n",
"$\\frac{\\mathrm{d} V_{nn}}{\\mathrm{d} \\xi}$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "a3d5fc57-e15d-41c5-bb39-20436d39dbdf",
"metadata": {
"jp-MarkdownHeadingCollapsed": true,
"tags": []
},
"source": [
"(sec:dft-gradients)=\n",
"#### DFT\n",
"\n",
"The [DFT](kohn-sham-equation) gradient can be derived in a similar way, (partially) replacing the exact exchange integrals with the corresponding exchange-correlation (xc) functional contributions.\n",
"Here, we only give equations for the simplest case of the [local density approximation](lda) (LDA),\n",
"but the approach is easily generalizable to other types of xc functionals.\n",
"Instead of using the Kohn--Sham (KS) matrix in the energy expression, we employ the core Hamiltonian,\n",
"\n",
"$$\n",
"E_{\\text{DFT}} = \\sum_{i} h_{ii} + \\sum_{i This has been written for ADC(1) -- has to be re-written for TDA. (Or should we call it CIS?)\n",
"\n",
"To illustrate this procedure, we will take the configuration interaction singles (CIS) method {cite}`Foresman1992`\n",
"as an example, which is equivalent to linear-response time-dependent Hartree--Fock (TDHF) theory within the Tamm--Dancoff approximation {cite}`Dreuw2005`.\n",
"Note that for excitation energies and excited-state properties, CIS also yields the same results as the ADC(1) scheme {cite}`Dreuw2015`.\n",
"The approach is then easily generalizable to TDHF and TDDFT.\n",
"%Furthermore, excited-state gradients of linear-response time-dependent density functional theory (TDDFT) within the TDA %{cite}`hirata99`\n",
"%can be derived in a completely analogous way, only the additional exchange-correlation terms need to be taken into account {cite}`Furche2002`.\n",
"%%%To illustrate this procedure, we will refer to ADC(1) (add link). Note that ADC(1) is equivalent to time-dependent Hartree-Fock (TDHF) in the Tamm--Dancoff approximation (TDA) (add link) and configuration interaction singles (CIS) (add link).\n",
"\n",
"In the CIS scheme, a Hermitian eigenvalue equation of the following form is solved,\n",
"```{math}\n",
"%:label: eq:cis_eigenvalue_eq\n",
" \\mathbf{AX}_n = \\omega_n \\mathbf{X}_n \\, ,\n",
"```\n",
"where $\\omega_n$ is the excitation energy for excited state $n$ with corresponding eigenvector $\\mathbf{X}_n$\n",
"(normalized according to $\\mathbf{X}_n^\\dagger \\mathbf{X}_n = 1$),\n",
"and $\\mathbf{A}$ is the CIS matrix given by the elements\n",
"\n",
"(eq:cis-matrix-elemts)=\n",
"$$\n",
"A_{ia,jb} = (\\epsilon_i - \\epsilon_a) \\delta_{ij} \\delta_{ab} - \\braket{ja||ib} \\, .\n",
"$$\n",
"\n",
"(The CIS matrix elements correspond to the sum of [the zeroth- and first-order terms](adcmat_phph) of the ADC(1) matrix elements.)\n",
"Besides the density matrices required for the HF reference state derived [above](sec:hf-gradients),\n",
"we need to identify additional one- and two-particle density matrices for the excitation energy.\n",
"We therefore formally represent the excitation energy $\\omega_n = \\mathbf{X}_n^\\dagger\\mathbf{A}\\mathbf{X}_n$\n",
"in terms of one- and two-particle density matrices:\n",
"```{math}\n",
"%:label: eq:excitation_energy_adc1\n",
"\\mathbf{X}_n^\\dagger\\mathbf{A}\\mathbf{X}_n = \\sum_{p,q}\\gamma^{(n)}_{pq}f_{pq} + \\frac{1}{4}\\sum_{p,q,r,s}\\Gamma^{(n)}_{pqrs}\\braket{pq||rs}\\, ,\n",
"```\n",
"%where $\\mathbf{M}$ is the TDA or CIS matrix, which can be decomposed into a zeroth- and first-order part, $\\mathbf{M}=\\mathbf{M}^{(0)}+\\mathbf{M}^{(1)}$,\n",
"%where $\\mathbf{M}^{(0)}$ contains Fock-matrix elements and $\\mathbf{M}^{(1)}$ the two-electron integrals,\n",
"%$\\mathbf{X}_n$ is the excitation vector corresponding to excited state $\\ket{n}$, and\n",
"%, which are correct through first order in perturbation theory.\n",
"where the superscript $(n)$ indicates that we are referring to the difference density matrices for the $n$-th excited state.\n",
"To identify these density matrices, we carry out explicitly the matrix-vector multiplication on the left hand side,\n",
"```{math}\n",
"%:label: eq:adc1_identify_dms\n",
"\\mathbf{X}_n^\\dagger\\mathbf{A}\\mathbf{X}_n &= \\sum_{i,a,j,b} x_{jb} \\left[(\\epsilon_a - \\epsilon_i)\\delta_{ab}\\delta_{ij}-\\braket{ja||ib}\\right] x_{ia}\\nonumber\\\\\n",
"&=\\sum_{i,a,j,b} \\left(x_{jb}x_{ia}f_{ab}\\delta_{ij} - x_{jb}x_{ia}f_{ij}\\delta_{ab}-x_{jb}x_{ia}\\braket{ja||ib}\\right)\\,, \n",
"```\n",
"where we have used the explicit form of the CIS matrix elements,\n",
"and $x_{ia}$ are the elements of a specific eigenvector $\\mathbf{X}_n$.\n",
"%and used the fact that the Fock matrix of the underlying reference state is diagonal\n",
"%with orbital eigenvalues on the diagonal, $f_{pq}=\\epsilon_{p}\\delta_{pq}$.\n",
"\n",
"(eq:cis-densities)=\n",
"From the above equation, we identify the excited-state density matrix contributions:\n",
"```{math}\n",
"%:label: eq:adc1_gamma_oo\n",
"\\gamma^{(n)}_{ij} = - \\sum_{a}x_{ja}x_{ia}\\,,\\label{eq:adc1_gamma_oo}\n",
"```\n",
"```{math}\n",
"%:label: eq:adc1_gamma_vv\n",
"\\gamma^{(n)}_{ab} = \\sum_{i} x_{ib}x_{ia}\\,,\\label{eq:adc1_gamma_vv}\n",
"```\n",
"```{math}\n",
"%:label: eq:adc1_Gamma_ovov\n",
"\\Gamma^{(n)}_{iajb} = -x_{ib}x_{ja}\\, ,\\label{eq:adc1_Gamma_ovov}\n",
"```\n",
"where the indices of the vectors in the two-particle density matrix have been renamed.\n",
"\n",
"To obtain the molecular gradient of the excited state ${n}$, the density matrices of the ground state must also be included.\n",
"For CIS, this is the [HF reference state](sec:hf-gradients), so the density matrices for the excited state are:\n",
"(eq:cis_density_matrices)=\n",
"\n",
"$$\n",
"\\gamma'_{ij} &= \\gamma_{ij}+ \\gamma^{(n)}_{ij}= \\delta_{ij} - \\sum_{a}x_{ja}x_{ia}\\,, \\\\\n",
"\\gamma'_{ab} &= \\gamma^{(n)}_{ab} =\\sum_{i} x_{ib}x_{ia}\\,, \\\\\n",
"\\Gamma'_{ijkl} &= \\Gamma_{ijkl} = -2\\delta_{ik}\\delta_{jl}\\, \\\\\n",
"\\Gamma'_{iajb} &= \\Gamma^{(n)}_{iajb} = -x_{ib}x_{ja}\\, .\n",
"\n",
"$$\n",
"\n",
"Since, when written in terms of one- and two-particle density matrices,\n",
"the excited state Lagrangian is virtually identical to the Lagrangians written for [HF](sec:hf-gradients) and [MP2](sec:mp2-gradients),\n",
"we leave the exercise of plugging in the expressions of the DMs to the reader.\n",
"The orbital response equations are also straightforwad to derive by using the above [density matrices](eq:cis_density_matrices)\n",
"in the general [orbital response](eq:OrbRspEq) [equations](eq:omega).\n",
"We leave the step-by-step derivation to the reader, with a note that care should be given\n",
"to the symmetry of the two-particle density matrix $\\Gamma_{pqrs}$.\n",
"Here, we provide the final expressions for $\\boldsymbol{\\Lambda}$ (to be determined iteratively):\n",
"```{math}\n",
"%:label: eq:lambda_adc1_ov\n",
"(\\epsilon_i-\\epsilon_a)\\lambda_{ia} - \\sum_{j,b}\\lambda_{jb}(\\braket{ji||ba} - \\braket{ja||ib}) &= \\sum_{j,k}\\gamma^{(n)}_{jk}\\braket{ji||ka} - \\sum_{b,c}\\gamma^{(n)}_{bc}\\braket{ib||ca} \\nonumber\\\\\n",
"&+ \\sum_{j,k,b}\\Gamma^{(n)}_{jakb}\\braket{ij||kb} + \\sum_{b,k,c}\\Gamma_{ibkc}^{(n)}\\braket{kc||ab} \\label{eq:lambda_adc1_ov}\\,,\n",
"```\n",
"and $\\boldsymbol{\\Omega}$:\n",
"\n",
"$$\n",
"\\omega_{ij} &= -(\\delta_{ij} +\\gamma^{(n)}_{ij})\\epsilon_{i}-\\sum_{k,l}\\gamma^{(n)}_{kl}\\braket{ki||lj}+\\sum_{k,a}\\lambda_{ka}\\left(\\braket{ki||ja}+\\braket{kj||ia}\\right) \\nonumber \\\\\n",
"&\\quad-\\sum_{a,b}\\gamma^{(n)}_{ab}\\braket{ia||jb} - \\sum_{a,k,b} \\Gamma^{(n)}_{jakb}\\braket{ia||kb}\\label{eq:omega_adc1_oo} \\\\\n",
"\\omega_{ia} &= - \\lambda_{ia}\\epsilon_{i} + \\sum_{k,b,j}\\Gamma^{(n)}_{kajb}\\braket{ik||jb}\\label{eq:omega_adc1_ov} \\\\\n",
"\\omega_{ab} &=-\\gamma^{(n)}_{ab}\\epsilon_a - \\sum_{k,c,j}\\Gamma^{(n)}_{kajc}\\braket{kb||jc}\\ \\label{eq:omega_adc1_vv} \\,.\n",
"$$\n",
"\n",
"Using the density matrices and Lagrange multipliers, the analytical CIS gradient can now be determined from the partial derivative of the Lagrangian with respect to $\\xi$.\n"
]
},
{
"cell_type": "markdown",
"id": "fbb17b75-1ae6-4208-a1ab-edc5b928734d",
"metadata": {},
"source": [
"(sec:tdhf-gradients)=\n",
"#### TDHF\n",
"\n",
"In linear-response time-dependent Hartree--Fock (TDHF), also known as the _random phase approximation_ (RPA) {cite}`Dreuw2005`, one solves a pseudo-eigenvalue equation of the form\n",
"\n",
"$$\n",
"\\begin{pmatrix} \\mathbf{A} & \\mathbf{B} \\\\ \\mathbf{B}^* & \\mathbf{A}^* \\end{pmatrix} \\begin{pmatrix} \\mathbf{X}_n \\\\ \\mathbf{Y}_n \\end{pmatrix} = \\omega_n \\begin{pmatrix} \\mathbf{1} & \\mathbf{0} \\\\ \\mathbf{0} & -\\mathbf{1} \\end{pmatrix} \\begin{pmatrix} \\mathbf{X}_n \\\\ \\mathbf{Y}_n \\end{pmatrix} \\, ,\n",
"$$\n",
"\n",
"where $\\mathbf{X}_n$ and $\\mathbf{Y}_n$ are referred to as the _excitation_ and _de-excitation_ parts of the response vectors, respectively, the matrix $\\mathbf{A}$ corresponds to the [one](eq:cis-matrix-elemts) from CIS,\n",
"and the $\\mathbf{B}$ matrix is given by the elements $B_{ia,jb} = -\\braket{ab||ij}$.\n",
"The vectors are usually normalized according to $\\mathbf{X}_n^\\dagger \\mathbf{X}_n - \\mathbf{Y}_n^\\dagger \\mathbf{Y}_n = 1$.\n",
"\n",
"The following procedure is analogous to the CIS case and only a few changes are required in the definition of the one- and two-particle density matrices. The non-vanishing contributions to $\\boldsymbol{\\gamma}^{(n)}$ and $\\boldsymbol{\\Gamma}^{(n)}$ for TDHF are given in the following in terms the elements $x_{ia}$ and $y_{ia}$ of $\\mathbf{X}_n$ and $\\mathbf{Y}_n$,\n",
"respectively, or rather their linear combinations $\\mathbf{X}_n \\pm \\mathbf{Y}_n$:\n",
"\n",
"$$\n",
" \\gamma_{ij}^{(n)} &= -\\frac12 \\sum_{a} \\Big[ (x_{ia} + y_{ia})(x_{ja} + y_{ja})\n",
" + (x_{ia} - y_{ia})(x_{ja} - y_{ja}) \\Big] \\\\\n",
" \\gamma_{ab}^{(n)} &= +\\frac12 \\sum_{i} \\Big[ (x_{ia} + y_{ia})(x_{ib} + y_{ib})\n",
" + (x_{ia} - y_{ia})(x_{ib} - y_{ib}) \\Big] \\\\\n",
" \\Gamma_{iajb}^{(n)} &= -\\frac{1}{2} \\big[ (x_{ib}+y_{ib})(x_{ja}+y_{ja})+(x_{ib}-y_{ib})(x_{ja}-y_{ja}) \\big]\\\\\n",
" \\Gamma_{ijab}^{(n)} &= - \\big[ (x_{jb}+y_{jb})(x_{ia}+y_{ia})-(x_{jb}-y_{jb})(x_{ia}-y_{ia}) \\big] \\, .\n",
"$$\n",
"\n",
"Note that there is an additional non-vanishing block in the two-particle density matrix as compared to CIS,\n",
"namely $\\Gamma_{ijab}^{(n)}$, that enters both the right-hand side of the $\\boldsymbol{\\Lambda}$ orbital response equations\n",
"and the $\\boldsymbol{\\Omega}$ multipliers."
]
},
{
"cell_type": "markdown",
"id": "435d1f97-5eaa-48c9-b5f6-129221d05668",
"metadata": {},
"source": [
"(sec:tddft-gradients)=\n",
"#### TDDFT\n",
"\n",
"Analogous to the [ground state](sec:dft-gradients), analytical gradients in the case of\n",
"linear-response [time-dependent density functional theory](sec:tddft) (TDDFT) are virtually identical to the TDHF ones {cite}`Furche2002`.\n",
"Only the exchange-correlation terms need to be considered additionally,\n",
"meaning the matrix elements of $\\mathbf{A}$ and $\\mathbf{B}$ need to be [modified accordingly](sec:tddft-rpa).\n",
"\n",
"The excitation energy $\\omega_n$ for a general hybrid functional can be written as\n",
"\n",
"$$\n",
"\\omega_n = \\sum_{pq}\\gamma_{pq}^{(n)} \\, F_{pq} + \\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}^{(n)} \\, \\big( \\langle pq | rs \\rangle - c_{\\text{x}} \\langle pq | sr \\rangle + f^{\\text{xc}}_{pqrs} \\big) \\, ,\n",
"$$\n",
"\n",
"where the KS matrix $\\mathbf{F}$ is given by\n",
"\n",
"$$\n",
"F_{pq} = h_{pq} + \\sum_{i} \\big( \\langle pi | qi \\rangle - c_{\\text{x}} \\langle pi | iq \\rangle \\big) + v_{pq}^{\\text{xc}} \\, ,\n",
"$$\n",
"\n",
"and the xc contributions $v_{pq}^{\\text{xc}}$ and $f^{\\text{xc}}_{pqrs}$ introduced above, sometimes referred to as the xc potential and xc kernel, respectively {cite}`Dreuw2005`, are given in the LDA\n",
"in a real MO basis as\n",
"(eq:xc-potential-kernel)=\n",
"\n",
"$$\n",
"v_{pq}^{\\text{xc}} &= \\int \\frac{\\partial e_{\\text{xc}}}{\\partial \\rho} \\phi_{p} \\phi_{q} \\mathrm{d} \\mathbf{r} \\\\\n",
"f_{pqrs}^{\\text{xc}} &= \\int \\frac{\\partial^2 e_{\\text{xc}}}{\\partial \\rho \\partial \\rho'} \\phi_p \\phi_q \\phi_r \\phi_s \\mathrm{d}\\mathbf{r} \\, .\n",
"$$\n",
"\n",
"For the orbital response, derivatives of these two quantities with respect to the MO coefficients $\\mathbf{C}$ are required.\n",
"Derivatives of the orbitals $\\phi_p$ behave exactly like for the normal integrals given above,\n",
"however, additional terms occur since $e_{\\text{xc}}$ depends on $\\mathbf{C}$ through the density.\n",
"The derivative of $v_{pq}^{\\text{xc}}$ thus gives a term analogous to $f_{pqrs}^{\\text{xc}}$,\n",
"while the derivative of the latter gives a term with a third-order functional derivative {cite}`Furche2002`,\n",
"\n",
"$$\n",
"g_{pqrstu}^{\\text{xc}} = \\int \\frac{\\partial^3 e_{\\text{xc}}}{\\partial \\rho \\partial \\rho' \\partial \\rho''} \\phi_p \\phi_q \\phi_r \\phi_s \\phi_t \\phi_u \\mathrm{d} \\mathbf{r} \\, .\n",
"$$\n",
"\n",
"For the nuclear gradient of the TDDFT excitation energy we thus need the partial derivatives of those two terms.\n",
"(eq:vxc-deriv)=\n",
"For the xc potential, this is given by\n",
"\n",
"$$\n",
"v_{pq}^{\\text{xc}(\\xi)} = \\frac{\\partial}{\\partial \\xi} \\int \\frac{\\partial e_{\\text{xc}}}{\\partial \\rho} \\phi_{p} \\phi_{q} \\mathrm{d} \\mathbf{r}\n",
"= \\int \\frac{\\partial e_{\\text{xc}}}{\\partial \\rho} \\frac{\\partial (\\phi_{p} \\phi_{q})}{\\partial \\xi} \\mathrm{d} \\mathbf{r}\n",
"+ \\int \\frac{\\partial^2 e_{\\text{xc}}}{\\partial \\rho \\partial \\rho'} \\rho^{(\\xi)} \\phi_p \\phi_q \\mathrm{d} \\mathbf{r} \\, ,\n",
"$$\n",
"\n",
"where the first term is completely analogous to the ground-state contribution,\n",
"except that it gets contracted with the relaxed one-particle density matrix $\\boldsymbol{\\gamma} + \\boldsymbol{\\Lambda}$ (instead of the ground-state density matrix $\\mathbf{P}$),\n",
"and the second term corresponds to an $f^{\\text{xc}}$-like term with the partial derivative of the [ground-state density](eq:density-partial-deriv).\n",
"The partial derivative of the xc kernel with respect to $\\xi$ gives\n",
"\n",
"$$\n",
"f_{pqrs}^{\\text{xc}(\\xi)} = \\frac{\\partial}{\\partial \\xi} \\int \\frac{\\partial^2 e_{\\text{xc}}}{\\partial \\rho \\partial \\rho'} \\phi_p \\phi_q \\phi_r \\phi_s \\mathrm{d}\\mathbf{r}\n",
"= \\int \\frac{\\partial^2 e_{\\text{xc}}}{\\partial \\rho \\partial \\rho'} \\frac{\\partial(\\phi_p \\phi_q \\phi_r \\phi_s)}{\\partial \\xi} \\mathrm{d}\\mathbf{r}\n",
"+ \\int \\frac{\\partial^3 e_{\\text{xc}}}{\\partial \\rho \\partial \\rho' \\partial \\rho''} \\rho^{(\\xi)} \\phi_p \\phi_q \\phi_r \\phi_s \\mathrm{d} \\mathbf{r} \\, ,\n",
"$$\n",
"\n",
"where the first term is again $f^\\text{xc}$-like with an orbital derivative,\n",
"and the second term is analogous to $g^\\text{xc}$ from the orbital response contributions,\n",
"including again the partial derivative of the ground-state density.\n",
"This term eventually gets contracted with the two-particle density matrix $\\boldsymbol{\\Gamma}$,\n",
"so with two excitation or response vectors.\n",
"\n",
"TDDFT within the Tamm--Dancoff approximation (TDA) {cite}`hirata99` is obtained by setting \n",
"$\\mathbf{B} = \\mathbf{0}$ (or $\\mathbf{Y}_n = \\mathbf{0}$).\n",
"All the above considerations are equally valid for TDDFT/TDA,\n",
"with the one- and two-particle density matrices simplifying to the ones from [CIS](eq:cis-densities).\n",
"Note that TDHF (and thus CIS, which is TDHF within the TDA) can be considered a special case\n",
"of general hybrid TDDFT with $c_{\\text{x}} = 1$ and $e_{\\text{xc}} = 0$ {cite}`Furche2002`."
]
},
{
"cell_type": "markdown",
"id": "beaf9fbf-c9fd-4150-a2ff-a9c0286281cd",
"metadata": {},
"source": [
"(sec:first-order-prop)=\n",
"### First-order properties\n",
"\n",
"Not only nuclear gradients, but many other time-independent (or \"static\") molecular properties can be calculated as derivatives of the energy {cite}`jensen2006`.\n",
"For instance, the electric dipole moment can be calculated as the derivative of the energy with respect to an external electric field,\n",
"and the magnetic dipole moment with respect to an external magnetic field.\n",
"Second derivatives with respect to the external field give the electric polarizability and magnetizability, respectively.\n",
"Properties that can be calculated as first derivatives of the energy are referred to as __first-order properties__.\n",
"\n",
"For exact wave functions $\\ket{\\Psi}$ and energies $E = \\langle \\Psi | \\hat{H} | \\Psi \\rangle$, the **Hellmann--Feynman theorem** holds {cite}`jensen2006`,\n",
"\n",
"(eq:hellmann_feynman)=\n",
"```{math}\n",
" \\frac{\\mathrm{d} E}{\\mathrm{d} \\chi} = \\frac{\\mathrm{d}}{\\mathrm{d} \\chi} \\langle \\Psi | \\hat{H} | \\Psi \\rangle\n",
" = \\langle \\Psi | \\frac{\\mathrm{d} \\hat{H}}{\\mathrm{d} \\chi} | \\Psi \\rangle \\, ,\n",
"```\n",
"stating that derivative of the energy with respect to an external perturbation $\\chi$ is identical\n",
"to the expectation value of the perturbed Hamiltonian with the unperturbed wave function.\n",
"The above derivative is to be taken at zero perturbation strength, $\\chi = 0$.\n",
"If the basis functions do not depend on the perturbation, above equation\n",
"also holds for fully variational methods like the SCF and MCSCF schemes.\n",
"\n",
"Consider the electric dipole moment $\\boldsymbol{\\mu}$ as an example.\n",
"Numerically, each component of $\\boldsymbol{\\mu}$ can be obtained by calculating the energy $E$\n",
"in presence of a static electric field $\\mathcal{\\boldsymbol{F}}$, once with a positive and once with a negative sign in one of its components,\n",
"and then applying the formula of the [symmetric difference quotient](eq:symmetric-difference-quotient).\n",
"Starting from a Lagrangian $L$ of [this form](eq:MP2_Lagrangian),\n",
"dipole moments are obtained analytically from a perturbed Hamiltonian, $\\hat{H}_{\\mathcal{\\boldsymbol{F}}} = \\hat{H} + \\mathcal{\\boldsymbol{F}} \\hat{\\mu}$,\n",
"where $\\hat{\\mu}$ is the dipole operator,\n",
"in the Lagrange expression followed by differentiation,\n",
"```{math}\n",
"%:label: eq:dipmom_derivative\n",
" \\boldsymbol{\\mu} = \\frac{\\partial L}{\\partial \\mathcal{\\boldsymbol{F}}} = \\frac{\\partial E}{\\partial \\mathcal{\\boldsymbol{F}}}\n",
" + \\boldsymbol{\\Lambda} \\frac{\\partial f_c(\\mathbf{C})}{\\partial \\mathcal{\\boldsymbol{F}}}\n",
" + \\mathbf{\\tilde{T}} \\frac{\\partial f_t(\\mathbf{T})}{\\partial \\mathcal{\\boldsymbol{F}}} \\, .\n",
"```\n",
"\n",
"The determination of the Lagrange multipliers $\\boldsymbol{\\Lambda}$ and $\\mathbf{\\tilde{T}}$ is identical to the procedures described above.\n",
"The dipole moment as a derivative of the energy can then be calculated as\n",
"```{math}\n",
"%:label: eq:dipmom_explicit\n",
" \\boldsymbol{\\mu} = \\sum_{pq} ( \\gamma_{pq}' + \\lambda_{pq} ) \\mu_{qp} \\, ,\n",
"```\n",
"where $\\boldsymbol{\\gamma}' = \\boldsymbol{\\gamma} + \\boldsymbol{\\gamma}^{\\text{A}}$ is the (orbital) __unrelaxed__ one-particle density matrix,\n",
"which includes contributions from the amplitude response, and $\\boldsymbol{\\gamma}' + \\boldsymbol{\\Lambda}$\n",
"is referred to as the (orbital) __relaxed density matrix__, and $\\mu_{qp}$ are elements of the dipole operator in the MO basis.\n",
"Taking all terms of the above equation into account is referred to as the _relaxed_ dipole moment,\n",
"whereas neglecting $\\boldsymbol{\\Lambda}$ yields so-called _unrelaxed_ dipole moments.\n",
"The latter do not correspond to a proper energy derivative, but the computation is somewhat simplified\n",
"since the iterative solution of the orbital response equations is avoided.\n",
"\n",
"Another approach to one-electron properties such as dipole moments is to set out from the expectation value\n",
"of the corresponding operator with the wave function, such as $\\boldsymbol{\\mu} = \\langle \\Psi | \\hat{\\mu} | \\Psi \\rangle$,\n",
"see the [Hellmann-Feynman theorem](eq:hellmann_feynman) above.\n",
"Depending on the wave-function model, this approach can be equivalent to the orbital-unrelaxed approach,\n",
"such as in configuration interaction, but in particular for schemes based on perturbation theory,\n",
"all three different approaches to dipole moments differ {cite}`Hodecker2019`.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "69e0ee53-2f95-413c-a479-f52b610c35a8",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}