{
"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 when taking the derivative with respect to external **electromagnetic field**, or NMR and EPR parameters 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 required, 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(x)$ with respect to some parameter $x$ is to use finite difference approximations.\n",
"Choosing a small change $x_0$ in $x$, a two-point estimation is given by\n",
"```{math}\n",
"%:label: eq:divided_difference\n",
"\\frac{\\mathrm{d} E}{\\mathrm{d} x} \\approx \\frac{E(x + x_0) - E(x)}{x_0} \n",
"```\n",
"which is known as a first-order **divided difference** and its error to the true derivative is approximately proportional to $x_0$.\n",
"The true derivative of $E$ at $x$ is given by the limit\n",
"```{math}\n",
"%:label: eq:true_derivative\n",
"\\frac{\\mathrm{d} E}{\\mathrm{d} x} = \\lim_{x_0 \\to 0} \\frac{E(x + x_0) - E(x)}{x_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} x} \\approx \\frac{E(x + x_0) - E(x - x_0)}{2x_0} \n",
"```\n",
"where it can be shown that the first-order error cancels, so it is approximately proportional to $x_0^2$,\n",
"which means that for small $x_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(x)$ need to be performed in order to obtain\n",
"the derivative with respect to one variable $x$.\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} x} \\approx \\frac{E(x - 2 x_0) - 8 E(x - x_0) + 8 E(x + x_0) - E(x + 2 x_0)}{12x_0} \n",
"```\n",
"where the error is approximately proportional to $x_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 analytically, 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-centered 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 $x$ 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}x}=\\frac{\\partial E}{\\partial x}+\\frac{\\partial E}{\\partial\\mathbf{C}}\\frac{\\mathrm{d}\\mathbf{C}}{\\mathrm{d}x} \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 x$, is the Hellmann--Feynman contribution which describes the explicit dependence of the energy on the nuclear coordinate $x$ 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 $x$ due to the fact that the molecular orbitals are expanded in a finite atomic-centered 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 SCF 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 $x$ 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",
"%\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 $\\mathbf{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}x}$ 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 $x$ can be computed as\n",
"\\begin{equation*}\n",
"\\frac{\\mathrm{d}E}{\\mathrm{d}x}=\\frac{\\mathrm{d}L}{\\mathrm{d}x}=\\frac{\\partial L}{\\partial x}\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} \\langle ij || ij \\rangle\n",
"```\n",
"where $F_{pq}$ are Fock matrix elements, $\\langle pq || rs \\rangle$ are anti-symmetrized two-electron integrals in physicists' (\"1212\") notation {cite}`Szabo2012`. As usual, the indices $i,j,\\ldots$ denote occupied molecular orbitals, $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, $\\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 $x$ 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 x}=\\frac{\\partial E_\\mathrm{HF}}{\\partial x}+\\sum_{p,q}\\lambda_{pq}\\frac{\\partial F_{pq}}{\\partial x}+\\sum_{p,q}\\omega_{pq}\\frac{\\partial S_{pq}}{\\partial x}\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}\\langle pq || rs \\rangle\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 x}&=\\sum_{p,q}(\\lambda_{pq}+\\gamma_{pq}) F^{(x)}_{pq}+\\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}\\braket{pq||rs}^{(x)}+\\sum_{p,q}\\omega_{pq}S^{(x)}_{pq}\\nonumber\\\\\n",
"&=\\sum_{p,q}(\\lambda_{pq}+\\gamma_{pq}) h^{(x)}_{pq}+\\sum_{p,q}(\\lambda_{pq}+\\gamma_{pq})\\sum_{i}\\braket{pi||qi}^{(x)}+\\frac{1}{4}\\sum_{pqrs}\\Gamma_{pqrs}\\braket{pq||rs}^{(x)}+\\sum_{p,q}\\omega_{pq}S^{(x)}_{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 $x$,\n",
"i.e., with fixed MO coefficients $\\mathbf{C}$. Explicitly, they are given by {cite}`Levchenko2005`\n",
"\n",
"(eq:hpq_Spq)=\n",
"```{math}\n",
"%:label: eq:hpq\n",
"h_{pq}^{(x)}&=\\frac{\\partial h_{pq}}{\\partial x}=\\sum_{\\mu,\\nu}C_{\\mu p}h^{x}_{\\mu\\nu}C_{\\nu q} \\\\\n",
"h^{x}_{\\mu\\nu}&=\\bra{\\phi_\\mu}\\frac{\\partial \\hat{h}}{\\partial x}\\ket{\\phi_\\nu}+\\bra{\\frac{\\partial\\phi_\\mu}{\\partial x}}\\hat{h}\\ket{\\phi_\\nu}+\\bra{\\phi_\\mu}\\hat{h}\\ket{\\frac{\\partial\\phi_\\nu}{\\partial x}} \\\\\n",
"\\braket{pq||rs}^{(x)}&=\\frac{\\partial \\braket{pq||rs}}{\\partial x}=\\sum_{\\mu,\\nu,\\theta,\\sigma}C_{\\mu p}C_{\\nu q}\\braket{\\phi_{\\mu}\\phi_{\\nu}||\\phi_{\\theta}\\phi_{\\sigma}}^{x}C_{\\theta r}C_{\\sigma s} \\\\\n",
"\\braket{\\phi_{\\mu}\\phi_{\\nu}||\\phi_{\\theta}\\phi_{\\sigma}}^{x}&=\\braket{\\frac{\\partial\\phi_{\\mu}}{\\partial x}\\phi_{\\nu}||\\phi_{\\theta}\\phi_{\\sigma}}+\\braket{\\phi_{\\mu}\\frac{\\partial\\phi_{\\nu}}{\\partial x}||\\phi_{\\theta}\\phi_{\\sigma}}+\\braket{\\phi_{\\mu}\\phi_{\\nu}||\\frac{\\partial\\phi_{\\theta}}{\\partial x}\\phi_{\\sigma}}\\nonumber + \\braket{\\phi_{\\mu}\\phi_{\\nu}||\\phi_{\\theta}\\frac{\\partial\\phi_{\\sigma}}{\\partial x}} \\\\\n",
"S_{pq}^{(x)} &=\\frac{\\partial S_{pq}}{\\partial x}=\\sum_{\\mu,\\nu}C_{\\mu p}S_{\\mu\\nu}^x C_{\\nu q} \\\\\n",
"S_{\\mu\\nu}^x &=\\braket{\\frac{\\partial\\phi_\\mu}{\\partial x}|\\phi_\\nu}+\\braket{\\phi_\\mu|\\frac{\\partial\\phi_\\nu}{\\partial x}} \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 $x$, 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 x$: (1) the derivatives of the core-Hamiltonian matrix, 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 [the HF energy expression](eq:HF_energy_fct), we can immediately identify the non-vanishing 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",
"\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",
"\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} x}&=\\sum_{i,j}\\delta_{ij} h^{(x)}_{ij}+\\sum_{i,j}\\delta_{ij}\\sum_{k}\\braket{ik||jk}^{(x)} - \\frac{1}{2}\\sum_{i,j,k,l}\\delta_{ik}\\delta_{jl}\\braket{ij||kl}^{(x)}-\\sum_{i,j}\\epsilon_{i}\\delta_{ij}S^{(x)}_{ij}\\nonumber \\\\\n",
"&=\\sum_{i}h^{(x)}_{ii}+\\frac{1}{2}\\sum_{i,j}\\braket{ij||ij}^{(x)}-\\sum_{i}\\epsilon_{i}S^{(x)}_{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` $\\mathrm{d} V_{nn}/\\mathrm{d} x$.\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