CoCalc -- hamiltonian.py

GitHub Repository: NVIDIA/cuda-q-academic
Path: blob/main/chemistry-simulations/aux_files/krylov/qchem/hamiltonian.py
⁵⁸³ views
1
import numpy as np
2
import itertools
3

4
from cudaq import spin
5

6

7
############################################################
8
def generate_molecular_spin_ham_restricted(h1e, h2e, ecore):
9

10
    # This function generates the molecular spin Hamiltonian
11
    # H= E_core+sum_{`pq`}  h_{`pq`} a_p^dagger a_q +
12
    #                          0.5 * h_{`pqrs`} a_p^dagger a_q^dagger a_r a_s
13
    # h1e: one body integrals h_{`pq`}
14
    # h2e: two body integrals h_{`pqrs`}
15
    # `ecore`: constant (nuclear repulsion or core energy in the active space Hamiltonian)
16

17
    # Total number of qubits equals the number of spin molecular orbitals
18
    nqubits = 2 * h1e.shape[0]
19

20
    # Initialization
21
    one_body_coeff = np.zeros((nqubits, nqubits))
22
    two_body_coeff = np.zeros((nqubits, nqubits, nqubits, nqubits))
23

24
    for p in range(nqubits // 2):
25
        for q in range(nqubits // 2):
26

27
            # p & q have the same spin <a|a>= <b|b>=1
28
            # <a|b>=<b|a>=0 (orthogonal)
29
            one_body_coeff[2 * p, 2 * q] = h1e[p, q]
30
            one_body_coeff[2 * p + 1, 2 * q + 1] = h1e[p, q]
31

32
            for r in range(nqubits // 2):
33
                for s in range(nqubits // 2):
34

35
                    # Same spin (`aaaa`, `bbbbb`) <a|a><a|a>, <b|b><b|b>
36
                    two_body_coeff[2 * p, 2 * q, 2 * r,
37
                                   2 * s] = 0.5 * h2e[p, q, r, s]
38
                    two_body_coeff[2 * p + 1, 2 * q + 1, 2 * r + 1,
39
                                   2 * s + 1] = 0.5 * h2e[p, q, r, s]
40

41
                    # Mixed spin(`abab`, `baba`) <a|a><b|b>, <b|b><a|a>
42
                    #<a|b>= 0 (orthogonal)
43
                    two_body_coeff[2 * p, 2 * q + 1, 2 * r + 1,
44
                                   2 * s] = 0.5 * h2e[p, q, r, s]
45
                    two_body_coeff[2 * p + 1, 2 * q, 2 * r,
46
                                   2 * s + 1] = 0.5 * h2e[p, q, r, s]
47

48
    return one_body_coeff, two_body_coeff, ecore
49

50

51
#######################################################################
52

53

54
def jordan_wigner_one_body(p, q, coef):
55

56
    # Diagonal term: 0.5 h_{pp} (I_p - Z_p)
57
    if p == q:
58
        spin_hamiltonian = 0.5 * coef * spin.i(p)
59
        spin_hamiltonian -= 0.5 * coef * spin.z(p)
60

61
    # h_`pq`(a_p^dagger a_q + a_q^dagger a_p) = R(h_`pq`) (a_p^dagger a_q + a_q^dagger a_p) +
62
    #                                          `imag` (h_`pq`) (a_p^dagger a_q - a_q^dagger a_p)
63
    # Off-diagonal real part: 0.5 * real(h_{`pq`}) [ X_p (Z_{p+1}^{q-1}) X_q + Y_p (Z_{p+1}^{q-1}) Y_q ]
64
    # Off-diagonal imaginary part: 0.5* `im`(h_`pq`) [y_p (Z_{p+1}^{q-1}) x_q - x_p (Z_{p+1}^{q-1}) y_q]
65

66
    else:
67
        if p > q:
68
            p, q = q, p
69
            coef = np.conj(coef)
70

71
        # Compute the parity string (Z_{p+1}^{q-1})
72
        z_indices = [i for i in range(p + 1, q)]
73
        parity_string = 1.0
74
        for i in z_indices:
75
            parity_string *= spin.z(i)
76

77
        spin_hamiltonian = 0.5 * coef.real * spin.x(p) * parity_string * spin.x(
78
            q)
79
        spin_hamiltonian += 0.5 * coef.real * spin.y(
80
            p) * parity_string * spin.y(q)
81
        spin_hamiltonian += 0.5 * coef.imag * spin.y(
82
            p) * parity_string * spin.x(q)
83
        spin_hamiltonian -= 0.5 * coef.imag * spin.x(
84
            p) * parity_string * spin.y(q)
85

86
    return spin_hamiltonian
87

88

89
#############
90

91

92
def jordan_wigner_two_body(p, q, r, s, coef):
93

94
    # Diagonal term: p=r, q=s or p=s,q=r --> (+ h_`pqpq` = + h_`qpqp` = - h_`qppq` = - h_`pqqp`)
95
    #
96
    # exchange operator:  h_`pqpq` (a_p^dagger a_q^dagger a_p a_q) + h_`qpqp` (a_q^dagger a_p^dagger a_q a_p)
97
    # p<q: -1/4 (I_p I_q - I_p Z_q - Z_p I_q+Z_p Z_q)
98
    #
99
    # coulomb operator: h_`qppq` (a_q^dagger a_p^dagger a_p a_q) + h_`pqqp` (a_p^dagger a_q^dagger a_q a_p)
100
    # p<q: 1/4 (I_p I_q - I_p Z_q - Z_p I_q + Z_p Z_q)
101

102
    if len(set([p, q, r, s])) == 2:
103

104
        if p == r:
105
            spin_hamiltonian = -0.25 * coef * spin.i(p) * spin.i(q)
106
            spin_hamiltonian += 0.25 * coef * spin.i(p) * spin.z(q)
107
            spin_hamiltonian += 0.25 * coef * spin.z(p) * spin.i(q)
108
            spin_hamiltonian -= 0.25 * coef * spin.z(p) * spin.z(q)
109

110
        elif q == r:
111
            spin_hamiltonian = 0.25 * coef * spin.i(p) * spin.i(q)
112
            spin_hamiltonian -= 0.25 * coef * spin.i(p) * spin.z(q)
113
            spin_hamiltonian -= 0.25 * coef * spin.z(p) * spin.i(q)
114
            spin_hamiltonian += 0.25 * coef * spin.z(p) * spin.z(q)
115

116
    # Off-diagonal term with three different sets of non-equal indices
117
    # Number with excitation operator
118
    # + h_`pqqs` = + h_`qpsq` = - h_`qpqs` = - h_`pqsq` and their hermitian conjugate
119
    # Real (h_`pqqs`) (a_p^dagger a_q^dagger a_q a_s + a_s^dagger a_q^dagger a_q a_p) +
120
    # `imag` (h_`pqqs`) (a_p^dagger a_q^dagger a_q a_s - a_s^dagger a_q^dagger a_q a_p)
121
    # p <q <s: (1/4)(Z_{p+1}^{s-1}) [ I_q {real (h_`pqqs`/4) (x_p x_s + y_p y_s) + {`imag` (h_`pqqs`/4) (y_p x_s - x_p y_s)}
122
    #                           - Z_q {real (h_`pqqs`/4) (x_p x_s + y_p y_s) + `imag`(h_`pqqs`) (y_p x_s -x_p y_s)}]
123

124
    if len(set([p, q, r, s])) == 3:
125

126
        if q == r:
127
            if p > r:
128
                a, b = s, p
129
                coef = np.conj(coef)
130
            else:
131
                a, b = p, s
132
            c = q
133

134
        elif q == s:
135
            if p > r:
136
                a, b = r, p
137
                coef = -1.0 * np.conj(coef)
138
            else:
139
                a, b = p, r
140
                coef *= -1.0
141
            c = q
142

143
        elif p == r:
144
            if q > s:
145
                a, b = s, q
146
                coef = -1.0 * np.conj(coef)
147
            else:
148
                a, b = q, s
149
                coef = -1.0 * coef
150
            c = p
151

152
        elif p == s:
153
            if q > r:
154
                a, b = r, q
155
                coef = np.conj(coef)
156
            else:
157
                a, b = q, r
158
            c = p
159

160
        parity_string = 1.0
161
        z_qubit = [i for i in range(a + 1, b)]
162
        for i in z_qubit:
163
            parity_string *= spin.z(i)
164

165
        spin_hamiltonian = 0.25 * coef.real * spin.x(
166
            a) * parity_string * spin.x(b) * spin.i(c)
167
        spin_hamiltonian += 0.25 * coef.real * spin.y(
168
            a) * parity_string * spin.y(b) * spin.i(c)
169
        spin_hamiltonian += 0.25 * coef.imag * spin.y(
170
            a) * parity_string * spin.x(b) * spin.i(c)
171
        spin_hamiltonian -= 0.25 * coef.imag * spin.x(
172
            a) * parity_string * spin.y(b) * spin.i(c)
173

174
        spin_hamiltonian -= 0.25 * coef.real * spin.x(
175
            a) * parity_string * spin.x(b) * spin.z(c)
176
        spin_hamiltonian -= 0.25 * coef.real * spin.y(
177
            a) * parity_string * spin.y(b) * spin.z(c)
178
        spin_hamiltonian -= 0.25 * coef.imag * spin.y(
179
            a) * parity_string * spin.x(b) * spin.z(c)
180
        spin_hamiltonian += 0.25 * coef.imag * spin.x(
181
            a) * parity_string * spin.y(b) * spin.z(c)
182

183
    # Off-diagonal term with four different sets of non-equal indices
184
    # h_`pqrs` = h_`qpsr` = - h_`qprs` = - h_`pqsr`
185
    # real {h_`pqrs`} (a_p^dagger a_q^dagger a_r a_s + a_s^dagger a_r^dagger a_q a_p) +
186
    # `imag` (h_`pqrs`) (a_p^dagger a_q^dagger a_r a_s - a_s^dagger a_r^dagger a_q a_p)
187
    # p<q<r<s real part: -1/8 (Z_{p+1}^{q-1}) (Z_{r+1}^{s-1}) (x_p x_q x_r x_s - x_p x_q y_r y_s + x_p y_q x_r y_s + x_p y_q y_r x_s
188
    #                         + y_p x_q x_r y_s + y_p x_q y_r x_s - y_p y_q x_r x_s + y_p y_q y_r y_s)
189
    # p<q<r<s `imag` part: -1/8 (x_p x_q x_r y_s + x_p x_q y_r x_s - x_p y_q x_r x_s + x_p y_q y_r y_s
190
    #                         - y_p x_q x_r x_s + y_p x_q y_r y_s - y_p y_q x_r y_s - y_p y_q y_r x_s)
191
    # also we need to compute p<r<q<s and p<r<s<q
192

193
    elif len(set([p, q, r, s])) == 4:
194

195
        if (p > q) ^ (r > s):
196
            coef *= -1.0
197

198
        if p < q < r < s:
199
            a, b, c, d = p, q, r, s
200

201
            parity_string_a = 1.0
202
            z_qubit_a = [i for i in range(a + 1, b)]
203
            for i in z_qubit_a:
204
                parity_string_a *= spin.z(i)
205

206
            parity_string_b = 1.0
207
            z_qubit_b = [i for i in range(c + 1, d)]
208
            for i in z_qubit_b:
209
                parity_string_b *= spin.z(i)
210

211
            spin_hamiltonian = -0.125 * coef.real * spin.x(
212
                a) * parity_string_a * spin.x(b) * spin.x(
213
                    c) * parity_string_b * spin.x(d)
214
            spin_hamiltonian -= -0.125 * coef.real * spin.x(
215
                a) * parity_string_a * spin.x(b) * spin.y(
216
                    c) * parity_string_b * spin.y(d)
217
            spin_hamiltonian += -0.125 * coef.real * spin.x(
218
                a) * parity_string_a * spin.y(b) * spin.x(
219
                    c) * parity_string_b * spin.y(d)
220
            spin_hamiltonian += -0.125 * coef.real * spin.x(
221
                a) * parity_string_a * spin.y(b) * spin.y(
222
                    c) * parity_string_b * spin.x(d)
223
            spin_hamiltonian += -0.125 * coef.real * spin.y(
224
                a) * parity_string_a * spin.x(b) * spin.x(
225
                    c) * parity_string_b * spin.y(d)
226
            spin_hamiltonian += -0.125 * coef.real * spin.y(
227
                a) * parity_string_a * spin.x(b) * spin.y(
228
                    c) * parity_string_b * spin.x(d)
229
            spin_hamiltonian -= -0.125 * coef.real * spin.y(
230
                a) * parity_string_a * spin.y(b) * spin.x(
231
                    c) * parity_string_b * spin.x(d)
232
            spin_hamiltonian += -0.125 * coef.real * spin.y(
233
                a) * parity_string_a * spin.y(b) * spin.y(
234
                    c) * parity_string_b * spin.y(d)
235

236
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
237
                a) * parity_string_a * spin.x(b) * spin.x(
238
                    c) * parity_string_b * spin.y(d)
239
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
240
                a) * parity_string_a * spin.x(b) * spin.y(
241
                    c) * parity_string_b * spin.x(d)
242
            spin_hamiltonian -= 0.125 * coef.imag * spin.x(
243
                a) * parity_string_a * spin.y(b) * spin.x(
244
                    c) * parity_string_b * spin.x(d)
245
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
246
                a) * parity_string_a * spin.y(b) * spin.y(
247
                    c) * parity_string_b * spin.y(d)
248
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
249
                a) * parity_string_a * spin.x(b) * spin.x(
250
                    c) * parity_string_b * spin.x(d)
