ABACUS/src/SCAN/General_Scan.cc

/**********************************************************

This software is part of J.-S. Caux's ABACUS library.

Copyright (c) J.-S. Caux.

-----------------------------------------------------------

File:  src/SCAN/General_Scan.cc

Purpose:  universal implementation of state scanning:
            functions to descend down hierarchy of intermediate states.

NOTE:  since templated functions have to be in the same file,
       we put all scanning functions here. The externally-accessible
       functions are defined at the end of this file.

***********************************************************/

#include <omp.h>
#include "ABACUS.h"


using namespace std;
using namespace ABACUS;


string LABEL_TO_CHECK = "bla";
//string LABEL_TO_CHECK = "6_0_";
//string LABEL_TO_CHECK = "6_2_22y32";


namespace ABACUS {

  // Types of descendents:
  // 14 == iK stepped up, leading exc one step further (can lead to ph recombination)
  // 13 == iK up, next exc (nr of ph preserved, not taking possible ph recombination into account)
  // 12 == iK up, next ext (nr ph increased, not taking possible ph recombination into account)
  // 11 == iK down, leading exc one step further (can lead to ph recombination)
  // 10 == iK down, next exc (nr of ph preserved, not taking possible ph recombination into account)
  //  9 == iK down, next exc (nr ph increased, not taking possible ph recombination into account)
  //  8 == iK preserved, 14 and 11 (up to 2 ph recombinations)
  //  7 == iK preserved, 14 and 10
  //  6 == iK preserved, 14 and 9
  //  5 == iK preserved, 13 and 11
  //  4 == iK preserved, 13 and 10
  //  3 == iK preserved, 13 and 9
  //  2 == iK preserved, 12 and 11
  //  1 == iK preserved, 12 and 10
  //  0 == iK preserved, 12 and 9

  // For scanning over symmetric states, the interpretation is slightly different.
  // Types 14, 13 and 12 step iK up using the Ix2 on the right only,
  //  and mirrors the change on the left Ix2.
  // Types 11, 10 and 9 step iK down using the Ix2 on the right only,
  //  and mirrors the change on the left Ix2.
  // There is then no need of scanning over types 0 - 8.
  // By convention, types 9, 10 and 11 can call types 9 - 14; types 12-14 can only call types 12-14.

  bool Expect_ph_Recombination_iK_Up (string ScanIx2_label, const Vect<Vect<int> >& OriginIx2,
				      const Vect<Vect<int> >& BaseScanIx2)
  {
    // This function returns true if descending further can lead to a particle-hole recombination.
    // The criteria which are used are:
    // - the active excitation has moved at least one step (so it has already created its p-h pair)
    // - there exists an OriginIx2 between the active Ix2 and the next Ix2
    // (to right or left depending on type of descendent)

    Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);

    // Determine the level and index of the bottom-most left-most right-moving quantum number sits:
    int exclevel = -1;
    int excindex = 0;
    bool excfound = false;
    do {
      exclevel++;
      if (exclevel == ScanIx2.size()) { // there is no right-moving quantum number in ScanIx2
	break;
      }
      for (int alpha = 0; alpha < ScanIx2[exclevel].size(); ++alpha)
	if (ScanIx2[exclevel][alpha] > BaseScanIx2[exclevel][alpha]) {
	  excindex = alpha;
	  excfound = true;
	  break;
	}
    } while (!excfound);
    // If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;

    if (excfound && !BaseScanIx2[exclevel].includes(ScanIx2[exclevel][excindex])) {
      // there exists an already dispersing excitation which isn't in Origin
      // Is there a possible recombination?
      if (excindex < ScanIx2[exclevel].size() - 1) {
	// a particle to the right of excitation has already move right, so there is a hole
	// check that there exists an occupied Ix2 in Origin sitting between the excitation
	// and the next Ix2 to its right in ScanIx2
	for (int alpha = BaseScanIx2[exclevel].size() - 1; alpha >= 0; --alpha)
	  if (BaseScanIx2[exclevel][alpha] > ScanIx2[exclevel][excindex]
	      && BaseScanIx2[exclevel][alpha] < ScanIx2[exclevel][excindex + 1]) {
	    return(true);
	  }
      }
    } // if (excfound)

