Fix models

src/hpc/contact_process_stat.py

@@ -1,3 +1,5 @@
+import argparse
+import itertools
 import logging
 import numpy as np
 import time
@@ -24,11 +26,14 @@
 OMEGAS = np.linspace(0.5, 10, 10)
 
 # TN algorithm parameters
-bond_dims = np.array([50, 70])
+bond_dims = np.array([35, 50])
 conv_eps = 1e-6
 
-### Stationary simulation
-print("Stationary simulation")
+### storing observables
+params = list(itertools.product(bond_dims, OMEGAS))
+indices = list(itertools.product(list(range(len(bond_dims))), list(range(len(OMEGAS)))))
+indices = dict(zip(params, indices))
+
 eval_0 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 eval_1 = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
 evp_residual = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
@@ -45,80 +50,90 @@
 number_mpo = [None]
 number_mpo = tt.TT(number_op)
 
-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=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)
-        mps = mps.ortho()
-        mps = (1 / mps.norm()**2) * mps
-        time1 = time.time()
-        eigenvalues, eigentensors, _ = als(lindblad_hermitian, mps, number_ev=2, repeats=10, conv_eps=conv_eps, sigma=0)
-        time2 = time.time()
-        print(f"Elapsed time: {time2 - time1} seconds")
-
-        evp_residual[j, i] = (lindblad_hermitian @ eigentensors[0] - eigenvalues[0] * eigentensors[0]).norm()**2
-        eval_0[j, i] = eigenvalues[0]
-        eval_1[j, i] = eigenvalues[1]
-        print(f"Eigensolver error: {evp_residual[j, i]}")
-
-        print(f"Ground state energy per site E: {eigenvalues/L}")
-        print(f"Norm of ground state: {eigentensors[0].norm()**2}")
-
-        gs_mps = canonicalize_mps(eigentensors[0])
-        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}")
-
-        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 observables
-        print("Compute particle numbers")
-        particle_nums = compute_site_expVal(hermit_mps, number_mpo)
-        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=}")
-            print("Compute spectral gap of L†L for largest bond dimension")
-            spectral_gaps[i] = abs(eigenvalues[1] - eigenvalues[0])
-
-            print("Compute purity of state for largest bond dimension")
-            purities[i] = compute_purity(gs_mps)
-            print(f"Purity: {purities[-1]}")
-
-            logger.info("Compute two-point correlation for largest bond dimension")
-            for k in range(L - 1):
-                correlations[i, k] = compute_correlation(gs_mps, number_mpo, r0=0, r1=k + 1)
-
-            print("Compute density-density correlation for largest bond dimension")
-            for k in range(L - 1):
-                dens_dens_corr[i, k] = compute_dens_dens_corr(gs_mps, number_mpo, r=k + 1)
-
-            print("Compute half-chain entanglement spectrum for largest bond dimension")
-            entanglement_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
-
-            basis_0 = np.array([1, 0])
-            dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
-            print(f"Overlap with dark state: {compute_overlap(dark_state, gs_mps)}")  #NOTE: 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"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"entanglement_spectrum_L_{L}_{time3}.txt", entanglement_spectrum, delimiter=',')
-np.savetxt(PATH + f"purities_L_{L}_{time3}.txt", purities, delimiter=',')
-np.savetxt(PATH + f"correlations_L_{L}_{time3}.txt", correlations, delimiter=',')
-np.savetxt(PATH + f"dens_dens_corr_L_{L}_{time3}.txt", dens_dens_corr, delimiter=',')
+
+def main():
+    # parse environment variables
+    parser = argparse.ArgumentParser()
+    parser.add_argument('slurmarrayid', type=int)
+    args = parser.parse_args()
+
+    OMEGA, bond_dim = params(args.slurmarrayid)
+
+    j, i = indices[(bond_dim, OMEGA)]
+
+    print("Stationary simulation")
+    print(f"Run ALS for {L=}, {OMEGA=} and {bond_dim=}")
+    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)
+    mps = mps.ortho()
+    mps = (1 / mps.norm()**2) * mps
+    time1 = time.time()
+    eigenvalues, eigentensors, _ = als(lindblad_hermitian, mps, number_ev=2, repeats=10, conv_eps=conv_eps, sigma=0)
+    time2 = time.time()
+    print(f"Elapsed time: {time2 - time1} seconds")
+
+    evp_residual[j, i] = (lindblad_hermitian @ eigentensors[0] - eigenvalues[0] * eigentensors[0]).norm()**2
+    eval_0[j, i] = eigenvalues[0]
+    eval_1[j, i] = eigenvalues[1]
+    print(f"Eigensolver error: {evp_residual[j, i]}")
+
+    print(f"Ground state energy per site E: {eigenvalues/L}")
+    print(f"Norm of ground state: {eigentensors[0].norm()**2}")
+
+    gs_mps = canonicalize_mps(eigentensors[0])
+    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}")
+
+    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 observables
+    print("Compute particle numbers")
+    particle_nums = compute_site_expVal(hermit_mps, number_mpo)
+    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=}")
+        print("Compute spectral gap of L†L for largest bond dimension")
+        spectral_gaps[i] = abs(eigenvalues[1] - eigenvalues[0])
+
+        print("Compute purity of state for largest bond dimension")
+        purities[i] = compute_purity(gs_mps)
+        print(f"Purity: {purities[-1]}")
+
+        logger.info("Compute two-point correlation for largest bond dimension")
+        for k in range(L - 1):
+            correlations[i, k] = compute_correlation(gs_mps, number_mpo, r0=0, r1=k + 1)
+
+        print("Compute density-density correlation for largest bond dimension")
+        for k in range(L - 1):
+            dens_dens_corr[i, k] = compute_dens_dens_corr(gs_mps, number_mpo, r=k + 1)
+
+        print("Compute half-chain entanglement spectrum for largest bond dimension")
+        entanglement_spectrum[i, :] = compute_entanglement_spectrum(gs_mps)
+
+        basis_0 = np.array([1, 0])
+        dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
+        print(f"Overlap with dark state: {compute_overlap(dark_state, gs_mps)}")  #NOTE: 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"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"entanglement_spectrum_L_{L}_{time3}.txt", entanglement_spectrum, delimiter=',')
+    np.savetxt(PATH + f"purities_L_{L}_{time3}.txt", purities, delimiter=',')
+    np.savetxt(PATH + f"correlations_L_{L}_{time3}.txt", correlations, delimiter=',')
+    np.savetxt(PATH + f"dens_dens_corr_L_{L}_{time3}.txt", dens_dens_corr, delimiter=',')

