Document Rayleigh-Ritz theory

docs/src/Rayleigh-Ritz.md

@@ -4,19 +4,94 @@ CurrentModule = TwoBody
 
 # Rayleigh-Ritz Method
 
-This method is one of the variational method. It solves the generalized eigenvalue problem,
+The Rayleigh-Ritz method is a variational method for approximating the eigenvalues and eigenfunctions of a Hamiltonian. The expectation value obtained from a trial wave function is an upper bound on the exact ground-state energy. Note that nonlinear parameters, such as the exponents of Gaussian basis functions, are not optimized by the linear procedure described below.
+
+## Theory
+
+In the Rayleigh-Ritz method, the trial wave function is written as a linear combination of basis functions ``\phi_1, \phi_2, \phi_3, \ldots``:
+
+```math
+\psi = \sum_i c_i \phi_i.
+```
+
+We then seek the expansion coefficients ``c_i`` that minimize the energy ``E[\psi]``. The coefficients are ultimately obtained by solving the generalized eigenvalue problem
+
+```math
+\boldsymbol{H}\boldsymbol{c} = E\boldsymbol{S}\boldsymbol{c}.
+```
+
+The derivation is given below. First, let us define the notation. ``\boldsymbol{H}`` is the Hamiltonian matrix, whose elements are
+
+```math
+H_{ij} = \langle \phi_i | \hat{H} | \phi_j \rangle,
+```
+
+and ``\boldsymbol{S}`` is the overlap matrix, whose elements are
+
+```math
+S_{ij} = \langle \phi_i | \phi_j \rangle.
+```
+
+Substituting ``\psi = \sum_i c_i \phi_i`` into ``E[\psi]`` gives
+
 ```math
-\pmb{H} \pmb{c} = E \pmb{S} \pmb{c}.
+\begin{aligned}
+E[\psi]
+&= \frac{\langle \psi | \hat{H} | \psi \rangle}
+        {\langle \psi | \psi \rangle} \\
+&= \frac{\left\langle \sum_i c_i \phi_i \middle| \hat{H} \middle| \sum_j c_j \phi_j \right\rangle}
+        {\left\langle \sum_i c_i \phi_i \middle| \sum_j c_j \phi_j \right\rangle} \\
+&= \frac{\sum_{i,j} c_i c_j \langle \phi_i | \hat{H} | \phi_j \rangle}
+        {\sum_{i,j} c_i c_j \langle \phi_i | \phi_j \rangle} \\
+&= \frac{\sum_{i,j} c_i c_j H_{ij}}
+        {\sum_{i,j} c_i c_j S_{ij}}.
+\end{aligned}
 ```
-The Hamiltonian matrix is defined as ``H_{ij} = \langle \phi_{i} | \hat{H} | \phi_{j} \rangle`` and the overlap matrix is defined as ``S_{ij} = \langle \phi_{i} | \phi_{j} \rangle``. The eigenvector ``\pmb{c}`` is the column of the optimal coefficients $c_i$ for the linear combination,
+
+Here and below, the basis functions and coefficients are assumed to be real. For a complex-valued basis, the coefficient on the bra side must be replaced by its complex conjugate ``c_i^*``. The expansion above follows directly from the distributive law and can be seen explicitly by writing ``\sum_i c_i \phi_i = c_1 \phi_1 + c_2 \phi_2 + \cdots``.
+
+We minimize this expression subject to the normalization constraint
+
 ```math
-\psi(r) = \sum_i c_i \phi_i(r),
+\sum_{i,j} c_i c_j S_{ij} = 1.
 ```
-to minmize the expectation value of the energy,
+
+Using the method of Lagrange multipliers, define
+
 ```math
-E = \frac{\langle\psi|\hat{H}|\psi\rangle}{\langle\psi|\psi\rangle}.
+L(c_i,E)
+= \sum_{i,j} c_i c_j H_{ij}
+- E \sum_{i,j} c_i c_j S_{ij}
++ E.
 ```
-Note that the nonlinear parameters (e.g., exponents of the Gaussian basis functions) are not optimized. The expectation by trial wavefunction is the upper bound for the exact energy.
+
+When ``L`` is differentiated with respect to ``c_k``, only terms for which ``i=k`` or ``j=k`` remain. Therefore,
+
+```math
+\begin{aligned}
+\frac{\partial L}{\partial c_k} &= 0, \\
+\sum_i c_i H_{ik} + \sum_j c_j H_{kj}
+- E \sum_i c_i S_{ik} - E \sum_j c_j S_{kj} &= 0, \\
+2 \sum_j H_{kj} c_j - 2E \sum_j S_{kj} c_j &= 0, \\
+\sum_j \left(H_{kj}c_j - E S_{kj}c_j\right) &= 0, \\
+\boldsymbol{H}\boldsymbol{c} - E\boldsymbol{S}\boldsymbol{c} &= 0, \\
+\boldsymbol{H}\boldsymbol{c} &= E\boldsymbol{S}\boldsymbol{c}.
+\end{aligned}
+```
+
+In the third line, we used ``H_{ij}=H_{ji}`` and ``S_{ij}=S_{ji}``, which hold for the real-valued case considered here. Similar derivations can be found in the following references:
+
+- J. M. Thijssen, *Computational Physics* (Japanese translation, Maruzen Publishing, 2012), Sec. 3.1, [Variational Calculations](https://www.maruzen-publishing.co.jp/item/b294596.html).
+- A. Szabo and N. S. Ostlund, *Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory*, Vol. 1 (Japanese translation, University of Tokyo Press, 1987), Sec. 1.3.2, [The Linear Variational Problem](http://www.utp.or.jp/book/b302128.html).
+
+In Julia, the generalized eigenvalue problem ``\boldsymbol{H}\boldsymbol{c} = E\boldsymbol{S}\boldsymbol{c}`` can be solved easily with the `LinearAlgebra` standard library:
+
+```julia
+using LinearAlgebra
+E, c = eigen(H, S)
+```
+
+It is therefore not necessary to know the details of the eigensolver here. The practical challenge is instead to evaluate the matrix elements ``H_{ij} = \langle \phi_i | \hat{H} | \phi_j \rangle`` and ``S_{ij} = \langle \phi_i | \phi_j \rangle``, and to choose basis functions for which these elements can be calculated efficiently.
 
 ## Usage