ABACUS/src/HEIS/Heis.cc

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

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

Copyright (c) J.-S. Caux.

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

File:  src/HEIS/Heis.cc

Purpose:  defines functions in all HEIS classes.

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

#include "ABACUS.h"

using namespace std;

namespace ABACUS {

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

  // Function definitions:  class Heis_Chain

  Heis_Chain::Heis_Chain()
    : J(0.0), Delta(0.0), anis(0.0), hz(0.0), Nsites(0), Nstrings(0), Str_L(0), par(0),
      si_n_anis_over_2(new DP[1]), co_n_anis_over_2(new DP[1]), ta_n_anis_over_2(new DP[1]), prec(ITER_REQ_PREC) {}

  Heis_Chain::Heis_Chain (DP JJ, DP DD, DP hh, int NN)
    : J(JJ), Delta(DD), anis(0.0), hz(hh), Nsites(NN), Nstrings(MAXSTRINGS), Str_L(new int[MAXSTRINGS]), par(new int[MAXSTRINGS]),
      si_n_anis_over_2(new DP[10*MAXSTRINGS]), co_n_anis_over_2(new DP[10*MAXSTRINGS]), ta_n_anis_over_2(new DP[10*MAXSTRINGS]), prec(ITER_REQ_PREC)
  {
    // We restrict to even chains everywhere

    if (NN % 2) ABACUSerror("Please use an even-length chain.");

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

      // gapless XXZ case

      anis = acos(DD);

      // Set the Str_L and par vectors:

      DP gammaoverpi = acos(DD)/PI;
      //cout << "gammaoverpi = " << gammaoverpi << endl;

      Vect<int> Nu(MAXSTRINGS);

      DP gammaoverpi_eff = gammaoverpi;
      DP gammaoverpi_reached = 0.0;

      // First, define the series Nu[i], like in TakahashiBOOK eqn 9.35
      // CAREFUL:  the labelling is different, Nu[0] is \nu_1, etc.
      int l = 0;
      int ml_temp = 0;  // checks the sum of all Nu's

      while (fabs(gammaoverpi - gammaoverpi_reached) > 1.0e-8){

	Nu[l] = int(1.0/gammaoverpi_eff);
	ml_temp += Nu[l];
	if (Nu[l] == 0) gammaoverpi_eff = 0.0;
	else gammaoverpi_eff = 1.0/gammaoverpi_eff - DP(Nu[l]);

	// compute gammaoverpi_reached:
	gammaoverpi_reached = Nu[l];
	for (int p = 0; p < l; ++p) gammaoverpi_reached = Nu[l - p - 1] + 1.0/gammaoverpi_reached;
	gammaoverpi_reached = 1.0/gammaoverpi_reached;

	//cout << "gammaoverpi_reached = " << gammaoverpi_reached << "\tdiff = " << fabs(gammaoverpi - gammaoverpi_reached) << endl;
	//cout << "ml_temp = " << ml_temp << "\tMAXSTRINGS = " << MAXSTRINGS << endl;

	l++;

	if (ml_temp > MAXSTRINGS) break;  // we defined Str_L and par as arrays of at most MAXSTRINGS elements, so we cut off here...

      }
      //cout << "l = " << l << endl;

      // Check:  make sure the last Nu is greater than one:  if one, add 1 to previous Nu

      if (l > 1) {
	if (Nu[l-1] == 1) {
	  Nu[l-1] = 0;
	  Nu[l-2] += 1;
	  l -= 1;
	}
      }

      // Length of continued fraction is l-1, which is denoted l in Takahashi

      // Second, define the series y[i] and m[i] as in TakahashiBOOK eqn 9.36
      // y_{-1} is treated separately, and here y[i] = y_i, m[i] = m_i, i = 0, ..., l
      Vect<int> y(0, l+1);
      Vect<int> m(0, l+1);
      y[0] = 1;
      m[0] = 0;
      y[1] = Nu[0];
      m[1] = Nu[0];
      for (int i = 2; i <= l; ++i) {
	y[i] = y[i - 2] + Nu[i-1]*y[i-1];
	m[i] = 0;
	for (int k = 0; k < i; ++k) m[i] += Nu[k];
      }
      // Now determine the lengths and parity, following TakahashiBOOK eqn 9.37
      //      Nstrings = ABACUS::min(m[l] + 1, MAXSTRINGS);  // number of different strings that are possible for this chain
      // Always set to MAXSTRINGS

      Str_L[0] = 1;
      par[0] = 1;
      Str_L[1] = 1;
      par[1] = -1;

      int max_j = 0;

      for (int i = 0; i < l; ++i) {
	for (int j = ABACUS::max(1, m[i]) + 1; j < ABACUS::min(m[i+1] + 1, MAXSTRINGS); ++j) {
	  if (i == 0) Str_L[j] = (j - m[0])* y[0];
	  else if (i >= 1) Str_L[j] = y[i-1] + (j - m[i])*y[i];
	  par[j] = (int(floor((Str_L[j] - 1)*gammaoverpi)) % 2) ? -1 : 1;
	  max_j = j;
	}
      }

      // Set the rest of Str_L and par vector elements to zero

      for (int i = max_j + 1; i < MAXSTRINGS; ++i) {
	Str_L[i] = 0;
	par[i] = 0;
      }

      // Set the sin, cos and tan_n_zeta_over_2 vectors:

      for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = sin(i * anis/2.0);
      for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = cos(i * anis/2.0);
      for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = tan(i * anis/2.0);

    } // if XXZ gapless

    else if (Delta == 1.0) {

      // Isotropic antiferromagnet

      anis = 0.0;

      // Set Str_L and par:

      for (int i = 0; i < MAXSTRINGS; ++i) {
	Str_L[i] = i + 1;
	par[i] = 1;
      }

    } // if XXX AFM