    return(false);
  }

  // Specialization for Lieb-Liniger:
  bool Expect_ph_Recombination_iK_Up (string ScanIx2_label, const LiebLin_Bethe_State& OriginState)
  {
    Vect<Vect<int> > OriginIx2here(1);
    OriginIx2here[0] = OriginState.Ix2;
    Vect<Vect<int> > BaseScanIx2here(1);
    BaseScanIx2here[0] = OriginState.Ix2;
    return(Expect_ph_Recombination_iK_Up (ScanIx2_label, OriginIx2here, BaseScanIx2here));
  }

  // Specialization for Heis
  bool Expect_ph_Recombination_iK_Up (string ScanIx2_label, const Heis_Bethe_State& OriginState)
  {
    return(Expect_ph_Recombination_iK_Up (ScanIx2_label, OriginState.Ix2, OriginState.Ix2));
  }



  bool Expect_ph_Recombination_iK_Down (string ScanIx2_label, const Vect<Vect<int> >& OriginIx2,
					const Vect<Vect<int> >& BaseScanIx2)
  {
    // This function returns true if descending further can lead to a particle-hole recombination.
    // The criteria which are used are:
    // - the active excitation has moved at least one step (so it has already created its p-h pair)
    // - there exists an OriginIx2 between the active Ix2 and the next Ix2
    //   (to right or left depending on type of descendent)

    Vect<Vect<int> > ScanIx2 = Return_Ix2_from_Label (ScanIx2_label, OriginIx2);

    // Determine the level and index of the bottom-most right-most left-moving quantum number sits:
    int exclevel = -1;
    int excindex = 0;
    bool excfound = false;

    do {
      exclevel++;
      if (exclevel == ScanIx2.size()) { // there isn't a single left-moving quantum number in ScanIx2
	break;
      }
      for (int alpha = ScanIx2[exclevel].size() - 1; alpha >= 0; --alpha) {
	if (ScanIx2[exclevel][alpha] < BaseScanIx2[exclevel][alpha]) {
	  excindex = alpha;
	  excfound = true;
	  break;
	}
      }
    } while (!excfound);
    // If we haven't found an excitation, then exclevel == ScanIx2.size() and excindex = 0;
    if (!excfound) excindex = ScanIx2[exclevel].size() - 1;

    if (excfound && !BaseScanIx2[exclevel].includes(ScanIx2[exclevel][excindex])) {
      // there exists an already dispersing excitation which isn't in Origin
      // Is there a possible recombination?
      if (excindex > 0) {
	// a particle to the left of excitation has already moved left, so there is a hole
	// check that there exists an occupied Ix2 in Origin sitting between the excitation
	// and the next Ix2 to its left in ScanIx2
	for (int alpha = 0; alpha < BaseScanIx2[exclevel].size(); ++alpha)
	  if (BaseScanIx2[exclevel][alpha] > ScanIx2[exclevel][excindex - 1]
	      && BaseScanIx2[exclevel][alpha] < ScanIx2[exclevel][excindex]) {
	    return(true);
	  }
      }
    } // if (excfound)

    return(false);
  }

  // Specialization for Lieb-Liniger:
  bool Expect_ph_Recombination_iK_Down (string ScanIx2_label, const LiebLin_Bethe_State& OriginState)
  {
    Vect<Vect<int> > OriginIx2here(1);
    OriginIx2here[0] = OriginState.Ix2;
    Vect<Vect<int> > BaseScanIx2here(1);
    BaseScanIx2here[0] = OriginState.Ix2;
    return(Expect_ph_Recombination_iK_Down (ScanIx2_label, OriginIx2here, BaseScanIx2here));
  }

  // Specialization for Heis
  bool Expect_ph_Recombination_iK_Down (string ScanIx2_label, const Heis_Bethe_State& OriginState)
  {
    return(Expect_ph_Recombination_iK_Down (ScanIx2_label, OriginState.Ix2, OriginState.Ix2));
  }




  template<class Tstate>
  Scan_Info General_Scan (char whichDSF, int iKmin, int iKmax, int iKmod, DP kBT,
			  Tstate& AveragingState, Tstate& SeedScanState, string defaultScanStatename,
			  int Max_Secs, DP target_sumrule, bool refine,
			  int paralevel, Vect<int> rank, Vect<int> nr_processors)
  {
    // Performs the scan over excited states, writing data to file.

    // AveragingState is the state on which the correlations are calculated.
    // SeedScanState is the originator of all scan states.
    // This distinction is kept to allow for quenches and finite temperatures.

    // This function is also called by the parallel implementation of ABACUS.
    // In this case, file names carry a rank and nr_processors suffix.
    // In fact, the parallelization can be done in incremental levels.
    // If paralevel == 0, the run is serial.
    // If paralevel == n, the run is parallelized in a tree with n levels of branching.
    // A paralevel == 1 branching's files have a suffix of the form "_3_8", meaning that this
    //  is the rank 3 out of 8 processors.
    // A paralevel == 2 branching's files have a suffix of the form "_3_8_2_8", meaning that this
    //  is the rank 2 out of 8 subscan of the _3_8 scan.

    bool in_parallel = (paralevel > 0);
    if (in_parallel && (rank.size() != paralevel || nr_processors.size() != paralevel)) {
      cout << "paralevel = " << paralevel << "\trank.size() = " << rank.size()
	   << "\tnr_processors.size() = " << nr_processors.size() << endl;
      cout << "rank = " << rank << endl;
      cout << "nr_processors = " << nr_processors << endl;
      ABACUSerror("Inconsistent paralevel, rank or nr_processors in General_Scan.");
    }

    if (in_parallel && !refine) ABACUSerror("Must refine when using parallel ABACUS");

    // expected cost on data_value of adding a particle-hole excitation.
    DP ph_cost = Particle_Hole_Excitation_Cost (whichDSF, AveragingState);

    int Max_Secs_used = int(0.9 * Max_Secs); // we don't start any new ithread loop beyond this point
    int Max_Secs_alert = int(0.95 * Max_Secs); // we break any ongoing ithread loop beyond this point

    /* TO_STR
    stringstream filenameprefix;
    Data_File_Name (filenameprefix, whichDSF, iKmin, iKmax, kBT,
		    AveragingState, SeedScanState, defaultScanStatename);

    if (in_parallel)
      for (int r = 0; r < paralevel; ++r)
	filenameprefix << "_" << rank[r] << "_" << nr_processors[r];

    string prefix = filenameprefix.str();
    */
    string prefix = Data_Filename(whichDSF, iKmin, iKmax, kBT,
				  AveragingState, SeedScanState, defaultScanStatename);
    if (in_parallel)
      for (int r = 0; r < paralevel; ++r)
	prefix += "_" + to_string(rank[r]) + "_" + to_string(nr_processors[r]);

    /* TO_STR
    stringstream filenameprefix_prevparalevel;
    // without the rank and nr_processors of the highest paralevel

    Data_File_Name (filenameprefix_prevparalevel, whichDSF, iKmin, iKmax, kBT,
		    AveragingState, SeedScanState, defaultScanStatename);
    if (in_parallel) for (int r = 0; r < paralevel - 1; ++r)
		       filenameprefix << "_" << rank[r] << "_" << nr_processors[r];

    string prefix_prevparalevel = filenameprefix_prevparalevel.str();
    */
    // without the rank and nr_processors of the highest paralevel
    string prefix_prevparalevel = Data_Filename(whichDSF, iKmin, iKmax, kBT,
						AveragingState, SeedScanState, defaultScanStatename);
    if (in_parallel)
      for (int r = 0; r < paralevel - 1; ++r)
	prefix_prevparalevel += "_" + to_string(rank[r]) + "_" + to_string(nr_processors[r]);

    /* TO_STR
    stringstream RAW_stringstream;    string RAW_string;
    stringstream INADM_stringstream;    string INADM_string;
    stringstream CONV0_stringstream;    string CONV0_string;
    stringstream STAT_stringstream;    string STAT_string;
    stringstream LOG_stringstream;    string LOG_string;
    stringstream THR_stringstream;    string THR_string;
    stringstream THRDIR_stringstream;    string THRDIR_string;
    stringstream SRC_stringstream;    string SRC_string;
    stringstream SUM_stringstream;    string SUM_string;

    RAW_stringstream << prefix << ".raw";
    INADM_stringstream << prefix << ".inadm";
    CONV0_stringstream << prefix << ".conv0";
    STAT_stringstream << prefix << ".stat";
    LOG_stringstream << prefix << ".log";
    THR_stringstream << prefix << ".thr";
    THRDIR_stringstream << prefix << "_thrdir";
    SRC_stringstream << prefix << ".src";
    SUM_stringstream << prefix << ".sum";

    RAW_string = RAW_stringstream.str();    const char* RAW_Cstr = RAW_string.c_str();
    INADM_string = INADM_stringstream.str();    const char* INADM_Cstr = INADM_string.c_str();
    CONV0_string = CONV0_stringstream.str();    const char* CONV0_Cstr = CONV0_string.c_str();
    STAT_string = STAT_stringstream.str();    const char* STAT_Cstr = STAT_string.c_str();
    LOG_string = LOG_stringstream.str();    const char* LOG_Cstr = LOG_string.c_str();
    THR_string = THR_stringstream.str();    const char* THR_Cstr = THR_string.c_str();
    SRC_string = SRC_stringstream.str();    const char* SRC_Cstr = SRC_string.c_str();
    SUM_string = SUM_stringstream.str();    const char* SUM_Cstr = SUM_string.c_str();

    THRDIR_string = THRDIR_stringstream.str();
    */
    string RAW_str = prefix + ".raw";
    string INADM_str = prefix + ".inadm";
    string CONV0_str = prefix + ".conv0";
    string STAT_str = prefix + ".stat";
    string LOG_str = prefix + ".log";
    string THR_str = prefix + ".thr";
    string THRDIR_str = prefix + "_thrdir";
    string SRC_str = prefix + ".src";
    string SUM_str = prefix + ".sum";


    fstream RAW_outfile;
    if (!refine || in_parallel) RAW_outfile.open(RAW_str, fstream::out | fstream::trunc);
    else RAW_outfile.open(RAW_str, fstream::out | fstream::app);
    if (RAW_outfile.fail()) {
      cout << RAW_str << endl;
      ABACUSerror("Could not open RAW_outfile... ");
    }
    RAW_outfile.precision(16);

    fstream INADM_outfile;
    if (!refine || in_parallel) INADM_outfile.open(INADM_str, fstream::out | fstream::trunc);
    else INADM_outfile.open(INADM_str, fstream::out | fstream::app);
    if (INADM_outfile.fail()) ABACUSerror("Could not open INADM_outfile... ");
    INADM_outfile.precision(16);

    fstream CONV0_outfile;
    if (!refine || in_parallel) CONV0_outfile.open(CONV0_str, fstream::out | fstream::trunc);
    else CONV0_outfile.open(CONV0_str, fstream::out | fstream::app);
    if (CONV0_outfile.fail()) ABACUSerror("Could not open CONV0_outfile... ");
    CONV0_outfile.precision(16);

    fstream STAT_outfile;
    if (!refine || in_parallel) STAT_outfile.open(STAT_str, fstream::out | fstream::trunc);
    else STAT_outfile.open(STAT_str, fstream::out | fstream::app);
    if (STAT_outfile.fail()) ABACUSerror("Could not open STAT_outfile... ");
    STAT_outfile.precision(8);

    ofstream LOG_outfile;
    if (!in_parallel) {
      if (!refine) LOG_outfile.open(LOG_str, fstream::out | fstream::trunc);
      else LOG_outfile.open(LOG_str, fstream::out | fstream::app);
      if (LOG_outfile.fail()) ABACUSerror("Could not open LOG_outfile... ");
      LOG_outfile.precision(16);
    }
    else { // in_parallel
      LOG_outfile.open(LOG_str, fstream::out | fstream::trunc);
      if (LOG_outfile.fail()) ABACUSerror("Could not open LOG_outfile... ");
      LOG_outfile.precision(16);
    }

    Scan_Info scan_info;

    if (!refine) mkdir(THRDIR_str.c_str(), S_IRWXU | S_IRWXG | S_IRWXO);
    Scan_Thread_Data paused_thread_data (THRDIR_str, refine);

    if (refine) {
      paused_thread_data.Load();
      if (!in_parallel) scan_info.Load(SRC_str);
    }

    Scan_Info scan_info_before = scan_info;  // for LOG file
    Scan_Info scan_info_before_descent;
    Scan_Info scan_info_obtained_in_descent;

    Scan_State_List<Tstate> ScanStateList (whichDSF, SeedScanState);
    ScanStateList.Populate_List(whichDSF, SeedScanState);

    if (refine && !in_parallel) ScanStateList.Load_Info (SUM_str);
    else if (in_parallel && rank.sum() == 0) {}; // do nothing, keep info in the higher .sum file!

    DP Chem_Pot = Chemical_Potential (AveragingState);
    DP sumrule_factor = Sumrule_Factor (whichDSF, AveragingState, Chem_Pot, iKmin, iKmax);


    // Now go for it !

    bool at_least_one_new_flag_raised = false;

    int Ndata_previous_cycle = 0;

    int ninadm = 0; // nr of inadm states for which we save info in .inadm file. Save first 1000.
    int nconv0 = 0; // nr of unconv states for which we save info in .conv0 file. Save first 1000.

    double start_time_omp = omp_get_wtime();
    double current_time_omp = omp_get_wtime();

#pragma omp parallel
    do {

      int omp_thread_nr = omp_get_thread_num();

      if ((paused_thread_data.lowest_il_with_nthreads_neq_0 == paused_thread_data.nlists - 1)
	  && omp_thread_nr > 0) {
	double start_time_wait = omp_get_wtime();
	double stop_time_wait;
	do {
	  for (int i = 0; i < 100000; ++i) { }
	  stop_time_wait = omp_get_wtime();
	} while (stop_time_wait - start_time_wait < 5.0);
      }

      double start_time_cycle_omp = omp_get_wtime();

      at_least_one_new_flag_raised = false;


#pragma omp master
      {
	double start_time_flags = omp_get_wtime();

	// First flag the new base/type 's that we need to include:
	ScanStateList.Raise_Scanning_Flags
	  (exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0));


	// Get these base/type started:
	for (int i = 0; i < ScanStateList.ndef; ++i) {
	  if (ScanStateList.flag_for_scan[i] && ScanStateList.info[i].Ndata == 0
	      && !ScanStateList.scan_attempted[i]) {

	    Scan_Info scan_info_flags;

	    at_least_one_new_flag_raised = true;

	    ScanStateList.scan_attempted[i] = true;

	    Tstate ScanState;

	    ScanState = ScanStateList.State[i];

	    DP data_value = -1.0;

	    bool admissible = ScanState.Check_Admissibility(whichDSF);

	    if (admissible) {

	      ScanState.Compute_All(true);

	      if (ScanState.conv) {

		// Put momentum in fundamental window, if possible:
		int iKexc = ScanState.iK - AveragingState.iK;
		while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
		while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;

		stringstream rawfile_entry;
		data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState, Chem_Pot, rawfile_entry);
		{
#pragma omp critical
		  RAW_outfile << rawfile_entry.str();
		}

		{
#pragma omp critical
		  if (iKexc >= iKmin && iKexc <= iKmax) {
		    scan_info_flags.Ndata++;
		    scan_info_flags.Ndata_conv++;
		    scan_info_flags.sumrule_obtained += data_value*sumrule_factor;
		  }
		}

		// If we force descent:  modify data_value by hand so that descent is forced on next scanning pass
		for (int itype = 0; itype < 15; ++itype) {
		  DP data_value_used = 0.1* exp(-paused_thread_data.logscale * ABACUS::min(0, paused_thread_data.lowest_il_with_nthreads_neq_0));
		  if (Force_Descent(whichDSF, ScanState, AveragingState, itype, iKmod, Chem_Pot))
		    data_value = data_value_used;
		}

		Vect<bool> allowed(false, 15);
		if (whichDSF == 'B') { // symmetric state scanning
		  allowed[9] = true; allowed[10] = true; allowed[11] = true;
		  allowed[12] = true; allowed[13] = true; allowed[14] = true;
		}
		else {
		  allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
		  allowed[1] = allowed[0]; allowed[2] = allowed[0];
		  allowed[3] = allowed[0]; allowed[4] = allowed[0];
		  allowed[5] = allowed[0]; allowed[6] = allowed[0];
		  allowed[7] = allowed[0]; allowed[8] = allowed[0];
		  allowed[9] = (iKexc > iKmin);
		  allowed[10] = allowed[9]; allowed[11] = allowed[9];
		  allowed[12] = (iKexc < iKmax);
		  allowed[13] = allowed[12]; allowed[14] = allowed[12];
		}
		for (int type_required_here = 0; type_required_here < 15; ++type_required_here) {
		  if (!allowed[type_required_here]) continue;
		  // All cases here are such that the ScanState hasn't been descended yet,
		  // so we simply use data_value as expected data value:
		  {
#pragma omp critical
		    paused_thread_data.Include_Thread (abs(data_value), ScanState.label, type_required_here);
		  }
		}
	      } // if (ScanState.conv)

	      else {
		if (nconv0++ < 1000)
		  CONV0_outfile << setw(25) << ScanState.label << setw(25) << ScanState.diffsq
				<< setw(5) << ScanState.Check_Rapidities()
				<< setw(25) << ScanState.String_delta() << endl;
		scan_info_flags.Ndata++;
		scan_info_flags.Ndata_conv0++;
	      }
	    } // if admissible

	    else { // if inadmissible, modify data_value by hand so that descent is forced on next scanning pass

	      if (ninadm++ < 10000000) INADM_outfile << ScanState.label << endl;
	      scan_info_flags.Ndata++;
	      scan_info_flags.Ninadm++;

	      // Put momentum in fundamental window, if possible:
	      int iKexc = ScanState.iK - AveragingState.iK;
	      while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
	      while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;

	      DP data_value = 1.0e-32;
	      for (int itype = 0; itype < 15; ++itype)
		if (Force_Descent(whichDSF, ScanState, AveragingState, itype, iKmod, Chem_Pot))
		  data_value = 0.1* exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0);

	      Vect<bool> allowed(false, 15);
	      if (whichDSF == 'B') {
		// We scan over symmetric states. Only types 14 down to 9 are allowed.
		allowed[9] = true; allowed[10] = true; allowed[11] = true;
		allowed[12] = true; allowed[13] = true; allowed[14] = true;
	      }
	      else {
		allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
		allowed[1] = allowed[0]; allowed[2] = allowed[0];
		allowed[3] = allowed[0]; allowed[4] = allowed[0];
		allowed[5] = allowed[0]; allowed[6] = allowed[0];
		allowed[7] = allowed[0]; allowed[8] = allowed[0];
		allowed[9] = (iKexc > iKmin);
		allowed[10] = allowed[9]; allowed[11] = allowed[9];
		allowed[12] = (iKexc < iKmax);
		allowed[13] = allowed[12]; allowed[14] = allowed[12];
	      }
	      for (int type_required_here = 0; type_required_here < 15; ++type_required_here) {
		if (!allowed[type_required_here]) continue;
		{
#pragma omp critical
		  paused_thread_data.Include_Thread (abs(data_value), ScanState.label, type_required_here);
		}
	      }
	    } // inadmissible

	    scan_info_flags.TT += omp_get_wtime() - start_time_flags;

	    // Put this info into the appropriate ScanStateList.info
	    {
#pragma omp critical
	      ScanStateList.Include_Info(scan_info_flags, ScanStateList.base_label[i]);
	      scan_info += scan_info_flags;
	    }

	  } // if flag_for_scan
	} // for i

      } // #pragma omp master


      // Now we deal with the previously existing paused threads:

      Vect<Scan_Thread> threads_to_do;
      int il_to_do = paused_thread_data.lowest_il_with_nthreads_neq_0; // for resaving threads in case we're out of time
      DP current_threshold = exp(-paused_thread_data.logscale * il_to_do);
      {
#pragma omp critical
	threads_to_do = paused_thread_data.Extract_Next_Scan_Threads();
      }

      int ithread;

      {
	for (ithread = 0; ithread < threads_to_do.size(); ++ithread) {

	  Scan_Info scan_info_this_ithread;
	  double start_time_this_ithread = omp_get_wtime();

	  // If we don't have time anymore, resave the threads instead of computing them:
	  if (start_time_this_ithread - start_time_omp > Max_Secs_alert) {
	    for (int ith = ithread; ith < threads_to_do.size(); ++ith) {
#pragma omp critical
	      paused_thread_data.Include_Thread (il_to_do, threads_to_do[ith].label, threads_to_do[ith].type);
	    }
	    break; // jump out of ithread loop
	  }

	  Tstate ScanState;
	  {
#pragma omp critical
	    ScanState = ScanStateList.Return_State(Extract_Base_Label(threads_to_do[ithread].label));
	  }
	  Tstate BaseScanState;  BaseScanState = ScanState;

	  ScanState.Set_to_Label(threads_to_do[ithread].label, BaseScanState.Ix2);

	  Tstate ScanStateBeingDescended = ScanState;

	  // STARTING Descend_and_Compute block:
	  int type_required = threads_to_do[ithread].type;

	  ScanState.Compute_Momentum();
	  Vect<string> desc_label;

	  bool disperse_only_current_exc_up = false;
	  if (type_required == 14 || type_required == 8 || type_required == 7 || type_required == 6)
	    disperse_only_current_exc_up = true;
	  bool preserve_nexc_up = false;
	  if (type_required == 13 || type_required == 5 || type_required == 4 || type_required == 3)
	    preserve_nexc_up = true;
	  bool disperse_only_current_exc_down = false;
	  if (type_required == 11 || type_required == 8 || type_required == 5 || type_required == 2)
	    disperse_only_current_exc_down = true;
	  bool preserve_nexc_down = false;
	  if (type_required == 10 || type_required == 7 || type_required == 4 || type_required == 1)
	    preserve_nexc_down = true;

	  if (whichDSF == 'B') { // symmetric state scanning
	    if (type_required >= 9 && type_required <= 11)
	      desc_label = Descendent_States_with_iK_Stepped_Down_rightIx2only
		(ScanState.label, BaseScanState, disperse_only_current_exc_down, preserve_nexc_down);
	    else if (type_required >= 12 && type_required <= 14)
	      desc_label = Descendent_States_with_iK_Stepped_Up_rightIx2only
		(ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up);
	  }
	  else {
	    if (type_required >= 0 && type_required <= 8) {
	      desc_label = Descendent_States_with_iK_Preserved
		(ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up,
		 disperse_only_current_exc_down, preserve_nexc_down);
	    }
	    else if (type_required >= 9 && type_required <= 11)
	      desc_label = Descendent_States_with_iK_Stepped_Down
		(ScanState.label, BaseScanState, disperse_only_current_exc_down, preserve_nexc_down);
	    else if (type_required >= 12 && type_required <= 14)
	      desc_label = Descendent_States_with_iK_Stepped_Up
		(ScanState.label, BaseScanState, disperse_only_current_exc_up, preserve_nexc_up);
	  }

	  string label_here = ScanState.label;

	  for (int idesc = 0; idesc < desc_label.size(); ++idesc) {

	    ScanState.Set_to_Label (desc_label[idesc], BaseScanState.Ix2);

	    bool admissible = ScanState.Check_Admissibility(whichDSF);

	    DP data_value = 0.0;

	    ScanState.conv = false;
	    ScanState.Compute_Momentum(); // since momentum is used as forced descent criterion


	    if (admissible) {

	      ScanState.Compute_All (idesc == 0);

	      if (ScanState.conv) {

		// Put momentum in fundamental window, if possible:
		int iKexc = ScanState.iK - AveragingState.iK;
		while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
		while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;

		stringstream rawfile_entry;
		data_value = Compute_Matrix_Element_Contrib (whichDSF, iKmin, iKmax, ScanState, AveragingState,
							     Chem_Pot, rawfile_entry);

		{
#pragma omp critical
		  RAW_outfile << rawfile_entry.str();
		  if (iKexc >= iKmin && iKexc <= iKmax) {
		    scan_info_this_ithread.Ndata++;
		    scan_info_this_ithread.Ndata_conv++;
		    scan_info_this_ithread.sumrule_obtained += data_value*sumrule_factor;
		  }
		}

		// Uncomment line below if .stat file is desired:
		// STAT_outfile << setw(20) << label_here << "\t" << setw(5) << type_required
		// 	     << "\t" << setw(16) << std::scientific
		// 	     << exp(-paused_thread_data.logscale * il_to_do)
		// 	     << "\t" << setw(20) << ScanState.label << "\t" << setw(16) << data_value
		// 	     << "\t" << setw(16) << std::fixed << setprecision(8)
		// 	     << data_value/exp(-paused_thread_data.logscale * il_to_do) << endl;

		// Uncomment below if alerts for unexpectedly high data_value (as compared to threshold) are desired
		// if (fabs(data_value) > 10.0* current_threshold) {
		//   cout << "\nAlert: data_value > 10* threshold, " << data_value << "\t" << current_threshold << endl;
		//   cout << " for state " << ScanState.label << " descendent of type " << type_required
		//        << " of state " << ScanStateBeingDescended.label << endl;
		//   cout << AveragingState.Ix2 << endl;
		//   cout << ScanStateBeingDescended.Ix2 << endl;
		//   cout << ScanState.Ix2 << endl;
		// }

	      } // if (ScanState.conv)
	      else {
		if (nconv0++ < 1000)
		  CONV0_outfile << setw(25) << ScanState.label << setw(25)
				<< ScanState.diffsq << setw(5) << ScanState.Check_Rapidities()
				<< setw(25) << ScanState.String_delta() << endl;
		scan_info_this_ithread.Ndata++;
		scan_info_this_ithread.Ndata_conv0++;
	      }
	    } // if (admissible)

	    else {
	      if (ninadm++ < 1000000) INADM_outfile << ScanState.label << endl;
	      scan_info_this_ithread.Ndata++;
	      scan_info_this_ithread.Ninadm++;
	    }

	    Tstate state_to_descend;  state_to_descend = ScanState; // for checking

	    ScanState.Compute_Momentum();
	    // Put momentum in fundamental window, if possible:
	    int iKexc = ScanState.iK - AveragingState.iK;
	    while (iKexc > iKmax && iKexc - iKmod >= iKmin) iKexc -= iKmod;
	    while (iKexc < iKmin && iKexc + iKmod <= iKmax) iKexc += iKmod;

	    // Momentum-preserving are only descended to momentum-preserving.
	    // Momentum-increasing are only descended to momentum-preserving and momentum-increasing.
	    // Momentum-decreasing are only descended to momentum-preserving and momentum-decreasing.
	    Vect<bool> allowed(false, 15);
	    if (whichDSF == 'B') {
	      // We scan over symmetric states. Only types 14 down to 9 are allowed.
	      if (type_required >= 9 && type_required <= 11) { // iK stepped down on rightIx2; step further up or down
		allowed[9] = true; allowed[10] = true; allowed[11] = true;
		allowed[12] = true; allowed[13] = true; allowed[14] = true;
	      }
	      else if (type_required >= 12 && type_required <= 14) { // iK stepped up on rightIx2; only step further up
		allowed[12] = true; allowed[13] = true; allowed[14] = true;
	      }
	    }
	    else {
	      if (type_required >= 0 && type_required <= 8) { // momentum-preserving
		allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
		allowed[9] = false;
		allowed[12] = false;
	      }
	      if (type_required >= 9 && type_required <= 11) { // momentum-decreasing
		allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
		allowed[9] = (iKexc > iKmin);
		allowed[12] = false;
	      }
	      if (type_required >= 12 && type_required <= 14) { // momentum-increasing
		allowed[0] = (iKexc >= iKmin && iKexc <= iKmax);
		allowed[9] = false;
		allowed[12] = (iKexc < iKmax);
	      }
	      // The others are just copies of the ones above:
	      allowed[1] = allowed[0]; allowed[2] = allowed[0]; allowed[3] = allowed[0]; allowed[4] = allowed[0]; allowed[5] = allowed[0]; allowed[6] = allowed[0]; allowed[7] = allowed[0]; allowed[8] = allowed[0];
	      allowed[10] = allowed[9]; allowed[11] = allowed[9];
	      allowed[13] = allowed[12]; allowed[14] = allowed[12];
	    }


	    for (int type_required_here = 0; type_required_here < 15; ++type_required_here) {

	      if (!allowed[type_required_here]) continue;

	      // Reset ScanState to what it was, if change on first pass
	      if (type_required_here > 0) ScanState = state_to_descend;

	      // We determine if we carry on scanning based on the data_value obtained, or forcing conditions:
	      DP expected_abs_data_value = abs(data_value);

	      //++G_7 logic
	      if ((type_required_here == 14 || type_required_here == 8
		   || type_required_here == 7 || type_required_here == 6)
		  && Expect_ph_Recombination_iK_Up (ScanState.label, BaseScanState))
		expected_abs_data_value /= ph_cost;
	      if (type_required_here == 12 || type_required_here == 2
		  || type_required_here == 1 || type_required_here == 0)
		expected_abs_data_value *= ph_cost;
	      if ((type_required_here == 11 || type_required_here == 8
		   || type_required_here == 5 || type_required_here == 2)
		  && Expect_ph_Recombination_iK_Down (ScanState.label, BaseScanState))
		expected_abs_data_value /= ph_cost;
	      if (type_required_here == 9 || type_required_here == 6
		  || type_required_here == 3 || type_required_here == 0)
		expected_abs_data_value *= ph_cost;

	      {
#pragma omp critical
		paused_thread_data.Include_Thread (expected_abs_data_value, ScanState.label, type_required_here);
	      }
	    } // for type_required_here

	  } // for idesc

	  // FINISHED Descend_and_Compute block


	  scan_info_this_ithread.TT += omp_get_wtime() - start_time_this_ithread;

#pragma omp critical
	  {
	    scan_info += scan_info_this_ithread;
	    ScanStateList.Include_Info(scan_info_this_ithread, Extract_Base_Label(threads_to_do[ithread].label));
	  }
	} // for ithread

      } // omp parallel region


