Add scripts and functions for catalogue over k range

scripts/build_LiebLin_c_scan_k_fixed.sh

@@ -78,7 +78,7 @@ basedir="$(pwd)"
 for c in $clist
 do
     echo '** Starting run for c =' $c', N = '$N
-    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/N_'$N'/data/store/c_'$c
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/data/N_'$N'/c_'$c
     mkdir -p $dir
     cd $dir
     LiebLin_DSF $whichDSF $c $N $N $iK $iK $kBT $Max_Secs $target_sumrule $refine

scripts/build_LiebLin_c_scan_k_fixed_dsfs.sh

@@ -84,7 +84,7 @@ basedir="$(pwd)"
 for c in $clist
 do
     echo '** Starting computation of DSFs for c =' $c', N = '$N
-    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/N_'$N'/data/store/c_'$c
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/data/N_'$N'/c_'$c
     cd $dir
     Smoothen_LiebLin_DSF_Scaled $whichDSF $c $N $N $iK $iK $kBT 0 $ommin $ommax $Nom $width
     cd $basedir

scripts/build_LiebLin_c_scan_k_fixed_dsfs_plot.sh

@@ -82,9 +82,9 @@ done
 basedir="$(pwd)"
 
 echo '** Starting run for N = '$N
-dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/N_'$N'/data/plots'
+dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/plots'
 mkdir -p $dir
 cd $dir
-plot_c_scan_k_fixed_dsfs.py $width
+plot_c_scan_k_fixed_dsfs.py $N $width
 cd $basedir
 echo ' Successfully computed DSFs for N = '$N'.\n'

scripts/build_LiebLin_c_scan_k_range.sh

