ABACUS/src/TBA/TBA_LiebLin.cc

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

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

Copyright (c) J.-S. Caux.

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

File:  TBA_LiebLin.cc

Purpose:  TBA for Lieb-Liniger

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

#include "ABACUS.h"

using namespace std;
using namespace ABACUS;

namespace ABACUS {

  // T = 0 functions

  DP LiebLin_a2_kernel (DP c_int, DP lambdadiff)
  {
    return(c_int/(PI * (lambdadiff * lambdadiff + c_int * c_int)));
  }

  void Iterate_LiebLin_rho_GS (DP c_int, DP k_F, Root_Density& rho_GS)
  {
    // Performs an iteration of the TBA equation for the dressed charge:

    // First, copy the actual values into the previous ones:
    for (int i = 0; i < rho_GS.Npts; ++i) rho_GS.prev_value[i] = rho_GS.value[i];

    // Calculate the new values:
    DP conv = 0.0;
    for (int i = 0; i < rho_GS.Npts; ++i) {
      // First, calculate the convolution
      conv = 0.0;
      for (int j = 0; j < rho_GS.Npts; ++j)
	conv += fabs(rho_GS.lambda[j]) > k_F ? 0.0 :
	  LiebLin_a2_kernel (c_int, rho_GS.lambda[i] - rho_GS.lambda[j])
	  * rho_GS.prev_value[j] * rho_GS.dlambda[j];
      rho_GS.value[i] = fabs(rho_GS.lambda[i]) < k_F ? 1.0/twoPI + conv : 0.0;
    }

    // Calculate the sum of differences:
    rho_GS.diff = 0.0;
    for (int i = 0; i < rho_GS.Npts; ++i)
      rho_GS.diff += fabs(rho_GS.value[i] - rho_GS.prev_value[i]);

    return;
  }

  Root_Density LiebLin_rho_GS (DP c_int, DP k_F, DP lambdamax, int Npts, DP req_prec)
  {

    // This function returns a Root_Density object corresponding to the ground state
    // density function of the Lieb-Liniger model.

    Root_Density rho_GS(Npts, lambdamax);

    // We now attempt to find a solution...

    do {
      Iterate_LiebLin_rho_GS (c_int, k_F, rho_GS);
    } while (rho_GS.diff > req_prec);

    return(rho_GS);
  }

  DP Density_GS (Root_Density& rho_GS)
  {
    DP n = 0.0;
    for (int i = 0; i < rho_GS.Npts; ++i) n += rho_GS.value[i] * rho_GS.dlambda[i];

    return(n);
  }

  DP k_F_given_n (DP c_int, DP n, DP lambdamax, int Npts, DP req_prec)
  {
    DP k_F = 1.0;
    DP dk_F = k_F;
    Root_Density rho_GS = LiebLin_rho_GS (c_int, k_F, lambdamax, Npts, req_prec);
    DP n_found = Density_GS (rho_GS);
    do {
      dk_F *= 0.7;
      k_F += (n_found < n) ? dk_F : -dk_F;
      rho_GS = LiebLin_rho_GS (c_int, k_F, lambdamax, Npts, req_prec);
      n_found = Density_GS (rho_GS);

    } while (fabs(dk_F) > req_prec && dk_F > 1.0/Npts
	     && fabs(n - n_found) > req_prec && fabs(n - n_found) > 1.0/Npts);

    return(k_F);
  }

  void Iterate_LiebLin_Z_GS (DP c_int, DP k_F, Root_Density& Z_GS)
  {
    // Performs an iteration of the TBA equation for the dressed charge:

    // First, copy the actual values into the previous ones:
    for (int i = 0; i < Z_GS.Npts; ++i) Z_GS.prev_value[i] = Z_GS.value[i];

    // Calculate the new values:
    DP conv = 0.0;
    for (int i = 0; i < Z_GS.Npts; ++i) {
      // First, calculate the convolution
      conv = 0.0;
      for (int j = 0; j < Z_GS.Npts; ++j)
	conv += fabs(Z_GS.lambda[j]) > k_F ? 0.0 :
	  LiebLin_a2_kernel (c_int, Z_GS.lambda[i] - Z_GS.lambda[j])
	  * Z_GS.prev_value[j] * Z_GS.dlambda[j];
      Z_GS.value[i] = 1.0 + conv;
    }

    // Calculate the sum of differences:
    Z_GS.diff = 0.0;
    for (int i = 0; i < Z_GS.Npts; ++i)
      Z_GS.diff += fabs(Z_GS.value[i] - Z_GS.prev_value[i]);

    return;
  }

