Update hpc script

src/hpc/contact_process_stat_exc.py

@@ -1,11 +1,10 @@
 """
 Do static simulations for first excited state
 """
-
 import logging
 import numpy as np
+import sys
 import time
-from datetime import datetime
 
 import scikit_tt.tensor_train as tt
 from scikit_tt.solvers.evp import als
@@ -15,114 +14,124 @@
                                  compute_site_expVal, 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")
-logging.basicConfig(filename="/scratch/nguyed99/qcp-1d/logging/" + log_filename, level=logging.INFO)
-
-# path for results
-PATH = "/scratch/nguyed99/qcp-1d/results/"
+OUTPUT_PATH = "/scratch/nguyed99/qcp-1d/results/"
+LOG_PATH = "/scratch/nguyed99/qcp-1d/logging/"
 
 # system parameters
 L = 50
 d = 2
-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
-logger.info("Stationary simulation")
-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]))
-spectral_gaps = np.zeros(len(OMEGAS))
-n_s = np.zeros((bond_dims.shape[0], OMEGAS.shape[0]))
-entanglement_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))
-dens_dens_corr = np.zeros((len(OMEGAS), L - 1))
-
 ### observable operator
 number_op = np.array([[0, 0], [0, 1]])
 number_op.reshape((1, 2, 2, 1))
 number_mpo = [None]
 number_mpo = tt.TT(number_op)
 