@@ -0,0 +1,87 @@
+#! /bin/zsh
+
+# For given N (L), this
+# computes the required DSF (up to required sumrule)
+# over a range of values of c_int.
+
+if [[ $# -ne 5 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, N, nkmax (max momentum in units of kF/4)."
+    exit 1
+fi
+
+if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
+    echo "Only the d, g and o scanning options are implemented."
+    exit 1
+fi
+
+if [[ $2 -lt 0 ]]; then
+    echo "Temperature kBT must be > 0."
+    exit 1
+fi
+
+if [[ $3 -lt 0 || $3 -gt 1 ]]; then
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
+fi
+
+if [[ $4 -lt 0 ]]; then
+    echo "N must be > 0."
+    exit 1
+fi
+
+if [[ $5 -lt 0 ]]; then
+    echo "nkmax must be > 0."
+    exit 1
+fi
+
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+N=$4
+nkmax=$5
+
+correlator='rho-rho'
+if [[ $whichDSF == 'o' ]]; then
+    correlator='psidag-psi'
+elif [[ $whichDSF == 'g' ]]; then
+    correlator='psi-psidag'
+fi
+
+iKmax=$(($nkmax * $N/8))
+
+Max_Secs=3600
+refine=0
+
+#clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
+
+# List of interactions: fractional powers of 2
+clist_raw=()
+for nc in {-128..128}
+do
+    clist_raw=($clist_raw $(( 4 * 2**($nc/16.) )))
+done
+# Now cast the integer values to true integers
+zmodload zsh/mathfunc
+clist=()
+for c in $clist_raw
+do
+    if [[ $((floor($c))) == $((ceil($c))) ]]; then
+	clist=($clist $((int($c))))
+    else
+	clist=($clist $c)
+    fi
+done
+
+basedir="$(pwd)"
+
+for c in $clist
+do
+    echo '** Starting run for c =' $c', N = '$N
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_range/k_0_to_'${(l:2::0:)nkmax}'kFo4/sr_'$target_sumrule'/store/data/N_'$N'/c_'$c
+    mkdir -p $dir
+    cd $dir
+    LiebLin_DSF $whichDSF $c $N $N 0 $iKmax $kBT $Max_Secs $target_sumrule $refine
+    cd $basedir
+    echo ' Successfully computed DSFs for c =' $c', N = '$N'.\n'
+done

scripts/build_LiebLin_c_scan_k_range_dsfs.sh

@@ -0,0 +1,92 @@
+#! /bin/zsh
+
+# For given N (L) this
+# computes the smoothened dsfs for a `build_LiebLin_c_scan_k_range` run.
+
+if [[ $# -ne 6 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, N, nkmax (max momentum in units of kF/4), width."
+    exit 1
+fi
+
+if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
+    echo "Only the d, g and o scanning options are implemented."
+    exit 1
+fi
+
+if [[ $2 -lt 0 ]]; then
+    echo "Temperature kBT must be > 0."
+    exit 1
+fi
+
+if [[ $3 -lt 0 || $3 -gt 1 ]]; then
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
+fi
+
+if [[ $4 -lt 0 ]]; then
+    echo "N must be > 0."
+    exit 1
+fi
+
+if [[ $5 -lt 0 ]]; then
+    echo "nkmax must be > 0."
+    exit 1
+fi
+
+if [[ $6 -lt 0 ]]; then
+    echo "width must be > 0."
+    exit 1
+fi
+
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+N=$4
+nkmax=$5
+width=$6
+
+correlator='rho-rho'
+if [[ $whichDSF == 'o' ]]; then
+    correlator='psidag-psi'
+elif [[ $whichDSF == 'g' ]]; then
+    correlator='psi-psidag'
+fi
+
+iKmax=$(($nkmax * $N/8))
+
+ommin=0
+ommax=$(($nkmax*($nkmax - 8))) # use lower branch of 1ph continuum
+Nom=1000
+
+#clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
+
+# List of interactions: fractional powers of 2
+clist_raw=()
+for nc in {-128..128}
+do
+    clist_raw=($clist_raw $(( 4 * 2**($nc/16.) )))
+done
+# Now cast the integer values to true integers
+zmodload zsh/mathfunc
+clist=()
+for c in $clist_raw
+do
+    if [[ $((floor($c))) == $((ceil($c))) ]]; then
+	clist=($clist $((int($c))))
+    else
+	clist=($clist $c)
+    fi
+done
+
+basedir="$(pwd)"
+
+for c in $clist
+do
+    echo '** Starting computation of DSFs for c =' $c', N = '$N
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_range/k_0_to_'${(l:2::0:)nkmax}'kFo4/sr_'$target_sumrule'/store/data/N_'$N'/c_'$c
+    cd $dir
+    Smoothen_LiebLin_DSF_Scaled $whichDSF $c $N $N 0 $iKmax $kBT 0 $ommin $ommax $Nom $width
+    cd $basedir
+    echo ' Successfully computed DSFs for c =' $c', N = '$N'.\n'
+done

scripts/build_LiebLin_c_scan_k_range_dsfs_plot.sh

@@ -0,0 +1,90 @@
+#! /bin/zsh
+
+# For given N (L) this
+# produces an interactive (tunable c_int) plot for S(k, \omega)
+# from data produced from `build_LiebLin_c_scan_k_range[_dsfs]` scripts.
+
+if [[ $# -ne 6 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, N, nkmax (max momentum in units of kF/4), width."
+    exit 1
+fi
+
+if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
+    echo "Only the d, g and o scanning options are implemented."
+    exit 1
+fi
+
+if [[ $2 -lt 0 ]]; then
+    echo "Temperature kBT must be > 0."
+    exit 1
+fi
+
+if [[ $3 -lt 0 || $3 -gt 1 ]]; then
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
+fi
+
+if [[ $4 -lt 0 ]]; then
+    echo "N must be > 0."
+    exit 1
+fi
+
+if [[ $5 -lt 0 ]]; then
+    echo "nkmax must be > 0."
+    exit 1
+fi
+
+if [[ $6 -lt 0 ]]; then
+    echo "width must be > 0."
+    exit 1
+fi
+
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+N=$4
+nkmax=$5
+width=$6
+
+correlator='rho-rho'
+if [[ $whichDSF == 'o' ]]; then
+    correlator='psidag-psi'
+elif [[ $whichDSF == 'g' ]]; then
+    correlator='psi-psidag'
+fi
+
+iKmax=$(($nkmax * $N/8))
+
+Max_Secs=3600
+refine=0
+
+#clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
+
+# List of interactions: fractional powers of 2
+clist_raw=()
+for nc in {-128..128}
+do
+    clist_raw=($clist_raw $(( 4 * 2**($nc/16.) )))
+done
+# Now cast the integer values to true integers
+zmodload zsh/mathfunc
+clist=()
+for c in $clist_raw
+do
+    if [[ $((floor($c))) == $((ceil($c))) ]]; then
+	clist=($clist $((int($c))))
+    else
+	clist=($clist $c)
+    fi
+done
+
+basedir="$(pwd)"
+
+echo '** Starting run for N = '$N
+dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_scan/k_range/k_0_to_'${(l:2::0:)nkmax}'kFo4/sr_'$target_sumrule'/store/plots'
+mkdir -p $dir
+cd $dir
+plot_c_scan_k_range_dsfs.py $N $width
+cd $basedir
+echo ' Successfully computed DSFs for N = '$N'.\n'

scripts/build_LiebLin_catalogue_k_fixed.sh

@@ -50,7 +50,7 @@ elif [[ $whichDSF == 'g' ]]; then
 fi
 
 mkdir -p logs
-logfile='logs/run_'$whichDSF'_kBT_'$kBT'_sr_'$target_sumrule'_'$(date '+%Y-%m-%d-%Hh%M')'.log'
+logfile='logs/run_'$whichDSF'_kBT_'$kBT'_k_fixed_sr_'$target_sumrule'_'$(date '+%Y-%m-%d-%Hh%M')'.log'
 touch $logfile
 
 clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
@@ -63,7 +63,7 @@ do
     for nk in {1..16}
     do
 	echo '** Starting run for c =' $c', kfact = '$nk | tee -a $logfile
-	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/prelim'
+	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_fixed/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/prelim'
 	mkdir -p $dir
 	cd $dir
 	LiebLin_Catalogue_Fixed_c_k_Nscaling $whichDSF $c $nk $kBT $target_sumrule 0 $nr_minutes | tee -a $basedir/$logfile

scripts/build_LiebLin_catalogue_k_fixed_dsfs.sh

@@ -61,7 +61,7 @@ do
 	# Number of points:
 	Nom=1000
 	echo 'Starting computation of DSFs for c =' $c', kfact = '$nk
-	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/data'
+	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_fixed/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/data'
 	cd $dir
 	for datadir in *
 	do

scripts/build_LiebLin_catalogue_k_fixed_dsfs_plots.sh

@@ -49,28 +49,35 @@ do
     for nk in {1..16}
     do
 	echo '** Starting plotting of DSFs for c =' $c', k = '$nk' kFo4.'
-	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/plots'
+	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_fixed/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/plots'
 	mkdir -p $dir
 	cd $dir
-	# Move all K_, Omega_ and .dsfs files to the plots directory
-	mv ../data/*/K_* .
-	mv ../data/*/Omega_* .
-	mv ../data/*/*dsfs .
-	dirlist="$(ls)"
-	# Find the Omega file
-	for file in *
+
+	for datadir in ../data/*
 	do
-	    if [[ $file = Omega* ]]; then
-		omegafile=$file
-	    fi
-	done
-	echo 'Omega file: ' $omegafile
-	# For each .dsfs file, plot
-	for file in *.dsfs
-	do
-	    echo 'Found file ' $file
-	    # invoke the python script
-	    plot_dsf_k_fixed.py $omegafile $file
+	    N=${datadir#"../data/N_"}
+	    mkdir -p 'N_'$N
+	    cd 'N_'$N
+	    # Move all K_, Omega_ and .dsfs files here
+	    mv ../../data/N_$N/K_* .
+	    mv ../../data/N_$N/Omega_* .
+	    mv ../../data/N_$N/*dsfs .
+
+	    # Find the K and Omega file
+	    for file in *
+	    do
+		if [[ $file = Omega* ]]; then
+		    omegafile=$file
+		fi
+	    done
+	    # For each .dsfs file, plot
+	    for file in *.dsfs
+	    do
+		echo 'Found file ' $file
+		# invoke the python script
+		plot_dsf_k_fixed.py $omegafile $file
+	    done
+	    cd ..
 	done
 
 	cd $basedir

scripts/build_LiebLin_catalogue_k_range.sh

@@ -0,0 +1,76 @@
+#! /bin/zsh
+
+# We use zsh here to support floats (not supported in bash)
+
+# This script produces a whole stack of subdirectories,
+# for a range of values of c centered on 4 with 8 up/down factors of 2,
+# and momenta in range from 0 to 4kF.
+
+# For each c, the LiebLin_Catalogue_Fixed_c_Nscaling
+# executable is invoked.
+
+
+if [[ $# -ne 5 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, nkmax (max momentum in units of kF/4), nr_minutes (for each c value)."
+    exit 1
+fi
+
+if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
+    echo "Only the d, g and o scanning options are implemented."
+    exit 1
+fi
+
+if [[ $2 -lt 0 ]]; then
+    echo "Temperature kBT must be > 0."
+    exit 1
+fi
+
+if [[ $3 -lt 0 || $3 -gt 1 ]]; then
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
+fi
+
+if [[ $4 -lt 0 ]]; then
+    echo "nkmax must be > 0."
+    exit 1
+fi
+
+if [[ $5 -lt 1 ]]; then
+    echo "Please give a positive nr_minutes."
+    exit 1
+fi
+
+zmodload zsh/datetime
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+nkmax=$4
+nr_minutes=$5
+
+correlator='rho-rho'
+if [[ $whichDSF == 'o' ]]; then
+    correlator='psidag-psi'
+elif [[ $whichDSF == 'g' ]]; then
+    correlator='psi-psidag'
+fi
+
+mkdir -p logs
+logfile='logs/run_'$whichDSF'_kBT_'$kBT'_k_range_sr_'$target_sumrule'_'$(date '+%Y-%m-%d-%Hh%M')'.log'
+touch $logfile
+
+clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
+
+basedir="$(pwd)"
+
+
+for c in $clist
+do
+    echo '** Starting run for c =' $c | tee -a $logfile
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_fixed/c_'$c'/k_range/k_0_to_'$nkmax'kFo4/sr_'$target_sumrule'/prelim'
+    mkdir -p $dir
+    cd $dir
+    LiebLin_Catalogue_Fixed_c_Nscaling $whichDSF $c $nkmax $kBT $target_sumrule 0 $nr_minutes | tee -a $basedir/$logfile
+    cd $basedir
+    echo '** Successfully completed run for c =' $c', kfact = '$nk'.\n' | tee -a $logfile
+done

scripts/build_LiebLin_catalogue_k_range_dsfs.sh

@@ -0,0 +1,74 @@
+#! /bin/zsh
+
+# We use zsh here to support floats (not supported in bash)
+
+# This script produces the .dsfs files for the runs from `build_LiebLin_catalogue_k_range`
+# runs which achieved the required sumrule saturation.
+
+if [[ $# -ne 5 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, nkmax (max momentum in units of kF/4), width."
+    exit 1
+fi
+
+if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
+    echo "Only the d, g and o scanning options are implemented."
+    exit 1
+fi
+
+if [[ $2 -lt 0 ]]; then
+    echo "Temperature kBT must be > 0."
+    exit 1
+fi
+
+if [[ $3 -lt 0 || $3 -gt 1 ]]; then
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
+fi
+
+if [[ $4 -lt 0 ]]; then
+    echo "nkmax must be > 0."
+    exit 1
+fi
+
+if [[ $5 -lt 0 ]]; then
+    echo "width must be > 0."
+    exit 1
+fi
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+nkmax=$4
+width=$5
+
+correlator='rho-rho'
+if [[ $whichDSF == 'o' ]]; then
+    correlator='psidag-psi'
+elif [[ $whichDSF == 'g' ]]; then
+    correlator='psi-psidag'
+fi
+
+clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
+
+basedir="$(pwd)"
+
+ommin=0
+ommax=$(($nkmax*($nkmax - 8))) # use lower branch of 1ph continuum
+Nom=1000
+
+for c in $clist
+do
+    echo 'Starting computation of DSFs for c =' $c
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_fixed/c_'$c'/k_range/k_0_to_'${(l:2::0:)nkmax}'kFo4/sr_'$target_sumrule'/store/data'
+    cd $dir
+    for datadir in *
+    do
+	cd $datadir
+	N=${datadir#"N_"}
+	iKmax=$(($nkmax * $N/8))
+	Smoothen_LiebLin_DSF_Scaled $whichDSF $c $N $N 0 $iKmax $kBT 0 $ommin $ommax $Nom $width
+	cd ..
+    done
+    cd $basedir
+    echo ' Successfully computed DSFs for c =' $c
+done

scripts/build_LiebLin_catalogue_k_range_dsfs_plots.sh

@@ -0,0 +1,95 @@
+#! /bin/zsh
+
+# We use zsh here to support floats (not supported in bash)
+
+# This script produces the plots for the runs from `build_LiebLin_catalogue_k_range`
+# runs which achieved the required sumrule saturation.
+
+if [[ $# -ne 5 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, nkmax (max momentum in units of kF/4), width."
+    exit 1
+fi
+
+if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
+    echo "Only the d, g and o scanning options are implemented."
+    exit 1
+fi
+
+if [[ $2 -lt 0 ]]; then
+    echo "Temperature kBT must be > 0."
+    exit 1
+fi
+
+if [[ $3 -lt 0 || $3 -gt 1 ]]; then
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
+fi
+
+if [[ $4 -lt 0 ]]; then
+    echo "nkmax must be > 0."
+    exit 1
+fi
+
+if [[ $5 -lt 0 ]]; then
+    echo "width must be > 0."
+    exit 1
+fi
+
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+nkmax=$4
+width=$5
+
+correlator='rho-rho'
+if [[ $whichDSF == 'o' ]]; then
+    correlator='psidag-psi'
+elif [[ $whichDSF == 'g' ]]; then
+    correlator='psi-psidag'
+fi
+
+clist=(1024 512 256 128 64 32 16 8 4 2 1 0.5 0.25 0.125 0.0625 0.03125 0.015625)
+
+basedir="$(pwd)"
+
+
+for c in $clist
+do
+    echo '** Starting plotting of DSFs for c =' $c
+    dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_fixed/c_'$c'/k_range/k_0_to_'${(l:2::0:)nkmax}'kFo4/sr_'$target_sumrule'/store/plots'
+    mkdir -p $dir
+    cd $dir
+    for datadir in ../data/*
+    do
+	N=${datadir#"../data/N_"}
+	mkdir -p 'N_'$N
+	cd 'N_'$N
+	# Move all K_, Omega_ and .dsfs files here
+	mv ../../data/N_$N/K_* .
+	mv ../../data/N_$N/Omega_* .
+	mv ../../data/N_$N/*dsfs .
+
+	# Find the K and Omega file
+	for file in *
+	do
+	    if [[ $file = Omega* ]]; then
+		omegafile=$file
+	    elif [[ $file = K_* ]]; then
+		kfile=$file
+	    fi
+	done
+	# For each .dsfs file, plot
+	for file in *.dsfs
+	do
+	    echo 'Found file ' $file
+	    # invoke the python script
+	    plot_dsf_k_range.py $kfile $omegafile $file
+	done
+	cd ..
+    done
+
+    cd $basedir
+    echo '** Successfully plotted DSFs for c =' $c'.\n'
+
+done

scripts/plot_c_scan_k_fixed_dsfs.py

@@ -13,15 +13,16 @@ import numpy as np
 import os
 import sys
 
-width = str(sys.argv[1])
+N = str(sys.argv[1])
+width = str(sys.argv[2])
 
 
 # Get the list of interactions which have been computed
-dirlist = os.listdir('../store/')
+dirlist = os.listdir('../data/N_%s' % N)
 clist = sorted([float(c.lstrip('c_')) for c in dirlist])
 
 # Get the Omega file
-omegafile = glob.glob('../store/c_%s*/Omega*' % str(clist[0]).rstrip('.0')[:12])[0]
+omegafile = glob.glob('../data/N_%s/c_%s*/Omega*' % (N, str(clist[0]).rstrip('.0')[:12]))[0]
 omega = np.loadtxt(omegafile)
 
 # Load all the available dsfs from the data store
@@ -31,21 +32,32 @@ for c in clist:
     # first try for exact match, stripping '.0' to treat integer values, e.g. 4.0 into 4
     # and then (if it doesn't work) use only the first 12 characters, plus wildcard, to cover rounding errors
     try:
-        dsffile = glob.glob('../store/c_%s/*_w_%s.dsfs' % (str(c).rstrip('.0'), width))[0]
+        dsffile = glob.glob('../data/N_%s/c_%s/*_w_%s.dsfs' % (N, str(c).rstrip('.0'), width))[0]
     except IndexError:
-        dsffile = glob.glob('../store/c_%s*/*_w_%s.dsfs' % (str(c).rstrip('.0')[:12], width))[0]
+        dsffile = glob.glob('../data/N_%s/c_%s*/*_w_%s.dsfs' % (N, str(c).rstrip('.0')[:12], width))[0]
     dsfs[str(c)] = np.loadtxt(dsffile)
 
+# Read some useful parameters from (last) file name:
+elements = dsffile.rpartition('/')[2].split('_')
+L = elements[5]
+#N = elements[7]
+iK = elements[12]
+width = elements[20].rpartition('.')[0]
+# which leads to
+rho = int(N)/int(L)
+kokF = int(iK)*0.5/int(L)
+
+
 fig, ax = plt.subplots()
 
 ax.set_xlim(omega[0], omega[-1])
 # To determine the y axis limits, we look at the lowest value of c (sharpest peak)
 dsfsmax = max(dsfs[str(clist[0])])
-print('dsfsmax = ', dsfsmax)
 
-ymax = 1.2 * dsfsmax
+ymax = 1.01 * dsfsmax
+ymax=4
 xtext = 0.6 * omega[-1]
-ytext = 1.1 * dsfsmax
+ytext = 0.9 * ymax
 
 ax.set_ylim([0, ymax])
 ax.text(xtext, ytext, f'c = {clist[0]}')
@@ -54,6 +66,7 @@ def animate(i):
     ax.clear()
     ax.set_xlim(omega[0], omega[-1])
     ax.set_ylim([0, ymax])
+    ax.set_title(f'c scan, rho={rho} (N={N}), k={kokF}k_F, w={width}')
     ax.text(xtext, ytext, f'c = {clist[i]:10.6f}')
     ax.plot(omega, dsfs[str(clist[i])])
 

scripts/plot_c_scan_k_range_dsfs.py

@@ -0,0 +1,83 @@
+#! /usr/bin/env python
+
+"""
+Plot DSF.
+
+Usage: python plot_dsf_k_range.py [omega file] [dsf file]
+"""
+
+import glob
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+import numpy as np
+import os
+import sys
+
+N = str(sys.argv[1])
+width = str(sys.argv[2])
+
+
+# Get the list of interactions which have been computed
+dirlist = os.listdir('../data/N_%s' % N)
+clist = sorted([float(c.lstrip('c_')) for c in dirlist])
+
+# Get the K file
+kfile = glob.glob('../data/N_%s/c_%s*/K_*' % (N, str(clist[0]).rstrip('.0')[:12]))[0]
+k = np.loadtxt(kfile)
+
+# Get the Omega file
+omegafile = glob.glob('../data/N_%s/c_%s*/Omega*' % (N, str(clist[0]).rstrip('.0')[:12]))[0]
+omega = np.loadtxt(omegafile)
+
+# Load all the available dsfs from the data store
+dsfs = {}
+for c in clist:
+    # Do some black magic here: for matching the interaction,
+    # first try for exact match, stripping '.0' to treat integer values, e.g. 4.0 into 4
+    # and then (if it doesn't work) use only the first 12 characters, plus wildcard, to cover rounding errors
+    try:
+        dsffile = glob.glob('../data/N_%s/c_%s/*_w_%s.dsfs' % (N, str(c).rstrip('.0'), width))[0]
+    except IndexError:
+        dsffile = glob.glob('../data/N_%s/c_%s*/*_w_%s.dsfs' % (N, str(c).rstrip('.0')[:12], width))[0]
+    dsfs[str(c)] = np.loadtxt(dsffile)
+
+# Read some useful parameters from (last) file name:
+elements = dsffile.rpartition('/')[2].split('_')
+L = elements[5]
+#N = elements[7]
+iKmax = elements[14]
+width = elements[22].rpartition('.')[0]
+# which leads to
+rho = int(N)/int(L)
+#kokF = int(iK)*0.5/int(L)
+
+
+fig, ax = plt.subplots()
+
+# To determine the y axis limits, we look at the lowest value of c (sharpest peak)
+dsfsmax = 0.5
+
+ymax = omega[-1]
+xtext = 0.05 * k[-1]
+ytext = 0.9 * ymax
+
+ax.text(xtext, ytext, f'c = {clist[0]:10.6f}', color='white', fontsize='large')
+pcm = ax.pcolormesh(k, omega, dsfs[str(clist[0])], vmax=dsfsmax)
+fig.colorbar(pcm, ax=ax)
+def animate(i):
+    ax.clear()
+    ax.set_title(f'c scan, rho={rho} (N={N}), w={width}', fontsize='x-large')
+    ax.text(xtext, ytext, f'c = {clist[i]:10.6f}', color='white', fontsize='x-large')
+    ax.pcolormesh(k, omega, dsfs[str(clist[i])], vmax=dsfsmax)
+
+ani = animation.FuncAnimation(fig, animate, frames=len(clist), interval=100, repeat=False)
+
+outfilename = (dsffile.partition('_c_')[0].rpartition('/')[2] + '_c_scan_' +
+               dsffile.partition('_c_')[2].partition('_')[2].partition('.dsfs')[0])
+
+ani.save(outfilename + '.mp4')
+
+with open(outfilename + '.html', 'w') as file:
+    file.write(ani.to_html5_video())
+
+#plt.show()

scripts/plot_dsf_k_range.py

@@ -0,0 +1,39 @@
+#! /usr/bin/env python
+
+"""
+Plot DSF.
+
+Usage: python plot_dsf_k_range.py [k file] [omega file] [dsf file]
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+kfile = str(sys.argv[1])
+omegafile = str(sys.argv[2])
+dsffile = str(sys.argv[3])
+
+k = np.loadtxt(kfile)
+omega = np.loadtxt(omegafile)
+dsf = np.loadtxt(dsffile)
+
+dsfsmax = 0.5
+
+plt.pcolormesh(k, omega, dsf, vmax=dsfsmax)
+plt.colorbar()
+plt.xlabel('$k$')
+plt.ylabel('$\omega$')
+
+elements = dsffile.split('_')
+
+c_int = elements[3]
+L = elements[5]
+N = elements[7]
+width = elements[22].partition('.')[0]
+
+rho = int(N)/int(L)
+
+plt.title(f'c={c_int}, rho={rho} (N={N}), w={width}')
+
+plt.savefig(dsffile.replace('.', '_') + '.png')

src/EXECS/LiebLin_Catalogue_Fixed_c_Nscaling.cc

@@ -0,0 +1,142 @@
+/**********************************************************
+
+This software is part of J.-S. Caux's ABACUS library.
+
+Copyright (c) J.-S. Caux.
+
+-----------------------------------------------------------
+
+File:  LiebLin_Catalogue_Fixed_c_Nscaling.cc
+
+Purpose:  Produces sets of data files for correlations, increasing system size at fixed c.
+
+***********************************************************/
+
+#include <omp.h>
+#include "ABACUS.h"
+
+using namespace std;
+using namespace ABACUS;
+
+
+int main(int argc, char* argv[])
+{
+
+  if (argc != 8) { // provide some info
+
+    cout << endl << "Welcome to ABACUS\t(copyright J.-S. Caux)." << endl;
+    cout << endl << "Usage of LiebLin_Catalogue_Fixed_c_k_Nscaling executable: " << endl;
+    cout << endl << "Provide the following arguments:" << endl << endl;
+    cout << "char whichDSF \t\t Which structure factor should be calculated ?  Options are:  "
+      "d for rho rho, g for psi psi{dagger}, o for psi{dagger} psi" << endl;
+    cout << "DP c_int \t\t Value of the interaction parameter:  use positive real values only" << endl;
+    cout << "int kfactmax \t\t max momentum factor: max momentum will be set to kfact * kF/4" << endl;
+    cout << "DP kBT \t\t Temperature (positive only of course)" << endl;
+    cout << "DP target_sumrule \t sumrule saturation you're satisfied with" << endl;
+    cout << "int Hrs \t\t Allowed computational time (hours)" << endl;
+    cout << "int Mns \t\t Allowed computational time (minutes)" << endl;
+  }
+
+  else { // correct nr of arguments
+    int ia = 1;
+    char whichDSF = *argv[ia++];
+    DP c_int = atof(argv[ia++]);
+    int kfact = atoi(argv[ia++]);
+    DP kBT = atof(argv[ia++]);
+    DP target_sumrule = atof(argv[ia++]);
+    int Max_Secs = 3600*atoi(argv[ia++]) + 60*atoi(argv[ia++]);
+
+    double StartTime = omp_get_wtime();
+
+    double ActualTime = omp_get_wtime();
+
+    int Secs_left = Max_Secs;
+
+    int iN = 0;
+
+    int nN = 16;
+    Vect<int> Nv(nN);
+    // Multiples of 32 up to 256
+    for (int i = 1; i <= 8; ++i) Nv[i-1] = 32*i;
+    // Then steps of 64 up to 512
+    for (int i = 1; i <= 4; ++i) Nv[7+i] = 256 + 64*i;
+    // Then steps of 128 up to 1024
+    for (int i = 1; i <= 4; ++i) Nv[11+i] = 512 + 128*i;
+
+    for (int iN = 0; iN < nN; ++iN) {
+
+      int N = Nv[iN];
+      DP L = N;
+      int iKmin = 0;
+      int iKmax = (kfact * N)/8;
+      DP srsat = 0.0;
+      bool refine = false;
+
+      stringstream SRC_stringstream;  string SRC_string;
+      Data_File_Name (SRC_stringstream, whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "");
+      SRC_stringstream << ".src";
+      SRC_string = SRC_stringstream.str();    const char* SRC_Cstr = SRC_string.c_str();
+
+      fstream srcfile;
+      srcfile.open(SRC_Cstr, fstream::in);
+      if (srcfile.fail()) {
+	srsat = 0.0;
+	refine = false;
+      }
+      else {
+	srcfile >> srsat;
+	refine = true;
+      }
+      srcfile.close();
+
+      ActualTime = omp_get_wtime();
+      Secs_left = int(Max_Secs - (ActualTime - StartTime));
+
+      Scan_Info resulting_info;
+      if (srsat < target_sumrule && Secs_left > Max_Secs/2) {
+	// Improve the icmin calculation by one chunk:
+	cout << "---\nTime left = " << Secs_left << " seconds." << endl;
+	if (srsat > 0) {
+	  cout << "Continue with N = " << N << ". Sumrule previously achieved: " << srsat << endl;
+	} else {
+	  cout << "Start with N = " << N << "." << endl;
+	}
+	resulting_info = Scan_LiebLin (whichDSF, c_int, L, N, iKmin, iKmax, kBT,
+				       Secs_left, target_sumrule, refine);
+	cout << "Done with N = " << N
+	     << ". Sumrule obtained: " << resulting_info.sumrule_obtained
+	     << endl;
+      }
+
+      if (resulting_info.sumrule_obtained > target_sumrule) {
+	// Move files to storage, keeping a copy of the .src file in the current directory
+	string command1 = "mkdir -p ../store/data/N_" + to_string(N);
+	system(command1.c_str());
+	string command2 = "mv *_N_" + to_string(N) + "* ../store/data/N_" + to_string(N) + "/";
+	system(command2.c_str());
+	string command3 = "cp ../store/data/N_" + to_string(N) + "/*src .";
+	system(command3.c_str());
+      }
+
+      ActualTime = omp_get_wtime();
+
+      Secs_left = int(Max_Secs - (ActualTime - StartTime));
+
+      if (Secs_left < 30) {
+	if (resulting_info.sumrule_obtained > target_sumrule) {
+	  cout << "---\nBreaking out after completing N = " << N
+	       << " since time left = " << Secs_left << " seconds." << endl;
+	}
+	else {
+	  cout << "---\nBreaking out while working on N = " << N
+	       << " since allocated time is exhausted." << endl;
+	}
+	break;
+      }
+
+    } // while there is time
+
+  } // else if arguments given OK
+
+  return(0);
+}

src/EXECS/LiebLin_Catalogue_Fixed_c_k_Nscaling.cc

@@ -30,7 +30,7 @@ int main(int argc, char* argv[])
     cout << "char whichDSF \t\t Which structure factor should be calculated ?  Options are:  "
       "d for rho rho, g for psi psi{dagger}, o for psi{dagger} psi" << endl;
     cout << "DP c_int \t\t Value of the interaction parameter:  use positive real values only" << endl;
-    cout << "int kfact \t\t momentum factor: momemntum will be set to kfact * kF/4" << endl;
+    cout << "int kfact \t\t momentum factor: momentum will be set to kfact * kF/4" << endl;
     cout << "DP kBT \t\t Temperature (positive only of course)" << endl;
     cout << "DP target_sumrule \t sumrule saturation you're satisfied with" << endl;
     cout << "int Hrs \t\t Allowed computational time (hours)" << endl;