Improve logging

ABACUS_Doxyfile

@@ -38,7 +38,7 @@ PROJECT_NAME           = ABACUS
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         =
+PROJECT_NUMBER         = 1.0.0
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -790,7 +790,7 @@ WARN_LOGFILE           =
 # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
 # Note: If this tag is empty the current directory is searched.
 
-INPUT                  = /Users/jscaux/WORK/_ABACUS/ABACUS
+INPUT                  =
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
@@ -889,7 +889,7 @@ EXCLUDE_SYMLINKS       = NO
 # Note that the wildcards are matched against the file with absolute path, so to
 # exclude all test directories for example use the pattern */test/*
 
-EXCLUDE_PATTERNS       =
+EXCLUDE_PATTERNS       = */test/*
 
 # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
 # (namespaces, classes, functions, etc.) that should be excluded from the

include/ABACUS.h

@@ -17,7 +17,7 @@ Purpose:  Core header file, includes all descendents.
 
 // This core header file includes all the others
 
-const char ABACUS_VERSION[20] = "ABACUS_0a";
+const char ABACUS_VERSION[20] = "1.0.0";
 
 // Standard includes
 #include <cmath>

include/ABACUS_Heis.h

@@ -21,7 +21,7 @@ namespace ABACUS {
 
   // First, some global constants...
 
-  const long long int ID_UPPER_LIMIT = 10000000LL;  // max size of vectors we can define without seg fault
+  const long long int ID_UPPER_LIMIT = 10000000LL;  // max size we can define without seg fault
   const int INTERVALS_SIZE = 100000;                // size of Scan_Intervals arrays
   const int NBASESMAX = 1000; // max number of bases kept
 
@@ -30,7 +30,7 @@ namespace ABACUS {
   // Cutoffs on particle numbers
   const int MAXSTRINGS = 20;      // maximal number of particle types we allow in bases
 
-  const int NEXC_MAX_HEIS = 16; // maximal number of excitations (string binding/unbinding, particle-hole) considered
+  const int NEXC_MAX_HEIS = 16; // max nr of excitations (string binding/unbinding, particle-hole) considered
 
 
   //***********************************************************************
@@ -66,7 +66,8 @@ namespace ABACUS {
 
   //****************************************************************************
 
-  // Objects in class Heis_Base are a checked vector containing the number of rapidities of allowable types for a given state
+  // Objects in class Heis_Base are a checked vector containing the number of rapidities
+  // of allowable types for a given state
 
   class Heis_Base {
 

include/ABACUS_Integ.h

@@ -116,7 +116,8 @@ namespace ABACUS {
 	if (i % 2 && bdry[i] < xmax_ref) xmax_reg++;
       }
 
-      Vect<T> new_bdry(bdry.size() + 2 * (((xmin_reg + 1) % 2 && (xmax_reg + 1) % 2) - (xmax_reg - xmin_reg)/2));
+      Vect<T> new_bdry(bdry.size()
+		       + 2 * (((xmin_reg + 1) % 2 && (xmax_reg + 1) % 2) - (xmax_reg - xmin_reg)/2));
 
       int ishift = 0;
       for (int i = 0; i <= xmin_reg; ++i) new_bdry[i] = bdry[i];
@@ -201,18 +202,40 @@ namespace ABACUS {
 
   // ******************************** Recursive integration functions ******************************
 
-  DP Integrate_Riemann (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax, int Npts);
+  DP Integrate_Riemann (DP (*function) (Vect_DP),
-  DP Integrate_Riemann_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
+			Vect_DP& args, int arg_to_integ,
+			DP xmin, DP xmax,
+			int Npts);
+
+  DP Integrate_Riemann_using_table (DP (*function) (Vect_DP, I_table),
+				    Vect_DP& args, int arg_to_integ,
+				    I_table Itable,
 				    DP xmin, DP xmax, int Npts);
 
-  DP Integrate_rec (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax, DP req_prec, int max_rec_level);
+  DP Integrate_rec (DP (*function) (Vect_DP),
-  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
+		    Vect_DP& args, int arg_to_integ,
-				DP xmin, DP xmax, DP req_prec, int max_rec_level);
+		    DP xmin, DP xmax,
-  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable,
+		    DP req_prec, int max_rec_level);
-				DP xmin, DP xmax, DP req_prec, int max_rec_level, std::ofstream& outfile);
+
-  DP Integrate_rec_using_table_and_file (DP (*function) (Vect_DP, I_table, std::ofstream&), Vect_DP& args,
+  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table),
-					 int arg_to_integ, I_table Itable,
+				Vect_DP& args, int arg_to_integ,
-					 DP xmin, DP xmax, DP req_prec, int max_rec_level, std::ofstream& outfile);
+				I_table Itable,
+				DP xmin, DP xmax,
+				DP req_prec, int max_rec_level);
+
+  DP Integrate_rec_using_table (DP (*function) (Vect_DP, I_table),
+				Vect_DP& args, int arg_to_integ,
+				I_table Itable,
+				DP xmin, DP xmax,
+				DP req_prec, int max_rec_level,
+				std::ofstream& outfile);
+
+  DP Integrate_rec_using_table_and_file (DP (*function) (Vect_DP, I_table, std::ofstream&),
+					 Vect_DP& args, int arg_to_integ,
+					 I_table Itable,
+					 DP xmin, DP xmax,
+					 DP req_prec, int max_rec_level,
+					 std::ofstream& outfile);
 
 
 
@@ -248,34 +271,70 @@ namespace ABACUS {
     DP xmax;
 
   public:
-    Integral_data (DP (*function_ref) (Vect_DP), Vect_DP& args, int arg_to_integ_ref, DP xmin_ref, DP xmax_ref);
+    Integral_data (DP (*function_ref) (Vect_DP),
-    Integral_data (DP (*function_ref) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ_ref,
+		   Vect_DP& args, int arg_to_integ_ref,
-		   I_table Itable, DP xmin_ref, DP xmax_ref);
+		   DP xmin_ref, DP xmax_ref);
-    Integral_data (DP (*function_ref) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ_ref,
+
-		   Integral_table Itable, DP xmin_ref, DP xmax_ref);
+    Integral_data (DP (*function_ref) (Vect_DP, I_table),
+		   Vect_DP& args, int arg_to_integ_ref,
+		   I_table Itable,
+		   DP xmin_ref, DP xmax_ref);
+
+    Integral_data (DP (*function_ref) (Vect_DP, Integral_table),
+		   Vect_DP& args, int arg_to_integ_ref,
+		   Integral_table Itable,
+		   DP xmin_ref, DP xmax_ref);
+
     void Save (std::ofstream& outfile);
-    void Improve_estimate (DP (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, int Npts_max);
+
-    void Improve_estimate (DP (*function) (Vect_DP, I_table), Vect_DP& args, int arg_to_integ, I_table Itable, int Npts_max);
+    void Improve_estimate (DP (*function) (Vect_DP),
-    void Improve_estimate (DP (*function) (Vect_DP, Integral_table), Vect_DP& args, int arg_to_integ,
+			   Vect_DP& args, int arg_to_integ,
-			   Integral_table Itable, int Npts_max);
+			   int Npts_max);
+
+    void Improve_estimate (DP (*function) (Vect_DP, I_table),
+			   Vect_DP& args, int arg_to_integ,
+			   I_table Itable,
+			   int Npts_max);
+
+    void Improve_estimate (DP (*function) (Vect_DP, Integral_table),
+			   Vect_DP& args, int arg_to_integ,
+			   Integral_table Itable,
+			   int Npts_max);
+
     ~Integral_data ();
 
   };
 
-  Integral_result Integrate_optimal (DP (*function) (Vect_DP), Vect_DP& args,
+  Integral_result Integrate_optimal (DP (*function) (Vect_DP),
-				     int arg_to_integ, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
+				     Vect_DP& args, int arg_to_integ,
-  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
+				     DP xmin, DP xmax,
-						 I_table Itable, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
+				     DP req_rel_prec, DP req_abs_prec,
-  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, Integral_table Itable), Vect_DP& args, int arg_to_integ,
+				     int max_nr_pts);
-						 Integral_table Itable, DP xmin, DP xmax, DP req_rel_prec,
-						 DP req_abs_prec, int max_nr_pts);
 
-  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
+  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable),
-						 I_table Itable, DP xmin, DP xmax, DP req_rel_prec,
+						 Vect_DP& args, int arg_to_integ,
-						 DP req_abs_prec, int max_nr_pts, std::ofstream& outfile);
+						 I_table Itable,
+						 DP xmin, DP xmax,
+						 DP req_rel_prec, DP req_abs_prec,
+						 int max_nr_pts);
+
+  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, Integral_table Itable),
+						 Vect_DP& args, int arg_to_integ,
+						 Integral_table Itable,
+						 DP xmin, DP xmax,
+						 DP req_rel_prec, DP req_abs_prec,
+						 int max_nr_pts);
+
+  Integral_result Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable),
+						 Vect_DP& args, int arg_to_integ,
+						 I_table Itable,
+						 DP xmin, DP xmax,
+						 DP req_rel_prec, DP req_abs_prec,
+						 int max_nr_pts,
+						 std::ofstream& outfile);
 
 
-  // ******************************** Recursive version:  optimal, complex implementation ******************************
+  // ********************** Recursive version:  optimal, complex implementation ********************
 
   // NB:  function returns complex values but takes real arguments
 
@@ -307,23 +366,25 @@ namespace ABACUS {
     DP xmax;
 
   public:
-    Integral_data_CX (std::complex<DP> (*function_ref) (Vect_DP), Vect_DP& args, int arg_to_integ_ref, DP xmin_ref, DP xmax_ref);
+    Integral_data_CX (std::complex<DP> (*function_ref) (Vect_DP),
+		      Vect_DP& args, int arg_to_integ_ref,
+		      DP xmin_ref, DP xmax_ref);
+
     void Save (std::ofstream& outfile);
-    void Improve_estimate (std::complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, int Npts_max);
+
+    void Improve_estimate (std::complex<DP> (*function) (Vect_DP),
+			   Vect_DP& args, int arg_to_integ,
+			   int Npts_max);
+
     ~Integral_data_CX ();
 
   };
 
-  Integral_result_CX Integrate_optimal (std::complex<DP> (*function) (Vect_DP), Vect_DP& args, int arg_to_integ, DP xmin, DP xmax,
+  Integral_result_CX Integrate_optimal (std::complex<DP> (*function) (Vect_DP),
-					DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
+					Vect_DP& args, int arg_to_integ,
-  //Integral_result_CX Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
+					DP xmin, DP xmax,
-  //					 I_table Itable, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts);
+					DP req_rel_prec, DP req_abs_prec,
-  //Integral_result_CX Integrate_optimal_using_table (DP (*function) (Vect_DP, I_table Itable), Vect_DP& args, int arg_to_integ,
+					int max_nr_pts);
-  //					 I_table Itable, DP xmin, DP xmax, DP req_rel_prec, DP req_abs_prec, int max_nr_pts, std::ofstream& outfile);
-
-
-
-
 
 } // namespace ABACUS
 

include/ABACUS_Matrix.h

@@ -8,7 +8,7 @@ Copyright (c) J.-S. Caux.
 
 File:  ABACUS_Matrix.h
 
-Purpose:  Declares square matrix class.
+Purpose:  Declares square and rectangular matrix classes.
 
 ***********************************************************/
 

