Add further plotting scripts

scripts/DSF_iK.py

@@ -1,3 +1,4 @@
+#! /usr/bin/env python
 import argparse
 import plotly.graph_objects as go
 import numpy
@@ -16,5 +17,3 @@ y = [d for d in dsf]
 fig = go.Figure(data=go.Scatter(x=x, y=y))
 fig.update_layout(title=args.dsffile.rpartition('/')[2])
 fig.show()
-
-

scripts/build_LiebLin_catalogue_k_fixed.sh

@@ -11,28 +11,28 @@
 
 
 if [[ $# -ne 4 ]]; then
-	echo "Arguments needed: whichDSF, kBT, target_sumrule, nr_minutes (for each c, k combination)."
-	exit 1
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, nr_minutes (for each c, k combination)."
+    exit 1
 fi
 
 if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
-	echo "Only the d, g and o scanning options are implemented."
-	exit 1
+    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
+    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
+    echo "Requirement: 0 < target_sumrule < 1."
+    exit 1
 fi
 
 if [[ $4 -lt 1 ]]; then
-	echo "Please give a positive nr_minutes."
-	exit 1
+    echo "Please give a positive nr_minutes."
+    exit 1
 fi
 
 zmodload zsh/datetime
@@ -52,15 +52,14 @@ basedir="$(pwd)"
 
 for c in $clist
 do
-	for nk in {1..16}
-	do
-		echo 'Starting run for c =' $c', kfact = '$nk | tee -a $logfile
-		dir='c_'$c'/kBT_'$kBT'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule
-		echo $dir
-		mkdir -p $dir
-		cd $dir
-		LiebLin_Catalogue_Fixed_c_k_Nscaling $whichDSF $c $nk $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
+    for nk in {1..16}
+    do
+	echo 'Starting run for c =' $c', kfact = '$nk | tee -a $logfile
+	dir='c_'$c'/kBT_'$kBT'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule
+	mkdir -p $dir
+	cd $dir
+	LiebLin_Catalogue_Fixed_c_k_Nscaling $whichDSF $c $nk $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
 done

scripts/build_LiebLin_catalogue_k_fixed_dsf.sh

@@ -0,0 +1,68 @@
+#! /bin/zsh
+
+# We use zsh here to support floats (not supported in bash)
+
+# This script produces the .dsf files for the runs from `build_LiebLin_catalogue_k_fixed`
+# runs which achieved the required sumrule saturation.
+
+if [[ $# -ne 4 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, 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 "width must be > 0."
+    exit 1
+fi
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+width=$4
+
+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
+    for nk in {1..16}
+    do
+	# ommin is set to zero for all
+	ommin=0
+	# ommax is set to about 2x the Tonks-Girardeau dispersion at this momentum:
+	# w_TG = k^2 + 2\pi n |k|  with k=2pi*iK/L and iK=nk*N/8, so (with N=L) k=pi/4 nk and
+	# w_TG = pi^2/16 nk^2 + pi^2/2 nk = pi^2 nk (nk/16 + 1/2)
+	# We take pi^2 -> 8 for convenience, so 2*w_TG is about
+	ommax=$(($nk*($nk + 8)))
+	# Number of points:
+	Nom=1000
+	echo 'Starting computation of DSFs for c =' $c', kfact = '$nk | tee -a $logfile
+	dir='c_'$c'/kBT_'$kBT'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule
+	cd $dir
+	echo $dir
+	mkdir -p plots/
+	LiebLin_Catalogue_Fixed_c_k_Nscaling_Smoothen_DSF $whichDSF $c $nk $kBT $target_sumrule $ommin $ommax $Nom $width
+	mv *dsfs plots/
+	mv K_* plots/
+	mv Omega_* plots/
+	cd $basedir
+	echo ' Successfully computed DSFs for c =' $c', kfact = '$nk'\n' | tee -a $logfile
+    done
+done

scripts/build_LiebLin_catalogue_k_fixed_plot_dsf.sh

@@ -0,0 +1,69 @@
+#! /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_fixed`
+# runs which achieved the required sumrule saturation.
+
+if [[ $# -ne 3 ]]; then
+    echo "Arguments needed: whichDSF, kBT, target_sumrule."
+    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
+
+
+
+whichDSF=$1
+kBT=$2
+target_sumrule=$3
+width=$4
+
+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
+    for nk in {1..16}
+    do
+	echo 'Starting plotting of DSFs for c =' $c', kfact = '$nk | tee -a $logfile
+	dir='c_'$c'/kBT_'$kBT'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/plots'
+	cd $dir
+	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 *
+	do
+	    if [[ $file = LiebLin*dsfs ]]; then
+		echo 'Found file ' $file
+		# invoke the python script
+		plot_dsf_k_fixed.py $omegafile $file
+	    fi
+	done
+
+	cd $basedir
+	echo ' Successfully plotted DSFs for c =' $c', kfact = '$nk'\n' | tee -a $logfile
+    done
+done

scripts/plot_benchmarks.py

@@ -1,5 +1,5 @@
 """
-Plots benchmark files produced from Benchmark_RAW_File.
+Plot benchmark files produced from Benchmark_RAW_File.
 
 Usage: python plot_benchmarks.py [benchmark file name].
 

scripts/plot_dsf_k_fixed.py

@@ -0,0 +1,36 @@
+#! /usr/bin/env python
+
+"""
+Plot fixed momentum DSF.
+
+Usage: python plot_dsf_k_fixed.py [omega file] [dsf file]
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+omegafile = str(sys.argv[1])
+dsffile = str(sys.argv[2])
+
+omega = np.loadtxt(omegafile)
+dsf = np.loadtxt(dsffile)
+
+plt.plot(omega, dsf)
+
+plt.xlabel('$\omega$')
+plt.ylabel('$S (k, \omega)$')
+
+elements = dsffile.split('_')
+
+c_int = elements[4]
+L = elements[6]
+N = elements[8]
+iK = elements[13]
+width = elements[23].partition('.')[0]
+
+rho = int(N)/int(L)
+kokF = int(iK)*0.5/int(L)
+plt.title(f'c={c_int}, rho={rho} (N={N}), k={kokF}k_F, w={width}')
+
+plt.savefig(dsffile.replace('.', '_') + '.png')

src/EXECS/LiebLin_Catalogue_Fixed_c_k_Nscaling_Smoothen_DSF.cc

@@ -0,0 +1,109 @@
+/**********************************************************
+
+This software is part of J.-S. Caux's ABACUS library.
+
+Copyright (c) J.-S. Caux.
+
+-----------------------------------------------------------
+
+File:  LiebLin_Catalogue_Fixed_c_k_Nscaling_Smoothen_DSF.cc
+
+Purpose:  Produces smoothened DSF for output from LiebLin_Catalogue_Fixed_c_k_Nscaling.cc
+
+***********************************************************/
+
+#include <omp.h>
+#include "ABACUS.h"
+
+using namespace std;
+using namespace ABACUS;
+
+
+int main(int argc, char* argv[])
+{
+
+  if (argc != 10) { // provide some info
+
+    cout << endl << "Welcome to ABACUS\t(copyright J.-S. Caux)." << endl;
+    cout << endl << "Usage of LiebLin_Catalogue_Fixed_c_k_Nscaling_Smoothen_DSF 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 kfact \t\t momentum factor: momemntum 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 << "DP ommin" << endl << "DP ommax \t\t Min and max frequencies to cover in smoothened DSF" << endl;
+    cout << "Nom \t\t\t Number of frequency points used for discretization" << endl;
+    cout << "DP width \t\t Gaussian width used in smoothing, in units of two-particle level spacing" << 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 DiK = 1; // only do fixed momentum
+    DP ommin = atof(argv[ia++]);
+    DP ommax = atof(argv[ia++]);
+    int Nom = atoi(argv[ia++]);
+    DP width = atof(argv[ia++]);
+
+    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 = (kfact * N)/8;
+      int iKmax = iKmin;
+      DP srsat = 0.0;
+
+      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;
+      }
+      else {
+	srcfile >> srsat;
+      }
+      srcfile.close();
+
+      if (srsat > target_sumrule) {
+	// Compute the DSF
+	stringstream filenameprefix;
+	Data_File_Name (filenameprefix, whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "");
+	string prefix = filenameprefix.str();
+
+	DP normalization = twoPI * L;
+
+	Write_K_File (L, iKmin, iKmax);
+	Write_Omega_File (Nom, ommin, ommax);
+
+	// We use the scaled width function as default:
+	DP sumcheck;
+	sumcheck = Smoothen_RAW_into_SF_LiebLin_Scaled (prefix, L, N, iKmin, iKmax, DiK, ommin, ommax, Nom, width, normalization);
+      }
+
+    } // for N
+
+  } // else if arguments given OK
+
+  return(0);
+}