  Root_Density LiebLin_Z_GS (DP c_int, DP k_F, DP lambdamax, int Npts, DP req_prec)
  {

    // This function returns a Root_Density object corresponding to the ground state
    // dressed charge function of the Lieb-Liniger model.

    Root_Density Z_GS(Npts, lambdamax);

    // We now attempt to find a solution...

    do {
      Iterate_LiebLin_Z_GS (c_int, k_F, Z_GS);
    } while (Z_GS.diff > req_prec);

    return(Z_GS);
  }

  void Iterate_LiebLin_Fbackflow_GS (DP c_int, DP k_F, DP lambda, Root_Density& F_GS)
  {
    // Performs an iteration of the TBA equation for the backflow:

    // First, copy the actual values into the previous ones:
    for (int i = 0; i < F_GS.Npts; ++i) F_GS.prev_value[i] = F_GS.value[i];

    // Calculate the new values:
    DP conv = 0.0;
    for (int i = 0; i < F_GS.Npts; ++i) {
      // First, calculate the convolution
      conv = 0.0;
      for (int j = 0; j < F_GS.Npts; ++j)
	conv += fabs(F_GS.lambda[j]) > k_F ? 0.0 :
	  LiebLin_a2_kernel (c_int, F_GS.lambda[i] - F_GS.lambda[j])
	  * F_GS.prev_value[j] * F_GS.dlambda[j];
      F_GS.value[i] = 0.5 + atan((F_GS.lambda[i] - lambda)/c_int)/PI + conv;
    }

    // Calculate the sum of differences:
    F_GS.diff = 0.0;
    for (int i = 0; i < F_GS.Npts; ++i)
      F_GS.diff += fabs(F_GS.value[i] - F_GS.prev_value[i]);

    return;
  }

  Root_Density LiebLin_Fbackflow_GS (DP c_int, DP k_F, DP lambdamax, DP lambda, int Npts, DP req_prec)
  {

    // This function returns a Root_Density object corresponding to the ground state
    // F backflow for parameter \lambda.

    Root_Density F_GS(Npts, lambdamax);

    // We now attempt to find a solution...

    do {
      Iterate_LiebLin_Fbackflow_GS (c_int, k_F, lambda, F_GS);
    } while (F_GS.diff > req_prec);

    return(F_GS);
  }



  // Finite T functions:

  LiebLin_TBA_Solution::LiebLin_TBA_Solution (DP c_int_ref, DP mu_ref, DP kBT_ref, DP req_diff, int Max_Secs)
    : c_int(c_int_ref), mu(mu_ref), kBT(kBT_ref)
  {
    epsilon = LiebLin_epsilon_TBA (c_int, mu, kBT, req_diff, Max_Secs);
    depsilon_dmu = LiebLin_depsilon_dmu_TBA (c_int, mu, kBT, req_diff, Max_Secs, epsilon);
    rho = LiebLin_rho_TBA (kBT, epsilon, depsilon_dmu);
    rhoh = LiebLin_rhoh_TBA (kBT, epsilon, depsilon_dmu);
    nbar = LiebLin_nbar_TBA (rho);
    ebar = LiebLin_ebar_TBA (rho);
    sbar = LiebLin_sbar_TBA (rho, rhoh);
  }