src/models/contact_process_model.py

@@ -14,8 +14,8 @@ def construct_lindblad(gamma: float, omega: float, L: int) -> tt.TT:
     number_op = np.array([[0, 0], [0, 1]])
 
     # core components
-    S = gamma * np.kron(annihilation_op,
-                        annihilation_op) - 0.5 * (np.kron(number_op, identity) + np.kron(identity, number_op))
+    S = gamma * np.kron(annihilation_op, annihilation_op) - 0.5 * (
+        np.kron(number_op, identity) + np.kron(identity, number_op))  # on-site
     L_1 = np.kron(sigmax, identity)
     L_2 = np.kron(identity, sigmax)
     L_3 = np.kron(number_op, identity)

src/models/diss_ising_model.py

@@ -24,18 +24,19 @@ 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 = creation_op * annihilation_op
+    an_op = annihilation_op @ creation_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)
-    L_1 = sigmaz_L
-    L_2 = sigmaz_R
+        creation_op, creation_op) - (1 / 2) * gamma * (an_op_L + an_op_R)
+    L_1 = sigmaz_R
+    L_2 = sigmaz_L
     Id = np.eye(4)
-    M_1 = -1j * V / 4 * sigmaz_L
-    M_2 = -1j * V / 4 * -sigmaz_R
+    M_1 = 1j * V / 4 * sigmaz_R
+    M_2 = -1j * V / 4 * sigmaz_L
 
     # construct core
     op_cores = [None] * L
@@ -79,18 +80,19 @@ 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 = creation_op * annihilation_op
+    an_op = annihilation_op @ creation_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)
-    L_1 = sigmaz_L
-    L_2 = sigmaz_R
+        annihilation_op, annihilation_op) - (1 / 2) * gamma * (an_op_L + an_op_R)
+    L_1 = sigmaz_R
+    L_2 = sigmaz_L
     Id = np.eye(4)
-    M_1 = 1j * V / 4 * sigmaz_L
-    M_2 = -1j * V / 4 * sigmaz_R
+    M_1 = -1j * V / 4 * sigmaz_R
+    M_2 = 1j * V / 4 * sigmaz_L
 
     # construct core
     op_cores = [None] * L
@@ -123,43 +125,70 @@ def commpute_half_chain_corr(mps: tt.TT) -> np.ndarray:
     corr = []
 
     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)
+    an_op.cores[0] = an_op.cores[0].reshape(1, 2, 2, 1)
 
     for i in range(r0 + 1, L):
         corr.append(abs(compute_correlation(mps, an_op, r0=r0, r1=i)))
     return np.array(corr)
 
 
-def compute_staggered_mag(mps: tt.TT) -> float:
+def compute_staggered_mag_1(mps: tt.TT) -> float:
 
     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)
-
-    staggered_mag = np.tensordot(an_op.cores[0], an_op.cores[0], axes=([2, 3], [0, 1]))
-    staggered_mag = tt.TT(staggered_mag)
-    staggered_mag.cores[0] = staggered_mag.cores[0].reshape(2, 2, 2, 2)
+    sigmaz_sq = sigmaz @ sigmaz
+    staggered_mag = tt.TT(sigmaz_sq)
+    staggered_mag.cores[0] = staggered_mag.cores[0].reshape(1, 2, 2, 1)
 
     site_vals = np.zeros(mps.order, dtype=float)
     for i in range(mps.order):