251
            spin_hamiltonian += 0.125 * coef.imag * spin.y(
252
                a) * parity_string_a * spin.x(b) * spin.y(
253
                    c) * parity_string_b * spin.y(d)
254
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
255
                a) * parity_string_a * spin.y(b) * spin.x(
256
                    c) * parity_string_b * spin.y(d)
257
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
258
                a) * parity_string_a * spin.y(b) * spin.y(
259
                    c) * parity_string_b * spin.x(d)
260

261
        elif p < r < q < s:
262
            a, b, c, d = p, r, q, s
263

264
            parity_string_a = 1.0
265
            z_qubit_a = [i for i in range(a + 1, b)]
266
            for i in z_qubit_a:
267
                parity_string_a *= spin.z(i)
268

269
            parity_string_b = 1.0
270
            z_qubit_b = [i for i in range(c + 1, d)]
271
            for i in z_qubit_b:
272
                parity_string_b *= spin.z(i)
273

274
            spin_hamiltonian = -0.125 * coef.real * spin.x(
275
                a) * parity_string_a * spin.x(b) * spin.x(
276
                    c) * parity_string_b * spin.x(d)
277
            spin_hamiltonian += -0.125 * coef.real * spin.x(
278
                a) * parity_string_a * spin.x(b) * spin.y(
279
                    c) * parity_string_b * spin.y(d)
