add tests

README.md

@@ -16,8 +16,16 @@
 - Finite system: check for ||Ax - \lambda x|| < eps -> smallest required ranks
 - Thermodynamic limits for critical exponents-> VUMPS?
 - OBSERVATION: convergence for E doesn't require large bond dimensions, but for other observables alr!
-- Result from time-propagation is strange.
-- Negative values for site expectation value.
+- Negative values for site expectation value
+- Simulations for first excited state
+- Reproduce results from 2 Banuls' papers
+- Add TN diagrams to thesis
+- Open issues on TensorKit
+
+## For discussions
+- DMRG for dissipative Ising chain: GS is not dark state
+- Result from time-propagation is strange
+
 
 ## Technical details 
 Required tools: `poetry`

src/hpc/plotting_results.py

@@ -2,10 +2,11 @@
 import numpy as np
 
 PATH = "/home/psireal42/study/quantum-contact-process-1D/hpc/results/"
-time_string = "2024_04_26_01_02_06"
+time_string = "2024_05_01_16_12_13"
 
 # system parameters
-L = 10
+L = 50
+d = 2
 GAMMA = 1
 OMEGAS = np.linspace(0, 10, 10)
 
@@ -17,7 +18,7 @@
 n_s = np.loadtxt(PATH + f"n_s_L_{L}_{time_string}.txt", delimiter=',')
 eval_0 = np.loadtxt(PATH + f"eval_0_L_{L}_{time_string}.txt", delimiter=',')
 eval_1 = np.loadtxt(PATH + f"eval_1_L_{L}_{time_string}.txt", delimiter=',')
-eval_spectrum = np.loadtxt(PATH + f"eval_spectrum_L_{L}_{time_string}.txt", delimiter=',')
+ent_ent_spectrum = np.loadtxt(PATH + f"ent_ent_spectrum_L_{L}_{time_string}.txt", delimiter=',')
 purities = np.loadtxt(PATH + f"purities_L_{L}_{time_string}.txt", delimiter=',')
 correlations = np.loadtxt(PATH + f"correlations_L_{L}_{time_string}.txt", delimiter=',')
 evp_residual = np.loadtxt(PATH + f"evp_residual_L_{L}_{time_string}.txt", delimiter=',')
@@ -57,13 +58,13 @@
 # plot eigenvalue spectrum
 plt.figure()
 for i, omega in enumerate(OMEGAS):
-    plt.plot(list(range(bond_dims[-1])), eval_spectrum[i, :], 'o-', label=f"$\Omega=${omega}")
+    plt.plot(list(range(bond_dims[-1] * d**2)), ent_ent_spectrum[i, :], 'o-', label=f"$\Omega=${omega}")
 plt.xlabel(r"$\lambda_i$")
 plt.legend()
 plt.title(f"{L=}, $\chi=${bond_dims[-1]}")
 plt.grid()
 plt.tight_layout()
-plt.savefig(PATH + f"eval_spectrum_L_{L}_{time_string}.png")
+plt.savefig(PATH + f"ent_ent_spectrum_L_{L}_{time_string}.png")
 
 # plot stationary densities
 plt.figure()

src/simulations/contact_process_dyn.py

@@ -11,14 +11,13 @@
 # # system parameters
 L = 21
 OMEGA = 6.0
-GAMMA = 1.0
 
 # # dynamics parameter
 step_size = 1
 number_of_steps = 7
 max_rank = 50
 
