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;

      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;

	l++;

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

      }

      // 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;
  }


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

  // 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())
      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];

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

    //State_Label_Data labeldata = Read_State_Label (label_ref, dummyOriginIx2);
    State_Label_Data labeldata = Read_Base_Label (label_ref);

    // 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));

	  sum2 = Theta_XXZ(1.0, RefChain.Str_L[j], RefChain.Str_L[k], RefChain.par[j],
			   RefChain.par[k], RefChain.ta_n_anis_over_2);

	  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;

	// 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 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 - 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;
	}

	Ix2_min[j] = -Ix2_max[j];

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

      // New attempt
      // for (int j = 0; j < RefChain.Nstrings; ++j) {
      // 	sum1 = 0;
      // 	for (int k = 0; k < RefChain.Nstrings; ++k) {
      // 	  sum1 += (j == k ? 0 : Nrap[k] * 2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]));
      // 	}
      // 	Ix2_max[j] = RefChain.Nsites - (2 * RefChain.Str_L[j] - 1) * (Nrap[j] - 1) - sum1;
      // 	Ix2_min[j] = -Ix2_max[j] + 2;
      // Ix2_max[j] += 2;
      // Ix2_min[j] -= 2;
      //}

    } // 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[0] = new 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[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) {
      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()), 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), label("")
  {
  };

  Heis_Bethe_State::Heis_Bethe_State (const Heis_Bethe_State& RefState) // copy constructor
    : chain(RefState.chain), base(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),
      label(RefState.label)
  {
    // copy arrays into new ones
    for (int j = 0; j < RefState.chain.Nstrings; ++j) {
      for (int alpha = 0; alpha < RefState.base[j]; ++j) {
	lambda[j][alpha] = RefState.lambda[j][alpha];
      }
    }
  }

  Heis_Bethe_State::Heis_Bethe_State (const Heis_Chain& RefChain, int M)
    : chain(RefChain), base(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)
  {
    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), 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)
  {
    // Check that the number of rapidities is consistent with Mdown
    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)
  {
    State_Label_Data labeldata = Read_State_Label (label_ref, OriginIx2);

    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]) {
	    Ix2[labeldata.type[it] ][i] = labeldata.Ix2exc[it][iexc];
	  }

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

  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 = -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;

    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;
    }

    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

    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();
    (*this).Compute_BE();

    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 (!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 > 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 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:

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

    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);
    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;
      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];

      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 {
	return; // this procedure has failed
      }

      prec_obtained = BEmid * BEmid;
      niter_here++;

    }

    return;
  }

  void Heis_Bethe_State::Iterate_BAE (DP iter_factor)
  {
    for (int j = 0; j < chain.Nstrings; ++j) {
      for (int alpha = 0; alpha < base[j]; ++alpha)
	{
	  lambda[j][alpha] += iter_factor * (Iterate_BAE (j, alpha) - lambda[j][alpha]);
	}
    }

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

      (*this).Iterate_BAE(iter_factor);
      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;
    }
  }

  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;

      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];

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

	// 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]);
	  }
	}

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

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

    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::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;
    }

    // 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]);
	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.
      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);
      complex<DP> lnnorm_check = 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;
  }


  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.

    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

  }


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