Update models and tests

README.md

@@ -16,11 +16,8 @@
 - 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!
-- scikit-tt norm needs to be squared
-- Automate testing pipeline
-- reproducibility test
-- results from time-propagation is strange
-- check for unphysicality of solution [it doesn't make sense to have purity = 0?]
+- Result from time-propagation is strange.
+- Negative values for site expectation value.
 
 ## Technical details 
 Required tools: `poetry`

src/hpc/contact_process_stat.py

@@ -8,7 +8,7 @@
 
 from src.models.contact_process_model import (construct_lindblad, construct_num_op)
 from src.utilities.utils import (canonicalize_mps, compute_correlation_vMPO, compute_purity, compute_site_expVal_vMPO,
-                                 compute_expVal, compute_eigenvalue_spectrum)
+                                 compute_expVal, compute_entanglement_spectrum, construct_basis_mps, compute_overlap)
 
 logger = logging.getLogger(__name__)
 log_filename = datetime.now().strftime("%Y-%m-%d_%H-%M-%S.log")
@@ -18,8 +18,7 @@
 PATH = "/scratch/nguyed99/qcp-1d/results/"
 
 # system parameters
-L = 10
-GAMMA = 1
+L = 50
 OMEGAS = np.linspace(0, 10, 10)
 
 # TN algorithm parameters
@@ -28,19 +27,20 @@
 
 ### Stationary simulation
 logger.info("Stationary simulation")
-spectral_gaps = np.zeros(len(OMEGAS))
-n_s = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
-evp_residual = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
+d = 2
 eval_0 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 eval_1 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