-lindblad = construct_lindblad(gamma=GAMMA, omega=OMEGA, L=L)
+lindblad = construct_lindblad(gamma=1.0, omega=OMEGA, L=L)
 lindblad_hermitian = lindblad.transpose(conjugate=True) @ lindblad
 num_op = construct_num_op(L)
 mps = tt.unit(L * [4], inds=L // 2 * [0] + [3] + L // 2 * [0])
@@ -33,7 +32,6 @@
 particle_num_t = np.zeros((len(t), L))
 
 for i in range(len(t)):
-    print(f"{i=}")
     hermit_mps = deepcopy(mps_t[i])
     hermit_mps = canonicalize_mps(hermit_mps)
     mps_dag = hermit_mps.transpose(conjugate=True)

src/simulations/contact_process_stat.py

@@ -15,11 +15,11 @@
 # system parameters
 L = 10
 # OMEGAS = np.linspace(0, 10, 10)
-OMEGAS = np.array([6.0])
+OMEGAS = np.array([0.0])
 
 # TN algorithm parameters
 # bond_dims = np.array([8, 16, 20])
-bond_dims = np.array([30])
+bond_dims = np.array([8])
 conv_eps = 1e-6
 
 ### Stationary simulation
@@ -84,7 +84,7 @@
             purities[i] = compute_purity(gs_mps)
             print(f"Purity: {purities[-1]}")
 
-            print("Compute density correlation for largest bond dimension")
+            print("Compute density-density correlation for largest bond dimension")
             an_op = construct_num_op(1)
             for k in range(L - 1):
                 correlations[i, k] = abs(compute_correlation_vMPO(gs_mps, an_op, r0=0,

tests/simple_test.py

@@ -0,0 +1,72 @@
+"""
+Checking the total polarization of the steady state
+"""
+
+import numpy as np
+
+import scikit_tt.tensor_train as tt
+from scikit_tt.solvers.evp import als
+
+import src.models.diss_ising_model as diss_ising_model
+import src.utilities.utils as utils
+
+L = 10
+GAMMA = 1.0
+OMEGA = 1.5
+V = 5.0
+# DELTAS = np.linspace(-4, 4, 9)
+DELTA = 0.0
+
+bond_dim = 16
+conv_eps = 1e-6
+
+print(f"{L=}")
+print(f"{DELTA=}")
+print(f"{bond_dim=}")
+
+ising_lindblad = diss_ising_model.construct_lindblad(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
+ising_lindblad_dag = diss_ising_model.construct_lindblad_dag(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
+ising_lindblad_hermitian = ising_lindblad_dag @ ising_lindblad
+
+mps = tt.ones(row_dims=L * [4], col_dims=L * [1], ranks=bond_dim)
+mps = mps.ortho()
+mps = (1 / mps.norm()**2) * mps
+eigenvalue, eigentensor, _ = als(ising_lindblad_hermitian, mps, number_ev=1, repeats=10, conv_eps=conv_eps, sigma=0)
+
+gs_mps = utils.canonicalize_mps(eigentensor)
+gs_mps_dag = gs_mps.transpose(conjugate=True)
+
+print(f"expVal_0, eigenvalue_0: {utils.compute_expVal(gs_mps, ising_lindblad_hermitian)}, {eigenvalue}")
+
+# compute non-Hermitian part of mps
+non_hermit_mps = (1 / 2) * (gs_mps - gs_mps_dag)
+print(f"The norm of the non-Hermitian part: {non_hermit_mps.norm()**2}")
+
+# compute Hermitian part of mps
+hermit_mps = (1 / 2) * (gs_mps + gs_mps_dag)
+
+print("Compute overlap with dark state")
+basis_0 = np.array([1, 0])
+dark_state = utils.construct_basis_mps(L, basis=[np.kron(basis_0, basis_0)] * L)
+print(f"{utils.compute_overlap(dark_state, gs_mps)=}")
+
+print("Compute local polarization")
+sigmaz = 0.5 * np.array([[1, 0], [0, -1]])
+sigmaz = np.kron(sigmaz, np.eye(2)) + np.kron(np.eye(2), sigmaz)
+an_op = tt.TT(sigmaz)
+an_op.cores[0] = an_op.cores[0].reshape(2, 2, 2, 2)
+polarizations = utils.compute_site_expVal_vMPO(hermit_mps, an_op)
+print(f"Local polarization of ground state: {polarizations}")
+print(f"Local polarization of dark state: {utils.compute_site_expVal_vMPO(dark_state, an_op)}")
+
+print("Compute two-site correlation")
+
+correlations = diss_ising_model.commpute_half_chain_corr(hermit_mps)
+print(f"Magnetization correlation of ground state: {diss_ising_model.commpute_half_chain_corr(hermit_mps)}")
+print(f"Magnetization correlation of dark state: {diss_ising_model.commpute_half_chain_corr(dark_state)}")
+
+# import matplotlib.pyplot as plt
+
+# plt.figure()
+# plt.plot(list(range(L//2+1, L)), correlations)
+# plt.show()

tests/test_diss_ising.py

@@ -0,0 +1,72 @@
+import numpy as np
+from datetime import datetime
+
+import scikit_tt.tensor_train as tt
+from scikit_tt.solvers.evp import als
+
+import src.models.diss_ising_model as diss_ising_model
+import src.utilities.utils as utils
+
+# path for results
+PATH = "/scratch/nguyed99/qcp-1d/results/"
+
+L = 10
+GAMMA = 1.0
+OMEGA = 1.5
+V = 5.0
+# DELTAS = np.linspace(-4, 4, 9)
+DELTAS = [0.0]
+
+bond_dim = 8
+conv_eps = 1e-6
+
+purities = np.zeros(len(DELTAS))
+staggered_mag = np.zeros(len(DELTAS))
+correlations = np.zeros((3, L // 2 - 1))
+polarizations = np.zeros((3, L))
+
+for i, DELTA in enumerate(DELTAS):
+    ising_lindblad = diss_ising_model.construct_lindblad(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
+    ising_lindblad_dag = diss_ising_model.construct_lindblad_dag(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
+    ising_lindblad_hermitian = ising_lindblad_dag @ ising_lindblad
+
+    mps = tt.ones(row_dims=L * [4], col_dims=L * [1], ranks=bond_dim)
+    mps = mps.ortho()
+    mps = (1 / mps.norm()**2) * mps
+    eigenvalue, eigentensor, _ = als(ising_lindblad_hermitian, mps, number_ev=1, repeats=10, conv_eps=conv_eps, sigma=0)
+
+    gs_mps = utils.canonicalize_mps(eigentensor)
+    gs_mps_dag = gs_mps.transpose(conjugate=True)
+
+    # compute non-Hermitian part of mps
+    non_hermit_mps = (1 / 2) * (gs_mps - gs_mps_dag)
+    print(f"The norm of the non-Hermitian part: {non_hermit_mps.norm()**2}")
+
+    # compute Hermitian part of mps
+    hermit_mps = (1 / 2) * (gs_mps + gs_mps_dag)
+
+    purities[i] = utils.compute_purity(gs_mps)
+
+    # compute observables
+    staggered_mag[i] = diss_ising_model.compute_staggered_mag(hermit_mps)
+
+    j = 0
+    if DELTA in [-4, 4, 0]:
+        correlations[j, :] = diss_ising_model.commpute_half_chain_corr(hermit_mps)
+
+        # compute local polarization
+        sigmaz = 0.5 * np.array([[1, 0], [0, -1]])
+        sigmaz = np.kron(sigmaz, np.eye(2)) + np.kron(np.eye(2), sigmaz)
+        an_op = tt.TT(sigmaz)
+        an_op.cores[0] = an_op.cores[0].reshape(2, 2, 2, 2)
+
+        polarizations[j, :] = utils.compute_site_expVal_vMPO(hermit_mps, an_op)
+
+        j += 1
+
+time3 = "{:%Y_%m_%d_%H_%M_%S}".format(datetime.now())
+
+np.savetxt(PATH + f"ising_purities_L_{L}_{time3}.txt", purities, delimiter=',')
+np.savetxt(PATH + f"ising_staggered_mag_L_{L}_{time3}.txt", staggered_mag, delimiter=',')
+np.savetxt(PATH + f"ising_correlations_L_{L}_{time3}.txt", correlations, delimiter=',')
+np.savetxt(PATH + f"ising_polarizations_L_{L}_{time3}.txt", polarizations, delimiter=',')