Skip to content

Commit 11fc3b6

Browse files
refraction-rayrefraction-ray
committed
add 2-otoc support for fgs
1 parent edc7da7 commit 11fc3b6

File tree

2 files changed

+189
-6
lines changed

2 files changed

+189
-6
lines changed

‎tensorcircuit/fgs.py

+170-6
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
"""
22
Fermion Gaussian state simulator
33
"""
4-
from typing import Any, List, Optional, Tuple
4+
from typing import Any, Dict, List, Optional, Tuple
55

66
import numpy as np
77

@@ -74,8 +74,10 @@ def __init__(
7474
_, _, self.alpha = self.fermion_diagonalization(hc, L)
7575
else:
7676
self.alpha = alpha
77+
self.alpha0 = self.alpha
7778
self.wtransform = self.wmatrix(L)
7879
self.cmatrix = cmatrix
80+
self.otcmatrix: Dict[Tuple[int, int], Tensor] = {}
7981

8082
@classmethod
8183
def fermion_diagonalization(
@@ -141,13 +143,35 @@ def init_alpha(filled: List[int], L: int) -> Tensor:
141143
def get_alpha(self) -> Tensor:
142144
return self.alpha
143145

144-
def get_cmatrix(self) -> Tensor:
145-
if self.cmatrix is not None:
146+
def get_cmatrix(self, now_i: bool = True, now_j: bool = True) -> Tensor:
147+
# otoc in FGS language, see: https://arxiv.org/pdf/1908.03292.pdf
148+
# https://journals.aps.org/prb/pdf/10.1103/PhysRevB.99.054205
149+
# otoc for non=hermitian system is more subtle due to the undefined normalization
150+
# of operator and not considered here, see: https://arxiv.org/pdf/2305.12054.pdf
151+
if self.cmatrix is not None and now_i is True and now_j is True:
146152
return self.cmatrix
147-
else:
153+
elif now_i is True and now_j is True:
148154
cmatrix = self.alpha @ backend.adjoint(self.alpha)
149155
self.cmatrix = cmatrix
150156
return cmatrix
157+
elif now_i is True: # new i old j
158+
if (1, 0) in self.otcmatrix:
159+
return self.otcmatrix[(1, 0)]
160+
cmatrix = self.alpha @ backend.adjoint(self.alpha0)
161+
self.otcmatrix[(1, 0)] = cmatrix
162+
return cmatrix
163+
elif now_j is True: # old i new j
164+
if (0, 1) in self.otcmatrix:
165+
return self.otcmatrix[(0, 1)]
166+
cmatrix = self.alpha0 @ backend.adjoint(self.alpha)
167+
self.otcmatrix[(0, 1)] = cmatrix
168+
return cmatrix
169+
else: # initial cmatrix
170+
if (0, 0) in self.otcmatrix:
171+
return self.otcmatrix[(0, 0)]
172+
cmatrix = self.alpha0 @ backend.adjoint(self.alpha0)
173+
self.otcmatrix[(0, 0)] = cmatrix
174+
return cmatrix
151175

152176
def get_reduced_cmatrix(self, subsystems_to_trace_out: List[int]) -> Tensor:
153177
m = self.get_cmatrix()
@@ -216,6 +240,7 @@ def evol_hamiltonian(self, h: Tensor) -> None:
216240
h = backend.cast(h, dtype=dtypestr)
217241
self.alpha = backend.expm(-1.0j * h) @ self.alpha
218242
self.cmatrix = None
243+
self.otcmatrix = {}
219244

220245
def evol_ihamiltonian(self, h: Tensor) -> None:
221246
r"""
@@ -229,6 +254,7 @@ def evol_ihamiltonian(self, h: Tensor) -> None:
229254
self.alpha = backend.expm(h) @ self.alpha
230255
self.orthogonal()
231256
self.cmatrix = None
257+
self.otcmatrix = {}
232258

233259
def evol_ghamiltonian(self, h: Tensor) -> None:
234260
r"""
@@ -242,6 +268,7 @@ def evol_ghamiltonian(self, h: Tensor) -> None:
242268
self.alpha = backend.expm(-1.0j * backend.adjoint(h)) @ self.alpha
243269
self.orthogonal()
244270
self.cmatrix = None
271+
self.otcmatrix = {}
245272

246273
def orthogonal(self) -> None:
247274
q, _ = backend.qr(self.alpha)
@@ -350,7 +377,9 @@ def get_covariance_matrix(self) -> Tensor:
350377
m = self.get_cmatrix_majorana()
351378
return -1.0j * (2 * m - backend.eye(self.L * 2))
352379

353-
def expectation_2body(self, i: int, j: int) -> Tensor:
380+
def expectation_2body(
381+
self, i: int, j: int, now_i: bool = True, now_j: bool = True
382+
) -> Tensor:
354383
"""
355384
expectation of two fermion terms
356385
convention: (c, c^\dagger)
@@ -363,7 +392,7 @@ def expectation_2body(self, i: int, j: int) -> Tensor:
363392
:return: _description_
364393
:rtype: Tensor
365394
"""
366-
return self.get_cmatrix()[i][(j + self.L) % (2 * self.L)]
395+
return self.get_cmatrix(now_i, now_j)[i][(j + self.L) % (2 * self.L)]
367396

368397
def expectation_4body(self, i: int, j: int, k: int, l: int) -> Tensor:
369398
"""
@@ -520,6 +549,7 @@ def __init__(
520549
self.state = self.fermion_diagonalization(hc, L)
521550
else:
522551
self.state = self.init_state(filled, L)
552+
self.state0 = self.state
523553

524554
@staticmethod
525555
def init_state(filled: List[int], L: int) -> Tensor:
@@ -621,6 +651,140 @@ def evol_sp(self, i: int, j: int, chi: Tensor = 0) -> None:
621651
def orthogonal(self) -> None:
622652
self.state /= backend.norm(self.state)
623653

654+
def get_ot_cmatrix(self, h: Tensor, now_i: bool = True) -> Tensor:
655+
alpha1_jw = self.state
656+
cmatrix = np.zeros([2 * self.L, 2 * self.L], dtype=complex)
657+
for i in range(self.L):
658+
for j in range(self.L):
659+
op1 = openfermion.FermionOperator(f"{str(i)}")
660+
op2 = openfermion.FermionOperator(f"{str(j)}^")
661+
662+
m1 = openfermion.get_sparse_operator(op1, n_qubits=self.L).todense()
663+
m2 = openfermion.get_sparse_operator(op2, n_qubits=self.L).todense()
664+
eh = npb.expm(-1 / 2 * 1.0j * h)
665+
eh1 = npb.expm(1 / 2 * 1.0j * h)
666+
if now_i is True:
667+
cmatrix[i, j] = backend.item(
668+
(
669+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
670+
@ eh1
671+
@ m1
672+
@ eh
673+
@ m2
674+
@ backend.reshape(alpha1_jw, [-1, 1])
675+
)[0, 0]
676+
)
677+
else:
678+
cmatrix[i, j] = backend.item(
679+
(
680+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
681+
@ m1
682+
@ eh1
683+
@ m2
684+
@ eh
685+
@ backend.reshape(alpha1_jw, [-1, 1])
686+
)[0, 0]
687+
)
688+
for i in range(self.L, 2 * self.L):
689+
for j in range(self.L):
690+
op1 = openfermion.FermionOperator(f"{str(i-self.L)}^")
691+
op2 = openfermion.FermionOperator(f"{str(j)}^")
692+
693+
m1 = openfermion.get_sparse_operator(op1, n_qubits=self.L).todense()
694+
m2 = openfermion.get_sparse_operator(op2, n_qubits=self.L).todense()
695+
eh = npb.expm(-1 / 2 * 1.0j * h)
696+
eh1 = npb.expm(1 / 2 * 1.0j * h)
697+
698+
if now_i is True:
699+
cmatrix[i, j] = backend.item(
700+
(
701+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
702+
@ eh1
703+
@ m1
704+
@ eh
705+
@ m2
706+
@ backend.reshape(alpha1_jw, [-1, 1])
707+
)[0, 0]
708+
)
709+
else:
710+
cmatrix[i, j] = backend.item(
711+
(
712+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
713+
@ m1
714+
@ eh1
715+
@ m2
716+
@ eh
717+
@ backend.reshape(alpha1_jw, [-1, 1])
718+
)[0, 0]
719+
)
720+
721+
for i in range(self.L):
722+
for j in range(self.L, 2 * self.L):
723+
op1 = openfermion.FermionOperator(f"{str(i)}")
724+
op2 = openfermion.FermionOperator(f"{str(j-self.L)}")
725+
726+
m1 = openfermion.get_sparse_operator(op1, n_qubits=self.L).todense()
727+
m2 = openfermion.get_sparse_operator(op2, n_qubits=self.L).todense()
728+
eh = npb.expm(-1 / 2 * 1.0j * h)
729+
eh1 = npb.expm(1 / 2 * 1.0j * h)
730+
731+
if now_i is True:
732+
cmatrix[i, j] = backend.item(
733+
(
734+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
735+
@ eh1
736+
@ m1
737+
@ eh
738+
@ m2
739+
@ backend.reshape(alpha1_jw, [-1, 1])
740+
)[0, 0]
741+
)
742+
else:
743+
cmatrix[i, j] = backend.item(
744+
(
745+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
746+
@ m1
747+
@ eh1
748+
@ m2
749+
@ eh
750+
@ backend.reshape(alpha1_jw, [-1, 1])
751+
)[0, 0]
752+
)
753+
for i in range(self.L, 2 * self.L):
754+
for j in range(self.L, 2 * self.L):
755+
op1 = openfermion.FermionOperator(f"{str(i-self.L)}^ ")
756+
op2 = openfermion.FermionOperator(f"{str(j-self.L)}")
757+
758+
m1 = openfermion.get_sparse_operator(op1, n_qubits=self.L).todense()
759+
m2 = openfermion.get_sparse_operator(op2, n_qubits=self.L).todense()
760+
eh = npb.expm(-1 / 2 * 1.0j * h)
761+
eh1 = npb.expm(1 / 2 * 1.0j * h)
762+
763+
if now_i is True:
764+
cmatrix[i, j] = backend.item(
765+
(
766+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
767+
@ eh1
768+
@ m1
769+
@ eh
770+
@ m2
771+
@ backend.reshape(alpha1_jw, [-1, 1])
772+
)[0, 0]
773+
)
774+
else:
775+
cmatrix[i, j] = backend.item(
776+
(
777+
backend.reshape(backend.adjoint(alpha1_jw), [1, -1])
778+
@ m1
779+
@ eh1
780+
@ m2
781+
@ eh
782+
@ backend.reshape(alpha1_jw, [-1, 1])
783+
)[0, 0]
784+
)
785+
786+
return cmatrix
787+
624788
def get_cmatrix(self) -> Tensor:
625789
alpha1_jw = self.state
626790
cmatrix = np.zeros([2 * self.L, 2 * self.L], dtype=complex)

‎tests/test_fgs.py

+19
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,25 @@ def asymmetric_hopping_jw(J, gamma, i, j, L):
5959
)
6060

6161

62+
@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")])
63+
def test_otoc(backend, highp):
64+
c = tc.FGSSimulator(4, [0, 1])
65+
c1 = tc.fgs.FGSTestSimulator(4, [0, 1])
66+
h = c.hopping(0.6, 1, 2, 4)
67+
h += c.hopping(-0.8, 3, 2, 4)
68+
h += c.chemical_potential(0.3, 2, 4)
69+
h1 = c1.hopping_jw(0.6, 1, 2, 4)
70+
h1 += c1.hopping_jw(-0.8, 3, 2, 4)
71+
h1 += c1.chemical_potential_jw(0.3, 2, 4)
72+
c.evol_hamiltonian(h)
73+
m = c.get_cmatrix(now_i=True, now_j=False)
74+
m0 = c1.get_ot_cmatrix(h1)
75+
np.testing.assert_allclose(m, m0, atol=1e-5)
76+
m = c.get_cmatrix(now_i=False, now_j=True)
77+
m0 = c1.get_ot_cmatrix(h1, now_i=False)
78+
np.testing.assert_allclose(m, m0, atol=1e-5)
79+
80+
6281
@pytest.mark.parametrize("backend", [lf("tfb"), lf("jaxb")])
6382
def test_fgs_ad(backend, highp):
6483
import optax

0 commit comments

Comments
 (0)