-        left_boundary = np.ones(1)
-        right_boundary = np.ones(1)
+        left_boundary = np.ones((1, 1))
+        right_boundary = np.ones((1, 1))
 
         for j in range(mps.order):
             if j == i:
-                contraction = np.tensordot(mps.cores[j], staggered_mag.cores[0], axes=([1, 2], [0, 1]))
-                contraction = np.einsum('ijlk->ij', contraction)  # tracing over the physical indices
+                contraction = np.tensordot(mps.cores[j], staggered_mag.cores[0], axes=([1, 2], [2, 1]))
             else:
-                contraction = np.einsum('ikjl->il', mps.cores[j])  # tracing over the physical indices
+                contraction = np.tensordot(mps.cores[j], np.eye(2).reshape(1, 2, 2, 1), axes=([1, 2], [2, 1]))
 
-            left_boundary = np.tensordot(left_boundary, contraction, axes=([0], [0]))
+            left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1], [0, 2]))
 
-        if (left_boundary @ right_boundary).imag < 1e-12:
-            site_vals[i] = (-1)**(i + 1) * (left_boundary @ right_boundary).real
+        if (left_boundary @ right_boundary).item().imag < 1e-12:
+            site_vals[i] = (-1)**(i + 1) * (left_boundary @ right_boundary).item().real
         else:
             raise ValueError("Complex expectation value is found.")
 
     return np.sqrt(np.mean(site_vals))
+
+
+def compute_staggered_mag_2(mps: tt.TT) -> float:
+
+    sigmaz = 0.5 * np.array([[1, 0], [0, -1]])
+    staggered_mag = tt.TT(sigmaz)
+    staggered_mag.cores[0] = staggered_mag.cores[0].reshape(1, 2, 2, 1)
+
+    left_boundary = np.ones((1, 1, 1))
+    right_boundary = np.ones((1, 1, 1))
+
+    site_vals = []
+    for i in range(mps.order):
+        for j in range(i, mps.order):
+            for k in range(mps.order):
+                if k == i == j:
+                    contraction = np.tensordot(staggered_mag.cores[0], staggered_mag.cores[0], axes=([1], [2]))
+                elif k != i and j != j:
+                    contraction = np.tensordot(np.eye(2).reshape(1, 2, 2, 1),
+                                               np.eye(2).reshape(1, 2, 2, 1),
+                                               axes=([1], [2]))
+                else:
+                    contraction = np.tensordot(staggered_mag.cores[0], np.eye(2).reshape(1, 2, 2, 1), axes=([1], [2]))
+
+                contraction = np.tensordot(mps.cores[k], contraction, axes=([1, 2], [1, 4]))
+
+                left_boundary = np.tensordot(left_boundary, contraction, axes=([0, 1, 2], [0, 2, 4]))
+
+            if (left_boundary @ right_boundary).item().imag < 1e-12:
+                site_vals.append((-1)**(i + j) * (left_boundary @ right_boundary).item().real)
+            else:
+                raise ValueError("Complex expectation value is found.")
+    return np.sqrt((1 / (mps.order)**2) * 2 * np.sum(site_vals))

src/simulations/contact_process_stat.py

@@ -17,7 +17,7 @@
 L = 10
 d = 2  # physical dimension
 # OMEGAS = np.linspace(0, 10, 10)
-OMEGAS = np.array([0.0])
+OMEGAS = np.array([2.0])
 
 # TN algorithm parameters
 # bond_dims = np.array([8, 16, 20])

src/simulations/diss_ising_chain.py

@@ -1,16 +1,15 @@
 import scikit_tt.tensor_train as tt
 from scikit_tt.solvers.evp import als
 
-from models.diss_ising_model import construct_lindblad, construct_lindblad_dag
-
+from src.models.diss_ising_model import construct_lindblad, construct_lindblad_dag
 # path for results
 # PATH = "/home/psireal42/study/quantum-contact-process-1D/results"
 
 # system parameters
 GAMMA = 1.0
 OMEGA = 1.5
 V = 5
-DELTAS = [0.0]
+DELTAS = [-4.0]
 
 # TN algorithm parameters
 conv_eps = 1e-6
@@ -24,7 +23,12 @@
     lindblad = construct_lindblad(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L)
     lindblad_hermitian = construct_lindblad_dag(gamma=GAMMA, V=V, omega=OMEGA, delta=DELTA, L=L) @ lindblad
 
-    mps = tt.ones(row_dims=L * [4], col_dims=L * [1], ranks=bond_dim)
+    # mps = [None]*L
+    # for i in range(L):
+    #     mps[i] = np.eye(2)
+    # mps = t
+
+    mps = tt.rand(row_dims=L * [4], col_dims=L * [1], ranks=bond_dim)
     mps = mps.ortho()
     mps = (1 / mps.norm()**2) * mps
     eigenvalues, eigentensors, _ = als(lindblad_hermitian, mps, number_ev=1, repeats=10, conv_eps=conv_eps, sigma=0)