Skip to content

Commit 924f56e

Browse files
author
Giovanni Benedetti Da Rosa
committed
#171 Reduded Hamiltonian and Energy Testing
1 parent 68d5223 commit 924f56e

File tree

5 files changed

+101
-229
lines changed

5 files changed

+101
-229
lines changed

‎examples/molkit.ipynb‎

Lines changed: 41 additions & 183 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@
3232
},
3333
{
3434
"cell_type": "code",
35-
"execution_count": 1,
35+
"execution_count": 2,
3636
"metadata": {},
3737
"outputs": [],
3838
"source": [
@@ -109,181 +109,22 @@
109109
},
110110
{
111111
"cell_type": "code",
112-
"execution_count": null,
112+
"execution_count": 3,
113113
"metadata": {},
114-
"outputs": [],
115-
"source": [
116-
"import pyscf\n",
117-
"from pyscf import gto, scf, fci\n",
118-
"from pyscf.tools import fcidump\n",
119-
"import os\n",
120-
"\n",
121-
"\n",
122-
"path_xyz = r'../data/hydrogenic/H4/Paldus/h4_0.456_2.xyz'\n",
123-
"mol = gto.M(atom=path_xyz, basis='sto-6g')\n",
124-
"hf = scf.RHF(mol)\n",
125-
"\n",
126-
"hf.conv_tol = 1.0e-12\n",
127-
"status = hf.run().converged\n",
128-
"print(status)\n",
129-
"\n",
130-
"\n",
131-
"mo_energy = hf.mo_energy\n",
132-
"mo_occ = hf.mo_occ\n",
133-
"\n",
134-
"print('starting writing')\n",
135-
"# Write FCIDUMP\n",
136-
"\n",
137-
"fcidump.from_mo(mol, 'tmp1.fcidump',\n",
138-
" hf.mo_coeff)\n",
139-
"\n",
140-
"# fcidump.from_scf(hf, \n",
141-
"# f\"tmp.fcidump\",\n",
142-
"# tol=0.0)\n",
143-
"\n",
144-
"print('done')\n",
145-
"\n",
146-
"data = fcidump.read(\"tmp1.fcidump\")\n",
147-
"e, ci_vec = fci.direct_spin1.kernel(data['H1'],\n",
148-
" data['H2'], \n",
149-
" data['NORB'], \n",
150-
" data['NELEC'], ecore=data['ECORE'], nroots=1)\n",
151-
"\n",
152-
"rdm1s, rdm2s = fci.direct_spin1.make_rdm12s(ci_vec, data['NORB'], data['NELEC'])\n",
153-
"rdm1 = np.block([[rdm1s[0], np.zeros((4, 4))],\n",
154-
" [np.zeros((4, 4)), rdm1s[1]]])\n",
155-
"\n",
156-
"rdm2 = np.zeros((8, 8, 8, 8))\n",
157-
"for i in range(4):\n",
158-
" for j in range(4):\n",
159-
" for k in range(4):\n",
160-
" for l in range(4):\n",
161-
" rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
162-
" rdm2[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l]\n",
163-
" rdm2[i, j+4, k, l+4] = rdm2s[1][i, k, j, l]\n",
164-
" # rdm2[i+4, j, k+4, l] = rdm2s[1][k, i, l, j]\n",
165-
" rdm2[i+4, j, k+4, l] = rdm2s[1][i, k, j, l]\n",
166-
" \n",
167-
"\n",
168-
"# converting rdm2 to physics notation\n",
169-
"# rdm2 = np.einsum('ijkl->ikjl', rdm2)\n",
170-
"print(\"trace of rdm1\", np.einsum('ii', rdm1))\n",
171-
"print(\"Trace of rdm2\", np.einsum('ijij', rdm2))\n",
172-
"\n",
173-
"# loading integrals\n",
174-
"# iodata stores integrals in physics notations\n",
175-
"# https://iodata.readthedocs.io/en/latest/_modules/iodata/formats/fcidump.html#load_one\n",
176-
"integrals = load_one(open('tmp1.fcidump', 'r'))\n",
177-
"os.remove('tmp1.fcidump')\n",
178-
"\n",
179-
"h1s = integrals['one_ints']['core_mo']\n",
180-
"h2s = integrals['two_ints']['two_mo']\n",
181-
"\n",
182-
"h1 = np.block([[h1s, np.zeros((4, 4))],\n",
183-
" [np.zeros((4, 4)), h1s]])\n",
184-
"\n",
185-
"h2 = np.zeros((8, 8, 8, 8))\n",
186-
"for i in range(4):\n",
187-
" for j in range(4):\n",
188-
" for k in range(4):\n",
189-
" for l in range(4):\n",
190-
" # aaaa\n",
191-
" h2[i, j, k, l] = h2s[i, j, k, l]\n",
192-
" # bbbb\n",
193-
" h2[i+4, j+4, k+4, l+4] = h2s[i, j, k, l]\n",
194-
" # abab\n",
195-
" h2[i, j+4, k, l+4] = h2s[i, j, k, l]\n",
196-
" # baba\n",
197-
"# h2[i+4, j, k+4, l] = h2s[k, l, i, j]\n",
198-
" h2[i+4, j, k+4, l] = h2s[i, j, k, l]\n",
199-
" \n",
200-
" \n",
201-
"res = 0\n",
202-
"for i in range(8):\n",
203-
" for j in range(8):\n",
204-
" assert rdm2[i, j, i, j] >=0\n",
205-
" res += rdm2[i, j, i, j]\n",
206-
" \n",
207-
"print(\"trace of 2dm with for loop: \", res)\n",
208-
"\n",
209-
"e_dm = np.einsum('pq, pq', h1, rdm1) + 0.5 * np.einsum('pqrs, pqrs', h2, rdm2)\n",
210-
"\n",
211-
"\n",
212-
"print(e_dm, e - integrals['core_energy'])\n",
213-
"\n",
214-
"\n",
215-
"#-----------\n",
216-
"eye = np.eye(8)\n",
217-
"N = 4\n",
218-
"N_alpha = N//2\n",
219-
"\n",
220-
"h1_ups = 0.5 * (np.kron(h1, eye) + np.kron(eye, h1)) / (N-1)\n",
221-
"h1_ups = h1_ups.reshape(8, 8, 8, 8)\n",
222-
"H = h1_ups + 0.5 * h2\n",
223-
"\n",
224-
"np.einsum('pqrs, pqrs', H, rdm2), e_dm\n",
225-
"\n",
226-
"# to geminal:\n",
227-
"# H_gem = to_geminal(H)\n",
228-
"# rdm_gem = to_geminal(rdm2)\n",
229-
"# np.einsum('pq, pq', H_gem, rdm_gem) == np.einsum('pqrs, pqrs', H, rdm2) == e_dm"
230-
]
231-
},
232-
{
233-
"cell_type": "code",
234-
"execution_count": null,
235-
"metadata": {},
236-
"outputs": [],
237-
"source": [
238-
"import numpy as np\n",
239-
"from moha.molkit.hamiltonians import MolHam \n",
240-
"h_sp = np.array([[1.0, 0.2],\n",
241-
" [0.2, 1.5]])\n",
242-
"\n",
243-
"V_sp = np.zeros((2, 2, 2, 2))\n",
244-
"V_sp[0, 0, 0, 0] = 0.7\n",
245-
"V_sp[1, 1, 1, 1] = 0.6\n",
246-
"\n",
247-
"N_ELEC = 2\n",
248-
"\n",
249-
"toy_ham = MolHam(one_body=h_sp, two_body=V_sp)\n",
250-
"\n",
251-
"h_spin, V_spin = toy_ham.spinize_H() \n",
252-
"n_spin = toy_ham.n_spin \n",
253-
"\n",
254-
"e_fci, ci_vec = fci.direct_spin1.kernel(\n",
255-
" h_spin, V_spin, n_spin, N_ELEC, ecore=0.0, nroots=1\n",
256-
")\n",
257-
"\n",
258-
"rdm1s, rdm2s = fci.direct_spin1.make_rdm12(ci_vec, n_spin, N_ELEC)\n",
259-
"rdm1 = np.block([[rdm1s[0], np.zeros((4, 4))],\n",
260-
" [np.zeros((4, 4)), rdm1s[1]]])\n",
261-
"\n",
262-
"rdm2 = np.zeros((8, 8, 8, 8))\n",
263-
"for i in range(4):\n",
264-
" for j in range(4):\n",
265-
" for k in range(4):\n",
266-
" for l in range(4):\n",
267-
" rdm2[i, j, k, l] = rdm2s[0][i, k, j, l]\n",
268-
" rdm2[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l]\n",
269-
" rdm2[i, j+4, k, l+4] = rdm2s[1][i, k, j, l]\n",
270-
" # rdm2[i+4, j, k+4, l] = rdm2s[1][k, i, l, j]\n",
271-
" rdm2[i+4, j, k+4, l] = rdm2s[1][i, k, j, l]"
272-
]
273-
},
274-
{
275-
"cell_type": "code",
276-
"execution_count": null,
277-
"metadata": {},
278-
"outputs": [],
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+
],
279124
"source": [
280125
"import numpy as np\n",
281126
"from pyscf import fci\n",
282-
"\n",
283127
"from moha.molkit.hamiltonians import MolHam\n",
284-
"from moha.utils.spinops import (\n",
285-
" to_geminal_rdm2\n",
286-
")\n",
287128
"\n",
288129
"\n",
289130
"h_sp = np.array([[1.0, 0.2],\n",
@@ -293,33 +134,50 @@
293134
"V_sp[0, 0, 0, 0] = 0.7\n",
294135
"V_sp[1, 1, 1, 1] = 0.6\n",
295136
"\n",
296-
"N_ELEC = 2 \n",
137+
"N_ELEC, n_spin = 2, 2 \n",
297138
"toy = MolHam(h_sp, V_sp)\n",
298139
"\n",
299140
"h_spin, V_spin = toy.spinize_H() \n",
300-
"\n",
301-
"e_fci, ci_vec = fci.direct_spin1.kernel(h_spin, V_spin, n_spin, N_ELEC, ecore=0.0)\n",
302-
"rdm1, rdm2 = fci.direct_spin1.make_rdm12(ci_vec, n_spin, N_ELEC)#is this correct?\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",
303164
"\n",
304165
"H = toy.to_reduced(n_elec=N_ELEC) # H = h↑ + 0.5 * V\n",
305-
"\n",
306166
"E_spin = np.einsum('pqrs,pqrs', H, rdm2)\n",
307167
"\n",
308-
"##teste\n",
309-
"H_gem = toy.to_geminal(H)\n",
310-
"rdm_gem = to_geminal_rdm2(rdm2)\n",
168+
"H_gem = toy.to_geminal(H, type='h2')\n",
169+
"rdm_gem = toy.to_geminal(rdm2, type='rdm2')\n",
311170
"E_gem = np.einsum('pq,pq', H_gem, rdm_gem)\n",
312171
"\n",
313172
"print(\"E_spin-orbital :\", E_spin)\n",
314173
"print(\"E_geminal :\", E_gem)\n",
315-
"\n",
316-
"assert np.allclose(E_spin, E_gem)"
174+
"assert np.allclose(E_spin, E_gem, atol=1e-10), \"Energies do not match!\""
317175
]
318176
}
319177
],
320178
"metadata": {
321179
"kernelspec": {
322-
"display_name": "base",
180+
"display_name": "moha",
323181
"language": "python",
324182
"name": "python3"
325183
},
@@ -333,7 +191,7 @@
333191
"name": "python",
334192
"nbconvert_exporter": "python",
335193
"pygments_lexer": "ipython3",
336-
"version": "3.12.3"
194+
"version": "3.11.13"
337195
}
338196
},
339197
"nbformat": 4,