-eval_spectrum = np.zeros((OMEGAS.shape[0], bond_dims[-1]))
+evp_residual = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
+spectral_gaps = np.zeros(len(OMEGAS))
+n_s = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
+ent_ent_spectrum = np.zeros((OMEGAS.shape[0], bond_dims[-1] * d**2))  # d² * bond_dim
 purities = np.zeros(len(OMEGAS))
 correlations = np.zeros((len(OMEGAS), L - 1))
 
 for i, OMEGA in enumerate(OMEGAS):
     for j, bond_dim in enumerate(bond_dims):
         logger.info(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
-        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
 
         mps = tt.ones(row_dims=L * [4], col_dims=L * [1], ranks=bond_dim)
@@ -61,19 +61,21 @@
 
         gs_mps = canonicalize_mps(eigentensors[0])
         gs_mps_dag = gs_mps.transpose(conjugate=True)
-        logger.info(f"expVal_0, eigenvalue_0: {compute_expVal(gs_mps, lindblad_hermitian)}, {eigenvalues[0]}")
-
-        # compute Hermitian part of mps
-        hermit_mps = (1 / 2) * (gs_mps + gs_mps_dag)
 
         # compute non-Hermitian part of mps
         non_hermit_mps = (1 / 2) * (gs_mps - gs_mps_dag)
         logger.info(f"The norm of the non-Hermitian part: {non_hermit_mps.norm()**2}")
 
+        logger.info(f"expVal_0, eigenvalue_0: {compute_expVal(gs_mps, lindblad_hermitian)}, {eigenvalues[0]}")
+
+        # compute Hermitian part of mps
+        hermit_mps = (1 / 2) * (gs_mps + gs_mps_dag)
+
         logger.info("Compute particle numbers")
         particle_nums = compute_site_expVal_vMPO(hermit_mps, construct_num_op(L))
         logger.info(f"Particle number/site: {particle_nums}")
         n_s[j, i] = np.mean(particle_nums)
+        logger.info(f"Mean Particle number: {n_s[j, i]}")
 
         if bond_dim == bond_dims[-1]:
             logger.info(f"{bond_dim=}")
@@ -89,17 +91,22 @@
             for k in range(L - 1):
                 correlations[i, k] = abs(compute_correlation_vMPO(gs_mps, an_op, r0=0, r1=k + 1))
 
-            logger.info("Compute half-chain eigenvalue spectrum for largest bond dimension")
-            eval_spectrum[i, :] = compute_eigenvalue_spectrum(eigentensors[0])
+            logger.info("Compute half-chain entanglement entropy spectrum for largest bond dimension")
+            ent_ent_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
+
+            logger.info("Compute overlap with dark state")
+            basis_0 = np.array([1, 0])
+            dark_state = construct_basis_mps(L, basis=[np.kron(basis_0, basis_0)] * L)
+            logger.info(f"{compute_overlap(dark_state, gs_mps)=}")  #TODO: negative?
 
 time3 = "{:%Y_%m_%d_%H_%M_%S}".format(datetime.now())
 
 # save result arrays
 np.savetxt(PATH + f"eval_0_L_{L}_{time3}.txt", eval_0, delimiter=',')
 np.savetxt(PATH + f"eval_1_L_{L}_{time3}.txt", eval_1, delimiter=',')
-np.savetxt(PATH + f"eval_spectrum_L_{L}_{time3}.txt", eval_spectrum, delimiter=',')
 np.savetxt(PATH + f"evp_residual_L_{L}_{time3}.txt", evp_residual, delimiter=',')
 np.savetxt(PATH + f"spectral_gaps_L_{L}_{time3}.txt", spectral_gaps, delimiter=',')
 np.savetxt(PATH + f"n_s_L_{L}_{time3}.txt", n_s, delimiter=',')
+np.savetxt(PATH + f"ent_ent_spectrum_L_{L}_{time3}.txt", ent_ent_spectrum, delimiter=',')
 np.savetxt(PATH + f"purities_L_{L}_{time3}.txt", purities, delimiter=',')
 np.savetxt(PATH + f"correlations_L_{L}_{time3}.txt", correlations, delimiter=',')

src/hpc/jobs_script.sh

@@ -45,7 +45,12 @@ export MKL_NUM_THREADS=$SLURM_CPUS_PER_TASK
 export VECLIB_MAXIMUM_THREADS=$SLURM_CPUS_PER_TASK
 export NUMEXPR_NUM_THREADS=$SLURM_CPUS_PER_TASK
 
-L=10
+# # set up environment
+# python3 -m venv /scratch/nguyed99/tensor # create venv
+# pip install -r requirements.txt
+
+
+L=50
 LOG_PATH="/scratch/nguyed99/qcp-1d/logging"
 timestamp=$(date +'%Y-%m-%d-%H-%M-%S')
 
@@ -72,9 +77,9 @@ log_time() {
 source /scratch/nguyed99/tensor/bin/activate
 
 # launch Python script
-log_memory_cpu_usage & log_time
+log_memory_cpu_usage & 
 LOG_PID=$!
 
-python3 contact_process_stat.py > "$LOG_PATH/contact_process_stat_L_${L}_${timestamp}.out"
+log_time python3 contact_process_stat.py > "$LOG_PATH/contact_process_stat_L_${L}_${timestamp}.out"
 
 kill $LOG_PID

src/hpc/requirements.txt

@@ -1,4 +1,4 @@
 numpy==1.26.3
 scipy==1.11.4
 matplotlib==3.8.2
-git+https://github.com/PGelss/scikit_tt.git@1f681e1
+git+https://github.com/PGelss/scikit_tt.git@616ff8b

src/models/contact_process_model.py

@@ -22,9 +22,9 @@ def construct_lindblad(gamma: float, omega: float, L: int) -> tt.TT:
     L_4 = np.kron(identity, number_op)
     Id = np.eye(4)
     M_1 = -1j * omega * np.kron(number_op, identity)
-    M_2 = -1j * omega * np.kron(identity, number_op)
+    M_2 = 1j * omega * np.kron(identity, number_op)
     M_3 = -1j * omega * np.kron(sigmax, identity)
-    M_4 = -1j * omega * np.kron(identity, sigmax)
+    M_4 = 1j * omega * np.kron(identity, sigmax)
 
     # construct core
     op_cores = [None] * L
@@ -78,9 +78,9 @@ def construct_lindblad_dag(gamma: float, omega: float, L: int) -> tt.TT:
     L_4 = np.kron(identity, number_op)
     Id = np.eye(4)
     M_1 = 1j * omega * np.kron(number_op, identity)
-    M_2 = 1j * omega * np.kron(identity, number_op)
+    M_2 = -1j * omega * np.kron(identity, number_op)
     M_3 = 1j * omega * np.kron(sigmax, identity)
-    M_4 = 1j * omega * np.kron(identity, sigmax)
+    M_4 = -1j * omega * np.kron(identity, sigmax)
 
     # construct core
     op_cores = [None] * L

src/models/ising_model.py

@@ -23,23 +23,23 @@ def construct_lindblad(gamma: float, V: float, omega: float, delta: float, L: in
 
     annihilation_op = np.array([[0, 1], [0, 0]])
     creation_op = np.array([[0, 0], [1, 0]])
-    an_op = annihilation_op * creation_op
+    an_op = creation_op * annihilation_op
     an_op_L = np.kron(an_op, identity)
     an_op_R = np.kron(identity, an_op)
 
     # core components
-    S = -1j * (omega / 2) * (sigmax_L + sigmax_R) + 1j * (V - delta) / 2 * (sigmaz_L + sigmaz_R) + gamma * np.kron(
-        creation_op, creation_op) - (1 / 2) * gamma * (an_op_L + an_op_R)
+    S = -1j * (omega / 2) * (sigmax_L - sigmax_R) + 1j * (V - delta) / 2 * (sigmaz_L - sigmaz_R) + gamma * np.kron(
+        annihilation_op, annihilation_op) - (1 / 2) * gamma * (an_op_L + an_op_R)
     L_1 = sigmaz_L
     L_2 = sigmaz_R
     Id = np.eye(4)
     M_1 = -1j * V / 4 * sigmaz_L
-    M_2 = -1j * V / 4 * sigmaz_R
+    M_2 = -1j * V / 4 * -sigmaz_R
 
     # construct core
     op_cores = [None] * L
     op_cores[0] = np.zeros([1, 4, 4, 4], dtype=complex)
-    op_cores[0][0, :, :, 0] = S + -1j * (V / 4) * (sigmaz_L + sigmaz_R)
+    op_cores[0][0, :, :, 0] = S + -1j * (V / 4) * (sigmaz_L - sigmaz_R)
     op_cores[0][0, :, :, 1] = L_1
     op_cores[0][0, :, :, 2] = L_2
     op_cores[0][0, :, :, 3] = Id
@@ -56,7 +56,7 @@ def construct_lindblad(gamma: float, V: float, omega: float, delta: float, L: in
     op_cores[-1][0, :, :, 0] = Id
     op_cores[-1][1, :, :, 0] = M_1
     op_cores[-1][2, :, :, 0] = M_2
-    op_cores[-1][3, :, :, 0] = S + -1j * (V / 4) * (np.kron(identity, sigmaz) + np.kron(sigmaz, identity))
+    op_cores[-1][3, :, :, 0] = S + -1j * (V / 4) * (sigmaz_L - sigmaz_R)
 
     return TT(op_cores)
 
@@ -78,23 +78,23 @@ def construct_lindblad_dag(gamma: float, V: float, omega: float, delta: float, L
 
     annihilation_op = np.array([[0, 1], [0, 0]])
     creation_op = np.array([[0, 0], [1, 0]])
-    an_op = annihilation_op * creation_op
+    an_op = creation_op * annihilation_op
     an_op_L = np.kron(an_op, identity)
     an_op_R = np.kron(identity, an_op)
 
     # core components
-    S = 1j * (omega / 2) * (sigmax_L + sigmax_R) - 1j * (V - delta) / 2 * (sigmaz_L + sigmaz_R) + gamma * np.kron(
-        annihilation_op, annihilation_op) - (1 / 2) * gamma * (an_op_L + an_op_R)
+    S = 1j * (omega / 2) * (sigmax_L - sigmax_R) - 1j * (V - delta) / 2 * (sigmaz_L - sigmaz_R) + gamma * np.kron(
+        creation_op, creation_op) - (1 / 2) * gamma * (an_op_L + an_op_R)
     L_1 = sigmaz_L
     L_2 = sigmaz_R
     Id = np.eye(4)
     M_1 = 1j * V / 4 * sigmaz_L
-    M_2 = 1j * V / 4 * sigmaz_R
+    M_2 = -1j * V / 4 * sigmaz_R
 
     # construct core
     op_cores = [None] * L
     op_cores[0] = np.zeros([1, 4, 4, 4], dtype=complex)
-    op_cores[0][0, :, :, 0] = S + 1j * (V / 4) * (sigmaz_L + sigmaz_R)
+    op_cores[0][0, :, :, 0] = S + 1j * (V / 4) * (sigmaz_L - sigmaz_R)
     op_cores[0][0, :, :, 1] = L_1
     op_cores[0][0, :, :, 2] = L_2
     op_cores[0][0, :, :, 3] = Id
@@ -111,6 +111,6 @@ def construct_lindblad_dag(gamma: float, V: float, omega: float, delta: float, L
     op_cores[-1][0, :, :, 0] = Id
     op_cores[-1][1, :, :, 0] = M_1
     op_cores[-1][2, :, :, 0] = M_2
-    op_cores[-1][3, :, :, 0] = S + 1j * (V / 4) * (np.kron(identity, sigmaz) + np.kron(sigmaz, identity))
+    op_cores[-1][3, :, :, 0] = S + 1j * (V / 4) * (sigmaz_L - sigmaz_R)
 
     return TT(op_cores)

src/simulations/contact_process_stat.py

@@ -14,31 +14,30 @@
 
 # system parameters
 L = 10
-GAMMA = 1
 # OMEGAS = np.linspace(0, 10, 10)
-OMEGAS = np.array([0.0])
+OMEGAS = np.array([6.0])
 
 # TN algorithm parameters
 # bond_dims = np.array([8, 16, 20])
-bond_dims = np.array([8])
+bond_dims = np.array([30])
 conv_eps = 1e-6
 
 ### Stationary simulation
 print("Stationary simulation")
 d = 2
-spectral_gaps = np.zeros(len(OMEGAS))
-n_s = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
-evp_residual = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 eval_0 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 eval_1 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
-entropy_spectrum = np.zeros((OMEGAS.shape[0], bond_dims[-1] * d**2))  # d² * bond_dim
+evp_residual = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
+spectral_gaps = np.zeros(len(OMEGAS))
+n_s = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
+ent_ent_spectrum = np.zeros((OMEGAS.shape[0], bond_dims[-1] * d**2))  # d² * bond_dim
 purities = np.zeros(len(OMEGAS))
 correlations = np.zeros((len(OMEGAS), L - 1))
 
 for i, OMEGA in enumerate(OMEGAS):
     for j, bond_dim in enumerate(bond_dims):
         print(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
-        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
 
         mps = tt.ones(row_dims=L * [4], col_dims=L * [1], ranks=bond_dim)
@@ -59,19 +58,21 @@
 
         gs_mps = canonicalize_mps(eigentensors[0])
         gs_mps_dag = gs_mps.transpose(conjugate=True)
-        print(f"expVal_0, eigenvalue_0: {compute_expVal(gs_mps, lindblad_hermitian)}, {eigenvalues[0]}")
-
-        # compute Hermitian part of mps
-        hermit_mps = (1 / 2) * (gs_mps + gs_mps_dag)
 
         # 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}")
 
+        print(f"expVal_0, eigenvalue_0: {compute_expVal(gs_mps, lindblad_hermitian)}, {eigenvalues[0]}")
+
+        # compute Hermitian part of mps
+        hermit_mps = (1 / 2) * (gs_mps + gs_mps_dag)
+
         print("Compute particle numbers")
         particle_nums = compute_site_expVal_vMPO(hermit_mps, construct_num_op(L))
         print(f"Particle number/site: {particle_nums}")
         n_s[j, i] = np.mean(particle_nums)
+        print(f"Mean Particle number: {n_s[j, i]}")
 
         if bond_dim == bond_dims[-1]:
             print(f"{bond_dim=}")
@@ -88,7 +89,7 @@
                 correlations[i, k] = abs(compute_correlation_vMPO(gs_mps, an_op, r0=0, r1=k + 1))
 
             print("Compute half-chain entanglement entropy spectrum for largest bond dimension")
-            entropy_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
+            ent_ent_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
 
             print("Compute overlap with dark state")
             basis_0 = np.array([1, 0])
@@ -100,9 +101,9 @@
 # save result arrays
 # np.savetxt(PATH + f"eval_0_L_{L}_{time3}.txt", eval_0, delimiter=',')
 # np.savetxt(PATH + f"eval_1_L_{L}_{time3}.txt", eval_1, delimiter=',')
-# np.savetxt(PATH + f"eval_spectrum_L_{L}_{time3}.txt", eval_spectrum, delimiter=',')
 # np.savetxt(PATH + f"evp_residual_L_{L}_{time3}.txt", evp_residual, delimiter=',')
 # np.savetxt(PATH + f"spectral_gaps_L_{L}_{time3}.txt", spectral_gaps, delimiter=',')
 # np.savetxt(PATH + f"n_s_L_{L}_{time3}.txt", n_s, delimiter=',')
+# np.savetxt(PATH + f"ent_ent_spectrum_L_{L}_{time3}.txt", ent_ent_spectrum, delimiter=',')
 # np.savetxt(PATH + f"purities_L_{L}_{time3}.txt", purities, delimiter=',')
 # np.savetxt(PATH + f"correlations_L_{L}_{time3}.txt", correlations, delimiter=',')

src/utilities/utils.py

@@ -138,7 +138,7 @@ def compute_expVal(mps: tt.TT, mpo: tt.TT) -> float:
         left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1, 2], [0, 2, 4]))
 
     exp_val = np.trace(left_boundary @ right_boundary).item()
-    if exp_val.imag < 1e-12:
+    if exp_val.imag < 1e-7:
         return exp_val.real
     else:
         raise ValueError("Complex expectation value is found.")

tests/test_model.py

@@ -37,10 +37,17 @@ def test_hermiticity(self):
         cp_lindblad_hermitian = cp_lindblad_dag @ cp_lindblad
         assert np.array_equal(np.conj(cp_lindblad_hermitian.matricize().T), cp_lindblad_hermitian.matricize())
 
-    def test_exact_diagonalization(self):
+    def test_exact_diagonalization(self):  #TODO: Add ising model too :)
         cp_lindblad = cp_model.construct_lindblad(gamma=GAMMA, omega=OMEGA, L=L)
         cp_lindblad_dag = cp_model.construct_lindblad_dag(gamma=GAMMA, omega=OMEGA, L=L)
         cp_lindblad_hermitian = cp_lindblad_dag @ cp_lindblad
 
         evals, _ = np.linalg.eig(cp_lindblad_hermitian.matricize())
-        assert np.min(evals.real) == 0.0
+        assert np.isclose(np.min(evals.real), 0.0)
+
+        ising_lindblad = ising_model.construct_lindblad(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
+        ising_lindblad_dag = ising_model.construct_lindblad_dag(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
+        ising_lindblad_hermitian = ising_lindblad_dag @ ising_lindblad
+
+        evals, _ = np.linalg.eig(ising_lindblad_hermitian.matricize())
+        assert np.isclose(np.min(evals.real), 0.0)