Skip to content

Commit a6f9eb6

Browse files
RichRick1RichRick1
authored
Merge pull request #172 from giovanni-br/molkit
Reduced Hamiltonian and Energy Testing
2 parents 1848a6f + 451094d commit a6f9eb6

File tree

5 files changed

+267
-110
lines changed

5 files changed

+267
-110
lines changed

‎examples/molkit.ipynb

Lines changed: 77 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -32,17 +32,9 @@
3232
},
3333
{
3434
"cell_type": "code",
35-
"execution_count": null,
35+
"execution_count": 2,
3636
"metadata": {},
37-
"outputs": [
38-
{
39-
"name": "stdout",
40-
"output_type": "stream",
41-
"text": [
42-
"Max error: 1.0\n"
43-
]
44-
}
45-
],
37+
"outputs": [],
4638
"source": [
4739
"import numpy as np\n",
4840
"from moha import HamHeisenberg\n",
@@ -105,17 +97,86 @@
10597
"h1_spin, h2_spin = Molecular_Hamiltonian.spinize_H()\n",
10698
"h2_spin = antisymmetrize_two_body(h2_spin)\n",
10799
"\n",
108-
"G = MolHam.to_geminal(h2_spin) \n",
109-
"n_orb = h2_spin.shape[0]\n",
100+
"reduced = Molecular_Hamiltonian.to_reduced(n_elec=2)"
101+
]
102+
},
103+
{
104+
"cell_type": "markdown",
105+
"metadata": {},
106+
"source": [
107+
"# Energies Tests"
108+
]
109+
},
110+
{
111+
"cell_type": "code",
112+
"execution_count": 4,
113+
"metadata": {},
114+
"outputs": [
115+
{
116+
"name": "stdout",
117+
"output_type": "stream",
118+
"text": [
119+
"E_spin-orbital : 1.004059262354666e-12\n",
120+
"E_geminal : 1.0040474662350293e-12\n"
121+
]
122+
}
123+
],
124+
"source": [
125+
"import numpy as np\n",
126+
"from pyscf import fci\n",
127+
"from moha.molkit.hamiltonians import MolHam\n",
128+
"\n",
110129
"\n",
111-
"V = from_geminal(G, n_orb = n_orb)\n",
112-
"print(\"Max error:\", np.max(np.abs(h2_spin - V)))\n"
130+
"h_sp = np.array([[1.0, 0.2],\n",
131+
" [0.2, 1.5]])\n",
132+
"\n",
133+
"V_sp = np.zeros((2, 2, 2, 2))\n",
134+
"V_sp[0, 0, 0, 0] = 0.7\n",
135+
"V_sp[1, 1, 1, 1] = 0.6\n",
136+
"\n",
137+
"N_ELEC, n_spin = 2, 2 \n",
138+
"toy = MolHam(h_sp, V_sp)\n",
139+
"\n",
140+
"h_spin, V_spin = toy.spinize_H() \n",
141+
"norb = h_spin.shape[0]\n",
142+
"e_fci, ci_vec = fci.direct_spin1.kernel(h_spin, V_spin, norb, N_ELEC, ecore=0.0)\n",
143+
"rdm1s, rdm2s = fci.direct_spin1.make_rdm12s(ci_vec, norb, N_ELEC)\n",
144+
"\n",
145+
"rdm1 = np.block([\n",
146+
" [rdm1s[0], np.zeros_like(rdm1s[0])],\n",
147+
" [np.zeros_like(rdm1s[1]), rdm1s[1]]\n",
148+
"])\n",
149+
"\n",
150+
"n = norb // 2 # spatial orbitals\n",
151+
"rdm2 = np.zeros((2*n, 2*n, 2*n, 2*n))\n",
152+
"for i in range(n):\n",
153+
" for j in range(n):\n",
154+
" for k in range(n):\n",
155+
" for l in range(n):\n",
156+
" # αααα\n",
157+
" rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
158+
" # ββββ\n",
159+
" rdm2[i+n, j+n, k+n, l+n] = rdm2s[2][i, k, j, l]\n",
160+
" # αβαβ\n",
161+
" rdm2[i, j+n, k, l+n] = rdm2s[1][i, k, j, l]\n",
162+
" # βαβα (hermitian transpose of αβαβ)\n",
163+
" rdm2[i+n, j, k+n, l] = rdm2s[1][i, k, j, l]\n",
164+
"\n",
165+
"H = toy.to_reduced(n_elec=N_ELEC) # H = h↑ + 0.5 * V\n",
166+
"E_spin = np.einsum('pqrs,pqrs', H, rdm2)\n",
167+
"H_gem = toy.to_geminal(H, type='h2')\n",
168+
"rdm_gem = toy.to_geminal(rdm2, type='rdm2')\n",
169+
"E_gem = np.einsum('pq,pq', H_gem, rdm_gem)\n",
170+
"\n",
171+
"print(\"E_spin-orbital :\", E_spin)\n",
172+
"print(\"E_geminal :\", E_gem)\n",
173+
"assert np.allclose(E_spin, E_gem, atol=1e-10), \"Energies do not match!\""
113174
]
114175
}
115176
],
116177
"metadata": {
117178
"kernelspec": {
118-
"display_name": "base",
179+
"display_name": "moha",
119180
"language": "python",
120181
"name": "python3"
121182
},
@@ -129,7 +190,7 @@
129190
"name": "python",
130191
"nbconvert_exporter": "python",
131192
"pygments_lexer": "ipython3",
132-
"version": "3.12.3"
193+
"version": "3.11.13"
133194
}
134195
},
135196
"nbformat": 4,

