Skip to content

Commit 39de6fd

Browse files
esraayantunaesraayantuna
authored
Merge branch 'main' into docs/compute-expval-guide
2 parents de88f40 + ec00d77 commit 39de6fd

37 files changed

+2754
-262
lines changed

‎docs/explanations/diag-coulomb-hamiltonian.ipynb‎

Lines changed: 307 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,307 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Diagonal Coulomb Hamiltonians\n",
8+
"\n",
9+
"This page explains the diagonal Coulomb Hamiltonian, a restricted but important class of fermionic Hamiltonians.\n",
10+
"\n",
11+
"## Definition\n",
12+
"\n",
13+
"A diagonal Coulomb Hamiltonian has the form\n",
14+
"\n",
15+
"$$\n",
16+
" H = \\sum_{\\sigma, pq} h_{pq} a^\\dagger_{\\sigma, p} a_{\\sigma, q}\n",
17+
" + \\frac12 \\sum_{\\sigma \\tau, pq} J^{\\sigma \\tau}_{pq} n_{\\sigma, p}\n",
18+
" n_{\\tau, q} + \\text{constant}\n",
19+
"$$\n",
20+
"\n",
21+
"where $n_{\\sigma, p} = a^\\dagger_{\\sigma, p} a_{\\sigma, p}$ is the number operator on orbital $p$ with spin $\\sigma$.\n",
22+
"\n",
23+
"The Hamiltonian is specified by the following data ($N$ is the number of spatial orbitals):\n",
24+
"\n",
25+
"- The one-body tensor $h_{pq}$, which is an $N \\times N$ Hermitian matrix.\n",
26+
"- The diagonal Coulomb matrices $J^{\\sigma \\tau}$, given as a pair of $N \\times N$ real symmetric matrices specifying the independent coefficients for alpha-alpha and alpha-beta interactions. We require that $J^{\\alpha\\alpha} = J^{\\beta\\beta}$ and $J^{\\alpha\\beta} = J^{\\beta\\alpha}$, so only two matrices are needed.\n",
27+
"- The constant, which is a real number.\n",
28+
"\n",
29+
"Compared to the general [molecular Hamiltonian](hamiltonians.ipynb), the two-body interaction is restricted to density-density form ($n_{\\sigma, p} n_{\\tau, q}$ terms). This restriction means the two-body part is specified by $O(N^2)$ parameters rather than $O(N^4)$, which enables more efficient simulation.\n",
30+
"\n",
31+
"## Data representation\n",
32+
"\n",
33+
"A diagonal Coulomb Hamiltonian is represented in ffsim using the [DiagonalCoulombHamiltonian](../api/ffsim.rst#ffsim.DiagonalCoulombHamiltonian) class. You can initialize it by passing the one-body tensor, the diagonal Coulomb matrices, and an optional constant."
34+
]
35+
},
36+
{
37+
"cell_type": "code",
38+
"execution_count": 1,
39+
"metadata": {},
40+
"outputs": [
41+
{
42+
"name": "stdout",
43+
"output_type": "stream",
44+
"text": [
45+
"Number of orbitals: 4\n",
46+
"One-body tensor shape: (4, 4)\n",
47+
"Diagonal Coulomb matrices shape: (2, 4, 4)\n",
48+
"Constant: -1.6404178369858733\n"
49+
]
50+
}
51+
],
52+
"source": [
53+
"import numpy as np\n",
54+
"\n",
55+
"import ffsim\n",
56+
"\n",
57+
"# Use 4 spatial orbitals, as an example.\n",
58+
"norb = 4\n",
59+
"\n",
60+
"rng = np.random.default_rng(12345)\n",
61+
"one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)\n",
62+
"diag_coulomb_mats = np.array(\n",
63+
" [\n",
64+
" ffsim.random.random_real_symmetric_matrix(norb, seed=rng),\n",
65+
" ffsim.random.random_real_symmetric_matrix(norb, seed=rng),\n",
66+
" ]\n",
67+
")\n",
68+
"constant = rng.standard_normal()\n",
69+
"\n",
70+
"hamiltonian = ffsim.DiagonalCoulombHamiltonian(\n",
71+
" one_body_tensor=one_body_tensor,\n",
72+
" diag_coulomb_mats=diag_coulomb_mats,\n",
73+
" constant=constant,\n",
74+
")\n",
75+
"\n",
76+
"print(f\"Number of orbitals: {hamiltonian.norb}\")\n",
77+
"print(f\"One-body tensor shape: {hamiltonian.one_body_tensor.shape}\")\n",
78+
"print(f\"Diagonal Coulomb matrices shape: {hamiltonian.diag_coulomb_mats.shape}\")\n",
79+
"print(f\"Constant: {hamiltonian.constant}\")"
80+
]
81+
},
82+
{
83+
"cell_type": "markdown",
84+
"metadata": {},
85+
"source": [
86+
"The `diag_coulomb_mats` attribute has shape `(2, N, N)`. The first matrix (`diag_coulomb_mats[0]`) contains the alpha-alpha (equivalently, beta-beta) interactions, and the second matrix (`diag_coulomb_mats[1]`) contains the alpha-beta (equivalently, beta-alpha) interactions.\n",
87+
"\n",
88+
"## Example: Fermi-Hubbard model\n",
89+
"\n",
90+
"The [Fermi-Hubbard model](https://en.wikipedia.org/wiki/Hubbard_model) is a widely studied example of a Hamiltonian with diagonal Coulomb form, since its two-body interaction is an onsite density-density interaction. The two-dimensional Fermi-Hubbard Hamiltonian on a square lattice is given by\n",
91+
"\n",
92+
"$$\n",
93+
" H = -t \\sum_{\\sigma, \\langle pq \\rangle}\n",
94+
" (a^\\dagger_{\\sigma, p} a_{\\sigma, q} + a^\\dagger_{\\sigma, q} a_{\\sigma, p})\n",
95+
" + U \\sum_p n_{\\alpha, p} n_{\\beta, p}\n",
96+
" - \\mu \\sum_p (n_{\\alpha, p} + n_{\\beta, p})\n",
97+
"$$\n",
98+
"\n",
99+
"where $t$ is the tunneling amplitude, $U$ is the onsite interaction strength, and $\\mu$ is the chemical potential. The index $\\langle pq \\rangle$ runs over pairs of neighboring orbitals on the lattice.\n",
100+
"\n",
101+
"The two-dimensional Fermi-Hubbard Hamiltonian can be constructed as a [FermionOperator](../api/ffsim.rst#ffsim.FermionOperator) using the [fermi_hubbard_2d](../api/ffsim.rst#ffsim.fermi_hubbard_2d) function, and then converted to a `DiagonalCoulombHamiltonian` using the `from_fermion_operator` method."
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": 2,
107+
"metadata": {},
108+
"outputs": [
109+
{
110+
"name": "stdout",
111+
"output_type": "stream",
112+
"text": [
113+
"One-body tensor:\n",
114+
"[[-2. -1. -1. 0.]\n",
115+
" [-1. -2. 0. -1.]\n",
116+
" [-1. 0. -2. -1.]\n",
117+
" [ 0. -1. -1. -2.]]\n",
118+
"\n",
119+
"Diagonal Coulomb matrix (alpha-alpha):\n",
120+
"[[0. 0. 0. 0.]\n",
121+
" [0. 0. 0. 0.]\n",
122+
" [0. 0. 0. 0.]\n",
123+
" [0. 0. 0. 0.]]\n",
124+
"\n",
125+
"Diagonal Coulomb matrix (alpha-beta):\n",
126+
"[[4. 0. 0. 0.]\n",
127+
" [0. 4. 0. 0.]\n",
128+
" [0. 0. 4. 0.]\n",
129+
" [0. 0. 0. 4.]]\n",
130+
"\n",
131+
"Constant: 0\n"
132+
]
133+
}
134+
],
135+
"source": [
136+
"norb_x = 2\n",
137+
"norb_y = 2\n",
138+
"norb = norb_x * norb_y\n",
139+
"nelec = (2, 2)\n",
140+
"\n",
141+
"# Build the Fermi-Hubbard Hamiltonian as a FermionOperator\n",
142+
"hubbard_op = ffsim.fermi_hubbard_2d(\n",
143+
" norb_x=norb_x,\n",
144+
" norb_y=norb_y,\n",
145+
" tunneling=1.0,\n",
146+
" interaction=4.0,\n",
147+
" chemical_potential=2.0,\n",
148+
")\n",
149+
"\n",
150+
"# Convert to a DiagonalCoulombHamiltonian\n",
151+
"hubbard_ham = ffsim.DiagonalCoulombHamiltonian.from_fermion_operator(hubbard_op)\n",
152+
"\n",
153+
"print(f\"One-body tensor:\\n{hubbard_ham.one_body_tensor.real}\")\n",
154+
"print(f\"\\nDiagonal Coulomb matrix (alpha-alpha):\\n{hubbard_ham.diag_coulomb_mats[0]}\")\n",
155+
"print(f\"\\nDiagonal Coulomb matrix (alpha-beta):\\n{hubbard_ham.diag_coulomb_mats[1]}\")\n",
156+
"print(f\"\\nConstant: {hubbard_ham.constant}\")"
157+
]
158+
},
159+
{
160+
"cell_type": "markdown",
161+
"metadata": {},
162+
"source": [
163+
"We can verify that the `DiagonalCoulombHamiltonian` representation yields the same ground state energy as the original `FermionOperator`."
164+
]
165+
},
166+
{
167+
"cell_type": "code",
168+
"execution_count": 3,
169+
"metadata": {},
170+
"outputs": [
171+
{
172+
"name": "stdout",
173+
"output_type": "stream",
174+
"text": [
175+
"Energy (DiagonalCoulombHamiltonian): -10.10274848346205\n",
176+
"Energy (FermionOperator): -10.102748483462063\n"
177+
]
178+
}
179+
],
180+
"source": [
181+
"import scipy.sparse.linalg\n",
182+
"\n",
183+
"# Compute the ground state energy using the DiagonalCoulombHamiltonian\n",
184+
"linop = ffsim.linear_operator(hubbard_ham, norb=norb, nelec=nelec)\n",
185+
"eigs, _ = scipy.sparse.linalg.eigsh(linop, k=1, which=\"SA\")\n",
186+
"energy_diag_coulomb = eigs[0]\n",
187+
"\n",
188+
"# Compute the ground state energy using the FermionOperator\n",
189+
"linop_op = ffsim.linear_operator(hubbard_op, norb=norb, nelec=nelec)\n",
190+
"eigs_op, _ = scipy.sparse.linalg.eigsh(linop_op, k=1, which=\"SA\")\n",
191+
"energy_op = eigs_op[0]\n",
192+
"\n",
193+
"print(f\"Energy (DiagonalCoulombHamiltonian): {energy_diag_coulomb}\")\n",
194+
"print(f\"Energy (FermionOperator): {energy_op}\")\n",
195+
"\n",
196+
"np.testing.assert_allclose(energy_diag_coulomb, energy_op)"
197+
]
198+
},
199+
{
200+
"cell_type": "markdown",
201+
"metadata": {},
202+
"source": [
203+
"## Time evolution via Trotter-Suzuki formulas\n",
204+
"\n",
205+
"The diagonal Coulomb Hamiltonian can be decomposed into two terms for the purpose of Trotter-Suzuki simulation:\n",
206+
"\n",
207+
"$$\n",
208+
" H = H_0 + H_1 + \\text{constant},\n",
209+
"$$\n",
210+
"\n",
211+
"where\n",
212+
"\n",
213+
"$$\n",
214+
" H_0 = \\sum_{\\sigma, pq} h_{pq} a^\\dagger_{\\sigma, p} a_{\\sigma, q}\n",
215+
"$$\n",
216+
"\n",
217+
"is a quadratic Hamiltonian and\n",
218+
"\n",
219+
"$$\n",
220+
" H_1 = \\frac12 \\sum_{\\sigma \\tau, pq} J^{\\sigma \\tau}_{pq} n_{\\sigma, p} n_{\\tau, q}\n",
221+
"$$\n",
222+
"\n",
223+
"is a diagonal Coulomb operator. Note,\n",
224+
"\n",
225+
"- $H_0$ can be simulated exactly using [orbital rotations](orbital-rotation.ipynb#Application-to-time-evolution-by-a-quadratic-Hamiltonian).\n",
226+
"- $H_1$ can be simulated exactly using [apply_diag_coulomb_evolution](../api/ffsim.rst#ffsim.apply_diag_coulomb_evolution).\n",
227+
"\n",
228+
"This split-operator approach is implemented in ffsim by the function [simulate_trotter_diag_coulomb_split_op](../api/ffsim.rst#ffsim.simulate_trotter_diag_coulomb_split_op). As with other Trotter functions in ffsim, the `order` parameter controls the order of the Trotter-Suzuki formula: `order=0` is the first-order asymmetric formula (the default), `order=1` is the first-order symmetric (second-order) formula, and so on.\n",
229+
"\n",
230+
"In the following example, we compare the Trotter-approximated time evolution against the exact time evolution (computed using `expm_multiply`)."
231+
]
232+
},
233+
{
234+
"cell_type": "code",
235+
"execution_count": 4,
236+
"metadata": {},
237+
"outputs": [
238+
{
239+
"name": "stdout",
240+
"output_type": "stream",
241+
"text": [
242+
"n_steps = 1, fidelity = 0.45702529\n",
243+
"n_steps = 2, fidelity = 0.95880093\n",
244+
"n_steps = 5, fidelity = 0.99915103\n",
245+
"n_steps = 10, fidelity = 0.99994861\n"
246+
]
247+
}
248+
],
249+
"source": [
250+
"# Initial state\n",
251+
"vec = ffsim.hartree_fock_state(norb, nelec)\n",
252+
"\n",
253+
"# Evolution time\n",
254+
"time = 1.0\n",
255+
"\n",
256+
"# Compute exact time evolution\n",
257+
"linop = ffsim.linear_operator(hubbard_ham, norb=norb, nelec=nelec)\n",
258+
"trace = ffsim.trace(hubbard_ham, norb=norb, nelec=nelec)\n",
259+
"exact = scipy.sparse.linalg.expm_multiply(\n",
260+
" -1j * time * linop, vec, traceA=-1j * time * trace\n",
261+
")\n",
262+
"\n",
263+
"# Compute Trotter-approximated time evolution with increasing number of steps\n",
264+
"for n_steps in [1, 2, 5, 10]:\n",
265+
" result = ffsim.simulate_trotter_diag_coulomb_split_op(\n",
266+
" vec,\n",
267+
" hubbard_ham,\n",
268+
" time,\n",
269+
" norb=norb,\n",
270+
" nelec=nelec,\n",
271+
" n_steps=n_steps,\n",
272+
" order=1,\n",
273+
" )\n",
274+
" fidelity = abs(np.vdot(result, exact))\n",
275+
" print(f\"n_steps = {n_steps:2d}, fidelity = {fidelity:.8f}\")"
276+
]
277+
},
278+
{
279+
"cell_type": "markdown",
280+
"metadata": {},
281+
"source": [
282+
"As expected, the fidelity increases with the number of Trotter steps."
283+
]
284+
}
285+
],
286+
"metadata": {
287+
"kernelspec": {
288+
"display_name": ".venv",
289+
"language": "python",
290+
"name": "python3"
291+
},
292+
"language_info": {
293+
"codemirror_mode": {
294+
"name": "ipython",
295+
"version": 3
296+
},
297+
"file_extension": ".py",
298+
"mimetype": "text/x-python",
299+
"name": "python",
300+
"nbconvert_exporter": "python",
301+
"pygments_lexer": "ipython3",
302+
"version": "3.12.11"
303+
}
304+
},
305+
"nbformat": 4,
306+
"nbformat_minor": 2
307+
}

0 commit comments

Comments
 (0)