include/ABACUS_Scan.h

@@ -101,35 +101,6 @@ namespace ABACUS {
   Vect<std::string> Descendent_States_with_iK_Stepped_Down_rightIx2only
     (std::string ScanIx2_label, const Heis_Bethe_State& OriginState, bool disperse_only_current_exc, bool preserve_nexc);
 
-  // Functions in src/SCAN/General_Scan.cc:
-  void Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
-		     int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
-  void Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
-		     int Max_Secs, DP target_sumrule, bool refine);
-  void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, std::string defaultScanStatename, int iKmin, int iKmax,
-		     int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
-  void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, std::string defaultname, int iKmin, int iKmax,
-		     int Max_Secs, DP target_sumrule, bool refine);
-  void Scan_LiebLin_Geometric_Quench (DP c_int, DP L_1, int type_id_1, long long int id_1, DP L_2, int N,
-				      int iK_UL, int Max_Secs, DP target_sumrule, bool refine);
-
-  void Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
-  void Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine);
-  void Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, std::string defaultScanStatename, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
-  void Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, std::string defaultScanStatename, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
-  void Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, std::string defaultScanStatename, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
-
-  void Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
-		   int Max_Secs, DP target_sumrule, bool refine, int rank, int nr_processors);
-  void Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax, int Max_Secs, bool refine);
-  void Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKneeded, int Max_Secs, bool refine);
-  void Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int Max_Secs, bool refine);
-
   // Functions to prepare and wrapup parallel scans:
   void Prepare_Parallel_Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
 				      std::string defaultname, int paralevel, Vect<int> rank_lower_paralevels,