  LiebLin_TBA_Solution LiebLin_TBA_Solution_fixed_nbar_ebar (DP c_int, DP nbar_required, DP ebar_required,
							     DP req_diff, int Max_Secs)
  {
    // This function finds the TBA solution for a required nbar and ebar (mean energy).
    // We here try to triangulate the temperature; the chemical potential is triangulated
    // using the function LiebLin_TBA_Solution_fixed_nbar.

    // FIRST VERSION
    DP lnkBT = 1.0;
    DP dlnkBT = 1.0;

    LiebLin_TBA_Solution tbasol = LiebLin_TBA_Solution_fixed_nbar (c_int, nbar_required, exp(lnkBT), req_diff, Max_Secs);
    DP lnkBT_prev = lnkBT;
    DP ebar_prev = LiebLin_ebar_TBA (tbasol.rho);

    DP dlnkBT_prev, debardlnkBT;

    lnkBT += dlnkBT;
    tbasol = LiebLin_TBA_Solution_fixed_nbar (c_int, nbar_required, exp(lnkBT), req_diff, Max_Secs);

    int niter = 1;
    while (fabs(tbasol.ebar - ebar_required) > req_diff) {

      dlnkBT_prev = dlnkBT;
      lnkBT_prev = lnkBT;
      // Calculate a new kBT based on earlier data:
      debardlnkBT = (tbasol.ebar - ebar_prev)/dlnkBT;
      ebar_prev = tbasol.ebar;
      dlnkBT = (ebar_required - tbasol.ebar)/debardlnkBT;
      if (fabs(dlnkBT) > fabs(dlnkBT_prev)) dlnkBT *= fabs(dlnkBT_prev/dlnkBT); // make sure we don't blow up.
      lnkBT += dlnkBT;
      tbasol = LiebLin_TBA_Solution_fixed_nbar(c_int, nbar_required, exp(lnkBT), req_diff, Max_Secs);
      niter++;
    }

    return(tbasol);
  }



  LiebLin_TBA_Solution LiebLin_TBA_Solution_fixed_nbar (DP c_int, DP nbar_required, DP kBT, DP req_diff, int Max_Secs)
  {
    // Find the required mu. Provide some initial guess.

    DP mu = 2.0;
    DP dmu = 1.0;

    LiebLin_TBA_Solution tbasol = LiebLin_TBA_Solution(c_int, mu, kBT, req_diff, Max_Secs);
    DP mu_prev = mu;
    DP nbar_prev = tbasol.nbar;

    DP dmu_prev, dnbardmu;

    mu += dmu;
    tbasol = LiebLin_TBA_Solution(c_int, mu, kBT, req_diff, Max_Secs);

    int niter = 1;
    while (fabs(tbasol.nbar - nbar_required) > req_diff) {

      dmu_prev = dmu;
      mu_prev = mu;
      // Calculate a new mu based on earlier data:
      dnbardmu = (tbasol.nbar - nbar_prev)/dmu;
      nbar_prev = tbasol.nbar;
      dmu = (nbar_required - tbasol.nbar)/dnbardmu;
      if (fabs(dmu) > 2.0 * fabs(dmu_prev)) dmu = 2.0*dmu * fabs(dmu_prev/dmu); // make sure we don't blow up.
      mu += dmu;
      tbasol = LiebLin_TBA_Solution(c_int, mu, kBT, req_diff, Max_Secs);
      niter++;
    }

    return(tbasol);
  }


