Wavefunctions/OrbitalRotation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "958b54cd-c89c-474a-a1a2-368fd69c0af4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "init_printing()\n",
    "import scipy.linalg\n",
    "import numpy as np\n",
    "import warnings\n",
    "from matplotlib import MatplotlibDeprecationWarning\n",
    "warnings.filterwarnings(\"ignore\",category=MatplotlibDeprecationWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "718952ee-5776-41f6-a3de-664e19a227b5",
   "metadata": {},
   "source": [
    "# Orbital Rotation\n",
    "\n",
    "One method to optimize the single particle orbitals in the Slater determinant is to mix occupied and unoccupied orbitals.  The orbitals are orthogonal and normalized, and any mixing transformation should preserve that.  A rotation matrix meets those conditions.  However, the entries in a rotation matrix are not independent.  A rotation matrix can be expressed as an exponential of a skew-symmetric matrix, and the entries in that matrix are independent.\n",
    "\n",
    "See Chapter 3 of \"Molecular Electronic Structure Theory\" by Trygve Helgaker, Poul Jorgensen, and Jeppe Olsen\n",
    "\n",
    "See also the Wikipedia page on rotation matrices https://en.wikipedia.org/wiki/Rotation_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d72a57d-7f14-45be-b315-0960f1c28289",
   "metadata": {},
   "source": [
    "### Size 2 example\n",
    "Write a 2x2 skew-symmetric matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "50e38434-f661-4b45-8866-f322a98b3e7d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & \\kappa\\\\- \\kappa & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   κ⎤\n",
       "⎢     ⎥\n",
       "⎣-κ  0⎦"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kappa = Symbol('kappa',real=True)\n",
    "m = Matrix([[0,kappa],[-kappa,0]])\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da647774-1ee1-4b03-ae45-ea93b1893b14",
   "metadata": {},
   "source": [
    "Apply the matrix exponential and we get the 2x2 rotation matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fc1b6833-b516-4041-97d4-9697405b01eb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\cos{\\left(\\kappa \\right)} & \\sin{\\left(\\kappa \\right)}\\\\- \\sin{\\left(\\kappa \\right)} & \\cos{\\left(\\kappa \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡cos(κ)   sin(κ)⎤\n",
       "⎢               ⎥\n",
       "⎣-sin(κ)  cos(κ)⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4241bc8a-f3d0-449b-bb7c-a77db8ef2d8f",
   "metadata": {},
   "source": [
    "We could do the 3x3 case, but the form looks different from how a 3x3 rotation matrix is usually expressed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84fa1b48-ace6-4f10-a607-b2856f6f1c76",
   "metadata": {},
   "source": [
    "# Combining rotation matrices\n",
    "Note that if X and Y do not commute, then $\\exp(X)\\exp(Y) \\ne \\exp(X+Y)$.\n",
    "We can demonstrate this with a 3x3 rotation matrix.  The symbolic form would get rather lengthy, so we will use an example with concrete values for the rotation parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ad98397b-29e9-47e0-b3de-b44c21446d1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "k1,k2,k3 = symbols('kappa_1 kappa_2 kappa_3',real=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "25897f8d-1f81-4490-a4a4-417635ceeb1e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & \\kappa_{1} & \\kappa_{2}\\\\- \\kappa_{1} & 0 & \\kappa_{3}\\\\- \\kappa_{2} & - \\kappa_{3} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ 0   κ₁   κ₂⎤\n",
       "⎢            ⎥\n",
       "⎢-κ₁   0   κ₃⎥\n",
       "⎢            ⎥\n",
       "⎣-κ₂  -κ₃  0 ⎦"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 3x3 skew-symmetric matrix.  Note the signs are different than is traditionally written.\n",
    "m1 = Matrix([[0, k1, k2],[-k1, 0, k3], [-k2, -k3, 0]])\n",
    "m1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "15eaf044-e4bf-4086-b09d-218251cee5ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}0 & 0.1 & 0.2\\\\-0.1 & 0 & 0\\\\-0.2 & 0 & 0\\end{matrix}\\right], \\  \\left[\\begin{matrix}0 & 0.15 & 0.5\\\\-0.15 & 0 & 0\\\\-0.5 & 0 & 0\\end{matrix}\\right]\\right)$"
      ],
      "text/plain": [
       "⎛⎡ 0    0.1  0.2⎤  ⎡  0    0.15  0.5⎤⎞\n",
       "⎜⎢              ⎥  ⎢                ⎥⎟\n",
       "⎜⎢-0.1   0    0 ⎥, ⎢-0.15   0     0 ⎥⎟\n",
       "⎜⎢              ⎥  ⎢                ⎥⎟\n",
       "⎝⎣-0.2   0    0 ⎦  ⎣-0.5    0     0 ⎦⎠"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m2 = m1.subs({k1:0.1, k2:0.2, k3:0})\n",
    "m3 = m1.subs({k1:0.15, k2:0.5, k3:0.0})\n",
    "m2,m3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "d1d07d0e-bb83-4dcf-99f7-28f63ab856b8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0\\\\0 & 0 & -0.02\\\\0 & 0.02 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0   0      0  ⎤\n",
       "⎢              ⎥\n",
       "⎢0   0    -0.02⎥\n",
       "⎢              ⎥\n",
       "⎣0  0.02    0  ⎦"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# m2 and m3 do not commute\n",
    "m2*m3 - m3*m2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "bfc62f57-0b1a-4bb1-88dd-cdfc2e8229b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.73630111157537 & 0.230520546096726 & 0.636176823627002\\\\-0.223771239225615 & 0.970234879766162 & -0.0925781323022542\\\\-0.638582105596693 & -0.074192694467908 & 0.765968888728051\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ 0.73630111157537   0.230520546096726    0.636176823627002 ⎤\n",
       "⎢                                                           ⎥\n",
       "⎢-0.223771239225615  0.970234879766162   -0.0925781323022542⎥\n",
       "⎢                                                           ⎥\n",
       "⎣-0.638582105596693  -0.074192694467908   0.765968888728051 ⎦"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The product of exponentials\n",
    "exp(m2) * exp(m3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "ff47753f-56ff-4a4e-9254-7a1fcff59c16",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0.736237065559526 & 0.227606815284241 & 0.637299082795875\\\\-0.227606815284241 & 0.970162563977322 & -0.0835448208634987\\\\-0.637299082795875 & -0.0835448208634986 & 0.766074501582204\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡0.736237065559526    0.227606815284241    0.637299082795875 ⎤\n",
       "⎢                                                            ⎥\n",
       "⎢-0.227606815284241   0.970162563977322   -0.0835448208634987⎥\n",
       "⎢                                                            ⎥\n",
       "⎣-0.637299082795875  -0.0835448208634986   0.766074501582204 ⎦"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The exponential of the sum.  Not the same as the product of exponentials\n",
    "exp(m2+m3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "922be045-2999-44c0-8531-c6aad6c02d88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0.25 & 0.7\\\\-0.25 & 0 & 0\\\\-0.7 & 0 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡  0    0.25  0.7⎤\n",
       "⎢                ⎥\n",
       "⎢-0.25   0     0 ⎥\n",
       "⎢                ⎥\n",
       "⎣-0.7    0     0 ⎦"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m2 + m3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "98036e55-8ea5-4ff9-b786-421113c1db9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.46454773e-16,  2.49492270e-01,  7.00084196e-01],\n",
       "       [-2.49492270e-01,  2.11772169e-16, -1.00970891e-02],\n",
       "       [-7.00084196e-01,  1.00970891e-02, -2.34862926e-16]])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Take the matrx log of the product.  Sympy doesn't have the matrix log, must use scipy instead\n",
    "# Read off the new values for k1,k2 and k3 from the skew-symmetric matrix form. \n",
    "# The rotation parameters are small, so the values are close to m2+m3, but not quite the same.\n",
    "m4 = exp(m2) * exp(m3)\n",
    "m4p = np.array(m4).astype(np.float64)\n",
    "scipy.linalg.logm(m4p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1facfa1b-88b6-41eb-b5fd-d56c3e150bb9",
   "metadata": {},
   "source": [
    "This is an issue with handling rotation parameters, because the code assumes that variational parameters can simply be summed.\n",
    "When applying multiple rotations, this is no longer the case, it requries more careful handling of the rotation parameters.\n",
    "To get an exact answer, we take the exponentials, multiple the matrices, then take the matrix log to recover the parameters.\n",
    "\n",
    "\n",
    "As an aside, the Baker-Campbell-Hausdorff formula could be used to approximate the matrix sum in terms of nested commutators."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40f1804f-8e60-4730-833e-8a2ec0a7dd66",
   "metadata": {},
   "source": [
    "# Fill-in\n",
    "This is an issue particular to QMC's usage of rotation matrices.  The only rotations of interest (that have non-zero derivatives) are between occupied and unoccupied orbitals.  The rotations between only occupied orbitals or only unoccupied orbitals do not change the energy.  However, when combining rotation matrices by the above procedure, the values corresponding to those rotations may become non-zero.  That means capturing the state of the rotation requires more than just the occupied-unoccupied rotational parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e3e3555-4941-40a5-9bd1-4f6cd288e0b5",
   "metadata": {},
   "source": [
    "## Derivatives\n",
    "\n",
    "To compute parameter derivatives of a wavefunction using a the rotation matrix, we can use Jacobi's formula (https://en.wikipedia.org/wiki/Jacobi%27s_formula)\n",
    "$$\n",
    "\\frac{1}{\\det(A)} \\frac{d}{dt} \\det(A) = \\mathrm{tr} (A^{-1} \\frac{dA}{dt})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84397a07-f98f-400e-ae5a-f64e197b2c0d",
   "metadata": {},
   "source": [
    "More detail about the argument of the determinant.  It is a matrix of single particle orbitals (SPO) evaluated at the different electron positions.  Without orbital rotation, exactly the same number of SPO's are used as electrons to yield a square matrix.  With the addition of a rotation matrix, more SPO's are used than electron positions, leaving us to deal with non-square matrices.\n",
    "$$\n",
    "A = E R \\Phi\n",
    "$$\n",
    "where $\\Phi$ is the matrix of SPOs, $R$ is the rotation matrix.\n",
    "Let $N$ be the number of electrons and $M$ be the number of SPOs.   The matrix $\\Phi$ is $M$ x $N$.  The rotation matrix, $R$, is $M$x$M$.  The product $R \\Phi$ is rectangular with dimensions $M$ x $N$.  The final matrix must be square, so formally we add a \"selection\" matrix, $E$, of dimensions $N$ x $M$.  The matrix $E$ is the identity in the square part, and zeros elsewhere."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beffc7c1-95c7-46e4-b337-69c32241168d",
   "metadata": {},
   "source": [
    "### Determinant of product\n",
    "\n",
    "For square matrices, $\\det(AB) = \\det(A)\\det(B)$.  For rectangular matrices, we need to use the Cauchy-Binet formulat (https://en.wikipedia.org/wiki/Cauchy%E2%80%93Binet_formula)\n",
    "The structure of E implies the only nonzero contribution ot the determinant comes from the \"square\" part."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9b3893c-7725-4964-87cd-85baefcb905f",
   "metadata": {},
   "source": [
    "### Derivative of rotation matrix\n",
    "We always evaluate the rotation matrix at zero angle.\n",
    "\n",
    "Expand the matrix exponential $R = \\exp(X)$ in the power series $1 + X + 1/2 X^2 + ...$\n",
    "\n",
    "The derivative is then\n",
    "$$\n",
    "\\frac{dR}{dp_i} = 0 + dX/dp_i + X (dX/dp_i)\n",
    "$$\n",
    "When evaluating at zero angle, all higher order terms invovling $X$ are zero.  So we are left with a single $dX/dp_i$ term.  That matrix has two entries corresponding to the location of the parameter $p_i$.  One is 1 and the other is -1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38e9ec35-1816-4de7-b5ab-9f8bcf2ead43",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}