@@ -228,10 +199,50 @@ namespace ABACUS {
 
   std::ostream& operator<< (std::ostream& s, const Scan_Info& info);
 
+  // Functions in src/SCAN/General_Scan.cc:
   template<class Tstate>
     Scan_Info General_Scan (char whichDSF, int iKmin, int iKmax, int iKmod, DP kBT, Tstate& AveragingState, Tstate& SeedScanState,
 			    std::string defaultScanStatename, int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors);
 
+  Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
+			  int Max_Secs, DP target_sumrule, bool refine,
+			  int paralevel, Vect<int> rank, Vect<int> nr_processors);
+  Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
+			  int Max_Secs, DP target_sumrule, bool refine);
+  Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, std::string defaultScanStatename,
+			  int iKmin, int iKmax,
+			  int Max_Secs, DP target_sumrule, bool refine,
+			  int paralevel, Vect<int> rank, Vect<int> nr_processors);
+  Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, std::string defaultname,
+			  int iKmin, int iKmax,
+			  int Max_Secs, DP target_sumrule, bool refine);
+  Scan_Info Scan_LiebLin_Geometric_Quench (DP c_int, DP L_1, int type_id_1, long long int id_1, DP L_2, int N,
+					   int iK_UL, int Max_Secs, DP target_sumrule, bool refine);
+
+  Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
+  Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine);
+  Scan_Info Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, std::string defaultScanStatename,
+		       int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
+  Scan_Info Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, std::string defaultScanStatename,
+		       int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
+  Scan_Info Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, std::string defaultScanStatename,
+		       int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors);
+
+  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
+			int Max_Secs, DP target_sumrule, bool refine, int rank, int nr_processors);
+  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax, int Max_Secs, bool refine);
+  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int iKneeded, int Max_Secs, bool refine);
+  Scan_Info Scan_ODSLF (char whichDSF, DP Delta, int N, int M, int Max_Secs, bool refine);
+
 
   //****************************************************************************
 

include/ABACUS_Spec_Fns.h