280
            spin_hamiltonian -= -0.125 * coef.real * spin.x(
281
                a) * parity_string_a * spin.y(b) * spin.x(
282
                    c) * parity_string_b * spin.y(d)
283
            spin_hamiltonian += -0.125 * coef.real * spin.x(
284
                a) * parity_string_a * spin.y(b) * spin.y(
285
                    c) * parity_string_b * spin.x(d)
286
            spin_hamiltonian += -0.125 * coef.real * spin.y(
287
                a) * parity_string_a * spin.x(b) * spin.x(
288
                    c) * parity_string_b * spin.y(d)
289
            spin_hamiltonian -= -0.125 * coef.real * spin.y(
290
                a) * parity_string_a * spin.x(b) * spin.y(
291
                    c) * parity_string_b * spin.x(d)
292
            spin_hamiltonian += -0.125 * coef.real * spin.y(
293
                a) * parity_string_a * spin.y(b) * spin.x(
294
                    c) * parity_string_b * spin.x(d)
295
            spin_hamiltonian += -0.125 * coef.real * spin.y(
296
                a) * parity_string_a * spin.y(b) * spin.y(
297
                    c) * parity_string_b * spin.y(d)
298

299
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
300
                a) * parity_string_a * spin.x(b) * spin.x(
301
                    c) * parity_string_b * spin.y(d)