  DP Refine_LiebLin_epsilon_TBA (Root_Density& epsilon, DP c_int, DP mu, DP kBT, DP refine_fraction)
  {
    // This function replaces Set by a new set with more points, where
    // Tln(...) needs to be evaluated more precisely.

    // The return value is the max of delta_tni found.

    // First, calculate the needed Tln...
    Vect<DP> Tln1plusemineps(epsilon.Npts);

    for (int i = 0; i < epsilon.Npts; ++i) {
      Tln1plusemineps[i] = epsilon.value[i] > 0.0 ?
	kBT * (epsilon.value[i] < 24.0 * kBT ? log(1.0 + exp(-epsilon.value[i]/kBT)) : exp(-epsilon.value[i]/kBT))
	: -epsilon.value[i] + kBT * log (1.0 + exp(epsilon.value[i]/kBT));
    }

    // Now find the achieved delta_tni
    DP max_delta_tni_dlambda = 0.0;
    DP sum_delta_tni_dlambda = 0.0;

    Vect<DP> tni(0.0, epsilon.Npts);
    Vect<DP> tni_ex(0.0, epsilon.Npts);

    DP measure_factor = 0.0;

    for (int i = 1; i < epsilon.Npts - 1; ++i) {
      measure_factor = (epsilon.value[i] > 0.0 ? exp(-epsilon.value[i]/kBT)/(1.0 + exp(-epsilon.value[i]/kBT))
			: 1.0/(1.0 + exp(epsilon.value[i]/kBT)));

      tni[i] = measure_factor * epsilon.value[i];
      tni_ex[i] = measure_factor * (epsilon.value[i-1] *  (epsilon.lambda[i+1] - epsilon.lambda[i])
				    + epsilon.value[i+1] * (epsilon.lambda[i] - epsilon.lambda[i-1]))
	/(epsilon.lambda[i+1] - epsilon.lambda[i-1]);

      max_delta_tni_dlambda = ABACUS::max(max_delta_tni_dlambda, fabs(tni[i] - tni_ex[i]) * epsilon.dlambda[i]);

      sum_delta_tni_dlambda += fabs(tni[i] - tni_ex[i]) * epsilon.dlambda[i];
    }

    // We now determine the locations where we need to add points
    Vect<bool> need_new_point_around(false, epsilon.Npts);
    int nr_new_points_needed = 0;

    for (int i = 1; i < epsilon.Npts - 1; ++i) {
      if (fabs(tni[i] - tni_ex[i]) * epsilon.dlambda[i] > (1.0 - refine_fraction) * max_delta_tni_dlambda) {
	need_new_point_around[i] = true;
	// Do also the symmetric ones...  Require need...[n][i] = need...[n][Npts - 1 - i]
	need_new_point_around[epsilon.Npts - 1 - i] = true;
      }
    }
    for (int i = 0; i < epsilon.Npts; ++i)
      if (need_new_point_around[i]) nr_new_points_needed++;

    // Now insert the new points between existing points:
    // Update all data via interpolation:
    Root_Density epsilon_before_update = epsilon;
    epsilon = Root_Density(epsilon_before_update.Npts + nr_new_points_needed, epsilon_before_update.lambdamax);

    int nr_pts_added = 0;
    for (int i = 0; i < epsilon_before_update.Npts; ++i) {
      if (!need_new_point_around[i]) {
	epsilon.lambda[i + nr_pts_added] = epsilon_before_update.lambda[i];
	epsilon.dlambda[i + nr_pts_added] = epsilon_before_update.dlambda[i];
	epsilon.value[i + nr_pts_added] = epsilon_before_update.value[i];
      }
      else if (need_new_point_around[i]) {
	epsilon.lambda[i + nr_pts_added] = epsilon_before_update.lambda[i] - 0.25 * epsilon_before_update.dlambda[i];
	epsilon.dlambda[i + nr_pts_added] = 0.5 * epsilon_before_update.dlambda[i];
	epsilon.value[i + nr_pts_added] = epsilon_before_update.Return_Value(epsilon.lambda[i + nr_pts_added]);
	nr_pts_added++;
	epsilon.lambda[i + nr_pts_added] = epsilon_before_update.lambda[i] + 0.25 * epsilon_before_update.dlambda[i];
	epsilon.dlambda[i + nr_pts_added] = 0.5 * epsilon_before_update.dlambda[i];
	epsilon.value[i + nr_pts_added] = epsilon_before_update.Return_Value(epsilon.lambda[i + nr_pts_added]);
      }
    }
    if (nr_pts_added != nr_new_points_needed) {
      cout << nr_pts_added << "\t" << nr_new_points_needed << endl;
      ABACUSerror("Wrong counting of new points in Insert_new_points.");
    }


    // Check boundary values;  if too different from value_infty, extend limits
    if (exp(-epsilon.value[0]/kBT) > 0.001 * max_delta_tni_dlambda) {

      int nr_pts_to_add_left = epsilon.Npts/10;
      int nr_pts_to_add_right = epsilon.Npts/10;
      Root_Density epsilon_before_update = epsilon;
      epsilon = Root_Density(epsilon_before_update.Npts + nr_pts_to_add_left + nr_pts_to_add_right,
			     epsilon_before_update.lambdamax + nr_pts_to_add_left * epsilon_before_update.dlambda[0]
			     + nr_pts_to_add_right * epsilon_before_update.dlambda[epsilon_before_update.Npts - 1]);

      // Initialize points added on left
      for (int i = 0; i < nr_pts_to_add_left; ++i) {
	epsilon.lambda[i] = epsilon_before_update.lambda[i] + (i - nr_pts_to_add_left) * epsilon_before_update.dlambda[0];
	epsilon.dlambda[i] = epsilon_before_update.dlambda[0];
	epsilon.value[i] = epsilon_before_update.value[0];
      }
      // Initialize middle (unchanged) points
      for (int i = 0; i < epsilon_before_update.Npts; ++i) {
	epsilon.lambda[nr_pts_to_add_left + i] = epsilon_before_update.lambda[i];
	epsilon.dlambda[nr_pts_to_add_left + i] = epsilon_before_update.dlambda[i];
	epsilon.value[nr_pts_to_add_left + i] = epsilon_before_update.value[i];
      }
      // Initialize points added on right
      for (int i = 0; i < nr_pts_to_add_right; ++i) {
	epsilon.lambda[nr_pts_to_add_left + epsilon_before_update.Npts + i] =
	  epsilon_before_update.lambda[epsilon_before_update.Npts - 1]
	  + (i+1) * epsilon_before_update.dlambda[epsilon_before_update.Npts - 1];
	epsilon.dlambda[nr_pts_to_add_left + epsilon_before_update.Npts + i] =
	  epsilon_before_update.dlambda[epsilon_before_update.Npts - 1];
	epsilon.value[nr_pts_to_add_left + epsilon_before_update.Npts + i] =
	  epsilon_before_update.value[epsilon_before_update.Npts - 1];
      }
    }

    return(sum_delta_tni_dlambda);
  }

