Add scripts and functions for catalogue over k range

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)

scripts/build_LiebLin_catalogue_k_fixed_dsfs_plots.sh

@@ -52,25 +52,32 @@ do
 	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 *
-	do
-	    if [[ $file = Omega* ]]; then
-		omegafile=$file
-	    fi
-	done
-	echo 'Omega file: ' $omegafile
-	# For each .dsfs file, plot
-	for file in *.dsfs
+
+	for datadir in ../data/*
 	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_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);
+}