302
            spin_hamiltonian -= 0.125 * coef.imag * spin.x(
303
                a) * parity_string_a * spin.x(b) * spin.y(
304
                    c) * parity_string_b * spin.x(d)
305
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
306
                a) * parity_string_a * spin.y(b) * spin.x(
307
                    c) * parity_string_b * spin.x(d)
308
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
309
                a) * parity_string_a * spin.y(b) * spin.y(
310
                    c) * parity_string_b * spin.y(d)
311
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
312
                a) * parity_string_a * spin.x(b) * spin.x(
313
                    c) * parity_string_b * spin.x(d)
314
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
315
                a) * parity_string_a * spin.x(b) * spin.y(
316
                    c) * parity_string_b * spin.y(d)
317
            spin_hamiltonian += 0.125 * coef.imag * spin.y(
318
                a) * parity_string_a * spin.y(b) * spin.x(
319
                    c) * parity_string_b * spin.y(d)
320
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
321
                a) * parity_string_a * spin.y(b) * spin.y(
322
                    c) * parity_string_b * spin.x(d)
323

324
        elif p < r < s < q:
325
            a, b, c, d = p, r, s, q
326

327
            parity_string_a = 1.0
328
            z_qubit_a = [i for i in range(a + 1, b)]
329
            for i in z_qubit_a:
330
                parity_string_a *= spin.z(i)
331

332
            parity_string_b = 1.0
333
            z_qubit_b = [i for i in range(c + 1, d)]
334
            for i in z_qubit_b:
335
                parity_string_b *= spin.z(i)
336

337
            spin_hamiltonian = -0.125 * coef.real * spin.x(
338
                a) * parity_string_a * spin.x(b) * spin.x(
339
                    c) * parity_string_b * spin.x(d)
340
            spin_hamiltonian += -0.125 * coef.real * spin.x(
341
                a) * parity_string_a * spin.x(b) * spin.y(
342
                    c) * parity_string_b * spin.y(d)
343
            spin_hamiltonian += -0.125 * coef.real * spin.x(
344
                a) * parity_string_a * spin.y(b) * spin.x(
345
                    c) * parity_string_b * spin.y(d)
346
            spin_hamiltonian -= -0.125 * coef.real * spin.x(
347
                a) * parity_string_a * spin.y(b) * spin.y(
348
                    c) * parity_string_b * spin.x(d)
349
            spin_hamiltonian -= -0.125 * coef.real * spin.y(
350
                a) * parity_string_a * spin.x(b) * spin.x(
351
                    c) * parity_string_b * spin.y(d)
352
            spin_hamiltonian += -0.125 * coef.real * spin.y(
353
                a) * parity_string_a * spin.x(b) * spin.y(
354
                    c) * parity_string_b * spin.x(d)
355
            spin_hamiltonian += -0.125 * coef.real * spin.y(
356
                a) * parity_string_a * spin.y(b) * spin.x(
357
                    c) * parity_string_b * spin.x(d)
358
            spin_hamiltonian += -0.125 * coef.real * spin.y(
359
                a) * parity_string_a * spin.y(b) * spin.y(
360
                    c) * parity_string_b * spin.y(d)
361

362
            spin_hamiltonian -= 0.125 * coef.imag * spin.x(
363
                a) * parity_string_a * spin.x(b) * spin.x(
364
                    c) * parity_string_b * spin.y(d)
365
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
366
                a) * parity_string_a * spin.x(b) * spin.y(
367
                    c) * parity_string_b * spin.x(d)
368
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
369
                a) * parity_string_a * spin.y(b) * spin.x(
370
                    c) * parity_string_b * spin.x(d)
371
            spin_hamiltonian += 0.125 * coef.imag * spin.x(
372
                a) * parity_string_a * spin.y(b) * spin.y(
373
                    c) * parity_string_b * spin.y(d)
374
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
375
                a) * parity_string_a * spin.x(b) * spin.x(
376
                    c) * parity_string_b * spin.x(d)
377
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
378
                a) * parity_string_a * spin.x(b) * spin.y(
379
                    c) * parity_string_b * spin.y(d)
380
            spin_hamiltonian -= 0.125 * coef.imag * spin.y(
381
                a) * parity_string_a * spin.y(b) * spin.x(
382
                    c) * parity_string_b * spin.y(d)
383
            spin_hamiltonian += 0.125 * coef.imag * spin.y(
384
                a) * parity_string_a * spin.y(b) * spin.y(
385
                    c) * parity_string_b * spin.x(d)
386

387
    return spin_hamiltonian
388

389

390
##########################################################
391
def jordan_wigner_fermion(h_pq, h_pqrs, ecore, tolerance=1e-12):
392

393
    # Compute the qubit `hamiltonian` using `jordan` `wigner`