  // epsilon for a given mu:
  Root_Density LiebLin_epsilon_TBA (DP c_int, DP mu, DP kBT, DP req_diff, int Max_Secs)
  {

    clock_t StartTime = clock();
    int Max_CPU_ticks = 98 * (Max_Secs - 0) * CLOCKS_PER_SEC/100;
    // give 30 seconds to wrap up, assume we time to 2% accuracy.

    // Set basic precision needed:
    DP running_prec = 1.0;

    DP refine_fraction = 0.5; // value fraction of points to be refined

    DP lambdamax_init = 10.0 * sqrt(ABACUS::max(1.0, kBT + mu));
    // such that exp(-(lambdamax^2 - mu - Omega)/T) <~ machine_eps

    int Npts = 50;

    Root_Density epsilon(Npts, lambdamax_init);

    // Initiate the function:
    for (int i = 0; i < epsilon.Npts; ++i) {
      epsilon.value[i] = epsilon.lambda[i] * epsilon.lambda[i] - mu;
      epsilon.prev_value[i] = epsilon.value[i];
    }

    clock_t StopTime = clock();

    int CPU_ticks = StopTime - StartTime;

    int ncycles = 0;
    int niter_tot = 0;
    DP previous_running_prec;

    DP oneoverpi = 1.0/PI;
    DP oneoverc = 1.0/c_int;

    do {

      StartTime = clock();

      // The running precision is an estimate of the accuracy of the free energy integral.
      // Refine... returns sum_delta_tni_dlambda, so running prec is estimated as...
      previous_running_prec = running_prec;
      running_prec = Refine_LiebLin_epsilon_TBA (epsilon, c_int, mu, kBT, refine_fraction);
      running_prec = ABACUS::min(running_prec, previous_running_prec);

      // Now iterate to convergence for given discretization
      int niter = 0;
      int niter_max = ncycles == 0 ? 300 : 300;

      do {

	StartTime = clock();

	// Iterate the TBA equation for epsilon:
	Vect<DP> Tln1plusemineps(epsilon.Npts);
	for (int i = 0; i < epsilon.Npts; ++i) {
	  Tln1plusemineps[i] = epsilon.value[i] > 0.0 ?
	    kBT * (epsilon.value[i] < 24.0 * kBT
		   ? log(1.0 + exp(-epsilon.value[i]/kBT))
		   : exp(-epsilon.value[i]/kBT) * (1.0 - 0.5 * exp(-epsilon.value[i]/kBT)))
	    :
	    -epsilon.value[i] + kBT * (-epsilon.value[i] < 24.0 * kBT
				       ? log (1.0 + exp(epsilon.value[i]/kBT))
				       : exp(epsilon.value[i]/kBT) * (1.0 - 0.5 * exp(epsilon.value[i]/kBT)));
	  // Keep previous rapidities:
	  epsilon.prev_value[i] = epsilon.value[i];
	}

	Vect<DP> a_2_Tln_conv(epsilon.Npts);

	for (int i = 0; i < epsilon.Npts; ++i) {
	  a_2_Tln_conv[i] = 0.0;

	  for (int j = 0; j < epsilon.Npts; ++j)
	    a_2_Tln_conv[i] +=
	      oneoverpi * (atan(oneoverc * (epsilon.lambda[i] - epsilon.lambda[j] + 0.5 * epsilon.dlambda[j]))
			   - atan(oneoverc * (epsilon.lambda[i] - epsilon.lambda[j]
					      - 0.5 * epsilon.dlambda[j]))) * Tln1plusemineps[j];

	} // a_2_Tln_conv is now calculated

	// Reconstruct the epsilon function:
	for (int i = 0; i < epsilon.Npts; ++i) {
	  epsilon.value[i] = pow(epsilon.lambda[i], 2.0) - mu;

	  // Add the convolution:
	  epsilon.value[i] -= a_2_Tln_conv[i];

	  // Include some damping:
	  epsilon.value[i] = 0.1 * epsilon.prev_value[i] + 0.9 * epsilon.value[i];
	}

	niter++;
	// epsilon is now fully iterated.

	// Calculate the diff:
	epsilon.diff = 0.0;
	for (int i = 0; i < epsilon.Npts; ++i)
	  epsilon.diff += epsilon.dlambda[i] *
	    (epsilon.value[i] > 0.0 ?
	     exp(-epsilon.value[i]/kBT)/(1.0 + exp(-epsilon.value[i]/kBT)) : 1.0/(1.0 + exp(epsilon.value[i]/kBT)))
	    * fabs(epsilon.value[i] - epsilon.prev_value[i]);

	StopTime = clock();
	CPU_ticks += StopTime - StartTime;

      } while (niter < 5 || niter < niter_max && CPU_ticks < Max_CPU_ticks && epsilon.diff > 0.1*running_prec);

      ncycles++;
      niter_tot += niter;

    } // do cycles
    while (ncycles < 5 || running_prec > req_diff && CPU_ticks < Max_CPU_ticks);

    return(epsilon);
  }