#pragma omp master
      {
	if (!in_parallel)
	  LOG_outfile << "Master cycling. Ndata_conv " << scan_info.Ndata_conv
		      << ". Threshold " << paused_thread_data.lowest_il_with_nthreads_neq_0 << " "
		      << setw(9) << setprecision(3)
		      << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0)
		      << ". " << setw(12) << scan_info.Ndata - Ndata_previous_cycle << " new data. Nr of threads: "
		      << setw(14) << paused_thread_data.nthreads_total.sum()
		      << ".  Saturation:  " << setprecision(12) << scan_info.sumrule_obtained << endl;

	Ndata_previous_cycle = scan_info.Ndata;
      }


      current_time_omp = omp_get_wtime();

    } while (current_time_omp - start_time_omp < Max_Secs_used
	     && scan_info.sumrule_obtained < target_sumrule
	     );
    // This closes the #pragram omp parallel block

    RAW_outfile.close();
    INADM_outfile.close();
    CONV0_outfile.close();
    STAT_outfile.close();

    scan_info.Save(SRC_str);

    Scan_Info scan_info_refine = scan_info;
    scan_info_refine -= scan_info_before;

    if (!in_parallel) {

      if (scan_info.sumrule_obtained >= target_sumrule)
	LOG_outfile << endl << "Achieved sumrule saturation of " << scan_info.sumrule_obtained
		    << "\t(target was " << target_sumrule << ")." << endl << endl;

      if (!refine) {
	LOG_outfile << "Main run info: " << scan_info << endl;
	LOG_outfile << "Latest threshold level " << paused_thread_data.lowest_il_with_nthreads_neq_0
		    << " " << std::scientific << setprecision(3)
		    << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
      }
      else if (refine) {
	LOG_outfile << "Refining info: " << scan_info_refine << endl;
	LOG_outfile << "Latest threshold level " << paused_thread_data.lowest_il_with_nthreads_neq_0
		    << " " << std::scientific << setprecision(3)
		    << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
	LOG_outfile << "Resulting info: " << scan_info << 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();
    }

    else { // in_parallel

      LOG_outfile << "rank " << rank << " out of " << nr_processors << " processors: "
		  << "run info: " << scan_info << endl << "Latest threshold = "
		  << exp(-paused_thread_data.logscale * paused_thread_data.lowest_il_with_nthreads_neq_0) << endl;
    }

    paused_thread_data.Save();

    ScanStateList.Order_in_SRC ();

    ScanStateList.Save_Info (SUM_str);


    // Evaluate f-sumrule:
    if (!in_parallel ) if (whichDSF != 'q')
			 Evaluate_F_Sumrule (prefix_prevparalevel, whichDSF, AveragingState, Chem_Pot, iKmin, iKmax);

    return(scan_info);
  }



  //******************************************************

  // Functions to initiate scans:


  // General version for equilibrium correlators at generic (possibly finite) temperature:

  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)
  {

    // This function scans the Hilbert space of the LiebLin gas,
    // for the function identified by whichDSF.

    // whichDSF == 'Z':  canonical partition function
    // whichDSF == 'd':  density-density correlation function
    // whichDSF == 'g':  Green's function < \Psi \Psi^{\dagger}>
    // whichDSF == 'o':  one-body function < \Psi^{\dagger} \Psi >

    // Delta is the number of sites involved in the smoothing of the entropy
    //int Delta = int(sqrt(N))/2;//6;//N/20;
    //DP epsilon = log(L)/L; // using Gaussian for density in entropy.
    //DP epsilon = 1.0/L; // using Lorentzian for density in entropy.

    // Construct the finite-size saddle-point state:
    // if we refine, read the quantum numbers of the saddle point state (and seed sps) from the sps file:

    /* TO_STR
    stringstream SPS_stringstream;  string SPS_string;
    Data_File_Name (SPS_stringstream, whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "");
    SPS_stringstream << ".sps";
    SPS_string = SPS_stringstream.str();    const char* SPS_Cstr = SPS_string.c_str();
    */
    string SPS_str = Data_Filename (whichDSF, c_int, L, N, iKmin, iKmax, kBT, 0.0, "") + ".sps";

    fstream spsfile;
    if (refine) spsfile.open(SPS_str, fstream::in);
    else spsfile.open(SPS_str, fstream::out | fstream::trunc);
    if (spsfile.fail()) {
      cout << SPS_str << endl; ABACUSerror("Could not open spsfile.");
    }

    LiebLin_Bethe_State spstate;

    if (!refine) { // obtain the sps from discretized TBA
      //spstate = Canonical_Saddle_Point_State (c_int, L, N, whichDSF == 'Z' ? 0.0 : kBT);
      spstate = Canonical_Saddle_Point_State (c_int, L, N, kBT);
    }
    else { // read it from the sps file
      // Check that the sps has the right number of Ix2:
      int Nspsread;
      spsfile >> Nspsread;
      if (Nspsread != N) {
	cout << Nspsread << "\t" << N << endl;
	ABACUSerror("Wrong number of Ix2 in saddle-point state.");
      }
      spstate = LiebLin_Bethe_State (c_int, L, N);
      for (int i = 0; i < N; ++i) spsfile >> spstate.Ix2[i];
    }
    spstate.Compute_All(true);

    int Nscan = N;
    if (whichDSF == 'o') Nscan = N - 1;
    if (whichDSF == 'g') Nscan = N + 1;

    // Now construct or read off the seed scan state:

    // TO MODIFY: this is not a good idea, since this might construct a state with many p-h w/r to the AveragingState.
    LiebLin_Bethe_State SeedScanState;

    if (whichDSF != 'o' && whichDSF != 'g') SeedScanState = spstate;

    else if (whichDSF == 'o' || whichDSF == 'g') {
      if (!refine) {
	if (whichDSF == 'o') SeedScanState = Remove_Particle_at_Center (spstate);
	else SeedScanState = Add_Particle_at_Center (spstate);
      }
      else { // read it from the sps file
	// Check that the sps has the right number of Ix2:
	int Nsspsread;
	spsfile >> Nsspsread;
	if (Nsspsread != Nscan) {
	  cout << Nsspsread << "\t" << Nscan << endl;
	  ABACUSerror("Wrong number of Ix2 in scan saddle-point state.");
	}
	SeedScanState = LiebLin_Bethe_State (c_int, L, Nscan);
	for (int i = 0; i < Nscan; ++i) spsfile >> SeedScanState.Ix2[i];
      }
    } // if one-body or Green's function

    SeedScanState.Compute_All(true);

    LiebLin_Bethe_State ScanState = SeedScanState;

    DP delta = sqrt(DP(N)) * (spstate.lambdaoc[N-1] - spstate.lambdaoc[0])/N;

    if (!refine) { // we write data to the sps file

      spsfile << N << endl;
      spsfile << spstate.Ix2 << endl;
      spsfile << Nscan << endl;
      spsfile << SeedScanState.Ix2 << endl;

      spsfile << endl << spstate << endl << endl;
      for (int i = 1; i < spstate.N - 2; ++i)
	spsfile << 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1])
		<< "\t" << 1.0/spstate.L * (0.25/(spstate.lambdaoc[i] - spstate.lambdaoc[i-1])
					    + 0.5/(spstate.lambdaoc[i+1] - spstate.lambdaoc[i])
					    + 0.25/(spstate.lambdaoc[i+2] - spstate.lambdaoc[i+1]))
		<< "\t" << rho_of_lambdaoc_1 (spstate, 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]), delta)
		<< "\t" << rho_of_lambdaoc_2 (spstate, 0.5 * (spstate.lambdaoc[i] + spstate.lambdaoc[i+1]), delta)
		<< endl;
    }
    spsfile.close();

    // Perform the scan:
    return General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT, spstate, SeedScanState, "",
			 Max_Secs, target_sumrule, refine, paralevel, rank, 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 = 0;
    Vect<int> rank(0,1);
    Vect<int> nr_processors(0,1);

    return Scan_LiebLin (whichDSF, c_int, L, N, iKmin, iKmax, kBT, Max_Secs, target_sumrule,
			 refine, paralevel, rank, nr_processors);
  }


  // Scanning on an excited state defined by a set of Ix2:
  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 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

    // PRECONDITIONS:
    // - the Ix2 of AveragingState are properly set.

    DP c_int = AveragingState.c_int;
    DP L = AveragingState.L;
    int N = AveragingState.N;

    // The label of the Averaging State is by definition the `empty' label
    AveragingState.Set_Label_from_Ix2 (AveragingState.Ix2);
    AveragingState.Compute_All(true);

    int Nscan = N;
    if (whichDSF == 'o') Nscan = N - 1;
    if (whichDSF == 'g') Nscan = N + 1;

    LiebLin_Bethe_State SeedScanState (c_int, L, Nscan);
    if (whichDSF == 'd' || whichDSF == 'B') SeedScanState.Ix2 = AveragingState.Ix2;
    // If 'o', remove midmost and shift quantum numbers by half-integer towards removed one:
    if (whichDSF == 'o') {
      for (int i = 0; i < N-1; ++i)
	SeedScanState.Ix2[i] = AveragingState.Ix2[i + (i >= N/2)] + 1 - 2*(i >= N/2);
    }
    // If 'g', add a quantum number in middle (explicitly: to right of index N/2)
    // and shift quantum numbers by half-integer away from added one:
    if (whichDSF == 'g') {
      SeedScanState.Ix2[N/2] = AveragingState.Ix2[N/2] - 1;
      for (int i = 0; i < N+1; ++i)
	SeedScanState.Ix2[i + (i >= N/2)] = AveragingState.Ix2[i] - 1 + 2*(i >= N/2);
    }
    SeedScanState.Compute_All(true);

    SeedScanState.Set_Label_from_Ix2 (SeedScanState.Ix2);

    DP kBT = 0.0;

    // Perform the scan:
    return General_Scan (whichDSF, iKmin, iKmax, 100000000, kBT,
			 AveragingState, SeedScanState, defaultScanStatename,
			 Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
  }

  // Simplified function call of the above:
  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 paralevel = 0;
    Vect<int> rank(0,1);
    Vect<int> nr_processors(0,1);

    return Scan_LiebLin (whichDSF, AveragingState, defaultScanStatename, iKmin, iKmax, Max_Secs,
			 target_sumrule, refine, paralevel, rank, nr_processors);
  }

  // Scanning on a previously-defined AveragingState
  Scan_Info Scan_Heis (char whichDSF, XXZ_Bethe_State& AveragingState, string defaultScanStatename,
		       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

    // PRECONDITIONS:
    // - the Ix2 of AveragingState are properly set.

    // Prepare the AveragingState:
    AveragingState.Compute_All(true);

    XXZ_Bethe_State SeedScanState;
    if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState;
    else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState);
    else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState);
    else ABACUSerror("Unknown whichDSF in Scan_Heis.");

    // Now the scan itself
    General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
		  defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);

  }

  // Scanning on a previously-defined AveragingState
  Scan_Info Scan_Heis (char whichDSF, XXX_Bethe_State& AveragingState, string defaultScanStatename,
		       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

    // PRECONDITIONS:
    // - the Ix2 of AveragingState are properly set.

    // Prepare the AveragingState:
    AveragingState.Compute_All(true);

    XXX_Bethe_State SeedScanState;
    if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState;
    else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState);
    else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState);
    else ABACUSerror("Unknown whichDSF in Scan_Heis.");

    // Now the scan itself
    return General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
			 defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);

  }

  // Scanning on a previously-defined AveragingState
  Scan_Info Scan_Heis (char whichDSF, XXZ_gpd_Bethe_State& AveragingState, string defaultScanStatename,
		       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

    // PRECONDITIONS:
    // - the Ix2 of AveragingState are properly set.

    // Prepare the AveragingState:
    AveragingState.Compute_All(true);

    XXZ_gpd_Bethe_State SeedScanState;
    if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = AveragingState;
    else if (whichDSF == 'm') SeedScanState = Remove_Particle_at_Center (AveragingState);
    else if (whichDSF == 'p') SeedScanState = Add_Particle_at_Center (AveragingState);
    else ABACUSerror("Unknown whichDSF in Scan_Heis.");

    // Now the scan itself
    return General_Scan (whichDSF, iKmin, iKmax, AveragingState.chain.Nsites, 0.0, AveragingState, SeedScanState,
			 defaultScanStatename, Max_Secs, target_sumrule, refine, paralevel, rank, 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,
		       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.

    // whichDSF == 'Z':  canonical partition function
    // whichDSF == 'm':  S^{-+}
    // whichDSF == 'z':  S^{zz}
    // whichDSF == 'p':  S^{+-}
    // whichDSF == 'a':  < S^z_j S^z_{j+1} S^z_l S^z_{l+1} > for RIXS
    // whichDSF == 'b':  < S^z_j S^-_{j+1} S^-_l S^z_{l+1} > + (m <-> z) for RIXS
    // whichDSF == 'c':  < S^-_j S^-_{j+1} S^-_l S^-_{l+1} > for RIXS

    Heis_Chain BD1(1.0, Delta, 0.0, N);

    Vect_INT Nrapidities_groundstate(0, BD1.Nstrings);

    Nrapidities_groundstate[0] = M;

    Heis_Base baseconfig_groundstate(BD1, Nrapidities_groundstate);


    if ((Delta > 0.0) && (Delta < 1.0)) {

      XXZ_Bethe_State GroundState(BD1, baseconfig_groundstate);
      GroundState.Compute_All(true);
      // The ground state is now fully defined.

      XXZ_Bethe_State SeedScanState;
      if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = GroundState;
      else if (whichDSF == 'm') SeedScanState = XXZ_Bethe_State(GroundState.chain, M - 1);
      else if (whichDSF == 'p') SeedScanState = XXZ_Bethe_State(GroundState.chain, M + 1);
      else ABACUSerror("Unknown whichDSF in Scan_Heis.");

      // Now the scan itself
      General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
		    Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);

    }

    else if (Delta == 1.0) {

      XXX_Bethe_State GroundState(BD1, baseconfig_groundstate);
      GroundState.Compute_All(true);
      // The ground state is now fully defined.

      XXX_Bethe_State SeedScanState;
      if (whichDSF == 'Z' || whichDSF == 'z' || whichDSF == 'a' || whichDSF == 'q') SeedScanState = GroundState;
      else if (whichDSF == 'm') SeedScanState = XXX_Bethe_State(GroundState.chain, M - 1);
      else if (whichDSF == 'p') SeedScanState = XXX_Bethe_State(GroundState.chain, M + 1);
      else if (whichDSF == 'c') SeedScanState = XXX_Bethe_State(GroundState.chain, M - 2);
      else ABACUSerror("Unknown whichDSF in Scan_Heis.");

      // Now the scan itself
      General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
		    Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
    }

    else if (Delta > 1.0) {

      XXZ_gpd_Bethe_State GroundState(BD1, baseconfig_groundstate);
      GroundState.Compute_All(true);
      // The ground state is now fully defined.

      XXZ_gpd_Bethe_State SeedScanState;
      if (whichDSF == 'Z' || whichDSF == 'z') SeedScanState = GroundState;
      else if (whichDSF == 'm') SeedScanState = XXZ_gpd_Bethe_State(GroundState.chain, M - 1);
      else if (whichDSF == 'p') SeedScanState = XXZ_gpd_Bethe_State(GroundState.chain, M + 1);
      else ABACUSerror("Unknown whichDSF in Scan_Heis.");

      // Now the scan itself
      return General_Scan (whichDSF, iKmin, iKmax, N, 0.0, GroundState, SeedScanState, "",
			   Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
    }

    else ABACUSerror("Delta out of range in Heis_Structure_Factor");

    return Scan_Info();
  }

  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 = 0;
    Vect<int> rank(0,1);
    Vect<int> nr_processors(0,1);

    return Scan_Heis (whichDSF, Delta, N, M, iKmin, iKmax,
		      Max_Secs, target_sumrule, refine, paralevel, rank, nr_processors);
  }


} // namespace ABACUS