394
    # one-body and two-body integrals could be real or complex.
395
    #
396
    # For one-body integrals: there are two terms (diagonal term (number operator) h_pp and
397
    #                                    off-diagonal term (excitation operator) h_`pq`)
398
    # H_1= sum_`pq` h_{`pq`} a_p^dagger a_q + h.c.
399
    #
400
    # For two-body integrals: There are three different terms
401
    #                 (diagonal term (coulomb and exchange operator):h_`pqqp`, h_`pqpq`,
402
    #                   off diagonal terms: (number with excitation operator)h_`pqqr`, h_`pqrp`,
403
    #                   and (double excitation operator) h_`pqrs`
404
    # H_2 = sum_{`pqrs`} h_`pqrs` a_p^dagger a_q^dagger a_r a_s + h.c.
405
    #
406
    # Jordan Wigner transformation
407
    # a_j^dagger = (Z_{k=1} ^{j-1}) [0.5(X_j - i Y_j)]
408
    # a_j = (Z_{k=1} ^{j-1}) [0.5(X_j + i Y_j)]
409
    #
410
    # Some combination of indices are not allowed because of the `pauli` principles and
411
    # the anti-commutation relations of the creation and annihilation operators.
412

413
    spin_hamiltonian = ecore
414

415
    nqubit = h_pq.shape[0]
416

417
    for p in range(nqubit):
418

419
        # Diagonal one-body term (number operator): sum_p h_{pp} a_p^dagger a_p
420
        coef = h_pq[p, p]
421
        if np.abs(coef) > tolerance:
422
            spin_hamiltonian += jordan_wigner_one_body(p, p, coef)
423

424
    for p, q in itertools.combinations(range(nqubit), 2):
425

426
        # Off-diagonal one-body term (excitation operator): sum_p<q (h_{`pq`} a_p^dagger a_q + h_{`qp`} a_q^dagger a_p)
427
        coef = 0.5 * (h_pq[p, q] + np.conj(h_pq[q, p]))
428
        if np.abs(coef) > tolerance:
429
            spin_hamiltonian += jordan_wigner_one_body(p, q, coef)
430

431
        # Diagonal term two-body (coulomb and exchange operators)
432
        # Diagonal term: p=r, q=s or p=s,q=r --> (+ h_`pqpq` = + h_`qpqp` = - h_`qppq` = - h_`pqqp`)
433

434
        # exchange operator
435
        coef = h_pqrs[p, q, p, q] + h_pqrs[q, p, q, p]
436
        if np.abs(coef) > tolerance:
437
            spin_hamiltonian += jordan_wigner_two_body(p, q, p, q, coef)
438

439
        # coulomb operator
440
        coef = h_pqrs[p, q, q, p] + h_pqrs[q, p, p, q]
441
        if np.abs(coef) > tolerance:
442
            spin_hamiltonian += jordan_wigner_two_body(p, q, q, p, coef)
443

444
    for (p, q), (r, s) in itertools.combinations(
445
            itertools.combinations(range(nqubit), 2), 2):
446

447
        # h_`pqrs` = - h_`qprs` = -h_`pqsr` = h_`qpsr`
448
        # Four point symmetry if integrals are complex: `pqrs` = `srqp` = `qpsr` = `rspq`
449
        # Eight point symmetry if integrals are real: `pqrs` = `rqps` = `psrq` = `srqp` = `qpsr` = `rspq` = `spqr` = `qrsp`
450

451

452
        coef=0.5*(h_pqrs[p,q,r,s] + np.conj(h_pqrs[s,r,q,p]) - h_pqrs[q,p,r,s] - np.conj(h_pqrs[s,r,p,q]) \
453
             - h_pqrs[p,q,s,r] - np.conj(h_pqrs[r,s,q,p]) + h_pqrs[q,p,s,r] + np.conj(h_pqrs[r,s,p,q]))
454

455
        # Compute number with excitation operator and double excitation operator
456
        if np.abs(coef) > tolerance:
457
            spin_hamiltonian += jordan_wigner_two_body(p, q, r, s, coef)
458

459
    # Remove term with zero coefficient.
460
    spin_hamiltonian = spin_hamiltonian.canonicalize().trim(tolerance)
461

462
    return spin_hamiltonian
463

464

465
############
466

467
Product

Resources

Company