  // depsilon/dmu for a given mu
  Root_Density LiebLin_depsilon_dmu_TBA (DP c_int, DP mu, DP kBT, DP req_diff, int Max_Secs, const Root_Density& epsilon)
  {

    clock_t StartTime = clock();
    int Max_CPU_ticks = 98 * (Max_Secs - 0) * CLOCKS_PER_SEC/100;
    // give 30 seconds to wrap up, assume we time to 2% accuracy.

    Root_Density depsilon_dmu = epsilon;

    // Initiate the functions:
    for (int i = 0; i < depsilon_dmu.Npts; ++i) {
      depsilon_dmu.value[i] = -1.0;
      depsilon_dmu.prev_value[i] = depsilon_dmu.value[i];
    }

    clock_t StopTime = clock();

    int CPU_ticks = StopTime - StartTime;

    int niter = 0;
    int niter_max = 1000;

    DP oneoverpi = 1.0/PI;
    DP oneoverc = 1.0/c_int;

    do {

      StartTime = clock();

      // Iterate the TBA equation for depsilon_dmu:
      Vect<DP> depsover1plusepluseps(depsilon_dmu.Npts);
      for (int i = 0; i < depsilon_dmu.Npts; ++i) {
	depsover1plusepluseps[i] = epsilon.value[i] > 0.0 ?
	  depsilon_dmu.value[i] * exp(-epsilon.value[i]/kBT)/(1.0 + exp(-epsilon.value[i]/kBT)) :
	  depsilon_dmu.value[i]/(1.0 + exp(epsilon.value[i]/kBT));

	// Keep previous rapidities:
	depsilon_dmu.prev_value[i] = depsilon_dmu.value[i];
      }

      Vect<DP> a_2_depsover1plusepluseps_conv(depsilon_dmu.Npts);

      for (int i = 0; i < depsilon_dmu.Npts; ++i) {
	a_2_depsover1plusepluseps_conv[i] = 0.0;

	for (int j = 0; j < depsilon_dmu.Npts; ++j)
	  a_2_depsover1plusepluseps_conv[i] +=
	    oneoverpi * (atan(oneoverc * (epsilon.lambda[i] - epsilon.lambda[j] + 0.5 * epsilon.dlambda[j]))
			 - atan(oneoverc * (epsilon.lambda[i] - epsilon.lambda[j] - 0.5 * epsilon.dlambda[j])))
	    * depsover1plusepluseps[j];
      }

      // Reconstruct the depsilon_dmu function:
      for (int i = 0; i < epsilon.Npts; ++i) {
	depsilon_dmu.value[i] = -1.0;

	// Add the convolution:
	depsilon_dmu.value[i] += a_2_depsover1plusepluseps_conv[i];

	// Include some damping:
	depsilon_dmu.value[i] = 0.1 * depsilon_dmu.prev_value[i] + 0.9 * depsilon_dmu.value[i];
      }
      niter++;
      // depsilon_dmu is now fully iterated.

      // Calculate the diff:
      depsilon_dmu.diff = 0.0;
      for (int i = 0; i < depsilon_dmu.Npts; ++i)
	depsilon_dmu.diff += fabs(depsilon_dmu.value[i] - depsilon_dmu.prev_value[i]);

      StopTime = clock();
      CPU_ticks += StopTime - StartTime;

    } while (niter < 5 || niter < niter_max && CPU_ticks < Max_CPU_ticks && depsilon_dmu.diff > req_diff);

    return(depsilon_dmu);
  }