‎moha/molkit/hamiltonians.py‎

Lines changed: 22 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -32,8 +32,6 @@ def __init__(self, one_body, two_body):
3232
self.n_spin = 2 * self.n_spatial
3333
self.reduced_ham = None
3434

35-
36-
3735
def antisymmetrize(self):
3836
"""Apply proper antisymmetrization to two-electron integrals."""
3937
self.two_body = antisymmetrize_two_body(self.two_body, inplace=True)
@@ -53,13 +51,18 @@ def get_spin_blocks(self):
5351

5452
return get_spin_blocks(self.two_body_spin, self.n_spatial)
5553

56-
def to_geminal(two_body):
57-
r"""Convert the two-body term to the geminal basis.
54+
def to_geminal(self, two_body=None, type='h2'):
55+
r"""
56+
Convert the two-body term to the geminal basis.
5857
5958
Parameters
6059
----------
6160
two_body : np.ndarray
6261
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
6366
6467
Returns
6568
-------
@@ -88,11 +91,6 @@ def to_geminal(two_body):
8891
V_{A B}
8992
=\frac{1}{2}\left(V_{p q r s}-V_{q p r s}-V_{p q r s}+V_{qprs}\right)
9093
91-
This is coming from the fact, that the Hamiltonian object
92-
retured from the fcidump (converted to the physics notation)
93-
doesn't obbey the permutation rules.
94-
Thus, the full formula has to be used.
95-
9694
"""
9795
n = two_body.shape[0]
9896
two_body_gem = []
@@ -102,13 +100,19 @@ def to_geminal(two_body):
102100
for q in range(p + 1, n):
103101
for r in range(n):
104102
for s in range(r + 1, n):
105-
term = 0.5 * (
106-
two_body[p, q, r, s]
107-
- two_body[q, p, r, s]
108-
- two_body[p, q, s, r]
109-
+ two_body[q, p, s, r]
110-
)
111-
two_body_gem.append(term)
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+
)
112116

113117
n_gem = n * (n - 1) // 2
114118
return np.array(two_body_gem).reshape(n_gem, n_gem)
@@ -160,7 +164,7 @@ def spinize_H(self) -> tuple[np.ndarray, np.ndarray]:
160164
one_body_spin[n:, n:] = one_body # ββ block
161165

162166
two_body_spin = np.zeros((2 * n, 2 * n, 2 * n, 2 * n),
163-
dtype=two_body.dtype)
167+
dtype=two_body.dtype)
164168
# αααα
165169
two_body_spin[:n, :n, :n, :n] = two_body
166170
# ββββ
@@ -203,9 +207,8 @@ def to_reduced(self, n_elec):
203207
``MolHam`` instance.
204208
"""
205209
h_spin, V_spin = self.spinize_H()
206-
210+
207211
h_upscaled = upscale_one_body(h_spin, n_elec)
208212

209213
k = h_upscaled + 0.5 * V_spin
210214
return k
211-

0 commit comments

Comments
 (0)