-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=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()
-        logger.info(f"Elapsed time: {time2 - time1} seconds")
-
-        evp_residual[j, i] = (lindblad_hermitian @ eigentensors[1] - eigenvalues[1] * eigentensors[1]).norm()**2
-        eval_0[j, i] = eigenvalues[0]
-        eval_1[j, i] = eigenvalues[1]
-        logger.info(f"Eigensolver error for first excited state: {evp_residual[j, i]}")
-
-        logger.info(f"First excited state energy per site: {eigenvalues[1]/L}")
-        logger.info(f"Norm of first excited state: {eigentensors[1].norm()**2}")
-
-        mps = canonicalize_mps(eigentensors[1])
-        mps_dag = mps.transpose(conjugate=True)
-
-        # compute non-Hermitian part of mps
-        non_hermit_mps = (1 / 2) * (mps - 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(mps, lindblad_hermitian)}, {eigenvalues[1]}")
-
-        # compute Hermitian part of mps
-        hermit_mps = (1 / 2) * (mps + mps_dag)
-
-        # compute observables
-        logger.info("Compute particle numbers")
-        particle_nums = compute_site_expVal(hermit_mps, number_mpo)
-        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=}")
-            logger.info("Compute spectral gap of L†L for largest bond dimension")
-            spectral_gaps[i] = abs(eigenvalues[1] - eigenvalues[0])
-
-            logger.info("Compute purity of state for largest bond dimension")
-            purities[i] = compute_purity(mps)
-            logger.info(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(mps, number_mpo, r0=0, r1=k + 1)
-
-            logger.info("Compute density-density correlation for largest bond dimension")
-            for k in range(L - 1):
-                dens_dens_corr[i, k] = compute_dens_dens_corr(mps, number_mpo, r=k + 1)
-
-            logger.info("Compute half-chain entanglement spectrum for largest bond dimension")
-            entanglement_spectrum[i, :] = compute_entanglement_spectrum(mps)
-
-            basis_0 = np.array([1, 0])
-            dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
-            logger.info(f"Overlap with dark state: {compute_overlap(dark_state, mps)}")  #NOTE: negative?
-
-time3 = "{:%Y_%m_%d_%H_%M_%S}".format(datetime.now())
-
-# save result arrays
-np.savetxt(PATH + f"eval_0_exc_L_{L}_{time3}.txt", eval_0, delimiter=',')
-np.savetxt(PATH + f"eval_1_exc_L_{L}_{time3}.txt", eval_1, delimiter=',')
-np.savetxt(PATH + f"evp_residual_exc_L_{L}_{time3}.txt", evp_residual, delimiter=',')
-np.savetxt(PATH + f"spectral_gaps_exc_L_{L}_{time3}.txt", spectral_gaps, delimiter=',')
-np.savetxt(PATH + f"n_s_L_exc_{L}_{time3}.txt", n_s, delimiter=',')
-np.savetxt(PATH + f"entanglement_spectrum_exc_L_{L}_{time3}.txt", entanglement_spectrum, delimiter=',')
-np.savetxt(PATH + f"purities_exc_L_{L}_{time3}.txt", purities, delimiter=',')
-np.savetxt(PATH + f"correlations_exc_L_{L}_{time3}.txt", correlations, delimiter=',')
-np.savetxt(PATH + f"dens_dens_corr_exc_L_{L}_{time3}.txt", dens_dens_corr, delimiter=',')
+
+def main():
+    # parse environment variables
+    OMEGA, bond_dim, SLURM_ARRAY_JOB_ID = float(sys.argv[1]), int(sys.argv[2]), str(sys.argv[3])
+    logger = logging.getLogger(__name__)
+    log_filename = f"L_{L}_OMEGA_{OMEGA}_D_{bond_dim}_{SLURM_ARRAY_JOB_ID}.log"
+    logging.basicConfig(filename=LOG_PATH + log_filename, level=logging.INFO)
+
+    logger.info("Stationary simulation")
+    logger.info(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()
+    logger.info(f"Elapsed time: {time2 - time1} seconds")
+
+    np.savetxt(OUTPUT_PATH + f"eval_0_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+               eigenvalues[0],
+               delimiter=',')
+    np.savetxt(OUTPUT_PATH + f"eval_1_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+               eigenvalues[1],
+               delimiter=',')
+
+    evp_residual = (lindblad_hermitian @ eigentensors[1] - eigenvalues[1] * eigentensors[1]).norm()**2
+    logger.info(f"{eigenvalues=}")
+    logger.info(f"Eigensolver error for first excited state: {evp_residual}")
+
+    logger.info(f"First excited state energy per site: {eigenvalues[1]/L}")
+    logger.info(f"Norm of first excited state: {eigentensors[1].norm()**2}")
+
+    mps = canonicalize_mps(eigentensors[1])
+    mps_dag = mps.transpose(conjugate=True)
+
+    # compute non-Hermitian part of mps
+    non_hermit_mps = (1 / 2) * (mps - mps_dag)
+    logger.info(f"The norm of the non-Hermitian part: {non_hermit_mps.norm()**2}")
+    logger.info(
+        f"expVal_1, eval_1_als, eval_1: {compute_expVal(mps, lindblad_hermitian)}, {eigenvalues[1]}, {mps_dag@lindblad_hermitian@mps}"
+    )
+
+    # compute Hermitian part of mps
+    hermit_mps = (1 / 2) * (mps + mps_dag)
+
+    # compute observables
+    logger.info("Compute particle numbers")
+    particle_nums = compute_site_expVal(hermit_mps, number_mpo)
+    logger.info(f"Particle number/site: {particle_nums}")
+    n_s = np.mean(particle_nums)
+    logger.info(f"Mean Particle number: {n_s}")
+    np.savetxt(OUTPUT_PATH + f"n_s_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt", n_s, delimiter=',')
+
+    if bond_dim == bond_dims[-1]:
+        logger.info(f"{bond_dim=}")
+        logger.info("Compute spectral gap of L†L for largest bond dimension")
+        spectral_gaps = abs(eigenvalues[1] - eigenvalues[0])
+        logger.info(f"{spectral_gaps=}")
+        np.savetxt(OUTPUT_PATH + f"spectral_gaps_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+                   spectral_gaps,
+                   delimiter=',')
+
+        logger.info("Compute purity of state for largest bond dimension")
+        purities = compute_purity(mps)
+        logger.info(f"Purity: {purities}")
+        np.savetxt(OUTPUT_PATH + f"purities_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+                   purities,
+                   delimiter=',')
+
+        logger.info("Compute two-point correlation for largest bond dimension")
+        correlations = np.zeros(L - 1)
+        for k in range(L - 1):
+            correlations[k] = compute_correlation(mps, number_mpo, r0=0, r1=k + 1)
+        logger.info(f"{correlations=}")
+        np.savetxt(OUTPUT_PATH + f"correlations_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+                   correlations,
+                   delimiter=',')
+
+        logger.info("Compute density-density correlation for largest bond dimension")
+        dens_dens_corr = np.zeros(L - 1)
+        for k in range(L - 1):
+            dens_dens_corr[k] = compute_dens_dens_corr(mps, number_mpo, r=k + 1)
+        logger.info(f"{dens_dens_corr=}")
+        np.savetxt(OUTPUT_PATH + f"dens_dens_corr_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+                   dens_dens_corr,
+                   delimiter=',')
+
+        logger.info("Compute half-chain entanglement spectrum for largest bond dimension")
+        entanglement_spectrum = compute_entanglement_spectrum(mps)
+        logger.info(f"{entanglement_spectrum=}")
+        np.savetxt(OUTPUT_PATH + f"entanglement_spectrum_exc_L_{L}_D_{bond_dim}_O_{OMEGA}_{SLURM_ARRAY_JOB_ID}.txt",
+                   entanglement_spectrum,
+                   delimiter=',')
+
+        basis_0 = np.array([1, 0])
+        dark_state = construct_basis_mps(L, basis=[np.outer(basis_0, basis_0)] * L)
+        logger.info(f"Overlap with dark state: {compute_overlap(dark_state, mps)}")  #NOTE: negative?
+
+
+if __name__ == "__main__":
+    main()

src/hpc/jobs_script.sh

@@ -1,7 +1,7 @@
 #!/bin/bash
 
 # job name, will show up in squeue output
-#SBATCH --job-name=jobName
+#SBATCH --job-name=cp-stat-exc
 
 # mail to which notifications will be sent 
 #SBATCH [email protected]
@@ -22,22 +22,40 @@
 #SBATCH --mem-per-cpu=4096
 
 # runtime in HH:MM:SS format (DAYS-HH:MM:SS format)
-#SBATCH --time=0-03:00:00
+#SBATCH --time=3-00:00:00
 
 # file to which standard output will be written (%A --> jobID, %a --> arrayID)
-#SBATCH --output=logFiles/simulationName_%A_%a_$(date +\%d\%m\%y\%H\%M).out
+#SBATCH --output=/scratch/nguyed99/qcp-1d/logging/cp_stat_exc_%A_%a.out
 
 # file to which standard errors will be written (%A --> jobID, %a --> arrayID)
-#SBATCH --error=logFiles/simulationName_%A_%a_$(date +\%d\%m\%y\%H\%M).err
+#SBATCH --error=/scratch/nguyed99/qcp-1d/logging/cp_stat_exc_%A_%a.err
+
+# job arrays
+#SBATCH --array=0-19
 
 # select partition
 #SBATCH --partition=main
 
 # set specific queue/nodes
 # SBATCH --reservation=bqa
 
+# simulation parameter
+L=50
+OMEGAS=(0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5)
+BOND_DIMS=(35 50)
+OMEGA_INDEX=$((SLURM_ARRAY_TASK_ID / 2))
+BOND_DIM_INDEX=$((SLURM_ARRAY_TASK_ID % 2))
+OMEGA=${OMEGAS[OMEGA_INDEX]}
+BOND_DIM=${BOND_DIMS[BOND_DIM_INDEX]}
+
+# paths and file names
+timestamp=$(date +'%Y-%m-%d-%H-%M-%S')
+LOG_PATH="/scratch/nguyed99/qcp-1d/logging"
+
 # store job info in output file, if you want...
 scontrol show job $SLURM_JOBID
+echo "slurm task ID = $SLURM_ARRAY_TASK_ID"
+echo $OMEGA $BOND_DIM
 
 export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK
 export OPENBLAS_NUM_THREADS=$SLURM_CPUS_PER_TASK
@@ -49,37 +67,9 @@ export NUMEXPR_NUM_THREADS=$SLURM_CPUS_PER_TASK
 # 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')
-
-# functions for logging memory usage and time
-log_memory_cpu_usage() {
-        while true; do
-                echo -n "$(date +'%Y-%m-%d %H:%M:%S') " >> "$LOG_PATH/usage_log_L_${L}_${timestamp}.txt"
-                # only include the used and total memory
-                free -h | awk '/^Mem:/ { print $3 "/" $2 }' >> "$LOG_PATH/usage_log_L_${L}_${timestamp}.txt"
-                top -bn1 | awk '/^%Cpu/ { print "CPU: " $2 " us, " $4 " sy" }' >> "$LOG_PATH/usage_log_L_${L}_${timestamp}.txt"
-                sleep 900
-        done
-}
-
-log_time() {
-        local start=$(date +%s)
-        $@
-        local end=$(date +%s)
-        local runtime=$((end - start))
-        echo "$(date +'%Y-%m-%d %H:%M:%S') Time taken for $@: ${runtime}s" >> "$LOG_PATH/time_log_L_${L}_${timestamp}.txt"
-}
-
 # activate virtualenv
 source /scratch/nguyed99/tensor/bin/activate
 
- # launch Python script
-log_memory_cpu_usage & 
-LOG_PID=$!
-
-log_time python3 contact_process_stat_exc.py > "$LOG_PATH/contact_process_stat_exc_L_${L}_${timestamp}.out"
-
-kill $LOG_PID
+export PYTHONPATH=$PYTHONPATH:/scratch/nguyed99/qcp-1d
+echo "Output log" >> "$LOG_PATH/${timestamp}.log"
+python3 contact_process_stat_exc.py $OMEGA $BOND_DIM $SLURM_ARRAY_JOB_ID 2>&1 | awk -v task_id=$SLURM_ARRAY_TASK_ID '{print "array task " task_id, $0}' >> "$LOG_PATH/${timestamp}.log"