diff --git a/README.md b/README.md
index 35f788c..069c9d7 100644
--- a/README.md
+++ b/README.md
@@ -2,6 +2,24 @@
 
 TwoBody.jl: a Julia package for quantum mechanical two-body problems
 
+## Quick Start
+
+Install TwoBody.jl from the Julia REPL or a notebook:
+
+```julia
+import Pkg; Pkg.add(url="https://github.com/JuliaFewBody/TwoBody.jl.git")
+```
+
+Load the package and solve the hydrogen ground state:
+
+```julia
+using TwoBody
+
+H = Hamiltonian(Kinetic(hbar=1, m=1), Coulomb(coefficient=-1))
+BS = BasisSet(SimpleGaussianBasis(13.00773), SimpleGaussianBasis(1.962079))
+solve(H, BS; info=0).E[1]
+```
+
 ## Documentation
 
 - Home: https://juliafewbody.github.io/TwoBody.jl
diff --git a/docs/src/Free-Complement.md b/docs/src/Free-Complement.md
index 0946849..3283c54 100644
--- a/docs/src/Free-Complement.md
+++ b/docs/src/Free-Complement.md
@@ -153,17 +153,6 @@ for reference_energy in reference_energies
 end
 ```
 
-You can also generate a complement from one basis function:
-
-```@repl fc-single
-using TwoBody
-H = Hamiltonian(
-  Kinetic(hbar = 1, m = 1),
-  Coulomb(coefficient = -1),
-)
-FC(H, PowerSlaterBasis(0, 1.5))
-```
-
 ## Acknowledgments
 
 This work was developed on the basis of the fourth lecture in Section I of the
diff --git a/docs/src/GEM.md b/docs/src/GEM.md
index 6b9580d..8efc6b7 100644
--- a/docs/src/GEM.md
+++ b/docs/src/GEM.md
@@ -49,7 +49,12 @@ nothing # hide
 Define the Gaussian basis set.
 
 ```@example gem
-BS = GeometricBasisSet(GaussianBasis, 0.1, 10.0, 20)
+BS = BasisSet(
+  GaussianBasis(13.00773),
+  GaussianBasis(1.962079),
+  GaussianBasis(0.444529),
+  GaussianBasis(0.1219492),
+)
 nothing # hide
 ```
 
diff --git a/docs/src/Rayleigh-Ritz.md b/docs/src/Rayleigh-Ritz.md
index dd70dbb..68759e8 100644
--- a/docs/src/Rayleigh-Ritz.md
+++ b/docs/src/Rayleigh-Ritz.md
@@ -200,6 +200,8 @@ save("assets/RR_HA.svg", fig) # hide
 ```
 ![](assets/RR_HA.svg)
 
+The numerical solution agrees well with the analytical solution even in the excited state.
+
 ## Example of Spherical Oscillator
 
 Analytical solutions are implemented in [spherical oscillator](https://ohno.github.io/Antique.jl/stable/SphericalOscillator/).
@@ -254,7 +256,9 @@ save("assets/RR_SO.svg", fig) # hide
 ```
 ![](assets/RR_SO.svg)
 
-## STO-3G
+The numerical solution agrees well with the analytical solution even in the excited state.
+
+## Example of STO-3G
 
 This example reproduces the STO-3G calculation for hydrogen reported by [Pérez-Torres (2019)](https://doi.org/10.1021/acs.jchemed.8b00959). In the contracted calculation, the published coefficients are held fixed, and the resulting contracted function is supplied to the solver. In the uncontracted calculation, the three primitive functions are supplied separately, allowing the Rayleigh–Ritz solver to optimize their linear coefficients.
 
diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 2e637c5..cba7d1b 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -100,8 +100,7 @@ Because successive Markov-chain samples are correlated, the naive `standard_erro
    arXiv:1607.07892 (2016).
 5. J. M. Thijssen,
    [*Computational Physics*, 2nd ed.](https://doi.org/10.1017/CBO9781139171397),
-   Cambridge University Press (2007);
-   邦訳: 松田和典, 道廣嘉隆, 谷村吉隆, 高須昌子, 吉江友照 訳, 『計算物理学』, 丸善出版 (2012).
+   Cambridge University Press (2007).
 
 ## API reference
 
diff --git a/docs/src/VNN.md b/docs/src/VNN.md
index a7990f7..79f1555 100644
--- a/docs/src/VNN.md
+++ b/docs/src/VNN.md
@@ -17,7 +17,7 @@ E[\psi_\theta] =
 where the grid, Hamiltonian matrix ``\pmb{H}``, and radial Jacobian ``\pmb{J}``
 are provided by `FiniteDifferenceMethod`.
 
-## Standard model
+## Usage
 
 The two-argument `solve` method constructs a Lux network from `architecture`.
 
diff --git a/docs/src/developer.md b/docs/src/developer.md
index f43c214..34f9ce2 100644
--- a/docs/src/developer.md
+++ b/docs/src/developer.md
@@ -98,7 +98,7 @@ To register a release in the [General](https://github.com/JuliaRegistries/Genera
 
 ## Architecture
 
-`src/TwoBody.jl` defines the `TwoBody` module and includes the source files in dependency order. `Hamiltonian.jl` defines the shared problem representation. `Basis.jl` supports the Rayleigh–Ritz implementation, and `FDM.jl` supplies the discretization used by the variational neural-network method. The solver files extend `solve` for their respective method types.
+`src/TwoBody.jl` defines the `TwoBody` module and includes the source files in dependency order. `Hamiltonian.jl` defines the shared problem representation. `Basis.jl` supports Rayleigh–Ritz and GEM, and `FDM.jl` supplies the discretization used by the variational neural-network method. The solver files extend `solve` for their respective method types.
 
 ```mermaid
 ---
@@ -111,6 +111,7 @@ flowchart TD
   D["DB.jl"]
   B["Basis.jl"]
   R["Rayleigh-Ritz.jl"]
+  G["GEM.jl"]
   F["FDM.jl"]
   Q["QTT.jl"]
   N["VNN.jl"]
@@ -118,8 +119,9 @@ flowchart TD
   T["TwoBody.jl"]
 
   H --> D
-  H --> R & F & Q & N & V
-  B --> R
+  H --> R & G & F & Q & N & V
+  B --> R & G
+  R --> G
   F --> N
-  H & D & B & R & F & Q & N & V --> T
+  H & D & B & R & G & F & Q & N & V --> T
 ```
diff --git a/src/Rayleigh-Ritz.jl b/src/Rayleigh-Ritz.jl
index 4ac8c17..a1d2298 100644
--- a/src/Rayleigh-Ritz.jl
+++ b/src/Rayleigh-Ritz.jl
@@ -568,7 +568,7 @@ This function returns the overlap matrix $\pmb{S}$. The element is written as ``
 `matrix(operator::Operator, basisset::BasisSet)`
 
 !!! note
-  This function is used for the expectation values and is not used in computing the Hamiltonian matrix.
+  This function computes expectation-value matrices, not the Hamiltonian matrix.
 
 This function returns the matrix corresponding to the operator in the given basis set. The element is written as ``O_{ij} = \langle \phi_{i} | \hat{o} | \phi_{j} \rangle``.
 """ matrix(operator::Operator, basisset::BasisSet)