  Root_Density LiebLin_rho_TBA (DP kBT, const Root_Density& epsilon, const Root_Density& depsilon_dmu)
  {
    Root_Density rho = epsilon;
    for (int i = 0; i < epsilon.Npts; ++i)
      rho.value[i] = -(1.0/twoPI) * depsilon_dmu.value[i]
	* (epsilon.value[i] > 0.0 ?
	   exp(-epsilon.value[i]/kBT)/(1.0 + exp(-epsilon.value[i]/kBT)) :
	   1.0/(1.0 + exp(epsilon.value[i]/kBT)));

    return(rho);
  }

  Root_Density LiebLin_rhoh_TBA (DP kBT, const Root_Density& epsilon, const Root_Density& depsilon_dmu)
  {
    Root_Density rhoh = epsilon;
    for (int i = 0; i < epsilon.Npts; ++i)
      rhoh.value[i] = -(1.0/twoPI) * depsilon_dmu.value[i]
	* (epsilon.value[i] > 0.0 ?
	   1.0/(1.0 + exp(-epsilon.value[i]/kBT)) :
	   exp(epsilon.value[i]/kBT)/(1.0 + exp(epsilon.value[i]/kBT)));

    return(rhoh);
  }

  DP LiebLin_nbar_TBA (const Root_Density& rho)
  {
    DP nbar = 0.0;
    for (int i = 0; i < rho.Npts; ++i) nbar += rho.value[i] * rho.dlambda[i];

    return(nbar);
  }

  DP LiebLin_ebar_TBA (const Root_Density& rho)
  {
    DP ebar = 0.0;
    for (int i = 0; i < rho.Npts; ++i) ebar += rho.lambda[i] * rho.lambda[i] * rho.value[i] * rho.dlambda[i];

    return(ebar);
  }

  DP LiebLin_sbar_TBA (const Root_Density& rho, const Root_Density& rhoh)
  {
    DP sbar = 0.0;
    for (int i = 0; i < rho.Npts; ++i)
      sbar += ((rho.value[i] + rhoh.value[i]) * log(rho.value[i] + rhoh.value[i])
	       - rho.value[i] * log(rho.value[i]+1.0e-30) - rhoh.value[i] * log(rhoh.value[i]+1.0e-30)) * rho.dlambda[i];

    return(sbar);
  }


