ABACUS/src/LIEBLIN/LiebLin_Bethe_State.cc

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

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

Copyright (c) J.-S. Caux.

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

File:  LiebLin_Bethe_State.cc

Purpose:  Definitions for LiebLin_Bethe_State class.

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

#include "ABACUS.h"

using namespace std;

namespace ABACUS {

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

  // Function definitions:  class LiebLin_Bethe_State

  LiebLin_Bethe_State::LiebLin_Bethe_State ()
    : c_int (0.0), L(0.0), cxL(0.0), N(0),
      Ix2_available(Vect<int>(0, 1)), index_first_hole_to_right (Vect<int>(0,1)),
      displacement (Vect<int>(0,1)), Ix2(Vect<int>(0, 1)), lambdaoc(Vect<DP>(0.0, 1)),
      S(Vect<DP>(0.0, 1)), dSdlambdaoc(Vect<DP>(0.0, 1)), diffsq(0.0), prec(ITER_REQ_PREC_LIEBLIN),
      conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0)
  {
    stringstream Nout; Nout << N; label = Nout.str() + LABELSEP + ABACUScoding[0] + LABELSEP;
  }

  LiebLin_Bethe_State::LiebLin_Bethe_State (DP c_int_ref, DP L_ref, int N_ref)
    : c_int(c_int_ref), L(L_ref), cxL(c_int_ref * L_ref), N(N_ref),
      Ix2_available(Vect<int>(0, 2)), index_first_hole_to_right (Vect<int>(0,N)),
      displacement (Vect<int>(0,N)), Ix2(Vect<int>(0, N)), lambdaoc(Vect<DP>(0.0, N)),
      S(Vect<DP>(0.0, N)), dSdlambdaoc(Vect<DP>(0.0, N)),
      diffsq(0.0), prec(ABACUS::max(1.0, 1.0/(c_int * c_int)) * ITER_REQ_PREC_LIEBLIN),
      conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0)
  {
    if (c_int < 0.0) ABACUSerror("You must use a positive interaction parameter !");
    if (N < 0) ABACUSerror("Particle number must be strictly positive.");

    stringstream Nout; Nout << N; label = Nout.str() + LABELSEP + ABACUScoding[0] + LABELSEP;

    // Set quantum numbers to ground-state configuration:
    for (int i = 0; i < N; ++i) Ix2[i] = -(N-1) + 2*i;

    Vect<int> OriginIx2 = Ix2;

    (*this).Set_Label_from_Ix2 (OriginIx2);
  }


  LiebLin_Bethe_State& LiebLin_Bethe_State::operator= (const LiebLin_Bethe_State& RefState)
  {
    if (this != &RefState) {
      c_int = RefState.c_int;
      L = RefState.L;
      cxL = RefState.cxL;
      N = RefState.N;
      label = RefState.label;
      Ix2_available = RefState.Ix2_available;
      index_first_hole_to_right = RefState.index_first_hole_to_right;
      displacement = RefState.displacement;
      Ix2 = RefState.Ix2;
      lambdaoc = RefState.lambdaoc;
      S = RefState.S;
      dSdlambdaoc = RefState.dSdlambdaoc;
      diffsq = RefState.diffsq;
      prec = RefState.prec;
      conv = RefState.conv;
      iter_Newton = RefState.iter_Newton;
      E = RefState.E;
      iK = RefState.iK;
      K = RefState.K;
      lnnorm = RefState.lnnorm;
    }
    return(*this);
  }

  void LiebLin_Bethe_State::Set_to_Label (string label_ref, const Vect<int>& OriginStateIx2)
  {
    State_Label_Data labeldata = Read_State_Label (label_ref, OriginStateIx2);

    if (N != labeldata.M[0]) {
      cout << label_ref << endl;
      cout << labeldata.M << endl;
      ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: N != M[0].");
    }
    if (N != OriginStateIx2.size()) {
      cout << label_ref << endl;
      cout << labeldata.M << endl;
      ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: "
		  "N != OriginStateIx2.size().");
    }

    label = label_ref;

    Vect<int> OriginStateIx2ordered = OriginStateIx2;
    OriginStateIx2ordered.QuickSort();

    // Set all Ix2 to OriginState's Ix2:
    for (int i = 0; i < N; ++i) Ix2[i] = OriginStateIx2ordered[i];

    // Now set the excitations:
    for (int iexc = 0; iexc < labeldata.nexc[0]; ++iexc)
      for (int i = 0; i < N; ++i)
	if (Ix2[i] == labeldata.Ix2old[0][iexc]) Ix2[i] = labeldata.Ix2exc[0][iexc];

    // Now reorder the Ix2 to follow convention:
    Ix2.QuickSort();

    (*this).Set_Label_from_Ix2 (OriginStateIx2ordered);
  }

  void LiebLin_Bethe_State::Set_to_Label (string label_ref)
  {
    // This function assumes that OriginState is the ground state.

    Vect<int> OriginStateIx2(N);
    for (int i = 0; i < N; ++i) OriginStateIx2[i] = -N + 1 + 2*i;

    (*this).Set_to_Label(label_ref, OriginStateIx2);
  }


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

    if (N != OriginStateIx2.size()) ABACUSerror("N != OriginStateIx2.size() in Set_Label_from_Ix2.");


    // Set the state label:
    Vect<int> type_ref(0,1);
    Vect<int> M_ref(N, 1);
    Vect<int> nexc_ref(0, 1);
    // Count nr of particle-holes:
    for (int i = 0; i < N; ++i) if (!OriginStateIx2.includes(Ix2[i])) nexc_ref[0] += 1;
    Vect<Vect<int> > Ix2old_ref(1);
    Vect<Vect<int> > Ix2exc_ref(1);
    Ix2old_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
    Ix2exc_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
    int nexccheck = 0;
    for (int i = 0; i < N; ++i)
      if (!OriginStateIx2.includes(Ix2[i])) Ix2exc_ref[0][nexccheck++] = Ix2[i];
    if (nexccheck != nexc_ref[0])
      ABACUSerror("Counting excitations wrong (1) in LiebLin_Bethe_State::Set_Label_from_Ix2");
    nexccheck = 0;
    for (int i = 0; i < N; ++i)
      if (!Ix2.includes (OriginStateIx2[i])) Ix2old_ref[0][nexccheck++] = OriginStateIx2[i];
    if (nexccheck != nexc_ref[0]) {
      cout << "nexc_ref[0] = " << nexc_ref[0] << "\tnexccheck = " << nexccheck << endl;
      cout << OriginStateIx2 << endl;
      cout << Ix2 << endl;
      cout << nexc_ref[0] << endl;
      cout << Ix2old_ref[0] << endl;
      cout << Ix2exc_ref[0] << endl;
      ABACUSerror("Counting excitations wrong (2) in LiebLin_Bethe_State::Set_Label_from_Ix2");
    }
    // Now order the Ix2old_ref and Ix2exc_ref:
    Ix2old_ref[0].QuickSort();
    Ix2exc_ref[0].QuickSort();

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

    label = Return_State_Label (labeldata, OriginStateIx2);
  }

  void LiebLin_Bethe_State::Set_Label_Internals_from_Ix2 (const Vect<int>& OriginStateIx2)
  {
    if (N != OriginStateIx2.size())
      ABACUSerror("N != OriginStateIx2.size() in Set_Label_Internals_from_Ix2.");

    Vect<int> OriginStateIx2ordered = OriginStateIx2;
    OriginStateIx2ordered.QuickSort();

    // Set the state label:
    Vect<int> type_ref(0,1);
    Vect<int> M_ref(N, 1);
    Vect<int> nexc_ref(0, 1);
    // Count nr of particle-holes:
    for (int i = 0; i < N; ++i) if (!OriginStateIx2.includes(Ix2[i])) nexc_ref[0] += 1;
    Vect<Vect<int> > Ix2old_ref(1);
    Vect<Vect<int> > Ix2exc_ref(1);
    Ix2old_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
    Ix2exc_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
    int nexccheck = 0;
    for (int i = 0; i < N; ++i)
      if (Ix2[i] != OriginStateIx2ordered[i]) {
	Ix2old_ref[0][nexccheck] = OriginStateIx2ordered[i];
	Ix2exc_ref[0][nexccheck++] = Ix2[i];
      }

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

    label = Return_State_Label (labeldata, OriginStateIx2);

    // Construct the Ix2_available vector: we give one more quantum number on left and right:
    int navailable = 2 + (ABACUS::max(Ix2.max(), OriginStateIx2.max())
			  - ABACUS::min(Ix2.min(), OriginStateIx2.min()))/2 - N + 1;
    Ix2_available = Vect<int>(navailable);
    index_first_hole_to_right = Vect<int>(N);

    // First set Ix2_available to all holes from left
    for (int i = 0; i < Ix2_available.size(); ++i)
      Ix2_available[i] = ABACUS::min(Ix2.min(), OriginStateIx2.min()) - 2 + 2*i;

    // Now shift according to Ix2 of OriginState:
    for (int j = 0; j < N; ++j) {
      int i = 0;
      while (Ix2_available[i] < OriginStateIx2ordered[j]) i++;
      // We now have Ix2_available[i] == OriginStateIx2[j].
      // Shift all Ix2_available to the right of this by 2;
      for (int i1 = i; i1 < navailable; ++i1) Ix2_available[i1] += 2;
      index_first_hole_to_right[j] = i;
    }
    // Ix2_available and index_first_hole_to_right are now fully defined.

    // Now set displacement vector:
    displacement = Vect<int>(0, N);
    // Set displacement vector from the Ix2:
    for (int j = 0; j < N; ++j) {
      if (Ix2[j] < OriginStateIx2ordered[j]) {
	// Ix2[j] must be equal to some OriginState_Ix2_available[i]
	// for i < OriginState_index_first_hole_to_right[j]
	while (Ix2[j] != Ix2_available[index_first_hole_to_right[j] + displacement[j] ]) {
	  if (index_first_hole_to_right[j] + displacement[j] == 0) {
	    cout << label << endl << j << endl << OriginStateIx2 << endl
		 << Ix2 << endl << Ix2_available
		 << endl << index_first_hole_to_right << endl << displacement << endl;
	    ABACUSerror("Going down too far in Set_Label_Internals...");
	  }
	  displacement[j]--;
	}
      }
      if (Ix2[j] > OriginStateIx2ordered[j]) {
	// Ix2[j] must be equal to some Ix2_available[i] for i >= index_first_hole_to_right[j]
	displacement[j] = 1; // start with this value to prevent segfault
	while (Ix2[j] != Ix2_available[index_first_hole_to_right[j] - 1 + displacement[j] ]) {
	  if (index_first_hole_to_right[j] + displacement[j] == Ix2_available.size() - 1) {
	    cout << label << endl << j << endl << OriginStateIx2 << endl
		 << Ix2 << endl << Ix2_available
		 << endl << index_first_hole_to_right << endl << displacement << endl;
	    ABACUSerror("Going up too far in Set_Label_Internals...");
	  }
	  displacement[j]++;
	}
      }
    }
  }

  bool LiebLin_Bethe_State::Check_Admissibility (char whichDSF)
  {
    // For testing with restricted Hilbert space:
    //if (Ix2.min() < -13 || Ix2.max() > 13) return(false);
    return(true);
  }

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

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

    diffsq = 1.0;

    if (reset_rapidities) (*this).Set_Free_lambdaocs();

    iter_Newton = 0;

    DP damping = 1.0;
    DP diffsq_prev = 1.0e+6;

    Vect_DP RHSBAE (0.0, N);  // contains RHS of BAEs
    Vect_DP dlambdaoc (0.0, N);  // contains delta lambdaoc computed from Newton's method
    SQMat_DP Gaudin (0.0, N);
    Vect_INT indx (N);

    while (diffsq > prec && !is_nan(diffsq) && iter_Newton < 100) {
      (*this).Iterate_BAE_Newton(damping, RHSBAE, dlambdaoc, Gaudin, indx);
      if (diffsq > diffsq_prev && damping > 0.5) damping /= 2.0;
      else if (diffsq < diffsq_prev) damping = 1.0;
      diffsq_prev = diffsq;
    }

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

    if (!conv) {
      cout << "Alert! State " << label << " did not converge... diffsq " << diffsq
	   << "\titer_Newton " << iter_Newton << (*this) << endl;
    }

    return;
  }

  bool LiebLin_Bethe_State::Check_Rapidities()
  {
    bool nonan = true;

    for (int j = 0; j < N; ++j) nonan *= !is_nan(lambdaoc[j]);

    return nonan;
  }

  /**
     For the repulsive Lieb-Liniger model, all eigenstates have real rapidities
     (there are no strings). All string deviations are thus set to zero.
   */
  DP LiebLin_Bethe_State::String_delta()
  {
    return(0.0);
  }

  /**
     Checks whether the quantum numbers are symmetrically distributed: \f$\{ I \} = \{ -I \}\f$.
  */
  bool LiebLin_Bethe_State::Check_Symmetry ()
  {
    bool symmetric_state = true;

    Vect<int> Ix2check = Ix2;
    Ix2check.QuickSort();

    for (int alpha = 0; alpha <= N/2; ++alpha)
      symmetric_state = symmetric_state && (Ix2check[alpha] == -Ix2check[N - 1 - alpha]);

    return(symmetric_state);
  }

  /**
     This function computes the log of the norm of a `LiebLin_Bethe_State` instance.
     This is obtained from the Gaudin-Korepin norm formula
     \f{eqnarray*}{
     \langle 0 | \prod_{j=0}^{N-1} C(\lambda_j) \prod_{j=0}^{N-1} B(\lambda_j) | 0 \rangle
     = \prod_{j=0}^{N-2}\prod_{k=j+1}^{N-1}
     \frac{(\lambda_j - \lambda_k )^2 + c^2}{(\lambda_j - \lambda_k)^2}
     \times ~\mbox{Det}_N G
     \f}
     in which \f$G\f$ is the Gaudin matrix computed by
     LiebLin_Bethe_State::Build_Reduced_Gaudin_Matrix.
   */
  void LiebLin_Bethe_State::Compute_lnnorm ()
  {
    if (lnnorm == -100.0) {  // else Gaudin part already calculated by Newton method

      SQMat_DP Gaudin_Red(N);

      (*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);

      complex<DP> lnnorm_CX = lndet_LU_dstry(Gaudin_Red);
      if (abs(imag(lnnorm_CX)) > 1.0e-6) {
	cout << "lnnorm_CX = " << lnnorm_CX << endl;
	ABACUSerror("Gaudin norm of LiebLin_Bethe_State should be real and positive.");
      }

      lnnorm = real(lnnorm_CX);

      // Add the pieces outside of Gaudin determinant

      for (int j = 0; j < N - 1; ++j) for (int k = j+1; k < N; ++k)
					lnnorm += log(1.0 + 1.0/pow(lambdaoc[j] - lambdaoc[k], 2.0));

    }

    return;
  }

  /**
     This function solves the Bethe equations, computes the energy, momentum and state norm.
     It calls LiebLin_Bethe_State::Find_Rapidities, LiebLin_Bethe_State::Compute_Energy,
     LiebLin_Bethe_State::Compute_Momentum and LiebLin_Bethe_State::Compute_lnnorm.

     Assumptions:
     - the set of quantum numbers is admissible.
  */
  void LiebLin_Bethe_State::Compute_All (bool reset_rapidities)
  {
    (*this).Find_Rapidities (reset_rapidities);
    if (conv == 1) {
      (*this).Compute_Energy ();
      (*this).Compute_Momentum ();
      (*this).Compute_lnnorm ();
    }
    return;
  }

  void LiebLin_Bethe_State::Set_Free_lambdaocs()
  {
    if (cxL >= 1.0)
      for (int a = 0; a < N; ++a) lambdaoc[a] = PI * Ix2[a]/cxL;

    // For small values of c, use better approximation using approximate
    // zeroes of Hermite polynomials: see Gaudin eqn 4.71.
    if (cxL < 1.0) {
      DP oneoversqrtcLN = 1.0/pow(cxL * N, 0.5);
      for (int a = 0; a < N; ++a) lambdaoc[a] = oneoversqrtcLN * PI * Ix2[a];
    }

    return;
  }

  void LiebLin_Bethe_State::Iterate_BAE (DP damping)
  {
    // does one step of simple iterations

    DP sumtheta = 0.0;
    Vect_DP dlambdaoc (0.0, N);

    for (int j = 0; j < N; ++j) {

      sumtheta = 0.0;
      for (int k = 0; k < N; ++k) sumtheta += atan((lambdaoc[j] - lambdaoc[k]));
      sumtheta *= 2.0;

      dlambdaoc[j] = damping * ((PI*Ix2[j] - sumtheta)/cxL - lambdaoc[j]);

    }

    diffsq = 0.0;
    for (int i = 0; i < N; ++i) {
      lambdaoc[i] += dlambdaoc[i];
      // Normalize the diffsq by the typical value of the rapidities:
      if (cxL > 1.0) diffsq += dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
      else diffsq += cxL * cxL * dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
    }
    diffsq /= DP(N);

    return;
  }

  void LiebLin_Bethe_State::Iterate_BAE_S (DP damping)
  {
    // This is essentially Newton's method but only in one variable.
    // The logic is that the derivative of the LHS of the BE_j w/r to lambdaoc_j is much larger
    // than with respect to lambdaoc_l with l != j.

    Vect_DP dlambdaoc (0.0, N);

    // Start by calculating S and dSdlambdaoc:
    for (int j = 0; j < N; ++j) {

      S[j] = 0.0;
      for (int k = 0; k < N; ++k) S[j] += atan((lambdaoc[j] - lambdaoc[k]));
      S[j] *= 2.0/cxL;

      dSdlambdaoc[j] = 0.0;
      for (int k = 0; k < N; ++k)
	dSdlambdaoc[j] += 1.0/((lambdaoc[j] - lambdaoc[k]) * (lambdaoc[j] - lambdaoc[k]) + 1.0);
      dSdlambdaoc[j] *= 2.0/(PI * cxL);

      dlambdaoc[j] = (PI*Ix2[j]/cxL - S[j] + lambdaoc[j] * dSdlambdaoc[j])
	/(1.0 + dSdlambdaoc[j]) - lambdaoc[j];

    }

    diffsq = 0.0;
    for (int i = 0; i < N; ++i) {
      lambdaoc[i] += damping * dlambdaoc[i];
      // Normalize the diffsq by the typical value of the rapidities:
      if (cxL > 1.0) diffsq += dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
      else diffsq += cxL * cxL * dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
    }
    diffsq /= DP(N);

    iter_Newton++;

    return;
  }

  /**
     Performs one step of the matrix Newton method on the rapidities.
   */
  void LiebLin_Bethe_State::Iterate_BAE_Newton (DP damping, Vect_DP& RHSBAE, Vect_DP& dlambdaoc, SQMat_DP& Gaudin, Vect_INT& indx)
  {
    DP sumtheta = 0.0;
    // for large |lambda|, use atan (lambda) = sgn(lambda) pi/2 - atan(1/lambda)
    int atanintshift = 0;
    DP lambdahere = 0.0;

    // Compute the RHS of the BAEs:

    for (int j = 0; j < N; ++j) {

      sumtheta = 0.0;
      atanintshift = 0;
      for (int k = 0; k < N; ++k)
	if (j != k) { // otherwise 0
	  if (fabs(lambdahere = lambdaoc[j] - lambdaoc[k]) < 1.0) { // use straight atan
	    sumtheta += atan(lambdahere);
	  }
	  else { // for large rapidities, use dual form of atan, extracting pi/2 factors
	    atanintshift += sgn_DP(lambdahere);
	    sumtheta -= atan(1.0/lambdahere);
	  }
	}
      sumtheta *= 2.0;

      RHSBAE[j] = cxL * lambdaoc[j] + sumtheta - PI*(Ix2[j] - atanintshift);
    }

    (*this).Build_Reduced_Gaudin_Matrix (Gaudin);

    for (int j = 0; j < N; ++j) dlambdaoc[j] = - RHSBAE[j];

    DP d;
    ludcmp (Gaudin, indx, d);
    lubksb (Gaudin, indx, dlambdaoc);

    bool ordering_changed = false;
    for (int j = 0; j < N-1; ++j)
      if (lambdaoc[j] + dlambdaoc[j] > lambdaoc[j+1] + dlambdaoc[j+1]) ordering_changed = true;

    // To prevent Newton from diverging, we limit the size of the rapidity changes.
    // The leftmost and rightmost rapidities can grow by one order of magnitude per iteration step.
    if (ordering_changed) {
      // We explicitly ensure that the ordering remains correct after the iteration step.
      bool ordering_still_changed = false;
      DP maxdlambdaoc = 0.0;
      do {
	ordering_still_changed = false;
	if (dlambdaoc[0] < 0.0 && fabs(dlambdaoc[0])
	    > (maxdlambdaoc = 10.0*ABACUS::max(fabs(lambdaoc[0]), fabs(lambdaoc[N-1]))))
	  dlambdaoc[0] = -maxdlambdaoc;
	if (lambdaoc[0] + dlambdaoc[0] > lambdaoc[1] + dlambdaoc[1]) {
	  dlambdaoc[0] = 0.25 * (lambdaoc[1] + dlambdaoc[1] - lambdaoc[0] ); // max quarter distance
	  ordering_still_changed = true;
	}
	if (dlambdaoc[N-1] > 0.0 && fabs(dlambdaoc[N-1])
	    > (maxdlambdaoc = 10.0*ABACUS::max(fabs(lambdaoc[0]), fabs(lambdaoc[N-1]))))
	  dlambdaoc[N-1] = maxdlambdaoc;
	if (lambdaoc[N-1] + dlambdaoc[N-1] < lambdaoc[N-2] + dlambdaoc[N-2]) {
	  dlambdaoc[N-1] = 0.25 * (lambdaoc[N-2] + dlambdaoc[N-2] - lambdaoc[N-1]);
	  ordering_still_changed = true;
	}
	for (int j = 1; j < N-1; ++j) {
	  if (lambdaoc[j] + dlambdaoc[j] > lambdaoc[j+1] + dlambdaoc[j+1]) {
	    dlambdaoc[j] = 0.25 * (lambdaoc[j+1] + dlambdaoc[j+1] - lambdaoc[j]);
	    ordering_still_changed = true;
	  }
	  if (lambdaoc[j] + dlambdaoc[j] < lambdaoc[j-1] + dlambdaoc[j-1]) {
	    dlambdaoc[j] = 0.25 * (lambdaoc[j-1] + dlambdaoc[j-1] - lambdaoc[j]);
	    ordering_still_changed = true;
	  }
	}
      } while (ordering_still_changed);
    }

    diffsq = 0.0;
    for (int i = 0; i < N; ++i) {
      // Normalize the diffsq by the typical value of the rapidities:
      if (cxL > 1.0) diffsq += dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
      else diffsq += cxL * cxL * dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
    }
    diffsq /= DP(N);

    if (ordering_changed) diffsq = 1.0; // reset if Newton wanted to change ordering

    for (int j = 0; j < N; ++j) lambdaoc[j] += damping * dlambdaoc[j];

    iter_Newton++;

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

    if (diffsq < prec && !ordering_changed) {

      lnnorm = 0.0;
      for (int j = 0; j < N; j++) lnnorm += log(fabs(Gaudin[j][j]));

      // Add the pieces outside of Gaudin determinant
      for (int j = 0; j < N - 1; ++j)
	for (int k = j+1; k < N; ++k)
	  lnnorm += log(1.0 + 1.0/pow(lambdaoc[j] - lambdaoc[k], 2.0));
    }

    return;
  }

  void LiebLin_Bethe_State::Compute_Energy()
  {
    E = 0.0;
    for (int j = 0; j < N; ++j) E += lambdaoc[j] * lambdaoc[j];
    E *= c_int * c_int;
  }

  void LiebLin_Bethe_State::Compute_Momentum()
  {
    iK = 0;
    for (int j = 0; j < N; ++j) {
      iK += Ix2[j];
    }
    if (iK % 2) {
      cout << Ix2 << endl;
      cout << iK << "\t" << iK % 2 << endl;
      ABACUSerror("Sum of Ix2 is not even:  inconsistency.");
    }
    iK /= 2;  // sum of Ix2 is guaranteed even.

    K = 2.0 * iK * PI/L;
  }

  /**
     The fundamental scattering kernel for Lieb-Liniger,
     \f[
     K (\tilde{\lambda}) = \frac{ 2 }{\tilde{\lambda}^2 + 1}
     \f]
     given two indices for the rapidities as arguments.
   */
  DP LiebLin_Bethe_State::Kernel (int a, int b)
  {
    return(2.0/(pow(lambdaoc[a] - lambdaoc[b], 2.0) + 1.0));
  }

  /**
     The fundamental scattering kernel for Lieb-Liniger,
     \f[
     K (\tilde{\lambda}) = \frac{ 2 }{\tilde{\lambda}^2 + 1}
     \f]
     given the rapidity difference as argument.
   */
  DP LiebLin_Bethe_State::Kernel (DP lambdaoc_ref)
  {
    return(2.0/(lambdaoc_ref * lambdaoc_ref + 1.0));
  }

  /**
     This function constructs the Gaudin matrix \f$G\f$
     which is defined as the Hessian of the Yang-Yang action \f$S^{YY}\f$,
     \f[
     G_{jk} = \delta_{jk} \left\{ cL + \sum_{l} K (\tilde{\lambda}_j - \tilde{\lambda}_l) \right\}
     - K (\tilde{\lambda}_j - \tilde{\lambda}_k).
     \f]
   */
  void LiebLin_Bethe_State::Build_Reduced_Gaudin_Matrix (SQMat<DP>& Gaudin_Red)
  {

    if (Gaudin_Red.size() != N)
      ABACUSerror("Passing matrix of wrong size in Build_Reduced_Gaudin_Matrix.");

    DP sum_Kernel = 0.0;

    for (int j = 0; j < N; ++j)
      for (int k = 0; k < N; ++k) {

	if (j == k) {
	  sum_Kernel = 0.0;
	  for (int kp = 0; kp < N; ++kp)
	    if (j != kp) sum_Kernel += Kernel (lambdaoc[j] - lambdaoc[kp]);
	  Gaudin_Red[j][k] = cxL + sum_Kernel;
	}

	else Gaudin_Red[j][k] = - Kernel (lambdaoc[j] - lambdaoc[k]);

      }

    return;
  }

  /**
     For a summetric `LiebLin_Bethe_State`, this function computes the Gaudin matrix
     as an \f$N/2 \times N/2\f$ matrix to accelerate computations.
   */

  void LiebLin_Bethe_State::Build_Reduced_BEC_Quench_Gaudin_Matrix (SQMat<DP>& Gaudin_Red)
  {
    // Passing  a matrix of dimension N/2

    if (N % 2 != 0) ABACUSerror("Choose a state with even numer of particles please");

    // Check Parity invariant

    bool ck = true;

    for (int j = 0; j < N/2; ++j){  if(Ix2[j] != - Ix2[N-j-1])  ck = false;}

    if (!ck) ABACUSerror("Choose a parity invariant state please");

    if (Gaudin_Red.size() != N/2)
      ABACUSerror("Passing matrix of wrong size in Build_Reduced_Gaudin_Matrix.");

    DP sum_Kernel = 0.0;

    for (int j = 0; j < N/2; ++j)
      for (int k = 0; k < N/2; ++k) {

	if (j == k) {
	  sum_Kernel = 0.0;
	  for (int kp = N/2; kp < N; ++kp)
	    if (j + N/2 != kp)
	      sum_Kernel += Kernel (lambdaoc[j+N/2] - lambdaoc[kp])
		+ Kernel (lambdaoc[j+N/2] + lambdaoc[kp]);
	  Gaudin_Red[j][k] = cxL + sum_Kernel;
	}

	else Gaudin_Red[j][k] = - (Kernel (lambdaoc[j+ N/2] - lambdaoc[k+ N/2])
				   + Kernel (lambdaoc[j+ N/2] + lambdaoc[k+ N/2]) );

      }

    return;
  }


  /**
     This function changes the Ix2 of a given state by annihilating a particle and hole
     pair specified by ipart and ihole (counting from the left, starting with index 0).
   */
  void LiebLin_Bethe_State::Annihilate_ph_pair (int ipart, int ihole,
						const Vect<int>& OriginStateIx2)
  {

    State_Label_Data currentdata = Read_State_Label ((*this).label, OriginStateIx2);

    if (ipart >= currentdata.nexc[0])
      ABACUSerror("Particle label too large in LiebLin_Bethe_State::Annihilate_ph_pair.");
    if (ihole >= currentdata.nexc[0])
      ABACUSerror("Hole label too large in LiebLin_Bethe_State::Annihilate_ph_pair.");

    // Simply remove the given pair:
    Vect<int> type_new = currentdata.type;
    Vect<int> M_new = currentdata.M;
    Vect<int> nexc_new = currentdata.nexc;
    nexc_new[0] -= 1; // we drill one more particle-hole pair at level 0
    int ntypespresent = 1; // only one type for LiebLin
    Vect<Vect<int> > Ix2old_new(ntypespresent);
    Vect<Vect<int> > Ix2exc_new(ntypespresent);
    for (int it = 0; it < ntypespresent; ++it)
      Ix2old_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
    for (int it = 0; it < ntypespresent; ++it)
      Ix2exc_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));

    // Copy earlier data in, leaving out ipart and ihole:
    for (int it = 0; it < ntypespresent; ++it) {
      for (int i = 0; i < nexc_new[it]; ++i) {
	Ix2old_new[it][i] = currentdata.Ix2old[it][i + (i >= ihole)];
	Ix2exc_new[it][i] = currentdata.Ix2exc[it][i + (i >= ipart)];
      }
    }

    State_Label_Data newdata (type_new, M_new, nexc_new, Ix2old_new, Ix2exc_new);

    (*this).Set_to_Label (Return_State_Label(newdata, OriginStateIx2));
  }

  /**
     Performs a spatial inversion of a `LiebLin_Bethe_State`, mapping all quantum numbers
     and rapidities to minus themselves.
   */
  void LiebLin_Bethe_State::Parity_Flip ()
  {
    Vect_INT Ix2buff = Ix2;
    Vect_DP lambdaocbuff = lambdaoc;
    for (int i = 0; i < N; ++i) Ix2[i] = -Ix2buff[N - 1 - i];
    for (int i = 0; i < N; ++i) lambdaoc[i] = -lambdaocbuff[N - 1 - i];
    iK = -iK;
    K = -K;
  }

  /**
     Send information about a `LiebLin_Bethe_State` to the output stream.
   */
  std::ostream& operator<< (std::ostream& s, const LiebLin_Bethe_State& state)
  {
    s << endl << "******** State for c = " << state.c_int << " L = " << state.L
      << " N = " << state.N << " with label " << state.label << " ********" << endl;
    s << "Ix2:" << endl;
    for (int j = 0; j < state.N; ++j) s << state.Ix2[j] << " ";
    s << endl << "lambdaocs:" << endl;
    for (int j = 0; j < state.N; ++j) s << state.lambdaoc[j] << " ";
    s << endl << "conv = " << state.conv << " iter_Newton = " << state.iter_Newton << endl;
    s << "E = " << state.E << " iK = " << state.iK << " K = " << state.K
      << " lnnorm = " << state.lnnorm << endl;

    return(s);
  }


} // namespace ABACUS