    else if (Delta > 1.0) {

      // Gapped antiferromagnet

      anis = acosh(DD);

      // Set the Str_L and par vectors:

      for (int i = 0; i < MAXSTRINGS; ++i) {
	Str_L[i] = i + 1;
	par[i] = 1;
      }

      // Set the sinh, cosh and tanh_n_eta_over_2 vectors:

      for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = sinh(i * anis/2.0);
      for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = cosh(i * anis/2.0);
      for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = tanh(i * anis/2.0);

    } // if XXZ_gpd

  }

  Heis_Chain::Heis_Chain (const Heis_Chain& RefChain)  // copy constructor
  {
    J = RefChain.J;
    Delta = RefChain.Delta;
    anis = RefChain.anis;
    hz = RefChain.hz;
    Nsites = RefChain.Nsites;
    Nstrings = RefChain.Nstrings;
    Str_L = new int[RefChain.Nstrings];
    for (int i = 0; i < RefChain.Nstrings; ++i) Str_L[i] = RefChain.Str_L[i];
    par = new int[RefChain.Nstrings];
    for (int i = 0; i < RefChain.Nstrings; ++i) par[i] = RefChain.par[i];
    si_n_anis_over_2 = new DP[10*MAXSTRINGS];
    for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = RefChain.si_n_anis_over_2[i];
    co_n_anis_over_2 = new DP[10*MAXSTRINGS];
    for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = RefChain.co_n_anis_over_2[i];
    ta_n_anis_over_2 = new DP[10*MAXSTRINGS];
    for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = RefChain.ta_n_anis_over_2[i];
    prec = RefChain.prec;
  }

  Heis_Chain& Heis_Chain::operator= (const Heis_Chain& RefChain)  // assignment operator
  {
    if (this != &RefChain) {
      J = RefChain.J;
      Delta = RefChain.Delta;
      anis = RefChain.anis;
      hz = RefChain.hz;
      Nsites = RefChain.Nsites;
      Nstrings = RefChain.Nstrings;
      if (Str_L != 0) delete[] Str_L;
      Str_L = new int[RefChain.Nstrings];
      for (int i = 0; i < RefChain.Nstrings; ++i) Str_L[i] = RefChain.Str_L[i];
      if (par != 0) delete[] par;
      par = new int[RefChain.Nstrings];
      for (int i = 0; i < RefChain.Nstrings; ++i) par[i] = RefChain.par[i];
      if (si_n_anis_over_2 != 0) delete[] si_n_anis_over_2;
      si_n_anis_over_2 = new DP[10*MAXSTRINGS];
      for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = RefChain.si_n_anis_over_2[i];
      if (co_n_anis_over_2 != 0) delete[] co_n_anis_over_2;
      co_n_anis_over_2 = new DP[10*MAXSTRINGS];
      for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = RefChain.co_n_anis_over_2[i];
      if (ta_n_anis_over_2 != 0) delete[] ta_n_anis_over_2;
      ta_n_anis_over_2 = new DP[10*MAXSTRINGS];
      for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = RefChain.ta_n_anis_over_2[i];
      prec = RefChain.prec;

    }
    return *this;
  }

  bool Heis_Chain::operator== (const Heis_Chain& RefChain)
  {
    bool answer = true;

    if ((J != RefChain.J) || (Nsites != RefChain.Nsites) || (Delta != RefChain.Delta)) answer = false;

    return(answer);
  }

  bool Heis_Chain::operator!= (const Heis_Chain& RefChain)
  {
    bool answer = false;

    if ((J != RefChain.J) || (Nsites != RefChain.Nsites) || (Delta != RefChain.Delta)) answer = true;

    return(answer);
  }

  Heis_Chain::~Heis_Chain()
  {
    if (Str_L != 0) delete[] Str_L;
    if (par != 0) delete[] par;
    if (si_n_anis_over_2 != 0) delete[] si_n_anis_over_2;
    if (co_n_anis_over_2 != 0) delete[] co_n_anis_over_2;
    if (ta_n_anis_over_2 != 0) delete[] ta_n_anis_over_2;
  }

  /* Deactivated in ++G_8
  void Heis_Chain::Scan_for_Possible_Bases (int Mdown_remaining, Vect<string>& possible_base_label, int& nfound, int nexc_max_used,
					    int base_level_to_scan, Vect<int>& Nrapidities)
  {
    if (Mdown_remaining < 0) { ABACUSerror("Scan_for_Possible_Bases:  shouldn't be here..."); }  // reached inconsistent point

    //cout << "Mdown_remaining " << Mdown_remaining << "\t" << possible_base_id << "\tnfound " << nfound
    // << "\tnexc_max_used " << nexc_max_used << "\tbase_level_to_scan " << base_level_to_scan << "\tNrap " << Nrapidities << endl;

    if (base_level_to_scan == 0) {
      Nrapidities[0] = Mdown_remaining;

      // Set label:
      stringstream M0out;
      M0out << Nrapidities[0];
      possible_base_label[nfound] = M0out.str();
      for (int itype = 1; itype < Nrapidities.size(); ++itype)
	if (Nrapidities[itype] > 0) {
	  possible_base_label[nfound] += TYPESEP;
	  stringstream typeout;
	  typeout << itype;
	  possible_base_label[nfound] += typeout.str();
	  possible_base_label[nfound] += EXCSEP;
	  stringstream Mout;
	  Mout << Nrapidities[itype];
	  possible_base_label[nfound] += Mout.str();
	}
      nfound++;
    }

    else {
      // Remove all higher strings:
      //Nrapidities[base_level_to_scan] = 0;
      //Scan_for_Possible_Bases (Mdown_remaining, possible_base_id, nfound, nexc_max_used, base_level_to_scan - 1, Nrapidities);

      for (int i = 0; i <= (Str_L[base_level_to_scan] == 0 ? 0 : nexc_max_used/Str_L[base_level_to_scan]); ++i) {
	Nrapidities[base_level_to_scan] = i;
	Scan_for_Possible_Bases (Mdown_remaining - i*Str_L[base_level_to_scan], possible_base_label, nfound,
				 nexc_max_used - i*Str_L[base_level_to_scan], base_level_to_scan - 1, Nrapidities);
      }

    }
  }


  Vect<string> Heis_Chain::Possible_Bases (int Mdown)  // returns a vector of possible bases
  {
    // We partition Mdown into up to NEXC_MAX_HEIS excitations

    int nexc_max_used = ABACUS::min(NEXC_MAX_HEIS, 2*(Mdown/2));  // since each inner sector can contain at most N/2 holes.

    Vect<string> possible_base_label (1000);
    int nfound = 0;
    Vect<int> Nrapidities (0, Nstrings);

    //cout << "In Possible_Bases:  start scan for Mdown = " << Mdown << endl;

    (*this).Scan_for_Possible_Bases (Mdown, possible_base_label, nfound, nexc_max_used, Nstrings - 1, Nrapidities);

    // Copy results into a clean vector:
    Vect<string> possible_base_label_found (nfound);
    for (int i = 0; i < nfound; ++i) possible_base_label_found[i] = possible_base_label[i];

    //cout << "In Possible_Bases:  possible_base_label_found = " << possible_base_label_found << endl;

    return(possible_base_label_found);
  }
  */

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

  // Function definitions:  class Heis_Base

  Heis_Base::Heis_Base () : Mdown(0), Nrap(Vect<int>()), Nraptot(0),
			    Ix2_infty(Vect<DP>()), Ix2_min(Vect<int>()), Ix2_max(Vect<int>()), dimH(0.0), baselabel("") {}

  Heis_Base::Heis_Base (const Heis_Base& RefBase)  // copy constructor
    : Mdown(RefBase.Mdown), Nrap(Vect<int>(RefBase.Nrap.size())), Nraptot(RefBase.Nraptot),
      Ix2_infty(Vect<DP>(RefBase.Nrap.size())), Ix2_min(Vect<int>(RefBase.Nrap.size())), Ix2_max(Vect<int>(RefBase.Nrap.size())),
      baselabel(RefBase.baselabel)
  {
    for (int i = 0; i < Nrap.size(); ++i) {
      Nrap[i] = RefBase.Nrap[i];
      Ix2_infty[i] = RefBase.Ix2_infty[i];
      Ix2_min[i] = RefBase.Ix2_min[i];
      Ix2_max[i] = RefBase.Ix2_max[i];
      dimH = RefBase.dimH;
    }
  }

  Heis_Base::Heis_Base (const Heis_Chain& RefChain, int M)
    : Mdown(M), Nrap(Vect<int>(RefChain.Nstrings)), Nraptot(0),
      Ix2_infty(Vect<DP>(RefChain.Nstrings)), Ix2_min(Vect<int>(RefChain.Nstrings)), Ix2_max(Vect<int>(RefChain.Nstrings))
  {
    for (int i = 0; i < RefChain.Nstrings; ++i) Nrap[i] = 0;
    Nrap[0] = M;

    Nraptot = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];

    stringstream M0out;
    M0out << M;
    baselabel = M0out.str();

    // Now compute the Ix2_infty numbers
    (*this).Compute_Ix2_limits(RefChain);

    // Compute dimensionality of this sub-Hilbert space:
    complex<double> ln_dimH_cx = 0.0;
    for (int i = 0; i < RefChain.Nstrings; ++i)
      if (Nrap[i] > 0) ln_dimH_cx += ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2)) - ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2 - Nrap[i])) - ln_Gamma(complex<double>(Nrap[i] + 1));
    dimH = exp(real(ln_dimH_cx));
  }

  Heis_Base::Heis_Base (const Heis_Chain& RefChain, const Vect<int>& Nrapidities)
    : Mdown(0), Nrap(Nrapidities), Nraptot(0),
      Ix2_infty(Vect<DP>(RefChain.Nstrings)), Ix2_min(Vect<int>(RefChain.Nstrings)), Ix2_max(Vect<int>(RefChain.Nstrings))
  {

    // Check consistency of Nrapidities vector with RefChain

    //if (RefChain.Nstrings != Nrapidities.size()) cout << "error:  Nstrings = " << RefChain.Nstrings << "\tNrap.size = " << Nrapidities.size() << endl;

    if (RefChain.Nstrings != Nrapidities.size()) ABACUSerror("Incompatible Nrapidities vector used in Heis_Base constructor.");

    int Mcheck = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
    Mdown = Mcheck;

    Nraptot = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];

    /*
    // Compute id
    id += Nrapidities[0];
    long long int factor = 100000LL;
    for (int i = 1; i < RefChain.Nstrings; ++i) {
      id += factor * Nrapidities[i];
      factor *= 100LL;
    }
    */

    // Set label:
    /*
    stringstream baselabel_strstream;
    baselabel_strstream << Nrapidities[0];
    for (int itype = 1; itype < Nrapidities.size(); ++itype)
      if (Nrapidities[itype] > 0) {
	baselabel_strstream << TYPESEP;
	stringstream typeout;
	typeout << itype;
	baselabel += typeout.str();
	baselabel += 'EXCSEP';
	stringstream Mout;
	Mout << Nrapidities[itype];
	baselabel += Mout.str();
      }
    */

    stringstream M0out;
    M0out << Nrapidities[0];
    baselabel = M0out.str();
    for (int itype = 1; itype < Nrapidities.size(); ++itype)
      if (Nrapidities[itype] > 0) {
	baselabel += TYPESEP;
	stringstream typeout;
	typeout << itype;
	baselabel += typeout.str();
	baselabel += EXCSEP;
	stringstream Mout;
	Mout << Nrapidities[itype];
	baselabel += Mout.str();
      }

    // Now compute the Ix2_infty numbers
    (*this).Compute_Ix2_limits(RefChain);

    // Compute dimensionality of this sub-Hilbert space:
    complex<double> ln_dimH_cx = 0.0;
    for (int i = 0; i < RefChain.Nstrings; ++i)
      if (Nrap[i] > 0) ln_dimH_cx += ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2)) - ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2 - Nrap[i])) - ln_Gamma(complex<double>(Nrap[i] + 1));
    dimH = exp(real(ln_dimH_cx));
  }

  Heis_Base::Heis_Base (const Heis_Chain& RefChain, string baselabel_ref)
    : Mdown(0), Nrap(Vect<int>(0, RefChain.Nstrings)), Nraptot(0),
      Ix2_infty(Vect<DP>(RefChain.Nstrings)), Ix2_min(Vect<int>(RefChain.Nstrings)), Ix2_max(Vect<int>(RefChain.Nstrings)),
      baselabel (baselabel_ref)
  {
    // Build Nrapidities vector from baselabel_ref.
    // This is simply done by using the state label standard reading function after conveniently
    // making baselabel into a label (as for a state with no excitations):
    string label_ref = baselabel + "_0_";
    Vect<Vect<int> > dummyOriginIx2(1);

    //cout << "Trying to build base from baselabel_ref " << baselabel_ref << "\t and label_ref " << label_ref << endl;
    //State_Label_Data labeldata = Read_State_Label (label_ref, dummyOriginIx2);
    State_Label_Data labeldata = Read_Base_Label (label_ref);
    //cout << "Read data." << endl;

    // Initialize Nrap:
    for (int i = 0; i < labeldata.type.size(); ++i)
      Nrap[labeldata.type[i] ] = labeldata.M[i];

    int Mcheck = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
    Mdown = Mcheck;

    Nraptot = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];

    // Now compute the Ix2_infty numbers
    (*this).Compute_Ix2_limits(RefChain);

    // Compute dimensionality of this sub-Hilbert space:
    complex<double> ln_dimH_cx = 0.0;
    for (int i = 0; i < RefChain.Nstrings; ++i)
      if (Nrap[i] > 0) ln_dimH_cx += ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2)) - ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2 - Nrap[i])) - ln_Gamma(complex<double>(Nrap[i] + 1));
    dimH = exp(real(ln_dimH_cx));
  }

  Heis_Base& Heis_Base::operator= (const Heis_Base& RefBase)
  {
    if (this != & RefBase) {
      Mdown = RefBase.Mdown;
      Nrap = RefBase.Nrap;
      Nraptot = RefBase.Nraptot;
      Ix2_infty = RefBase.Ix2_infty;
      Ix2_min = RefBase.Ix2_min;
      Ix2_max = RefBase.Ix2_max;
      dimH = RefBase.dimH;
      baselabel = RefBase.baselabel;
    }
    return(*this);
  }

  bool Heis_Base::operator== (const Heis_Base& RefBase)
  {
    bool answer = (Nrap == RefBase.Nrap);

    return (answer);
  }

  bool Heis_Base::operator!= (const Heis_Base& RefBase)
  {
    bool answer = (Nrap != RefBase.Nrap);

    return (answer);
  }

  void Heis_Base::Compute_Ix2_limits (const Heis_Chain& RefChain)
  {

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

      // Compute the Ix2_infty numbers

      DP sum1 = 0.0;
      DP sum2 = 0.0;

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	sum1 = 0.0;

	for (int k = 0; k < RefChain.Nstrings; ++k) {

	  sum2 = 0.0;

	  sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 : 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
										  - 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis));
	  sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
				 - 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis));

	  for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a)
	    sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
					 - 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis));

	  sum1 += (Nrap[k] - ((j == k) ? 1 : 0)) * sum2;
	}

	Ix2_infty[j] = (1.0/PI) * fabs(RefChain.Nsites * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j])
								      - 0.5 * RefChain.Str_L[j] * RefChain.anis)) - sum1);

      }  // The Ix2_infty are now set.

      // Now compute the Ix2_max limits

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	Ix2_max[j] = int(floor(Ix2_infty[j]));  // sets basic integer

	// Reject formally infinite rapidities (i.e. if Delta is root of unity)

	//cout << "Ix2_infty - Ix2_max = " << Ix2_infty[j] - Ix2_max[j] << endl;
	//if (Ix2_infty[j] == Ix2_max[j]) {
	//Ix2_max[j] -= 2;
	//}
	// If Nrap is even, Ix2_max must be odd.  If odd, then even.

	if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;

	while (Ix2_max[j] > RefChain.Nsites) {
	  Ix2_max[j] -= 2;
	}
	Ix2_min[j] = -Ix2_max[j];
      }
    } // if XXZ gapless

    else if (RefChain.Delta == 1.0) {

      // Compute the Ix2_infty numbers

      int sum1 = 0;

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	sum1 = 0;

	for (int k = 0; k < RefChain.Nstrings; ++k) {

	  sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
	}

	//Ix2_infty[j] = (RefChain.Nsites - 1.0 + 2.0 * RefChain.Str_L[j] - sum1);
	Ix2_infty[j] = (RefChain.Nsites + 1.0 - sum1); // to get counting right...

      }  // The Ix2_infty are now set.

      // Now compute the Ix2_max limits

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	Ix2_max[j] = int(floor(Ix2_infty[j]));  // sets basic integer

	// Give the correct parity to Ix2_max

	// If Nrap is even, Ix2_max must be odd.  If odd, then even.

	if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;

	// If Ix2_max equals Ix2_infty, we reduce it by 2:

	if (Ix2_max[j] == int(Ix2_infty[j])) Ix2_max[j] -= 2;

	while (Ix2_max[j] > RefChain.Nsites) {
	  Ix2_max[j] -= 2;
	}
	Ix2_min[j] = -Ix2_max[j];
      }

    } // if XXX AFM

    else if (RefChain.Delta > 1.0) {

      // Compute the Ix2_infty numbers

      int sum1 = 0;

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	sum1 = 0;

	for (int k = 0; k < RefChain.Nstrings; ++k) {

	  sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
	}

	Ix2_infty[j] = (RefChain.Nsites - 1 + 2 * RefChain.Str_L[j] - sum1);

      }  // The Ix2_infty are now set.

      // Now compute the Ix2_max limits

      for (int j = 0; j < RefChain.Nstrings; ++j) {

	Ix2_max[j] = int(floor(Ix2_infty[j]));  // sets basic integer

	// Give the correct parity to Ix2_max

	// If Nrap is even, Ix2_max must be odd.  If odd, then even.

	if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;

	// If Ix2_max equals Ix2_infty, we reduce it by 2:

	//if (Ix2_max[j] == Ix2_infty[j]) Ix2_max[j] -= 2;

	while (Ix2_max[j] > RefChain.Nsites) {
	  Ix2_max[j] -= 2;
	}

	// Fudge, for strings:
	//if (RefChain.Str_L[j] >= 1) Ix2_max[j] += 2;
	//Ix2_max[j] += 2;

	Ix2_min[j] = -Ix2_max[j];
      }

    } // if XXZ_gpd

  }



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

  // Function definitions:  class Lambda

  Lambda::Lambda () : lambda(NULL) {}

  Lambda::Lambda (const Heis_Chain& RefChain, int M)
    : Nstrings(1), Nrap(Vect<int>(M,1)), Nraptot(M), lambda(new DP*[1]) // single type of string here
      //: lambda(Vect<Vect<DP> > (1))
  {
    lambda[0] = new DP[M];
    //lambda[0] = Vect<DP> (M);

    for (int j = 0; j < M; ++j) lambda[0][j] = 0.0;
  }

  Lambda::Lambda (const Heis_Chain& RefChain, const Heis_Base& base)
    : Nstrings(RefChain.Nstrings), Nrap(base.Nrap), Nraptot(base.Nraptot), lambda(new DP*[RefChain.Nstrings])
      //: lambda(Vect<Vect<DP> > (RefChain.Nstrings))
  {
    //lambda[0] = new DP[base.Mdown];
    lambda[0] = new DP[base.Nraptot];
    for (int i = 1; i < RefChain.Nstrings; ++i) lambda[i] = lambda[i-1] + base[i-1];
    //for (int i = 0; i < RefChain.Nstrings; ++i) lambda[i] = Vect<DP> (base[i]);

    for (int i = 0; i < RefChain.Nstrings; ++i) {
      for (int j = 0; j < base[i]; ++j) lambda[i][j] = 0.0;
    }
  }

  Lambda& Lambda::operator= (const Lambda& RefLambda)
  {
    if (this != &RefLambda) {

      Nstrings = RefLambda.Nstrings;
      Nrap = RefLambda.Nrap;
      Nraptot = RefLambda.Nraptot;
      if (lambda != 0) {
	delete[] (lambda[0]);
	delete[] (lambda);
      }
      lambda = new DP*[Nstrings];
      lambda[0] = new DP[Nraptot];
      for (int i = 1; i < Nstrings; ++i) lambda[i] = lambda[i-1] + Nrap[i-1];
      for (int i = 0; i < Nstrings; ++i)
	for (int j = 0; j < Nrap[i]; ++j) lambda[i][j] = RefLambda.lambda[i][j];

   }
    return(*this);
  }

  Lambda::~Lambda()
  {
    if (lambda != 0) {
      delete[] (lambda[0]);
      delete[] lambda;
    }
  }


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

  // Function definitions:  class Heis_Bethe_State

  Heis_Bethe_State::Heis_Bethe_State ()
    : chain(Heis_Chain()), base(Heis_Base()), //offsets(Ix2_Offsets()),
      //Ix2(Ix2_Config(chain, 1)),
      Ix2(Vect<Vect<int> > (1)),
      lambda(Lambda(chain, 1)), BE(Lambda(chain, 1)), diffsq(0.0), conv(0), dev(1.0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0), //base_id(0LL), type_id(0LL), id(0LL), maxid(0LL), nparticles(0)
      label("")
  {
  };

  Heis_Bethe_State::Heis_Bethe_State (const Heis_Bethe_State& RefState) // copy constructor
    //: chain(RefState.chain), base(RefState.base), offsets(RefState.offsets), Ix2(Ix2_Config(RefState.chain, RefState.base.Mdown)),
    //  lambda(Lambda(RefState.chain, RefState.base.Mdown)), BE(Lambda(RefState.chain, RefState.base.Mdown)),
    : chain(RefState.chain), base(RefState.base), //offsets(RefState.offsets),
      //Ix2(Ix2_Config(RefState.chain, RefState.base)),
      Ix2 (RefState.Ix2),
      lambda(Lambda(RefState.chain, RefState.base)), BE(Lambda(RefState.chain, RefState.base)),
      diffsq(RefState.diffsq), conv(RefState.conv), dev(RefState.dev), iter(RefState.iter), iter_Newton(RefState.iter_Newton),
      E(RefState.E), iK(RefState.iK), K(RefState.K), lnnorm(RefState.lnnorm),
      //id(RefState.id), maxid(RefState.maxid)
      //base_id(RefState.base_id), type_id(RefState.type_id), id(RefState.id), maxid(RefState.maxid), nparticles(RefState.nparticles)
      label(RefState.label)
  {
    // copy arrays into new ones

    /*
    cout << "Here in Heis constructor state" << endl;
    cout << "lambda " << lambda[0][0] << endl;
    cout << "lambda OK" << endl;

    cout << "RefConfig: " << endl << RefState.Ix2 << endl;
    cout << "(*this).Ix2: " << endl << Ix2 << endl;
    */
    for (int j = 0; j < RefState.chain.Nstrings; ++j) {
      for (int alpha = 0; alpha < RefState.base[j]; ++j) {
	//Ix2[j][alpha] = RefState.Ix2[j][alpha]; // not needed anymore since Ix2 is Vect<Vect<int> >
	lambda[j][alpha] = RefState.lambda[j][alpha];
      }
    }
    //cout << "Heis constructor state OK" << endl;
  }

  Heis_Bethe_State::Heis_Bethe_State (const Heis_Chain& RefChain, int M)
    : chain(RefChain), base(RefChain, M), //offsets(base, 0LL),
      //Ix2(Ix2_Config(RefChain, M)),
      Ix2 (Vect<Vect<int> > (1)),
      lambda(Lambda(RefChain, M)),
      BE(Lambda(RefChain, M)), diffsq(1.0), conv(0), dev(1.0), iter(0), iter_Newton(0),
      E(0.0), iK(0), K(0.0), lnnorm(-100.0)
      //id(0LL), maxid(0LL)
      //base_id(0LL), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
  {
    Ix2[0] = Vect<int> (M);
    for (int j = 0; j < M; ++j) Ix2[0][j] = -(M - 1) + 2*j;

    stringstream M0out;
    M0out << M;
    label = M0out.str() + "_0_";
  }

  Heis_Bethe_State::Heis_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& RefBase)
    : chain(RefChain), base(RefBase), //offsets(RefBase, 0LL),
      //Ix2(Ix2_Config(RefChain, RefBase)),
      Ix2 (Vect<Vect<int> > (RefChain.Nstrings)),
      lambda(Lambda(RefChain, RefBase)),
      BE(Lambda(RefChain, RefBase)), diffsq(1.0), conv(0), dev(1.0), iter(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0)
      //id(0LL), maxid(0LL)
      //base_id(RefBase.id), type_id(0LL), id(0LL), maxid(offsets.maxid), nparticles(0)
  {
    // Check that the number of rapidities is consistent with Mdown

    //cout << "Here in Heis constructor chain base" << endl;
    //cout << "lambda " << lambda[0][0] << endl;
    //cout << "lambda OK" << endl;

    int Mcheck = 0;
    for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += base[i] * RefChain.Str_L[i];
    if (RefBase.Mdown != Mcheck) ABACUSerror("Wrong M from Nrapidities input vector, in Heis_Bethe_State constructor.");

    for (int i = 0; i < RefChain.Nstrings; ++i) Ix2[i] = Vect<int> (base[i]);

    // And now for the other string types...
    for (int i = 0; i < RefChain.Nstrings; ++i) {
      for (int j = 0; j < base[i]; ++j) Ix2[i][j] = -(base[i] - 1) + 2*j;
    }

    label = Return_State_Label (Ix2, Ix2);
  }

  void Heis_Bethe_State::Set_to_Label (string label_ref, const Vect<Vect<int> >& OriginIx2)
  {
    //cout << "Called Set_to_Label with label_ref " << label_ref << " on state " << (*this) << endl;

    //cout << "Check labeldata of label " << label_ref << endl;

    State_Label_Data labeldata = Read_State_Label (label_ref, OriginIx2);

    //cout << "type: " << labeldata.type << endl << "M: " << labeldata.M << endl << "nexc: " << labeldata.nexc << endl;
    //for (int i = 0; i < labeldata.Ix2old.size(); ++i) cout << "Ix2old[" << i << "] = " << labeldata.Ix2old[i] << endl << "Ix2exc[" << i << "] = " << labeldata.Ix2exc[i] << endl;

    label = label_ref;

    Vect<Vect<int> > OriginIx2ordered = OriginIx2;
    for (int i = 0; i < OriginIx2.size(); ++i) OriginIx2ordered[i].QuickSort();

    // Set all Ix2 to OriginIx2:
    for (int il = 0; il < OriginIx2.size(); ++il)
      for (int i = 0; i < OriginIx2[il].size(); ++i) Ix2[il][i] = OriginIx2ordered[il][i];

    // Now set the excitations:
    for (int it = 0; it < labeldata.type.size(); ++it)
      for (int iexc = 0; iexc < labeldata.nexc[it]; ++iexc)
	for (int i = 0; i < labeldata.M[it]; ++i)
	  if (Ix2[labeldata.type[it] ][i] == labeldata.Ix2old[it][iexc]) {
	    //cout << it << "\t" << iexc << "\t" << i << endl;
	    //cout << "\tSetting Ix2[" << labeldata.type[it] << "][" << i << "] to " << labeldata.Ix2exc[it][iexc] << endl;
	    //cout << Ix2[labeldata.type[it] ][i] << "\t" << labeldata.Ix2old[it][iexc] << "\t" << labeldata.Ix2exc[it][iexc] << endl;
	    Ix2[labeldata.type[it] ][i] = labeldata.Ix2exc[it][iexc];
	    //cout << Ix2[labeldata.type[it] ][i] << "\t" << labeldata.Ix2old[it][iexc] << "\t" << labeldata.Ix2exc[it][iexc] << endl;
	  }

    //cout << "State obtained(1): " << (*this) << endl;

    // Now reorder the Ix2 to follow convention:
    for (int il = 0; il < Ix2.size(); ++il) Ix2[il].QuickSort();

    //cout << "Setting label:" << label_ref << endl << "Ix2old = " << labeldata.Ix2old[0] << endl << "Ix2exc = " << labeldata.Ix2exc[0] << endl;
    //cout << "on " << OriginStateIx2ordered << endl << "giving " << Ix2 << endl;

    //cout << "State obtained(2): " << (*this) << endl;
    //char a; cin >> a;

    //(*this).Set_Label_from_Ix2 (OriginStateIx2ordered);
    //(*this).Set_Label_Internals_from_Ix2 (OriginStateIx2ordered);
  }

  void Heis_Bethe_State::Set_Label_from_Ix2 (const Vect<Vect<int> >& OriginIx2)
  {
    // This function does not assume any ordering of the Ix2.

    // ASSUMPTIONS:
    // (*this) has a base already identical to base of OriginIx2

    // First check that bases are consistent
    if ((*this).chain.Nstrings != OriginIx2.size()) ABACUSerror("Inconsistent base sizes in Heis_Bethe_State::Set_Label_from_Ix2.");

    // Then check that the filling at each level is equal
    for (int il = 0; il < (*this).chain.Nstrings; ++il)
      if ((*this).base.Nrap[il] != OriginIx2[il].size()) ABACUSerror("Inconsistent base filling in Heis_Bethe_State::Set_Label_from_Ix2.");

    // Determine how many types of particles are present
    int ntypes = 1; // level 0 always assumed present
    for (int il = 1; il < chain.Nstrings; ++il) if (base.Nrap[il] > 0) ntypes++;

    // Set the state label
    Vect<int> type_ref(0, ntypes);
    Vect<int> M_ref(ntypes);
    M_ref[0] = OriginIx2[0].size();
    type_ref[0] = 0;
    int itype = 1;
    for (int il = 1; il < chain.Nstrings; ++il)
      if (base.Nrap[il] > 0) {
	type_ref[itype] = il;
	M_ref[itype++] = OriginIx2[il].size();
      }
    Vect<int> nexc_ref(0, ntypes);
    Vect<Vect<int> > Ix2old_ref(ntypes);
    Vect<Vect<int> > Ix2exc_ref(ntypes);

    // Count nr of particle-holes at each level:
    //itype = 0;
    itype = -1;
    for (int il = 0; il < chain.Nstrings; ++il) {
      if (il == 0 || base.Nrap[il] > 0) {
	itype++;
	for (int alpha = 0; alpha < base.Nrap[il]; ++alpha)
	  if (!OriginIx2[il].includes(Ix2[il][alpha])) nexc_ref[itype] += 1;
	Ix2old_ref[itype] = Vect<int>(ABACUS::max(nexc_ref[itype], 1));
	Ix2exc_ref[itype] = Vect<int>(ABACUS::max(nexc_ref[itype], 1));
	int nexccheck = 0;
	for (int alpha = 0; alpha < base.Nrap[il]; ++alpha)
	  if (!OriginIx2[il].includes(Ix2[il][alpha])) Ix2exc_ref[itype][nexccheck++] = Ix2[il][alpha];
	if (nexccheck != nexc_ref[itype]) {
	  cout << "il = " << il << "\titype = " << itype << "\tnexccheck = " << nexccheck << "\tnexc_ref[itype] = " << nexc_ref[itype] << endl;
	  cout << OriginIx2[il] << endl << Ix2[il] << endl;
	  ABACUSerror("Counting excitations wrong (1) in Heis_Bethe_State::Set_Label_from_Ix2.");
	}
	nexccheck = 0;
	for (int alpha = 0; alpha < base.Nrap[il]; ++alpha)
	  if (!Ix2[il].includes (OriginIx2[il][alpha])) Ix2old_ref[itype][nexccheck++] = OriginIx2[il][alpha];
	if (nexccheck != nexc_ref[itype]) {
	  ABACUSerror("Counting excitations wrong (2) in Heis_Bethe_State::Set_Label_from_Ix2.");
	}
	// Now order the Ix2old_ref and Ix2exc_ref:
	Ix2old_ref[itype].QuickSort();
	Ix2exc_ref[itype].QuickSort();
      } // if (il == 0 || ...
    } // for il

    State_Label_Data labeldata(type_ref, M_ref, nexc_ref, Ix2old_ref, Ix2exc_ref);

    label = Return_State_Label (labeldata, OriginIx2);
  }

  bool Heis_Bethe_State::Check_Symmetry ()
  {
    // Checks whether the I's are symmetrically distributed.

    bool symmetric_state = true;
    int arg, test1, test3;

    //if (chain.Delta <= 1.0) {

      for (int j = 0; j < chain.Nstrings; ++j) {
	test1 = 0;
	test3 = 0;
	for (int alpha = 0; alpha < base[j]; ++alpha) {
	  arg = (Ix2[j][alpha] != chain.Nsites) ? Ix2[j][alpha] : 0;  // since Ix2 = N is same as Ix2 = -N by periodicity, this is symmetric.
	  test1 += arg;
	  test3 += arg * arg * arg;  // to make sure that all I's are symmetrical...
	}
	if (!(symmetric_state && (test1 == 0) && (test3 == 0))) symmetric_state = false;
      }
      //}
      /*
    else if (chain.Delta > 1.0) {
      // For the gapped antiferromagnet, we check symmetry *excluding* the Ix2_max values
      // The state is then inadmissible is symmetric && -Ix2_max occupied
      for (int j = 0; j < chain.Nstrings; ++j) {
	test1 = 0;
	test3 = 0;
	for (int alpha = 0; alpha < base[j]; ++alpha) {
	  arg = (Ix2[j][alpha] != chain.Nsites
		 && abs(Ix2[j][alpha]) != base.Ix2_max[j]) ? Ix2[j][alpha] : 0;  // since Ix2 = N is same as Ix2 = -N by periodicity, this is symmetric.
	  test1 += arg;
	  test3 += arg * arg * arg;  // to make sure that all I's are symmetrical...
	}
	if (!(symmetric_state && (test1 == 0) && (test3 == 0))) symmetric_state = false;
      }
    }
      */

    return symmetric_state;
  }


  // SOLVING BETHE EQUATIONS:

  void Heis_Bethe_State::Compute_diffsq ()
  {
    DP maxterm = 0.0;
    int jmax, alphamax;

    diffsq = 0.0;
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	diffsq += pow(BE[j][alpha], 2.0);
	if (pow(BE[j][alpha], 2.0)/chain.Nsites > maxterm) {
	  jmax = j;
	  alphamax = alpha;
	  maxterm = pow(BE[j][alpha], 2.0)/chain.Nsites;
	}
      }
    diffsq /= DP(chain.Nsites);
  }

  void Heis_Bethe_State::Find_Rapidities (bool reset_rapidities)
  {
    // This function finds the rapidities of the eigenstate

    //cout << endl << "Find_Rapidities called: " << (*this) << endl;

    lnnorm = -100.0;  // sentinel value, recalculated if Newton method used in the last step of iteration.

    Lambda lambda_ref(chain, base);

    diffsq = 1.0;

    if (reset_rapidities)
      (*this).Set_Free_lambdas();
    //cout << endl << "After Set_Free_lambdas: " << (*this) << endl;
    (*this).Compute_BE();
    //cout << endl << "After Compute_BE: " << (*this) << endl;


    iter = 0;
    iter_Newton = 0;

    // Start with conventional iterations

    DP iter_prec = 1.0e-2/DP(chain.Nsites);
    DP iter_factor = 0.99;
    DP iter_factor_extrap = 0.99;

    clock_t extrap_start;
    clock_t extrap_stop;
    clock_t Newton_start;
    clock_t Newton_stop;
    int iter_extrap_start, iter_extrap_stop, iter_Newton_start, iter_Newton_stop;
    DP diffsq_start;
    DP diffsq_extrap_start, diffsq_extrap_stop, diffsq_Newton_start, diffsq_Newton_stop;

    bool straight_improving = true;
    bool extrap_improving = true;
    bool Newton_improving = true;

    DP diffsq_iter_aim = 1.0e-5;

    bool info_findrap = false;

    if (info_findrap) cout << "Find_Rapidities called for state with label " << (*this).label << endl;

    while (diffsq > chain.prec && (iter < 500 && iter_Newton < 100
				   || straight_improving || extrap_improving || Newton_improving)) {

      // If we haven't reset, first try a few Newton steps...
      //if (!reset_rapidities && iter_Newton == 0) (*this).Solve_BAE_Newton (chain.prec, 10);

      if (!Newton_improving) diffsq_iter_aim *= 1.0e-5;

      // Start with a few straight iterations:

      if (!straight_improving) iter_factor *= 2.0/3.0;
      else iter_factor = 0.99;

      diffsq_start = diffsq;
      do {
	extrap_start = clock();
	iter_extrap_start = iter;
	diffsq_extrap_start = diffsq;

	if (diffsq > iter_prec) (*this).Solve_BAE_straight_iter (iter_prec, 10, iter_factor);

	extrap_stop = clock();
	iter_extrap_stop = iter;
	diffsq_extrap_stop = diffsq;

	iter_prec = diffsq * 0.1;

	if (info_findrap) cout << "Straight iter:  iter " << iter << "\titer_factor " << iter_factor << "\tdiffsq " << diffsq << endl;

	//} while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_extrap_stop/diffsq_extrap_start < 0.01);
      } while (diffsq > diffsq_iter_aim && !is_nan(diffsq) && diffsq_extrap_stop/diffsq_extrap_start < 0.01);

      straight_improving = (diffsq < diffsq_start);

      // Now try to extrapolate to infinite nr of iterations:

      if (!extrap_improving) iter_factor_extrap *= 2.0/3.0;
      else iter_factor_extrap = 0.99;

      diffsq_start = diffsq;

      if (diffsq > chain.prec)
	do {
	  extrap_start = clock();
	  iter_extrap_start = iter;
	  diffsq_extrap_start = diffsq;

	  if (diffsq > iter_prec) (*this).Solve_BAE_extrap (iter_prec, 10, iter_factor_extrap);

	  extrap_stop = clock();
	  iter_extrap_stop = iter;
	  diffsq_extrap_stop = diffsq;

	  iter_prec = diffsq * 0.1;

	  if (info_findrap) cout << "Extrap:  iter " << iter << "\tdiffsq " << diffsq << endl;

	} while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_extrap_stop/diffsq_extrap_start < 0.01);

      extrap_improving = (diffsq < diffsq_extrap_start);

      // Now try Newton

      diffsq_Newton_start = diffsq;
      if (diffsq > chain.prec)
	do {
	  Newton_start = clock();
	  iter_Newton_start = iter_Newton;
	  diffsq_Newton_start = diffsq;

	  if (diffsq > iter_prec) (*this).Solve_BAE_Newton (chain.prec, 5);

	  Newton_stop = clock();
	  iter_Newton_stop = iter_Newton;
	  diffsq_Newton_stop = diffsq;

	  if (info_findrap) cout << "Newton:  iter_Newton " << iter_Newton << "\tdiffsq " << diffsq << endl;

	} while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_Newton_stop/diffsq_Newton_start < 0.01);

      Newton_improving = (diffsq < diffsq_Newton_start);

      //if (diffsq > iter_prec) (*this).Solve_BAE_smackdown (0.1 * diffsq, 1);

      //cout << "Before silk:  diffsq = " << diffsq << endl;

      // If none of the methods are improving the result, try the silk gloves...
      if (!(straight_improving || extrap_improving || Newton_improving)) {
	if (info_findrap) cout << "Before silk gloves: diffsq " << diffsq << endl;
	(*this).Solve_BAE_with_silk_gloves (chain.prec, chain.Nsites, 0.9);
	if (info_findrap) cout << "After silk gloves: diffsq " << diffsq << endl;
      }

      if (is_nan(diffsq)) {
	(*this).Set_Free_lambdas();  // good start if we've messed up with Newton
	diffsq = 1.0;
      }

      iter_prec *= 1.0e-4;
      iter_prec = ABACUS::max(iter_prec, chain.prec);

    } // while (diffsq > chain.prec && (iter < 300 && iter_Newton < 50...

    // Check convergence:

    //cout << "Check_Rapidities: " << (*this).Check_Rapidities() << endl;

    //conv = (diffsq < chain.prec && (*this).Check_Rapidities() && ((*this).String_delta() < HEIS_deltaprec)) ? 1 : 0;
    conv = (diffsq < chain.prec && (*this).Check_Rapidities()) ? 1 : 0;
    dev = (*this).String_delta();

    //cout << "String delta: " << (*this).String_delta() << "\tBoolean: " << ((*this).String_delta() < HEIS_deltaprec) << endl;

    return;
  }

  void Heis_Bethe_State::Solve_BAE_bisect (int j, int alpha, DP req_prec, int itermax)
  {
    // Finds the root lambda[j][alpha] to precision req_prec using bisection

    DP prec_obtained = 1.0;
    int niter_here = 0;

    DP lambdajalphastart = lambda[j][alpha];
    DP lambdamax = (fabs(chain.Delta) <= 1.0 ? pow(log(chain.Nsites), 1.1) : PI - 1.0e-4);

    // Find values of lambda such that BE changes sign:
    DP lambdaleft = -lambdamax;
    DP lambdamid = 0.0;
    DP lambdaright = lambdamax;

    DP BEleft, BEmid, BEright;
    lambda[j][alpha] = lambdaleft;
    (*this).Compute_BE (j, alpha);
    //cout << "lambda: " << lambda[j][alpha] << "\ttanhlambda: " << tanh(lambda[j][alpha]) << endl;
    BEleft= BE[j][alpha];
    lambda[j][alpha] = lambdaright;
    (*this).Compute_BE (j, alpha);
    BEright= BE[j][alpha];

    if (BEleft * BEright > 0.0) {
      lambda[j][alpha] = lambdajalphastart;
      //cout << lambdaleft << "\t" << lambdaright << "\t" << BEleft << "\t" << BEright << endl;
      //cout << "Could not bisect BE[" << j << "][" << alpha << "]" << endl;
      return;
    }

    while (prec_obtained > req_prec && niter_here < itermax) {

      lambdamid = 0.5 * (lambdaleft + lambdaright);
      lambda[j][alpha] = lambdamid;
      (*this).Compute_BE (j, alpha);
      BEmid = BE[j][alpha];

      //cout << "niter_here = " << niter_here << "\t" << lambdaleft << "\t" << lambdamid << "\t" << lambdaright << endl;
      //cout << BEleft << "\t" << BEmid << "\t" << BEright << endl;

      if (BEmid * BEleft < 0.0) { // root is to the left of mid
	lambdaright = lambdamid;
	BEright = BEmid;
      }

      else if (BEmid * BEright < 0.0) { // root is to the right of mid
	lambdaleft = lambdamid;
	BEleft = BEmid;
      }

      else if (BEmid * BEmid < req_prec) return;

      else {
	//cout << lambdaleft << "\t" << lambdamid << "\t" << lambdaright << endl;
	//cout << BEleft << "\t" << BEmid << "\t" << BEright << endl;
	//ABACUSerror("Problem in Solve_BAE_bisect.");
	return; // this procedure has failed
      }

      prec_obtained = BEmid * BEmid;
      niter_here++;

      //cout << "bisect: " << lambdaleft << "\t" << lambdaright << "\t" << BEleft << "\t" << BEright << "\t" << prec_obtained << "\t" << req_prec << endl;
    }

    return;
  }

  void Heis_Bethe_State::Iterate_BAE (DP iter_factor)
  {
    //DP lambda_old;
    for (int j = 0; j < chain.Nstrings; ++j) {
      for (int alpha = 0; alpha < base[j]; ++alpha)
	{
	  //lambda_old = lambda[j][alpha];
	  //lambda[j][alpha] = Iterate_BAE (j, alpha);
	  lambda[j][alpha] += iter_factor * (Iterate_BAE (j, alpha) - lambda[j][alpha]);
	  //cout << j << "\t" << alpha << "\t" << Ix2[j][alpha] << "\t" << lambda_old << "\t" << lambda[j][alpha] << "\t";
	  //if (j > 0) cout << j << "\t" << alpha << "\t" << Ix2[j][alpha] << "\t" << lambda_old << "\t" << lambda[j][alpha] << endl;
	}
    }

    iter++;

    (*this).Compute_BE();
    (*this).Compute_diffsq();
  }

  void Heis_Bethe_State::Solve_BAE_straight_iter (DP straight_prec, int max_iter, DP iter_factor)
  {
    // This function attempts to get convergence diffsq <= straight_prec in at most max_iter steps.

    // If this fails, the lambda's are reset to original values, defined as...

    Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    diffsq = 1.0;
    int iter_done_here = 0;

    while ((diffsq > straight_prec) && (max_iter > iter_done_here)) {

      //cout << "BEFORE ITERATION" << endl << (*this) << endl << endl;
      (*this).Iterate_BAE(iter_factor);
      //cout << "ITERATION " << iter_done_here << endl << (*this) << endl << endl;
      iter_done_here++;
    }

    if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
      // This procedure has failed.  We reset everything to begin values.
      //cout << "Straight iter failed: resetting." << endl;
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      (*this).Compute_BE();
      diffsq = diffsq_ref;
    }
  }

  void Heis_Bethe_State::Solve_BAE_extrap (DP extrap_prec, int max_iter_extrap, DP iter_factor)
  {
    // This function attempts to get convergence diffsq <= extrap_prec in at most max_iter_extrap steps.

    // If this fails, the lambda's are reset to original values, defined as...

    Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    diffsq = 1.0;
    DP diffsq1, diffsq2, diffsq3, diffsq4;

    int iter_done_here = 0;

    Lambda lambda1(chain, base);
    Lambda lambda2(chain, base);
    Lambda lambda3(chain, base);
    Lambda lambda4(chain, base);
    Lambda lambda5(chain, base);

    while ((diffsq > extrap_prec) && (max_iter_extrap > iter_done_here)) {

      (*this).Iterate_BAE(iter_factor);
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda1[j][alpha] = lambda[j][alpha];
      diffsq1 = diffsq;
      (*this).Iterate_BAE(iter_factor);
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda2[j][alpha] = lambda[j][alpha];
      diffsq2 = diffsq;
      (*this).Iterate_BAE(iter_factor);
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda3[j][alpha] = lambda[j][alpha];
      diffsq3 = diffsq;
      (*this).Iterate_BAE(iter_factor);
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda4[j][alpha] = lambda[j][alpha];
      diffsq4 = diffsq;
      //(*this).Iterate_BAE(iter_factor);
      //for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda5[j][alpha] = lambda[j][alpha];
      //diffsq5 = diffsq;

      iter_done_here += 4;

      // now extrapolate the result to infinite number of iterations...

      if (diffsq4 < diffsq3 && diffsq3 < diffsq2 && diffsq2 < diffsq1 && diffsq1 < diffsq_ref) {

	Vect_DP rap(0.0, 4);
	Vect_DP oneoverP(0.0, 4);
	DP deltalambda = 0.0;

	for (int i = 0; i < 4; ++i) oneoverP[i] = 1.0/(1.0 + i*i);

	for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) {
	    rap[0] = lambda1[j][alpha];
	    rap[1] = lambda2[j][alpha];
	    rap[2] = lambda3[j][alpha];
	    rap[3] = lambda4[j][alpha];
	    //rap[4] = lambda5[j][alpha];

	    polint (oneoverP, rap, 0.0, lambda[j][alpha], deltalambda);

	    //cout << j << "\t" << alpha << "\t" << rap << "\t" << lambda[j][alpha] << "\t" << deltalambda << endl;
	  }

	// Iterate once to stabilize result
	(*this).Iterate_BAE(iter_factor);
      }

    } // ((diffsq > extrap_prec) && (max_iter_extrap > iter_done_here))

    if ((diffsq >= diffsq_ref) || (is_nan(diffsq))) {
      // This procedure has failed.  We reset everything to begin values.
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      (*this).Compute_BE();
      diffsq = diffsq_ref;
    }

    return;
  }

  void Heis_Bethe_State::Solve_BAE_with_silk_gloves (DP silk_prec, int max_iter_silk, DP iter_factor)
  {
    // Logic: do iterations, but change only one rapidity at a time.
    // This is slow, and called only if straight_iter, extrap and Newton methods don't work.

    // If this fails, the lambda's are reset to original values, defined as...

    Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    diffsq = 1.0;
    int iter_done_here = 0;

    while ((diffsq > silk_prec) && (max_iter_silk > iter_done_here)) {

      // Find the highest `deviant' rapidity
      int jmax = 0;
      int alphamax = 0;
      DP absBEmax = 0.0;
      for (int j = 0; j < chain.Nstrings; ++j)
	for (int alpha = 0; alpha < base[j]; ++alpha) {
	  //cout << j << "\t" << alpha << "\t" << BE[j][alpha] << endl;
	  if (fabs(BE[j][alpha]) > absBEmax) {
	    jmax = j;
	    alphamax = alpha;
	    absBEmax = fabs(BE[j][alpha]);
	  }
	}

      //cout << "jmax = " << jmax << "\talphamax = " << alphamax << "\t" << lambda[jmax][alphamax] << "\tBE before " << BE[jmax][alphamax] << endl;
      // Now recalculate this max deviant rapidity,
      //cout <<  lambda[jmax][alphamax] << "\t" << Iterate_BAE (jmax, alphamax) << endl;

      /*
      DP dlambda = 0.0;
      DP prevBEmax = 0.0;
      do {
	dlambda = Iterate_BAE (jmax, alphamax) - lambda[jmax][alphamax];
	lambda[jmax][alphamax] += iter_factor * dlambda;
	prevBEmax = BE[jmax][alphamax];
	(*this).Compute_BE();
	iter_done_here++;
	cout << "jmax = " << jmax << "\talphamax = " << alphamax << "\t" << lambda[jmax][alphamax] << "\t" << dlambda << "\tBE during " << BE[jmax][alphamax] << endl;
      } while (dlambda * dlambda > silk_prec && fabs(BE[jmax][alphamax]) < fabs(prevBEmax) && max_iter_silk > iter_done_here);

      */

      Solve_BAE_bisect (jmax, alphamax, silk_prec, max_iter_silk);
      iter_done_here++;

      //cout << "jmax = " << jmax << "\talphamax = " << alphamax << "\t" << lambda[jmax][alphamax] << "\tBE after " << BE[jmax][alphamax] << endl;

      // and reset all important arrays.
      (*this).Compute_diffsq();
    }

    //cout << "Silk gloves: diffsq from " << diffsq_ref << "\tto " << diffsq << endl;

    if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
      // This procedure has failed.  We reset everything to begin values.
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      (*this).Compute_BE();
      diffsq = diffsq_ref;
    }

    return;
  }

  /*
  void Heis_Bethe_State::BAE_smackdown (DP max_allowed)
  {
    // Re-solves for all rapidities lambda[j][alpha] such that BE[j][alpha]^2/N > max_allowed.
    // Assumes that BE[][] is up-to-date.

    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha)
	if (pow(BE[j][alpha], 2.0)/chain.Nsites > max_allowed) (*this).Solve_BAE (j, alpha, max_allowed, 100);
  }

  void Heis_Bethe_State::Solve_BAE_smackdown (DP max_allowed, int maxruns)
  {
    int runs_done = 0;

    (*this).Compute_BE();
    (*this).Compute_diffsq();

    while (diffsq > chain.prec && diffsq > max_allowed && runs_done < maxruns) {
      (*this).BAE_smackdown (max_allowed);
      (*this).Compute_BE();
      (*this).Compute_diffsq();
      runs_done++;
    }
  }
  */

  void Heis_Bethe_State::Iterate_BAE_Newton ()
  {
    // does one step of a Newton method on the rapidities...
    // Assumes that BE[j][alpha] have been computed

    Vect_CX dlambda (0.0, base.Nraptot);  // contains delta lambda computed from Newton's method
    SQMat_CX Gaudin (0.0, base.Nraptot);
    Vect_INT indx (base.Nraptot);

    Lambda lambda_ref(chain, base);
    for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];

    (*this).Build_Reduced_Gaudin_Matrix (Gaudin);

    int index = 0;
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	dlambda[index] = - BE[j][alpha] * chain.Nsites;
	index++;
      }

    DP d;

    try {
      ludcmp_CX (Gaudin, indx, d);
      lubksb_CX (Gaudin, indx, dlambda);
    }

    catch (Divide_by_zero) {
      diffsq = log(-1.0);  // reset to nan to stop Newton method
      return;
    }

    //cout << iter_Newton << "\t" << dlambda << endl;

    // Regularize dlambda:  max step is +-1.0 to prevent rapidity flying off into outer space.
    for (int i = 0; i < base.Nraptot; ++i) if (fabs(real(dlambda[i])) > 1.0) dlambda[i] = 0.0;//(real(dlambda[i]) > 0) ? 1.0 : -1.0;

    index = 0;
    for (int j = 0; j < chain.Nstrings; ++j)
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	lambda[j][alpha] = lambda_ref[j][alpha] + real(dlambda[index]);
	//cout << j << "\t" << alpha << "\t" << dlambda[index] << "\t" << lambda_ref[j][alpha] << "\t" << lambda[j][alpha] << endl;
	index++;
      }

    (*this).Compute_BE();
    (*this).Iterate_BAE(1.0);
    (*this).Compute_diffsq();

    // if we've converged, calculate the norm here, since the work has been done...

    if (diffsq < chain.prec) {
      lnnorm = 0.0;
      for (int j = 0; j < base.Nraptot; j++) lnnorm += log(abs(Gaudin[j][j]));
    }

    iter_Newton++;

    return;
  }

  void Heis_Bethe_State::Solve_BAE_Newton (DP Newton_prec, int max_iter_Newton)
  {
    // This function attempts to get convergence diffsq <= Newton_prec in at most max_iter_Newton steps.

    // The results are accepted if diffsq has decreased, otherwise the lambda's are reset to original values, defined as...

    Lambda lambda_ref(chain, base);
    DP diffsq_ref = 1.0;

    for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
    diffsq_ref = diffsq;

    // Now begin solving...

    int iter_done_here = 0;

    while ((diffsq > Newton_prec) && (diffsq < 10.0) && (!is_nan(diffsq)) && (iter_done_here < max_iter_Newton)) {
      (*this).Iterate_BAE_Newton();
      iter_done_here++;
      (*this).Iterate_BAE(1.0);
    }

    if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
      // This procedure has failed.  We reset everything to begin values.
      //cout << "Newton: failed, resetting." << "\t" << diffsq << endl;
      for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
      (*this).Compute_BE();
      diffsq = diffsq_ref;
    }

    return;
  }

  void Heis_Bethe_State::Compute_lnnorm ()
  {
    if (true || lnnorm == -100.0) {  // else norm already calculated by Newton method
      // Actually, compute anyway to increase accuracy !

      SQMat_CX Gaudin_Red(base.Nraptot);

      (*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);
      /*
      cout << endl << "Gaudin matrix: " << endl;

      for (int j = 0; j < Gaudin_Red.size(); ++j) {
	for (int k = 0; k < Gaudin_Red.size(); ++k) cout << Gaudin_Red[j][k] << "  ";
	cout << endl;
      }
      cout << endl << endl;
      */
      complex<DP> lnnorm_check = lndet_LU_CX_dstry(Gaudin_Red);
      //cout << "Calculated lnnorm = " << lnnorm_check;

      //lnnorm = real(lndet_LU_CX_dstry(Gaudin_Red));
      lnnorm = real(lnnorm_check);
    }

    return;
  }

  void Heis_Bethe_State::Compute_Momentum ()
  {
    int sum_Ix2 = 0;
    DP sum_M = 0.0;

    for (int j = 0; j < chain.Nstrings; ++j) {
      sum_M += 0.5 * (1.0 + chain.par[j]) * base[j];
      for (int alpha = 0; alpha < base[j]; ++alpha) {
	sum_Ix2 += Ix2[j][alpha];
      }
    }

    iK = (chain.Nsites/2) * int(sum_M + 0.1) - (sum_Ix2/2);  // + 0.1:  for safety...

    while (iK >= chain.Nsites) iK -= chain.Nsites;
    while (iK < 0) iK += chain.Nsites;

    K = PI * sum_M - PI * sum_Ix2/chain.Nsites;

    while (K >= 2.0*PI) K -= 2.0*PI;
    while (K < 0.0) K += 2.0*PI;

    return;
  }

  void Heis_Bethe_State::Compute_All (bool reset_rapidities)     // solves BAE, computes E, K and lnnorm
  {
    (*this).Find_Rapidities (reset_rapidities);
    if (conv == 1) {
      (*this).Compute_Energy ();
      (*this).Compute_Momentum ();
      (*this).Compute_lnnorm ();
    }
    return;
  }

  /*
  bool Heis_Bethe_State::Boost_Momentum (int iKboost)
  {
    if (iKboost == 0) return(true);  // done

    Ix2_Offsets offsets_here = offsets;

    bool success = false;

    if (iKboost < 0)
      success = offsets_here.Add_Boxes_From_Lowest(-iKboost, 0);  // add boxes in even sectors to decrease iK

    else if (iKboost > 0)
      success = offsets_here.Add_Boxes_From_Lowest(iKboost, 1);  // add boxes in odd sectors to increase iK

    if (success) (*this).Set_Ix2_Offsets(offsets_here);

    return(success);
  }
  */


  void Heis_Bethe_State::Set_to_Closest_Matching_Ix2_fixed_Base (const Heis_Bethe_State& StateToMatch)
  {
    // Given a state with given Ix2 distribution, set the Ix2 to closest match.
    // The base of (*this) is fixed, and does not necessarily match that of StateToMatch.

    //cout << "Matching Ix2 for base " << (*this).base.baselabel << " from base " << StateToMatch.base.baselabel << endl;

    if ((*this).chain != StateToMatch.chain)
      ABACUSerror("Heis_Bethe_State::Find_Closest_Matching_Ix2_fixed_Base: trying to match Ix2 for two states with different chains.");

    // Check level by level, match quantum numbers from center up.
    for (int il = 0; il < chain.Nstrings; ++il) {
      // Careful: if parity of Nraps is different, the Ix2 live on different lattices!
      int Ix2shift = ((StateToMatch.base.Nrap[il] - (*this).base.Nrap[il]) % 2);
      if (StateToMatch.base.Nrap[il] >= (*this).base.Nrap[il]) {
	// We match the rapidities from the center towards the sides.
	int ashift = (StateToMatch.base.Nrap[il] - (*this).base.Nrap[il])/2;
	for (int a = 0; a < (*this).base.Nrap[il]; ++a)
	  (*this).Ix2[il][a] = StateToMatch.Ix2[il][a + ashift] + Ix2shift;
      }
      else { // There are less Ix2 in StateToMatch than in (*this)
	// We start by filling all the (*this) Ix2 symmetrically from the middle.
	// We thereafter start from the left and identify half of the StateToMatch Ix2 with the leftmost (*this)Ix2
	// and then do the same thing from the right.
	for (int a = 0; a < (*this).base.Nrap[il]; ++a)
	  (*this).Ix2[il][a] = -(*this).base.Nrap[il] + 1 + 2*a;
	int nleft = StateToMatch.base.Nrap[il]/2;
	for (int a = 0; a < nleft; ++a)
	  if (StateToMatch.Ix2[il][a] - Ix2shift < (*this).Ix2[il][a]) (*this).Ix2[il][a] = StateToMatch.Ix2[il][a] - Ix2shift;
	for (int a = 0; a < StateToMatch.base.Nrap[il] - 1 - nleft; ++a)
	  if (StateToMatch.Ix2[il][StateToMatch.base.Nrap[il] - 1 - a] - Ix2shift > (*this).Ix2[il][(*this).base.Nrap[il] - 1 - a])
	    (*this).Ix2[il][(*this).base.Nrap[il] - 1 - a] = StateToMatch.Ix2[il][StateToMatch.base.Nrap[il] - 1 - a] - Ix2shift;
      }
    } // for il

    //cout << "StateToMatch:" << endl << StateToMatch << endl << "MatchingState:" << endl << (*this) << endl;

  }


  std::ostream& operator<< (std::ostream& s, const Heis_Bethe_State& state)
  {
    // sends all the state data to output stream

    s << endl << "******** Chain with Delta = " << state.chain.Delta << " Nsites = " << state.chain.Nsites << " Mdown = " << state.base.Mdown
      //<< ":  eigenstate with base_id " << state.base_id << ", type_id " << state.type_id << " id " << state.id << " maxid " << state.maxid << endl
      << ":  eigenstate with label " << state.label << endl
      << "E = " << state.E << " K = " << state.K << " iK = " << state.iK << " lnnorm = " << state.lnnorm << endl
      << "conv = " << state.conv << " dev = " << state.dev << " iter = " << state.iter << " iter_Newton = " << state.iter_Newton << "\tdiffsq " << state.diffsq << endl;

    for (int j = 0; j < state.chain.Nstrings; ++j) {
      if (state.base.Nrap[j] > 0) {
	s << "Type " << j << " Str_L = " << state.chain.Str_L[j] << " par = " << state.chain.par[j] << " M_j = " << state.base.Nrap[j]
	  << " Ix2_infty = " << state.base.Ix2_infty[j] << " Ix2_min = " << state.base.Ix2_min[j] << " Ix2_max = " << state.base.Ix2_max[j] << endl;
	Vect_INT qnumbers(state.base.Nrap[j]);
	Vect_DP rapidities(state.base.Nrap[j]);
	for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) {
	  qnumbers[alpha] = state.Ix2[j][alpha];
	  rapidities[alpha] = state.lambda[j][alpha];
	}
	qnumbers.QuickSort();
	rapidities.QuickSort();

	s << "Ix2 quantum numbers: " << endl;
	for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << qnumbers[alpha] << " ";
	s << endl;
	s << "Rapidities: " << endl;
	for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << rapidities[alpha] << " ";
	s << endl;
      }
    }
    s << endl;

    return s;

  }


} // namespace ABACUS