  LiebLin_Bethe_State Discretized_LiebLin_Bethe_State (DP c_int, DP L, int N, const Root_Density& rho)
  {
    // This function returns the Bethe state at finite size which is
    // the closest approximation to the continuum density rho(lambda)

    // Check that the provided rho has the expected filling
    DP rho_check = 0.0;
    for (int i = 0; i < rho.Npts; ++i) rho_check += L * rho.value[i] * rho.dlambda[i];
    // The integral of rho should be between N - 0.5 and N + 0.5 for the algorithm to work
    if (fabs(rho_check - N) > 0.5)
      cout << "WARNING: integral of rho != N in Discretized_LiebLin_Bethe_State, "
	   << "consistent discretization is impossible. \int rho dlambda = "
	   << rho_check << ", N = " << N << ", L = " << L << ", N/L = " << N/L << endl;

    // Using the counting function
    // c(\lambda) \equiv L \int_{-\infty}^\lambda d\lambda' \rho(\lambda'),
    // we seek to find lambda[i] such that c(lambda[i]) = i + 0.5, i = 0...N-1
    DP integral = 0.0;
    DP integral_prev = 0.0;
    int Nfound = 0;
    Vect<DP> lambda_found(0.0, N);
    int nr_to_set = 0;

    for (int i = 0; i < rho.Npts; ++i) {
      // Note: integral_prev is the counting function
      // at rapidity rho.lambda[i-1] + 0.5* rho.dlambda[i-1]
      // which equals rho.lambda[i] - 0.5* rho.dlambda[i]
      integral_prev = integral;
      // Note: integral gives the value of the counting function
      // at rapidity rho.lambda[i] + 0.5* rho.dlambda[i]
      integral += L * rho.value[i] * rho.dlambda[i];
      // We already have filled the RHS of the counting equation up to Nfound - 0.5.
      // Therefore, the additional number found is
      nr_to_set = floor(integral + 0.5 - Nfound);
      // if (nr_to_set > 1)
      // 	cout << "WARNING: setting " << nr_to_set << " rapidities in one step in Discretized_LiebLin_Bethe_State" << endl;
      for (int n = 1; n <= nr_to_set; ++n) {
	// Solve c(lambda[Nfound]) = Nfound + 0.5 for lambda[Nfound],
	// namely (using linear interpolation)
	// integral_prev + (lambda - rho.lambda[i-1] - 0.5*rho.dlambda[i-1]) * (integral - integra_prev)/rho.dlambda[i] = Nfound + 0.5
	lambda_found[Nfound] = rho.lambda[i-1] + 0.5*rho.dlambda[i-1] + (Nfound + 0.5 - integral_prev) * rho.dlambda[i]/(integral - integral_prev);
	Nfound++;
      }
    }

    Vect<DP> lambda(N);
    // Fill up the found rapidities:
    for (int il = 0; il < abacus::min(N, Nfound); ++il) lambda[il] = lambda_found[il];
    // If there are missing ones, put them at the end; ideally, this should never be called
    for (int il = Nfound; il < N; ++il) lambda[il] = lambda_found[Nfound-1] + (il - Nfound + 1) * (lambda_found[Nfound-1] - lambda_found[Nfound-2]);

    // Calculate quantum numbers: 2\pi * (L lambda + \sum_j 2 atan((lambda - lambda_j)/c)) = I_j
    // Use rounding.
    Vect<int> Ix2(N);
    for (int i = 0; i < N; ++i) {
      DP sum = 0.0;
      for (int j = 0; j < N; ++j) sum += 2.0 * atan((lambda[i] - lambda[j])/c_int);
      // For N is even/odd, we want to round off to the nearest odd/even integer.
      Ix2[i] = 2.0 * floor((L* lambda[i] + sum)/twoPI + 0.5 * (N%2 ? 1 : 2)) + (N%2) - 1;
    }

    // Check that the Ix2 are all ordered
    for (int i = 0; i < N-1; ++i) if (Ix2[i] >= Ix2[i+1]) cout << "Alert: Ix2 not ordered around index i = " << i << ": Ix2[i] = " << Ix2[i] << "\tIx2[i+1] = " << Ix2[i+1] << endl;

    // Check that the quantum numbers are all distinct:
    bool allOK = false;
    while (!allOK) {
      for (int i = 0; i < N-1; ++i) if (Ix2[i] == Ix2[i+1] && Ix2[i] < 0) Ix2[i] -= 2;
      for (int i = 1; i < N; ++i) if (Ix2[i] == Ix2[i-1] && Ix2[i] > 0) Ix2[i] += 2;
      allOK = true;
      for (int i = 0; i < N-1; ++i) if (Ix2[i] == Ix2[i+1]) allOK = false;
    }

    LiebLin_Bethe_State rhostate(c_int, L, N);
    rhostate.Ix2 = Ix2;
    rhostate.Compute_All(true);

    return(rhostate);
  }


} // namespace ABACUS