‎moha/molkit/hamiltonians.py

Lines changed: 117 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,10 @@
11
r"""Molkit Module."""
22

3-
from .utils.spinops import antisymmetrize_two_body, get_spin_blocks
3+
from .utils.spinops import (
4+
antisymmetrize_two_body,
5+
get_spin_blocks,
6+
upscale_one_body
7+
)
48
import numpy as np
59

610

@@ -28,12 +32,100 @@ def __init__(self, one_body, two_body):
2832
self.n_spin = 2 * self.n_spatial
2933
self.reduced_ham = None
3034

35+
def antisymmetrize(self):
36+
"""Apply proper antisymmetrization to two-electron integrals."""
37+
self.two_body = antisymmetrize_two_body(self.two_body, inplace=True)
38+
39+
def get_spin_blocks(self):
40+
"""Return the main spin blocks of the two-body spin-orbital tensor.
41+
42+
Returns
43+
-------
44+
dict
45+
Dictionary with spin blocks: 'aaaa', 'bbbb', 'abab'
46+
47+
"""
48+
if not hasattr(self, "two_body_spin"):
49+
raise RuntimeError(
50+
"Call .spinize() first to compute spin-orbital form.")
51+
52+
return get_spin_blocks(self.two_body_spin, self.n_spatial)
53+
54+
def to_geminal(self, two_body=None, type='h2'):
55+
r"""
56+
Convert the two-body term to the geminal basis.
57+
58+
Parameters
59+
----------
60+
two_body : np.ndarray
61+
Two-body term in spin-orbital basis in physics notation.
62+
type : str
63+
['rdm2', 'h2']. Type of the two-body term.
64+
- 'rdm2' : 2 body reduced density matrix
65+
- 'h2' : 2 body Hamiltonian
66+
67+
Returns
68+
-------
69+
two_body_gem : np.ndarray
70+
Two-body term in the geminal basis
71+
72+
Notes
73+
-----
74+
Assuming that rdm2 obbey the following permutation rules:
75+
- :math:`\Gamma_{p q r s}=-\Gamma_{q p r s}=-\Gamma_{p q s r}
76+
=\Gamma_{q p s r}`
77+
we can convert the two-body term to the geminal basis
78+
by the following formula:
79+
80+
.. math::
81+
82+
\Gamma_{p q}=0.5 * 4 \Gamma_{p q r s}
83+
84+
where:
85+
- :math:`\Gamma_{p q}` is the two-body term in the geminal basis
86+
- :math:`\Gamma_{p q r s}` is the two-body term in the spin-orbital
87+
Hamiltonian in the geminal basis is obtained by the following formula:
88+
89+
.. math::
90+
91+
V_{A B}
92+
=\frac{1}{2}\left(V_{p q r s}-V_{q p r s}-V_{p q r s}+V_{qprs}\right)
93+
94+
"""
95+
n = two_body.shape[0]
96+
two_body_gem = []
97+
98+
# i,j,k,l -> pqrs
99+
for p in range(n):
100+
for q in range(p + 1, n):
101+
for r in range(n):
102+
for s in range(r + 1, n):
103+
if type == 'rdm2':
104+
two_body_gem.append(
105+
0.5 * 4 * two_body[p, q, r, s]
106+
)
107+
elif type == 'h2':
108+
two_body_gem.append(
109+
0.5 * (
110+
two_body[p, q, r, s]
111+
- two_body[q, p, r, s]
112+
- two_body[p, q, s, r]
113+
+ two_body[q, p, s, r]
114+
)
115+
)
116+
117+
n_gem = n * (n - 1) // 2
118+
return np.array(two_body_gem).reshape(n_gem, n_gem)
119+
31120
def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
32121
r"""Convert the one/two body terms from spatial to spin-orbital basis.
33122
34123
Parameters
35124
----------
36-
None
125+
one_body : np.ndarray
126+
One-body term in spatial orbital basis (shape (n, n)).
127+
two_body : np.ndarray
128+
Two-body term in spatial orbital basis (shape (n, n, n, n)).
37129
38130
Returns
39131
-------
@@ -56,8 +148,7 @@ def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
56148
V_{pqrs}^{\\text{spatial}}`
57149
58150
"""
59-
one_body = self.one_body
60-
two_body = self.two_body
151+
one_body, two_body = self.one_body, self.two_body
61152
one_body = np.asarray(one_body)
62153
two_body = np.asarray(two_body)
63154