@@ -26,7 +26,8 @@ namespace ABACUS {
     // Refer to GR[6] 8.23
 
     if (x <= 0.0) {
-      std::cout << "Cosine_Integral called with real argument " << x << " <= 0, which is ill-defined because of the branch cut." << std::endl;
+      std::cout << "Cosine_Integral called with real argument "
+		<< x << " <= 0, which is ill-defined because of the branch cut." << std::endl;
       ABACUSerror("");
     }
 
@@ -59,7 +60,8 @@ namespace ABACUS {
 
     else { // Use high x power series
 
-      // Ci (x) = \frac{\sin x}{x} \sum_{n=0}^\infty (-1)^n (2n)! x^{-2n} - \frac{\cos x}{x} \sum_{n=0}^\infty (-1)^n (2n+1)! x^{-2n-1}
+      // Ci (x) = \frac{\sin x}{x} \sum_{n=0}^\infty (-1)^n (2n)! x^{-2n}
+      //        - \frac{\cos x}{x} \sum_{n=0}^\infty (-1)^n (2n+1)! x^{-2n-1}
 
       int n = 0;
       DP minonetothen = 1.0;
@@ -167,7 +169,8 @@ namespace ABACUS {
     int max_nr_pts = 10000;
     Integral_result integ_ln_Gamma = Integrate_optimal (ln_Gamma_for_Barnes_G_RE, args, 0, 0.0, z - 1.0, req_rel_prec, req_abs_prec, max_nr_pts);
 
-    return(0.5 * (z - 1.0) * (2.0 - z + logtwoPI) + (z - 1.0) * real(ln_Gamma(std::complex<double>(z - 1.0))) - integ_ln_Gamma.integ_est);
+    return(0.5 * (z - 1.0) * (2.0 - z + logtwoPI)
+	   + (z - 1.0) * real(ln_Gamma(std::complex<double>(z - 1.0))) - integ_ln_Gamma.integ_est);
   }
 
 } // namespace ABACUS

include/ABACUS_Utils.h

@@ -39,7 +39,9 @@ namespace ABACUS {
 
   // Inexplicably missing string functions in standard library:
 
-  inline std::string replace(const std::string& str, const std::string& from, const std::string& to) {
+  inline std::string replace(const std::string& str,
+			     const std::string& from,
+			     const std::string& to) {
     std::string repl = str;
     size_t start_pos = repl.find(from);
     if(start_pos < std::string::npos)
@@ -47,7 +49,9 @@ namespace ABACUS {
     return repl;
   }
 
-  inline std::string replace_all(const std::string& str, const std::string& from, const std::string& to) {
+  inline std::string replace_all(const std::string& str,
+				 const std::string& from,
+				 const std::string& to) {
     std::string repl = str;
     if(from.empty())
       return repl;
@@ -125,7 +129,8 @@ namespace ABACUS {
     int ans = 0;
 
     if (N < 0) {
-      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
+      std::cerr << "Error:  factorial of negative number.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if ( N == 1  ||  N == 0) ans = 1;
@@ -139,7 +144,8 @@ namespace ABACUS {
     DP ans = 0.0;
 
     if (N < 0) {
-      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
+      std::cerr << "Error:  factorial of negative number.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if ( N == 1  ||  N == 0) ans = 0.0;
@@ -153,7 +159,8 @@ namespace ABACUS {
     long long int ans = 0;
 
     if (N < 0) {
-      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
+      std::cerr << "Error:  factorial of negative number.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if ( N == 1  ||  N == 0) ans = 1;
@@ -167,7 +174,8 @@ namespace ABACUS {
     unsigned long long int ans = 0;
 
     if (N < 0) {
-      std::cerr << "Error:  factorial of negative number.  Exited." << std::endl;
+      std::cerr << "Error:  factorial of negative number.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if ( N == 1  ||  N == 0) ans = 1;
@@ -182,7 +190,8 @@ namespace ABACUS {
 
     int ans = 0;
     if (N1 < N2) {
-      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
+      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if (N1 == N2) ans = 1;
@@ -203,7 +212,8 @@ namespace ABACUS {
 
     DP ans = 0.0;
     if (N1 < N2) {
-      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
+      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if (N1 == N2) ans = 0.0;
@@ -219,13 +229,14 @@ namespace ABACUS {
 
     long long int ans = 0;
     if (N1 < N2) {
-      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
+      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if (N1 == N2) ans = 1;
     else if (N1 < 12) ans = fact_lli(N1)/(fact_lli(N2) * fact_lli(N1 - N2));
     else {
-      // Make sure that N2 is less than or equal to N1/2;  if not, just switch...
+      // Make sure that N2 is less than or equal to N1/2;  if not, just switch
       int N2_min = min(N2, N1 - N2);
 
       ans = 1;
@@ -244,13 +255,14 @@ namespace ABACUS {
 
     unsigned long long int ans = 0;
     if (N1 < N2) {
-      std::cout << "Error:  N1 smaller than N2 in choose.  Exited." << std::endl;
+      std::cout << "Error:  N1 smaller than N2 in choose.  Exited."
+		<< std::endl;
       exit(1);
     }
     else if (N1 == N2) ans = 1;
     else if (N1 < 12) ans = fact_ulli(N1)/(fact_ulli(N2) * fact_ulli(N1 - N2));
     else {
-      // Make sure that N2 is less than or equal to N1/2;  if not, just switch...
+      // Make sure that N2 is less than or equal to N1/2;  if not, just switch
       int N2_min = min(N2, N1 - N2);
 
       ans = 1;
@@ -302,6 +314,7 @@ namespace ABACUS {
     return(-0.5 * II * log((1.0 + II* x)/(1.0 - II* x)));
   }
 
+
   /****************  Gamma function *******************/
 
   inline std::complex<double> ln_Gamma (std::complex<double> z)
@@ -325,7 +338,8 @@ namespace ABACUS {
 	+ 0.5384136432509564062961e-7 / (z + 8.0)
 	- 0.4023533141268236372067e-8 / (z + 9.0);
 
-      return(0.5 * logtwoPI + (z - 0.5) * log(z + 8.5) - (z + 8.5) + log(series));
+      return(0.5 * logtwoPI + (z - 0.5) * log(z + 8.5)
+	     - (z + 8.5) + log(series));
     }
 
     return(log(0.0));  // never called
@@ -359,7 +373,8 @@ namespace ABACUS {
       for (int i = 1; i < g+2; ++i)
 	series += p[i]/(z_min_1 + std::complex<double>(i));
 
-      return(0.5 * logtwoPI + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
+      return(0.5 * logtwoPI
+	     + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
 	     - (z_min_1 + std::complex<double>(g) + 0.5) + log(series));
     }
 
@@ -376,16 +391,25 @@ namespace ABACUS {
 
       int g = 7;
 
-      double p[9] = { 0.99999999999980993, 676.5203681218851, -1259.1392167224028,
+      double p[9] = {
-		      771.32342877765313, -176.61502916214059, 12.507343278686905,
+	0.99999999999980993,
-		      -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
+	676.5203681218851,
+	-1259.1392167224028,
+	771.32342877765313,
+	-176.61502916214059,
+	12.507343278686905,
+	-0.13857109526572012,
+	9.9843695780195716e-6,
+	1.5056327351493116e-7
+      };
 
       std::complex<double> z_min_1 = z - 1.0;
       std::complex<double> series = p[0];
       for (int i = 1; i < g+2; ++i)
 	series += p[i]/(z_min_1 + std::complex<double>(i));
 
-      return(0.5 * logtwoPI + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
+      return(0.5 * logtwoPI
+	     + (z_min_1 + 0.5) * log(z_min_1 + std::complex<double>(g) + 0.5)
 	     - (z_min_1 + std::complex<double>(g) + 0.5) + log(series));
     }
 
@@ -396,7 +420,8 @@ namespace ABACUS {
 
   inline long long int Partition_Function (int n)
   {
-    // Returns the value of the partition function p(n), giving the number of partitions of n into integers.
+    // Returns the value of the partition function p(n),
+    // giving the number of partitions of n into integers.
 
     if (n < 0) ABACUSerror("Calling Partition_Function for n < 0.");
     else if (n == 0 || n == 1) return(1LL);

include/ABACUS_Young.h

@@ -45,7 +45,8 @@ namespace ABACUS {
   public:
     Young_Tableau ();    // empty constructor, does nothing
     Young_Tableau (int Nr, int Nc);  // constructs empty tableau
-    Young_Tableau (int Nr, int Nc, long long int idnr);  // constructs the tableau corresponding to identification number idnr
+    // constructs the tableau with identification number idnr:
+    Young_Tableau (int Nr, int Nc, long long int idnr);
     Young_Tableau (const Young_Tableau& RefTableau);  // copy constructor
     Young_Tableau (int Nr, int Nc, const Young_Tableau& RefTableau);
     Young_Tableau& operator= (const Young_Tableau& RefTableau);  // assignment
@@ -55,11 +56,13 @@ namespace ABACUS {
     Young_Tableau& Compute_Map (long long int idnr_to_reach);  // fills the map vector
     Young_Tableau& Distribute_boxes (int nboxes_to_dist, int level);
     Young_Tableau& Set_to_id (long long int idnr); // sets the tableau to the one corresponding to idnr
-    Young_Tableau& Set_to_id (long long int idnr, int option); // sets the tableau to the one corresponding to idnr according to rule option
+    // sets the tableau to the one corresponding to idnr according to rule option:
+    Young_Tableau& Set_to_id (long long int idnr, int option);
     Young_Tableau& Set_Row_L (Vect<int>& Row_Lengths);  // set row lengths
     Young_Tableau& Set_Col_L_given_Row_L ();   // sets the Col_L array self-consistently
     Young_Tableau& Set_Row_L_given_Col_L ();   // sets the Col_L array self-consistently
-    long long int Compute_Descendent_id (int option, Vect<int>& Desc_Row_L, int Nrows_Desc, int Ncols_Desc,
+    long long int Compute_Descendent_id (int option, Vect<int>& Desc_Row_L,
+					 int Nrows_Desc, int Ncols_Desc,
 					 const Young_Tableau& RefTableau);
     Young_Tableau& Compute_id();       // computes the id number of tableau
     Young_Tableau& Compute_id(int option);       // computes the id number of tableau according to rule option
@@ -69,7 +72,8 @@ namespace ABACUS {
     bool Raise_Row (int i);
     bool Lower_Col (int i);
     bool Raise_Col (int i);
-    bool Raise_Lowest_Nonzero_Row(); // adds a box to the lowest nonzero length Row, recomputes id, returns true if tableau has changed
+    // adds a box to the lowest nonzero length Row, recomputes id, returns true if tableau has changed:
+    bool Raise_Lowest_Nonzero_Row();
     bool Raise_Next_to_Lowest_Nonzero_Row(); // same thing, but for Row under lowest nonzero length one.
     bool Move_Box_from_Col_to_Col (int ifrom, int ito);
 
@@ -83,10 +87,12 @@ namespace ABACUS {
 
   inline int Nr_Nonzero_Rows (const Vect<Young_Tableau>& Tableau_ref)
   {
-    // This function checks the number of rows containing at least one box
+    /*
-    // in the whole vector of Young tableaux given as argument.
+    This function checks the number of rows containing at least one box
-    // The usefulness is to force descent of states in which only a few
+    in the whole vector of Young tableaux given as argument.
-    // excitations have started dispersing.
+    The usefulness is to force descent of states in which only a few
+    excitations have started dispersing.
+    */
 
     int nr_nonzero_rows = 0;
     for (int i = 0; i < Tableau_ref.size(); ++i)

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)."
+    echo "Arguments needed: whichDSF, kBT, target_sumrule, nr_minutes (for each c, k combination)."
-	exit 1
+    exit 1
 fi
 
 if [[ $1 != 'd' && $1 != 'g' && $1 != 'o' ]]; then
-	echo "Only the d, g and o scanning options are implemented."
+    echo "Only the d, g and o scanning options are implemented."
-	exit 1
+    exit 1
 fi
 
 if [[ $2 -lt 0 ]]; then
-	echo "Temperature kBT must be > 0."
+    echo "Temperature kBT must be > 0."
-	exit 1
+    exit 1
 fi
 
 if [[ $3 -lt 0 || $3 -gt 1 ]]; then
-	echo "Requirement: 0 < target_sumrule < 1."
+    echo "Requirement: 0 < target_sumrule < 1."
-	exit 1
+    exit 1
 fi
 
 if [[ $4 -lt 1 ]]; then
-	echo "Please give a positive nr_minutes."
+    echo "Please give a positive nr_minutes."
-	exit 1
+    exit 1
 fi
 
 zmodload zsh/datetime
@@ -42,7 +42,15 @@ kBT=$2
 target_sumrule=$3
 nr_minutes=$4
 
-logfile='run_'$whichDSF'_kBT_'$kBT'_sr_'$target_sumrule'_'$(date '+%Y-%m-%d-%Hh%M')'.log'
+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'_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)
@@ -52,15 +60,14 @@ basedir="$(pwd)"
 
 for c in $clist
 do
-	for nk in {1..16}
+    for nk in {1..16}
-	do
+    do
-		echo 'Starting run for c =' $c', kfact = '$nk | tee -a $logfile
+	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
+	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/prelim'
-		echo $dir
+	mkdir -p $dir
-		mkdir -p $dir
+	cd $dir
-		cd $dir
+	LiebLin_Catalogue_Fixed_c_k_Nscaling $whichDSF $c $nk $kBT $target_sumrule 0 $nr_minutes | tee -a $basedir/$logfile
-		LiebLin_Catalogue_Fixed_c_k_Nscaling $whichDSF $c $nk $kBT $target_sumrule 0 $nr_minutes | tee -a $basedir/$logfile
+	cd $basedir
-		cd $basedir
+	echo '** Successfully completed run for c =' $c', kfact = '$nk'.\n' | tee -a $logfile
-		echo ' Successfully completed run for c =' $c', kfact = '$nk'\n' | tee -a $logfile
+    done
-	done
 done

scripts/build_LiebLin_catalogue_k_fixed_dsfs.sh

@@ -0,0 +1,77 @@
+#! /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
+
+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
+    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
+	dir='Lieb-Liniger/'$correlator'/T_'$kBT'/c_'$c'/k_fixed/k_'${(l:2::0:)nk}'kFo4/sr_'$target_sumrule'/store/data'
+	cd $dir
+	for datadir in *
+	do
+	    cd $datadir
+	    N=${datadir#"N_"}
+	    iK=$(($nk * $N/8))
+	    Smoothen_LiebLin_DSF_Scaled $whichDSF $c $N $N $iK $iK $kBT 0 $ommin $ommax $Nom $width
+	    cd ..
+	done
+	cd $basedir
+	echo ' Successfully computed DSFs for c =' $c', kfact = '$nk'\n'
+    done
+done

scripts/build_LiebLin_catalogue_k_fixed_dsfs_plots.sh

@@ -0,0 +1,79 @@
+#! /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
+
+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
+    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'
+	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
+	do
+	    echo 'Found file ' $file
+	    # invoke the python script
+	    plot_dsf_k_fixed.py $omegafile $file
+	done
+
+	cd $basedir
+	echo '** Successfully plotted DSFs for c =' $c', k = '$nk'kFo4.\n'
+    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[3]
+L = elements[5]
+N = elements[7]
+iK = elements[12]
+width = elements[20].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.cc

@@ -92,17 +92,45 @@ int main(int argc, char* argv[])
       ActualTime = omp_get_wtime();
       Secs_left = int(Max_Secs - (ActualTime - StartTime));
 
-      if (srsat < target_sumrule && Secs_left > Max_Secs/2)
+      Scan_Info resulting_info;
+      if (srsat < target_sumrule && Secs_left > Max_Secs/2) {
 	// Improve the icmin calculation by one chunk:
-	Scan_LiebLin (whichDSF, c_int, L, N, iKmin, iKmax, kBT, Secs_left, target_sumrule, refine);
+	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));
-      cout << "Done with N = " << N << ". Time left = " << Secs_left << " seconds." << endl;
 
-      if (Secs_left < 60) {
+      if (Secs_left < 30) {
-	cout << "Breaking out after N = " << N << " since time left = " << Secs_left << endl;
+	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;
       }
 

src/SCAN/General_Scan.cc

@@ -819,7 +819,11 @@ namespace ABACUS {
 		    << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
 	LOG_outfile << "Resulting info: " << scan_info << endl;
       }
-      LOG_outfile << "Code version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl << endl;
+      time_t current_time = time(nullptr);
+      char timestr[100];
+      strftime(timestr, sizeof(timestr), "%Y-%m-%d %H:%M:%S", gmtime(&current_time));
+      LOG_outfile << "Run completion timestamp: " << timestr << " UTC" << endl;
+      LOG_outfile << "ABACUS version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl << endl;
       LOG_outfile.close();
     }
 
@@ -853,9 +857,9 @@ namespace ABACUS {
 
   // General version for equilibrium correlators at generic (possibly finite) temperature:
 
-  void Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
+  Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
-		     int Max_Secs, DP target_sumrule, bool refine,
+			  int Max_Secs, DP target_sumrule, bool refine,
-		     int paralevel, Vect<int> rank, Vect<int> nr_processors)
+			  int paralevel, Vect<int> rank, Vect<int> nr_processors)
   {
 
     // This function scans the Hilbert space of the LiebLin gas,
@@ -960,30 +964,26 @@ namespace ABACUS {
     spsfile.close();
 
     // Perform the scan:
-    General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT, spstate, SeedScanState, "",
+    return General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT, spstate, SeedScanState, "",
-		  Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
+			 Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-
-    return;
   }
 
-  void Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
+  Scan_Info Scan_LiebLin (char whichDSF, DP c_int, DP L, int N, int iKmin, int iKmax, DP kBT,
-		     int Max_Secs, DP target_sumrule, bool refine)
+			  int Max_Secs, DP target_sumrule, bool refine)
   {
     int paralevel = 0;
     Vect<int> rank(0,1);
     Vect<int> nr_processors(0,1);
 
-    Scan_LiebLin (whichDSF, c_int, L, N, iKmin, iKmax, kBT, Max_Secs, target_sumrule,
+    return Scan_LiebLin (whichDSF, c_int, L, N, iKmin, iKmax, kBT, Max_Secs, target_sumrule,
-		  refine, paralevel, rank, nr_processors);
+			 refine, paralevel, rank, nr_processors);
-
-    return;
   }
 
 
   // Scanning on an excited state defined by a set of Ix2:
-  void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename,
+  Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename,
-		     int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine,
+			  int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine,
-		     int paralevel, Vect<int> rank, Vect<int> nr_processors)
+			  int paralevel, Vect<int> rank, Vect<int> nr_processors)
   {
     // This function is as Scan_LiebLin for generic T defined above, except that the
     // averaging is now done on a state defined by AveragingStateIx2
@@ -1024,30 +1024,28 @@ namespace ABACUS {
     DP kBT = 0.0;
 
     // Perform the scan:
-    General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT, AveragingState, SeedScanState, defaultScanStatename,
+    return General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT,
-		  Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
+			 AveragingState, SeedScanState, defaultScanStatename,
-
+			 Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-    return;
-
   }
 
   // Simplified function call of the above:
-  void Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename,
+  Scan_Info Scan_LiebLin (char whichDSF, LiebLin_Bethe_State AveragingState, string defaultScanStatename,
-		     int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine)
+			  int iKmin, int iKmax, int Max_Secs, DP target_sumrule, bool refine)
   {
     int paralevel = 0;
     Vect<int> rank(0,1);
     Vect<int> nr_processors(0,1);
 
-    Scan_LiebLin (whichDSF, AveragingState, defaultScanStatename, iKmin, iKmax, Max_Secs,
+    return Scan_LiebLin (whichDSF, AveragingState, defaultScanStatename, iKmin, iKmax, Max_Secs,
-		  target_sumrule, refine, paralevel, rank, nr_processors);
+			 target_sumrule, refine, paralevel, rank, nr_processors);
-
-    return;
   }
 
   // Scanning on a previously-defined AveragingState
-  void Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax,
+  Scan_Info Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, string defaultScanStatename,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors)
+		       int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors)
   {
     // General state scanning for Heisenberg chains
 
@@ -1070,8 +1068,10 @@ namespace ABACUS {
   }
 
   // Scanning on a previously-defined AveragingState
-  void Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax,
+  Scan_Info Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, string defaultScanStatename,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors)
+		       int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors)
   {
     // General state scanning for Heisenberg chains
 
@@ -1088,14 +1088,16 @@ namespace ABACUS {
     else ABACUSerror("Unknown whichDSF in Scan_Heis.");
 
     // Now the scan itself
-    General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
+    return General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
-		  defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
+			 defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
 
   }
 
   // Scanning on a previously-defined AveragingState
-  void Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, string defaultScanStatename, int iKmin, int iKmax,
+  Scan_Info Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, string defaultScanStatename,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors)
+		       int iKmin, int iKmax,
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors)
   {
     // General state scanning for Heisenberg chains
 
@@ -1112,14 +1114,15 @@ namespace ABACUS {
     else ABACUSerror("Unknown whichDSF in Scan_Heis.");
 
     // Now the scan itself
-    General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
+    return General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
-		  defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
+			 defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
 
   }
 
 
-  void Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
+  Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine, int paralevel, Vect<int> rank, Vect<int> nr_processors)
+		       int Max_Secs, DP target_sumrule, bool refine,
+		       int paralevel, Vect<int> rank, Vect<int> nr_processors)
   {
     // This function scans the Hilbert space of the Heisenberg spin-1/2 chain
     // for the function identified by whichDSF.
@@ -1190,25 +1193,24 @@ namespace ABACUS {
       else ABACUSerror("Unknown whichDSF in Scan_Heis.");
 
       // Now the scan itself
-      General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
+      return General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
-		    Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
+			   Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
     }
 
     else ABACUSerror("Delta out of range in Heis_Structure_Factor");
 
-    return;
+    return Scan_Info();
   }
 
-  void Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
+  Scan_Info Scan_Heis (char whichDSF, DP Delta, int N, int M, int iKmin, int iKmax,
-		  int Max_Secs, DP target_sumrule, bool refine)
+		       int Max_Secs, DP target_sumrule, bool refine)
   {
     int paralevel = 0;
     Vect<int> rank(0,1);
     Vect<int> nr_processors(0,1);
 
-    Scan_Heis (whichDSF, Delta, N, M, iKmin, iKmax, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
+    return Scan_Heis (whichDSF, Delta, N, M, iKmin, iKmax,
-
+		      Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
-    return;
   }
 
 

src/SCAN/General_Scan_Parallel.cc

@@ -408,7 +408,7 @@ namespace ABACUS {
     LOG_outfile << "Refining in parallel mode using " << nr_processors_at_newlevel << " processors."
 		<< endl << "Refining info: " << scan_info_refinement
 		<< endl << "Resulting info:  " << scan_info << endl;
-    LOG_outfile << "Code version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl;
+    LOG_outfile << "ABACUS version " << ABACUS_VERSION << ", copyright J.-S. Caux." << endl;
     LOG_outfile.close();
 
     return;

src/UTILS/Data_File_Name.cc

@@ -33,9 +33,9 @@ namespace ABACUS {
   {
     name << "LiebLin_";
     if (whichDSF == 'Z') name << "Z";
-    else if (whichDSF == 'd') name << "Rho_Rho";
+    else if (whichDSF == 'd') name << "rho-rho";
-    else if (whichDSF == 'g') name << "Psi_Psidag";
+    else if (whichDSF == 'g') name << "psi-psidag";
-    else if (whichDSF == 'o') name << "Psidag_Psi";
+    else if (whichDSF == 'o') name << "psidag-psi";
     else if (whichDSF == 'q') name << "GeomQuench";
     else if (whichDSF == '1') name << "Type_I_Exp_Data";
     else if (whichDSF == 'B') name << "BECg2";
@@ -58,9 +58,9 @@ namespace ABACUS {
   {
     name << "LiebLin_";
     if (whichDSF == 'Z') name << "Z";
-    else if (whichDSF == 'd') name << "Rho_Rho";
+    else if (whichDSF == 'd') name << "rho-rho";
-    else if (whichDSF == 'g') name << "Psi_Psidag";
+    else if (whichDSF == 'g') name << "psi-psidag";
-    else if (whichDSF == 'o') name << "Psidag_Psi";
+    else if (whichDSF == 'o') name << "psidag-psi";
     else if (whichDSF == 'q') name << "GeomQuench";
     else if (whichDSF == '1') name << "Type_I_Exp_Data";
     else if (whichDSF == 'B') name << "BECg2";