@@ -85,85 +176,39 @@ def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
85176

86177
return one_body_spin, two_body_spin
87178

88-
def antisymmetrize(self):
89-
"""Apply proper antisymmetrization to two-electron integrals."""
90-
self.two_body = antisymmetrize_two_body(self.two_body, inplace=True)
91-
92-
def get_spin_blocks(self):
93-
"""Return the main spin blocks of the two-body spin-orbital tensor.
94-
95-
Returns
96-
-------
97-
dict
98-
Dictionary with spin blocks: 'aaaa', 'bbbb', 'abab'
99-
100-
"""
101-
if not hasattr(self, "two_body_spin"):
102-
raise RuntimeError(
103-
"Call .spinize() first to compute spin-orbital form.")
179+
def to_reduced(self, n_elec):
180+
r"""Return the reduced 4-index Hamiltonian tensor.
104181
105-
return get_spin_blocks(self.two_body_spin, self.n_spatial)
182+
.. math::
106183
107-
def to_geminal(two_body):
108-
r"""Convert the two-body term to the geminal basis.
184+
k_{pqrs}=\frac{1}{2\,(N-1)}\bigl(h_{pq}\,\delta_{rs}
185+
+h_{rs}\,\delta_{pq}\bigr)+\tfrac12\,V_{pqrs}.
109186
110187
Parameters
111188
----------
112-
two_body : np.ndarray
113-
Two-body term in spin-orbital basis in physics notation.
189+
one_body_spatial : ndarray, shape (n, n)
190+
One‑electron integral matrix :math:`h_{pq}` in a spatial‑orbital
191+
basis.
192+
two_body_spatial : ndarray, shape (n, n, n, n)
193+
Two‑electron integral tensor :math:`V_{pqrs}` (Dirac convention,
194+
spatial orbitals).
195+
n_elec : int
196+
Number of electrons *N* in the system.
114197
115198
Returns
116199
-------
117-
two_body_gem : np.ndarray
118-
Two-body term in the geminal basis
200+
ndarray, shape (2n, 2n, 2n, 2n)
201+
Reduced Hamiltonian tensor :math:`k_{pqrs}` in the
202+
spin‑orbital basis.
119203
120204
Notes
121205
-----
122-
Assuming that rdm2 obbey the following permutation rules:
123-
- :math:`\Gamma_{p q r s}=-\Gamma_{q p r s}=-\Gamma_{p q s r}
124-
=\Gamma_{q p s r}`
125-
we can convert the two-body term to the geminal basis
126-
by the following formula:
127-
128-
.. math::
129-
130-
\Gamma_{p q}=0.5 * 4 \Gamma_{p q r s}
131-
132-
where:
133-
- :math:`\Gamma_{p q}` is the two-body term in the geminal basis
134-
- :math:`\Gamma_{p q r s}` is the two-body term in the spin-orbital
135-
Hamiltonian in the geminal basis is obtained by the following formula:
136-
137-
.. math::
138-
139-
V_{A B}
140-
=\frac{1}{2}\left(V_{p q r s}-V_{q p r s}-V_{p q r s}+V_{qprs}\right)
141-
142-
This is coming from the fact, that the Hamiltonian object
143-
retured from the fcidump (converted to the physics notation)
144-
doesn't obbey the permutation rules.
145-
Thus, the full formula has to be used.
146-
206+
The function is stateless; it does not modify the parent
207+
``MolHam`` instance.
147208
"""
148-
n = two_body.shape[0]
149-
two_body_gem = []
209+
h_spin, V_spin = self.spinize_H()
150210

151-
# i,j,k,l -> pqrs
152-
for p in range(n):
153-
for q in range(p + 1, n):
154-
for r in range(n):
155-
for s in range(r + 1, n):
156-
term = 0.5 * (
157-
two_body[p, q, r, s]
158-
- two_body[q, p, r, s]
159-
- two_body[p, q, s, r]
160-
+ two_body[q, p, s, r]
161-
)
162-
two_body_gem.append(term)
163-
164-
n_gem = n * (n - 1) // 2
165-
return np.array(two_body_gem).reshape(n_gem, n_gem)
211+
h_upscaled = upscale_one_body(h_spin, n_elec)
166212

167-
def build_reduced(self):
168-
"""Build the reduced hamiltonian form."""
169-
pass
213+
k = h_upscaled + 0.5 * V_spin
214+
return k

0 commit comments

Comments
 (0)