ABACUS/src/TBA/TBA_2CBG.cc

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

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

Copyright (c) J.-S. Caux.

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

File:  TBA_2CBG.cc

Purpose:  solves the TBA equations for the 2-component Bose gas

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

#include "ABACUS.h"

using namespace std;
using namespace ABACUS;

namespace ABACUS {


  /********************* 2CBG specific *******************/


  DP Asymptotic_2CBG_epsilon (int n, DP Omega, DP kBT)
  {
    return(2.0 * Omega * n + kBT * log(pow((1.0 - exp(-2.0 * (n + 1.0) * Omega/kBT))
					   /(1.0 - exp(-2.0 * Omega/kBT)), 2.0) - exp(-2.0 * n * Omega/kBT)));
  }

  void Set_2CBG_Asymptotics (Root_Density_Set& TBA_Set, DP mu, DP Omega, DP kBT)
  {
    TBA_Set.epsilon[0].Set_Asymptotics (pow(TBA_Set.epsilon[0].lambdamax, 2.0) - mu - Omega);
    for (int n = 1; n < TBA_Set.ntypes; ++n) {
      TBA_Set.epsilon[n].Set_Asymptotics (Asymptotic_2CBG_epsilon(n, Omega, kBT));
    }
  }

  void Set_2CBG_deps_dchempot_Asymptotics (int option, const Root_Density_Set& Set, Root_Density_Set& DSet,
					   DP mu, DP Omega, DP kBT)
  {
    // option == 0:  deps/dmu
    // option == 1:  deps/dOmega

    DP zeroasymptote = -1.0;
    DP em2OoT = exp(-2.0 * Omega/kBT);
    for (int n = 1; n < DSet.ntypes; ++n) {
      if (option == 0) DSet.epsilon[n].Set_Asymptotics (0.0);
      else if (option == 1)
	DSet.epsilon[n].Set_Asymptotics (2.0 * (1.0 - pow(em2OoT, n+1.0))
					 * (n * (1.0 - pow(em2OoT, n+2.0)) - (n + 2.0) * em2OoT * (1.0 - pow(em2OoT, DP(n))))
					 /((1.0 - em2OoT) * (1.0 - pow(em2OoT, DP(n))) * (1.0 - pow(em2OoT, n + 2.0))));
      zeroasymptote += DSet.epsilon[n].value_infty
	* exp(-Set.epsilon[n].value_infty/kBT)/(1.0 + exp(-Set.epsilon[n].value_infty/kBT));
    }
    // For n > nmax sum in RHS of BE for epsilon, assuming epsilon_n = epsilon_n^\infty in those cases:
    // Remember:  nmax in notes is Set.ntypes - 1.
    zeroasymptote -= option == 0 ? 0.0
      : 2.0 * ((Set.ntypes + 1.0) * exp(-2.0 * (Set.ntypes + 1.0) * Omega/kBT)
	       /(1.0 - exp(-2.0 * (Set.ntypes + 1.0) * Omega/kBT))
	       - Set.ntypes * exp(-2.0 * Set.ntypes * Omega/kBT)/(1.0 - exp(-2.0 * Set.ntypes * Omega/kBT)));

    DSet.epsilon[0].Set_Asymptotics (zeroasymptote);

    return;
  }

  void Initiate_2CBG_TBA_Functions (Root_Density_Set& TBA_Set, DP mu, DP Omega)
  {
    for (int i = 0; i < TBA_Set.epsilon[0].Npts; ++i) {
      TBA_Set.epsilon[0].value[i] = TBA_Set.epsilon[0].lambda[i] * TBA_Set.epsilon[0].lambda[i] - mu - Omega;
      TBA_Set.epsilon[0].prev_value[i] = TBA_Set.epsilon[0].value[i];
    }
    for (int n = 1; n < TBA_Set.ntypes; ++n) {
      for (int i = 0; i < TBA_Set.epsilon[n].Npts; ++i)
	TBA_Set.epsilon[n].value[i] = TBA_Set.epsilon[n].value_infty;
    }
  }

  void Initiate_2CBG_deps_dchempot_Functions (Root_Density_Set& DSet)
  {
    for (int n = 0; n < DSet.ntypes; ++n) {
      for (int i = 0; i < DSet.epsilon[n].Npts; ++i)
	DSet.epsilon[n].value[i] = DSet.epsilon[n].value_infty;
    }
  }

  void Iterate_2CBG_TBAE (Root_Density_Set& Set, Vect<Vect<Vect_DP> >& a_n_dlambda, Vect<Vect<Vect_DP> >& fmin_dlambda,
			  Vect<Vect<Vect_DP> >& fplus_dlambda, DP c_int, DP mu, DP Omega, DP kBT)
  {
    // Produces a new Root_Density_Set from a previous iteration.
    // Does NOT add types or change Npts, lambdamax values.

    // First define some useful functions:
    Vect<Vect_DP> Tln1plusemineps(Set.ntypes);
    Vect_DP Tln1pluseminepsinfty(Set.ntypes);

    for (int n = 0; n < Set.ntypes; ++n) {

      Tln1plusemineps[n] = Vect_DP (0.0, Set.epsilon[n].Npts);

      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	Tln1plusemineps[n][i] = Set.epsilon[n].value[i] > 0.0 ?
	  kBT * (Set.epsilon[n].value[i] < 24.0 * kBT ? log(1.0 + exp(-Set.epsilon[n].value[i]/kBT))
		 : exp(-Set.epsilon[n].value[i]/kBT))
	  :
	  -Set.epsilon[n].value[i] + kBT * (-Set.epsilon[n].value[i] < 24.0 * kBT
					    ? log (1.0 + exp(Set.epsilon[n].value[i]/kBT)) : exp(Set.epsilon[n].value[i]/kBT));
	// Keep previous rapidities:
	Set.epsilon[n].prev_value[i] = Set.epsilon[n].value[i];
      }
      Tln1pluseminepsinfty[n] = kBT * (Set.epsilon[n].value_infty < 24.0 * kBT
				       ? log(1.0 + exp(-Set.epsilon[n].value_infty/kBT)) : exp(-Set.epsilon[n].value_infty/kBT));
    }

    // Now do the necessary convolutions for epsilon == epsilon[0].
    // For each value of lambda, do the convolutions:
    // Careful:  the lambda's used for lambda (index i) are those of epsilon[0], the lambda' (index j) are for epsilon[n] !!
    Vect<Vect_DP> a_n_Tln_conv(Set.ntypes);
    for (int n = 0; n < Set.ntypes; ++n) {
      a_n_Tln_conv[n] = Vect_DP (0.0, Set.epsilon[0].Npts);
      Vect_DP f(0.0, Set.epsilon[n].Npts);

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

	for (int j = 0; j < Set.epsilon[n].Npts; ++j) a_n_Tln_conv[n][i] += Tln1plusemineps[n][j] * a_n_dlambda[i][n][j];
	// Add alpha curvature terms:  VERY COSTLY, remove for now
	/*
	  for (int j = 1; j < Set.epsilon[n].Npts - 1; ++j)
	  a_n_Tln_conv[n][i] += (1.0/12.0) * pow(Set.epsilon[n].dlambda[j], 3.0)
	  * ((Tln1plusemineps[n][j+1] * a_n_dlambda[i][n][j+1]
	  - Tln1plusemineps[n][j] * a_n_dlambda[i][n][j])/(Set.epsilon[n].lambda[j+1] - Set.epsilon[n].lambda[j])
	  - (Tln1plusemineps[n][j] * a_n_dlambda[i][n][j]
	  - Tln1plusemineps[n][j-1] * a_n_dlambda[i][n][j-1])/(Set.epsilon[n].lambda[j] - Set.epsilon[n].lambda[j-1]))
	  /(Set.epsilon[n].lambda[j+1] - Set.epsilon[n].lambda[j-1]);
	*/
      } // for (int i ...

    } // for (int n...    We now have all the a_n * Tln... at our disposal.

    // For n > nmax sum in RHS of BE for epsilon, assuming epsilon_n = epsilon_n^\infty in those cases:
    // Remember:  nmax = Set.ntypes - 1
    DP Smaxsum = kBT * log((1.0 - exp(-2.0 * (Set.ntypes + 1.0) * Omega/kBT))/(1.0 - exp(-2.0 * Set.ntypes * Omega/kBT)));

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

      // Add the convolutions:
      for (int n = 0; n < Set.ntypes; ++n)
	Set.epsilon[0].value[i] -= a_n_Tln_conv[n][i];

      // Add the asymptotic parts of convolutions:
      for (int n = 1; n < Set.ntypes; ++n)
	Set.epsilon[0].value[i] -= Tln1pluseminepsinfty[n]
	  * (1.0 - (atan((Set.epsilon[n].lambdamax - Set.epsilon[0].lambda[i])/(0.5 * n * c_int))
		    + atan((Set.epsilon[n].lambdamax + Set.epsilon[0].lambda[i])/(0.5 * n * c_int)))/PI);

      // Add the leftover summation for species n > nmax, assuming epsilon_n = epsilon_n^\infty in those cases:
      Set.epsilon[0].value[i] -= Smaxsum;

      // Include some damping:
      Set.epsilon[0].value[i] = 0.1 * Set.epsilon[0].prev_value[i] + 0.9 * Set.epsilon[0].value[i];
      // No need to force boundaries here, epsilon[0] is inherently stable.
    }
    // epsilon[0] is now fully iterated.


    // Now do the remaining epsilons:

    for (int n = 1; n < Set.ntypes; ++n) {

      Vect_DP f_Tln_conv_min (0.0, Set.epsilon[n].Npts);  // 'down' convolution
      Vect_DP f_Tln_conv_plus (0.0, Set.epsilon[n].Npts); // 'up' convolution

      Vect_DP fmin(0.0, Set.epsilon[n-1].Npts);
      Vect_DP fplus(0.0, Set.epsilon[ABACUS::min(n+1, Set.ntypes - 1)].Npts);

      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	f_Tln_conv_min[i] = 0.0;
	f_Tln_conv_plus[i] = 0.0;

	// 'down' convolutions
	if (n == 1)
	  for (int j = 0; j < Set.epsilon[0].Npts; ++j)
	    f_Tln_conv_min[i] += Tln1plusemineps[0][j] * fmin_dlambda[n][i][j];

	else for (int j = 0; j < Set.epsilon[n - 1].Npts; ++j)
	       f_Tln_conv_min[i] += (Set.epsilon[n-1].prev_value[j] - Set.epsilon[n-1].value_infty
				     + Tln1plusemineps[n-1][j] - Tln1pluseminepsinfty[n-1])
		 * fmin_dlambda[n][i][j];

	// 'up' convolutions
	if (n < Set.ntypes - 1)
	  for (int j = 0; j < Set.epsilon[n+1].Npts; ++j)
	    f_Tln_conv_plus[i] += (Set.epsilon[n+1].prev_value[j] - Set.epsilon[n+1].value_infty
				   + Tln1plusemineps[n+1][j] - Tln1pluseminepsinfty[n+1])
	      * fplus_dlambda[n][i][j];

	else f_Tln_conv_plus[i] = 0.0;

	// Do some damping:
	Set.epsilon[n].value[i] = 0.1 * Set.epsilon[n].prev_value[i]
	  + 0.9 * (Set.epsilon[n].value_infty + f_Tln_conv_min[i] + f_Tln_conv_plus[i]);
	// Force boundary values to asymptotes:  force boundary 10 points on each side
	if (i < 10)
	  Set.epsilon[n].value[i] = (1.0 - 0.1 * i) * Set.epsilon[n].value_infty + 0.1 * i * Set.epsilon[n].value[i];
	if (i > Set.epsilon[n].Npts - 11)
	  Set.epsilon[n].value[i] = (1.0 - 0.1 * (Set.epsilon[n].Npts-1 - i)) * Set.epsilon[n].value_infty
	    + 0.1 * (Set.epsilon[n].Npts-1 - i) * Set.epsilon[n].value[i];

      } // for (int i = 0...

    } // for (int n = 1...

    // All functions have now been iterated.

    // Now calculate diff:

    DP eps0i = 0.0;
    DP eps1i = 0.0;

    Set.diff = 0.0;

    for (int n = 0; n < Set.ntypes; ++n) {
      Set.epsilon[n].diff = 0.0;
      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	// Measure based on delta f/delta epsilon:
	if (n == 0)
	  Set.epsilon[n].diff += Set.epsilon[n].dlambda[i] *
	    (Set.epsilon[0].value[i] > 0.0 ?
	     exp(-Set.epsilon[0].value[i]/kBT)/(1.0 + exp(-Set.epsilon[0].value[i]/kBT)) : 1.0/(1.0 + exp(Set.epsilon[0].value[i]/kBT)))
	    * fabs(Set.epsilon[n].value[i] - Set.epsilon[n].prev_value[i]);
	else {
	  eps0i = Set.epsilon[0].Return_Value(Set.epsilon[n].lambda[i]);
	  eps1i = Set.epsilon[1].Return_Value(Set.epsilon[n].lambda[i]);

	  Set.epsilon[n].diff += Set.epsilon[n].dlambda[i] *
	    // Logic:  simple 1/2 cascade
	    (eps0i > 0.0 ? exp(-eps0i/kBT)/(1.0 + exp(-eps0i/kBT)) : 1.0/(1.0 + exp(eps0i/kBT)))
	    * pow(0.5, n) //* (exp(-eps1i/kBT)/(1.0 + exp(-eps1i/kBT)))
	    * fabs(Set.epsilon[n].value[i] - Set.epsilon[n].prev_value[i]);
	}
      }
      Set.diff += Set.epsilon[n].diff;
    }

    return;
  }

  void Iterate_and_Extrapolate_2CBG_TBAE (Root_Density_Set& Set, Vect<Root_Density_Set>& IterSet,
					  Vect<Vect<Vect_DP> >& a_n_dlambda, Vect<Vect<Vect_DP> >& fmin_dlambda,
					  Vect<Vect<Vect_DP> >& fplus_dlambda, DP c_int, DP mu, DP Omega, DP kBT)
  {
    int nfit = IterSet.size();

    for (int ifit = 0; ifit < nfit; ++ifit) {
      Iterate_2CBG_TBAE (Set, a_n_dlambda, fmin_dlambda, fplus_dlambda, c_int, mu, Omega, kBT);
      IterSet[ifit] = Set;
    }

    // Now extrapolate each value to infinite nr of iterations:
    Vect_DP density(nfit);
    Vect_DP oneoverP(nfit);
    DP deltalambda = 0.0;
    for (int ifit = 0; ifit < nfit; ++ifit) oneoverP[ifit] = 1.0/(1.0 + ifit*ifit);
    for (int n = 0; n < Set.ntypes; ++n)
      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	for (int ifit = 0; ifit < nfit; ++ifit) density[ifit] = IterSet[ifit].epsilon[n].value[i];
	polint (oneoverP, density, 0.0, Set.epsilon[n].value[i], deltalambda);
      }

    // Now iterate a few times to stabilize:
    for (int iint = 0; iint < 2; ++iint)
      Iterate_2CBG_TBAE(Set, a_n_dlambda, fmin_dlambda, fplus_dlambda, c_int, mu, Omega, kBT);

    return;
  }

  DP Refine_2CBG_Set (Root_Density_Set& Set, DP c_int, DP mu, DP Omega, 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<Vect_DP> Tln1plusemineps(Set.ntypes);
    Vect_DP Tln1pluseminepsinfty(Set.ntypes);

    for (int n = 0; n < Set.ntypes; ++n) {

      Tln1plusemineps[n] = Vect_DP (0.0, Set.epsilon[n].Npts);

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

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

    Vect<Vect_DP> tni(Set.ntypes);
    Vect<Vect_DP> tni_ex(Set.ntypes);

    DP measure_factor = 0.0;
    DP eps0i = 0.0;
    DP eps1i = 0.0;

    for (int n = 0; n < Set.ntypes; ++n) {

      tni[n] = Vect_DP (0.0, Set.epsilon[n].Npts);
      tni_ex[n] = Vect_DP (0.0, Set.epsilon[n].Npts);  // extrapolation from adjacent points, to compare to obtained value

      for (int i = 1; i < Set.epsilon[n].Npts - 1; ++i) {
	if (n == 0) {
	  // Measure based on delta f/delta epsilon:
	  measure_factor = (Set.epsilon[0].value[i] > 0.0
			    ? exp(-Set.epsilon[0].value[i]/kBT)/(1.0 + exp(-Set.epsilon[0].value[i]/kBT))
			    : 1.0/(1.0 + exp(Set.epsilon[0].value[i]/kBT)));
	}
	else {
	  // Measure based on delta f/delta epsilon:
	  eps0i = Set.epsilon[0].Return_Value(Set.epsilon[n].lambda[i]);
	  eps1i = Set.epsilon[1].Return_Value(Set.epsilon[n].lambda[i]);

	  measure_factor = (eps0i > 0.0 ? exp(-eps0i/kBT)/(1.0 + exp(-eps0i/kBT)) : 1.0/(1.0 + exp(eps0i/kBT)))
	    // Logic:  simple 1/2 per level cascade down
	    * pow(0.5, n);

	}

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

	max_delta_tni_dlambda = ABACUS::max(max_delta_tni_dlambda, fabs(tni[n][i] - tni_ex[n][i]) * Set.epsilon[n].dlambda[i]);
	if (n == Set.ntypes - 1)
	  max_delta_tni_dlambda_toplevel = ABACUS::max(max_delta_tni_dlambda_toplevel, fabs(tni[n][i] - tni_ex[n][i]) * Set.epsilon[n].dlambda[i]);
	sum_delta_tni_dlambda += fabs(tni[n][i] - tni_ex[n][i]) * Set.epsilon[n].dlambda[i];
      }
    }


    // We now determine the locations where we need to add points
    Vect<Vect<bool> > need_new_point_around(Set.ntypes);
    Vect<bool> need_to_extend_limit(false, Set.ntypes);

    for (int n = 0; n < Set.ntypes; ++n) {

      need_new_point_around[n] = Vect<bool> (false, Set.epsilon[n].Npts);

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

      // Check boundary values;  if too different from value_infty, extend limits
      if (n == 0) {
	// Measure based on delta f/delta epsilon:
	if (exp(-Set.epsilon[0].value[0]/kBT) > 0.001 * max_delta_tni_dlambda)
	  need_to_extend_limit[0] = true;
      }
      else
	// Measure deviation from asymptote for 10th element, since we smoothly put the i<10 ones to the asymptote when damping:
	if (fabs(Set.epsilon[n].value[10] - Set.epsilon[n].value_infty) * Set.epsilon[0].dlambda[10] > max_delta_tni_dlambda)
	  need_to_extend_limit[n] = true;
    }

    // Check if we need to add a level
    bool need_new_epsilon_n_function = false;

    // We add new levels if the integral a_n * Tln1plusemineps at the highest level differs too much from
    // the asymptotic value.  Since such integrals appear for each point of the epsilon[0] function, these
    // errors should be compared to the individual delta_tni factors.
    DP a_2_Tln_conv_0_integ = 0.0;
    DP oneoverpi = 1.0/PI;
    DP twoovernc = 2.0/((Set.ntypes - 1) * c_int);
    int i0 = Set.epsilon[0].Npts/2;
    for (int j = 0; j < Set.epsilon[Set.ntypes - 1].Npts; ++j)
      a_2_Tln_conv_0_integ += (Tln1plusemineps[Set.ntypes - 1][j] - Tln1pluseminepsinfty[Set.ntypes - 1])
	* oneoverpi * (atan(twoovernc * (Set.epsilon[0].lambda[i0] - Set.epsilon[Set.ntypes - 1].lambda[j]
					 + 0.5 * Set.epsilon[Set.ntypes - 1].dlambda[j]))
		       - atan(twoovernc * (Set.epsilon[0].lambda[i0] - Set.epsilon[Set.ntypes - 1].lambda[j]
					   - 0.5 * Set.epsilon[Set.ntypes - 1].dlambda[j])));

    // Compare to prediction for this integral based on value_infty, which is simply 0.
    // Count this difference Set.ntypes times over, since it cascades down all levels
    if (fabs(a_2_Tln_conv_0_integ) * Set.ntypes > max_delta_tni_dlambda) need_new_epsilon_n_function = true;

    // Additionally, if the highest level needs updating, we automatically add new functions:
    for (int i = 0; i < Set.epsilon[Set.ntypes - 1].Npts; ++i)
      if (need_new_point_around[Set.ntypes - 1][i] || need_to_extend_limit[Set.ntypes - 1]) need_new_epsilon_n_function = true;


    // Now insert the new points between existing points:
    Set.Insert_new_points (need_new_point_around);

    // Now extend the integration limits if needed:
    Set.Extend_limits (need_to_extend_limit);

    // If we add a level, we do it here:
    if (need_new_epsilon_n_function) {
      // Insert more than one function per cycle...
      Set.Insert_new_function(Asymptotic_2CBG_epsilon(Set.ntypes, Omega, kBT));
      Set.Insert_new_function(Asymptotic_2CBG_epsilon(Set.ntypes, Omega, kBT));  // CAREFUL !!  ntypes is already updated
      Set.Insert_new_function(Asymptotic_2CBG_epsilon(Set.ntypes, Omega, kBT));  // CAREFUL !!  ntypes is already updated
      Set.Insert_new_function(Asymptotic_2CBG_epsilon(Set.ntypes, Omega, kBT));  // CAREFUL !!  ntypes is already updated;
      Set.Insert_new_function(Asymptotic_2CBG_epsilon(Set.ntypes, Omega, kBT));  // CAREFUL !!  ntypes is already updated
    }

    //return(max_delta_tni_dlambda);
    return(sum_delta_tni_dlambda);
  }

  DP Calculate_Gibbs_Free_Energy (const Root_Density_Set& Set, DP c_int, DP mu, DP Omega, DP kBT)
  {
    // Computes the Gibbs free energy, assuming that epsilon[0] is symmetric.

    // WORKING VERSION
    DP sum_f = 0.0;
    Vect_DP f(0.0, Set.epsilon[0].Npts);
    DP sum_f_check = 0.0;
    Vect_DP fcheck(0.0, Set.epsilon[0].Npts);
    DP sum_g_check = 0.0;
    Vect_DP gcheck(0.0, Set.epsilon[0].Npts);
    for (int i = 0; i < Set.epsilon[0].Npts; ++i) {
      f[i] = (Set.epsilon[0].value[i] > 0.0 ? kBT* log(1.0 + exp(-Set.epsilon[0].value[i]/kBT))
	      : -Set.epsilon[0].value[i] + kBT * log(1.0 + exp(Set.epsilon[0].value[i]/kBT)));
      sum_f += Set.epsilon[0].dlambda[i] * f[i];
      fcheck[i] = kBT * log(1.0 + exp(-(Set.epsilon[0].lambda[i] * Set.epsilon[0].lambda[i] - mu - Omega)/kBT));
      sum_f_check += Set.epsilon[0].dlambda[i] * fcheck[i];
      gcheck[i] = exp(-(Set.epsilon[0].lambda[i] * Set.epsilon[0].lambda[i])/kBT);
      sum_g_check += Set.epsilon[0].dlambda[i] * gcheck[i];
    }

    // Now add alpha curvature terms:
    for (int i = 1; i < Set.epsilon[0].Npts - 1; ++i)
      sum_f += pow(Set.epsilon[0].dlambda[i], 3.0) * ((f[i+1] - f[i])/(Set.epsilon[0].lambda[i+1] - Set.epsilon[0].lambda[i])
						      - (f[i] - f[i-1])/(Set.epsilon[0].lambda[i] - Set.epsilon[0].lambda[i-1]))
	/(12.0 * (Set.epsilon[0].lambda[i+1] - Set.epsilon[0].lambda[i-1]));

    // Now add alpha curvature terms:
    DP sum_gcorralphacheck = 0.0;
    Vect_DP gcorr_alpha_check(0.0, Set.epsilon[0].Npts);
    for (int i = 1; i < Set.epsilon[0].Npts - 1; ++i) {
      gcorr_alpha_check[i] = (1.0/12.0) * pow(Set.epsilon[0].dlambda[i], 3.0)
	* ((gcheck[i+1] - gcheck[i]) * (Set.epsilon[0].lambda[i] - Set.epsilon[0].lambda[i-1])
	   - (gcheck[i] - gcheck[i-1]) * (Set.epsilon[0].lambda[i+1] - Set.epsilon[0].lambda[i]))
	/((Set.epsilon[0].lambda[i+1] - Set.epsilon[0].lambda[i-1]) * (Set.epsilon[0].lambda[i+1] - Set.epsilon[0].lambda[i])
	  * (Set.epsilon[0].lambda[i] - Set.epsilon[0].lambda[i-1]));
      sum_gcorralphacheck += gcorr_alpha_check[i];
    }

    // Testing:
    int Npts_test = Set.epsilon[0].Npts;
    DP lambdamax_test = Set.epsilon[0].lambdamax;
    DP testsum = 0.0;
    DP testgauss = 0.0;
    DP lambda;
    DP dlambda;
    for (int i = 0; i < Npts_test; ++i) {
      lambda = lambdamax_test * (-Npts_test + 1.0 + 2*i)/Npts_test;
      dlambda = lambdamax_test * 2.0/Npts_test;
      testsum += kBT * log(1.0 + exp(-(lambda*lambda - mu - Omega)/kBT)) * dlambda;
      testgauss += exp(-lambda*lambda/kBT) * dlambda;
    }

    return(-sum_f/twoPI);
  }

  DP Calculate_dGibbs_dchempot (const Root_Density_Set& DSet, const Root_Density_Set& Set, DP c_int, DP mu, DP Omega, DP kBT)
  {
    // This calculates the derivative of the Gibbs free energy with respect to either of the two chemical potential,
    // given the fundametal set Set for eps and DSet for either deps_du or deps_dOmega.

    DP sum_f = 0.0;
    Vect_DP f(0.0, Set.epsilon[0].Npts);

    for (int i = 0; i < Set.epsilon[0].Npts; ++i) {
      f[i] = DSet.epsilon[0].value[i] * (Set.epsilon[0].value[i] > 0.0
					 ? exp(-Set.epsilon[0].value[i]/kBT)/(1.0 + exp(-Set.epsilon[0].value[i]/kBT))
					 : 1.0/(1.0 + exp(Set.epsilon[0].value[i]/kBT)));
      sum_f += DSet.epsilon[0].dlambda[i] * f[i];
    }

    // Now add alpha curvature terms:
    for (int i = 1; i < DSet.epsilon[0].Npts - 1; ++i)
      sum_f += pow(DSet.epsilon[0].dlambda[i], 3.0)
	* ((f[i+1] - f[i])/(DSet.epsilon[0].lambda[i+1] - DSet.epsilon[0].lambda[i])
	   - (f[i] - f[i-1])/(DSet.epsilon[0].lambda[i] - DSet.epsilon[0].lambda[i-1]))
	/(12.0 * (DSet.epsilon[0].lambda[i+1] - DSet.epsilon[0].lambda[i-1]));

    return(sum_f/twoPI);
  }

  Vect<Vect<Vect_DP> > Build_a_n_dlambda (const Root_Density_Set& Set, DP c_int, DP mu, DP Omega, DP kBT)
  {
    DP oneoverpi = 1.0/PI;
    DP oneoverc = 1.0/c_int;
    DP twoovernc = 2.0/c_int;

    Vect<Vect<Vect_DP> > a_n_dlambda(Set.epsilon[0].Npts);
    for (int i = 0; i < Set.epsilon[0].Npts; ++i) {
      a_n_dlambda[i] = Vect<Vect_DP>(Set.ntypes);
      for (int n = 0; n < Set.ntypes; ++n) {
	a_n_dlambda[i][n] = Vect_DP(0.0, Set.epsilon[n].Npts);
      }
    }

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

      // Do n == 0 separately:
      for (int j = 0; j < Set.epsilon[0].Npts; ++j)
	a_n_dlambda[i][0][j] = oneoverpi
	  * (atan(oneoverc * (Set.epsilon[0].lambda[i] - Set.epsilon[0].lambda[j] + 0.5 * Set.epsilon[0].dlambda[j]))
	     - atan(oneoverc * (Set.epsilon[0].lambda[i] - Set.epsilon[0].lambda[j] - 0.5 * Set.epsilon[0].dlambda[j])));

      // Now do n > 0:
      for (int n = 1; n < Set.ntypes; ++n) {
	twoovernc = 2.0/(n * c_int);
	for (int j = 0; j < Set.epsilon[n].Npts; ++j) {
	  a_n_dlambda[i][n][j] = oneoverpi
	    * (atan(twoovernc * (Set.epsilon[0].lambda[i] - Set.epsilon[n].lambda[j] + 0.5 * Set.epsilon[n].dlambda[j]))
	       - atan(twoovernc * (Set.epsilon[0].lambda[i] - Set.epsilon[n].lambda[j] - 0.5 * Set.epsilon[n].dlambda[j])));
	}
      } // for (int n
    } // for (int i

    return(a_n_dlambda);
  }

  Vect<Vect<Vect_DP> >  Build_fmin_dlambda (const Root_Density_Set& Set, DP c_int, DP mu, DP Omega, DP kBT)
  {
    DP oneoverpi = 1.0/PI;
    DP pioverc = PI/c_int;
    DP twopioverc = 2.0*PI/c_int;
    DP piovertwoc = 0.5 * pioverc;

    Vect<Vect<Vect_DP> > fmin_dlambda(Set.ntypes);
    for (int n = 0; n < Set.ntypes; ++n) {
      fmin_dlambda[n] = Vect<Vect_DP> (Set.epsilon[n].Npts);
      for (int i = 0; i < Set.epsilon[n].Npts; ++i)
	fmin_dlambda[n][i] = Vect_DP (0.0, Set.epsilon[ABACUS::max(n-1, 0)].Npts);
    }

    for (int n = 1; n < Set.ntypes; ++n) {
      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {

	for (int j = 0; j < Set.epsilon[n-1].Npts; ++j)
	  fmin_dlambda[n][i][j] = oneoverpi
	    * atan(exp(-pioverc * fabs(Set.epsilon[n].lambda[i] - Set.epsilon[n-1].lambda[j]))
		   * 2.0 * sinh(piovertwoc * Set.epsilon[n-1].dlambda[j])
		   /(1.0 + exp(-twopioverc * fabs(Set.epsilon[n].lambda[i] - Set.epsilon[n-1].lambda[j]))));
      } // for i
    } // for n

    return(fmin_dlambda);
  }

  Vect<Vect<Vect_DP> > Build_fplus_dlambda (const Root_Density_Set& Set, DP c_int, DP mu, DP Omega, DP kBT)
  {
    DP oneoverpi = 1.0/PI;
    DP pioverc = PI/c_int;
    DP twopioverc = 2.0*PI/c_int;
    DP piovertwoc = 0.5 * pioverc;

    Vect<Vect<Vect_DP> > fplus_dlambda(Set.ntypes);
    for (int n = 0; n < Set.ntypes; ++n) {
      fplus_dlambda[n] = Vect<Vect_DP> (Set.epsilon[n].Npts);
      for (int i = 0; i < Set.epsilon[n].Npts; ++i)
	fplus_dlambda[n][i] = Vect_DP (0.0, Set.epsilon[ABACUS::min(n+1, Set.ntypes - 1)].Npts);
    }

    for (int n = 0; n < Set.ntypes - 1; ++n) {
      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	for (int j = 0; j < Set.epsilon[n+1].Npts; ++j)
	  fplus_dlambda[n][i][j] = oneoverpi
	    * atan(exp(-pioverc * fabs(Set.epsilon[n].lambda[i] - Set.epsilon[n+1].lambda[j]))
		   * 2.0 * sinh(piovertwoc * Set.epsilon[n+1].dlambda[j])
		   /(1.0 + exp(-twopioverc * fabs(Set.epsilon[n].lambda[i] - Set.epsilon[n+1].lambda[j]))));
      }
    }

    return(fplus_dlambda);
  }


  Root_Density_Set Solve_2CBG_TBAE_via_refinements (DP c_int, DP mu, DP Omega, DP kBT, int Max_Secs,
						    ofstream& LOG_outfile, bool Save_data)
  {
    // This solves the 2CBG TBAE as best as possible given the time constraint.

    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

    // Set basic number of types needed:
    int ntypes_needed = int(kBT * log(kBT/1.0e-14)/Omega);
    int ntypes = ABACUS::max(ntypes_needed, 1);
    ntypes = ABACUS::min(ntypes, 10);

    if (Save_data)
      if (ntypes >= 10) LOG_outfile << "WARNING:  ntypes needs to be quite high for c_int = " << c_int
				    << " mu = " << mu << " Omega = " << Omega
				    << " kBT = " << kBT << ".  Set to " << ntypes << ", ideally needed: "
				    << ntypes_needed << ".  Accuracy might be incorrectly evaluated." << endl;

    DP lambdamax = 10.0 + sqrt(ABACUS::max(1.0, kBT * 36.0 + mu + Omega));
    // such that exp(-(lambdamax^2 - mu - Omega)/T) <~ machine_eps
    int Npts = 50;
    Vect_INT Npts_init(Npts, ntypes);
    if (Save_data) LOG_outfile << "Npts (basic) set to " << Npts_init << endl;

    Vect_DP lambdamax_init(lambdamax, ntypes);  // such that exp(-pi *lambdamax/c) <~ machine_eps
    Npts_init[0] = 1 * Npts;  // give more precision to lowest level
    lambdamax_init[0] = 10.0 + sqrt(ABACUS::max(1.0, kBT * 36.0 + mu + Omega));
    // such that exp(-(lambdamax^2 - mu - Omega)/T) <~ machine_eps
    Root_Density_Set TBA_Set (ntypes, Npts_init, lambdamax_init);

    // Set the asymptotics of the TBA_fns:
    Set_2CBG_Asymptotics (TBA_Set, mu, Omega, kBT);

    // Initiate the functions:
    Initiate_2CBG_TBA_Functions (TBA_Set, mu, Omega);

    clock_t StopTime = clock();
    clock_t Cycle_StartTime, Cycle_StopTime;

    int CPU_ticks = StopTime - StartTime;

    int ncycles = 0;
    int niter_tot = 0;

    do {

      StartTime = clock();

      Cycle_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...
      running_prec = Refine_2CBG_Set (TBA_Set, c_int, mu, Omega, kBT, refine_fraction);

      Vect<Vect<Vect_DP> > a_n_dlambda = Build_a_n_dlambda (TBA_Set, c_int, mu, Omega, kBT);
      Vect<Vect<Vect_DP> > fmin_dlambda = Build_fmin_dlambda (TBA_Set, c_int, mu, Omega, kBT);
      Vect<Vect<Vect_DP> > fplus_dlambda = Build_fplus_dlambda (TBA_Set, c_int, mu, Omega, kBT);

      StopTime = clock();

      CPU_ticks += StopTime - StartTime;

      int niter = 0;
      int niter_max = ncycles == 0 ? 300 : 300;

      // For extrapolations:
      Vect<Root_Density_Set> IterSet(4);

      do {

	StartTime = clock();

	if (niter <= 10 || niter > 100) {
	  Iterate_2CBG_TBAE (TBA_Set, a_n_dlambda, fmin_dlambda, fplus_dlambda, c_int, mu, Omega, kBT);
	  niter++;
	}
	else {
	  Iterate_and_Extrapolate_2CBG_TBAE (TBA_Set, IterSet, a_n_dlambda, fmin_dlambda, fplus_dlambda, c_int, mu, Omega, kBT);
	  niter += 6;
	}

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

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

      if (Save_data) {
	LOG_outfile << "ncycles = " << ncycles << "\trunning_prec = " << running_prec << "\t niter = " << niter
		    << "\tntypes = " << TBA_Set.ntypes << "\tNpts ";
	for (int n = 0; n < TBA_Set.ntypes; ++n) LOG_outfile << TBA_Set.epsilon[n].Npts << " ";
	LOG_outfile << "\tNpts_total = " << TBA_Set.Npts_total << endl
		    << "\tdiff = " << TBA_Set.diff
		    << "\tGSE = " << Calculate_Gibbs_Free_Energy (TBA_Set, c_int, mu, Omega, kBT) << endl;
      }

      ncycles++;
      niter_tot += niter;

      if (niter == niter_max) {
	if (Save_data) LOG_outfile << "Not able to improve functions enough after " << niter_max << " iterations." << endl;
      }

      Cycle_StopTime = clock();

    } while (CPU_ticks < Max_CPU_ticks - 2.0*(Cycle_StopTime - Cycle_StartTime));
    // Allow a new cycle only if there is time, assuming new cycle time < 2* last one

    if (Save_data) {
      LOG_outfile << "c_int " << c_int << "\tmu " << mu << "\tOmega " << Omega << "\tkBT " << kBT
		  << "\tncycles = " << ncycles << "\trunning_prec = " << running_prec << "\t niter_tot = " << niter_tot
		  << "\tntypes = " << TBA_Set.ntypes << "\tdiff = " << TBA_Set.diff << endl << "\tNpts ";
      for (int n = 0; n < TBA_Set.ntypes; ++n) LOG_outfile << TBA_Set.epsilon[n].Npts << " ";
      LOG_outfile << "\tNpts_total = " << TBA_Set.Npts_total << endl;
    }

    return(TBA_Set);
    //return;
  }


  // Iterative procedures for deps/dmu or /dOmega:
  void Iterate_2CBG_deps_dchempot (int option, Root_Density_Set& DSet, const Root_Density_Set& Set,
				   Vect<Vect<Vect_DP> >& a_n_dlambda, Vect<Vect<Vect_DP> >& fmin_dlambda,
				   Vect<Vect<Vect_DP> >& fplus_dlambda, DP c_int, DP mu, DP Omega, DP kBT)
  {
    // Produces a new Root_Density_Set for depsilon/dmu (option == 0) or depsilon/dOmega (option == 1) from a previous iteration.
    // Does NOT add types or change Npts, lambdamax values.

    // First define some useful functions:
    Vect<Vect_DP> depsover1plusemineps(Set.ntypes);
    Vect<Vect_DP> depsover1plusepluseps(Set.ntypes);
    Vect_DP depsover1pluseminepsinfty(Set.ntypes);
    Vect_DP depsover1pluseplusepsinfty(Set.ntypes);

    for (int n = 0; n < Set.ntypes; ++n) {

      depsover1plusemineps[n] = Vect_DP (0.0, Set.epsilon[n].Npts);
      depsover1plusepluseps[n] = Vect_DP (0.0, Set.epsilon[n].Npts);

      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	depsover1plusemineps[n][i] = Set.epsilon[n].value[i] > 0.0 ?
	  DSet.epsilon[n].value[i]/(1.0 + exp(-Set.epsilon[n].value[i]/kBT)) :
	  DSet.epsilon[n].value[i] * exp(Set.epsilon[n].value[i]/kBT)/(1.0 + exp(Set.epsilon[n].value[i]/kBT));
	depsover1plusepluseps[n][i] = Set.epsilon[n].value[i] > 0.0 ?
	  DSet.epsilon[n].value[i] * exp(-Set.epsilon[n].value[i]/kBT)/(1.0 + exp(-Set.epsilon[n].value[i]/kBT)) :
	  DSet.epsilon[n].value[i]/(1.0 + exp(Set.epsilon[n].value[i]/kBT));

	// Keep previous rapidities:
	DSet.epsilon[n].prev_value[i] = DSet.epsilon[n].value[i];

      }
      depsover1pluseminepsinfty[n] = DSet.epsilon[n].value_infty/(1.0 + exp(-Set.epsilon[n].value_infty/kBT));
      depsover1pluseplusepsinfty[n] = DSet.epsilon[n].value_infty * exp(-Set.epsilon[n].value_infty/kBT)
	/(1.0 + exp(-Set.epsilon[n].value_infty/kBT));
    }


    // Now do the necessary convolutions for epsilon == epsilon[0].
    // For each value of lambda, do the convolutions:
    // Careful:  the lambda's used for lambda (index i) are those of epsilon[0], the lambda' (index j) are for epsilon[n] !!
    Vect<Vect_DP> a_n_depsover_conv(Set.ntypes);
    for (int n = 0; n < Set.ntypes; ++n) {
      a_n_depsover_conv[n] = Vect_DP (0.0, Set.epsilon[0].Npts);
      Vect_DP f(0.0, Set.epsilon[n].Npts);

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

	for (int j = 0; j < Set.epsilon[n].Npts; ++j) {
	  f[j] = depsover1plusepluseps[n][j] * a_n_dlambda[i][n][j];
	  a_n_depsover_conv[n][i] += f[j];
	}
      }
    } // for (int n...    We now have all the a_n * deps... at our disposal.

    // For n > nmax sum in RHS of BE for epsilon, assuming epsilon_n = epsilon_n^\infty in those cases:
    // Remember: nmax = Set.ntypes - 1
    DP Smaxsum = option == 0 ? 0.0 : 2.0 * ((Set.ntypes + 1.0) * exp(-2.0 * (Set.ntypes + 1.0) * Omega/kBT)
					    /(1.0 - exp(-2.0 * (Set.ntypes + 1.0) * Omega/kBT))
					    - Set.ntypes * exp(-2.0 * Set.ntypes * Omega/kBT)
					    /(1.0 - exp(-2.0 * Set.ntypes * Omega/kBT)));

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

      // Add the convolutions:
      for (int n = 0; n < Set.ntypes; ++n)
	DSet.epsilon[0].value[i] += a_n_depsover_conv[n][i];
      // Add the asymptotic parts of convolutions:  n == 0 part is zero because of 1 + exp[epsilon[0] ] in denominator
      for (int n = 1; n < Set.ntypes; ++n)
	DSet.epsilon[0].value[i] += depsover1pluseplusepsinfty[n]
	  * (1.0 - (atan((DSet.epsilon[n].lambdamax - DSet.epsilon[0].lambda[i])/(0.5 * n * c_int))
		    + atan((DSet.epsilon[n].lambdamax + DSet.epsilon[0].lambda[i])/(0.5 * n * c_int)))/PI);

      // Add the leftover summation for species n > nmax, assuming epsilon_n = epsilon_n^\infty in those cases:
      DSet.epsilon[0].value[i] -= Smaxsum;
      // Include some damping:
      DSet.epsilon[0].value[i] = 0.1 * DSet.epsilon[0].prev_value[i] + 0.9 * DSet.epsilon[0].value[i];
      // Force boundary values to asymptotes:  force boundary 10 points on each side
      if (i < 10)
	DSet.epsilon[0].value[i] = (1.0 - 0.1 * i) * DSet.epsilon[0].value_infty + 0.1 * i * DSet.epsilon[0].value[i];
      if (i > DSet.epsilon[0].Npts - 11)
	DSet.epsilon[0].value[i] = (1.0 - 0.1 * (DSet.epsilon[0].Npts-1 - i)) * DSet.epsilon[0].value_infty
	  + 0.1 * (DSet.epsilon[0].Npts-1 - i) * DSet.epsilon[0].value[i];

    }
    // epsilon[0] is now fully iterated.

    // Now do the remaining epsilons:

    for (int n = 1; n < Set.ntypes; ++n) {

      Vect_DP f_depsover_conv_min (0.0, Set.epsilon[n].Npts);  // 'down' convolution
      Vect_DP f_depsover_conv_plus (0.0, Set.epsilon[n].Npts); // 'up' convolution

      Vect_DP fmin(0.0, Set.epsilon[n-1].Npts);
      Vect_DP fplus(0.0, Set.epsilon[ABACUS::min(n+1, Set.ntypes - 1)].Npts);

      for (int i = 0; i < Set.epsilon[n].Npts; ++i) {
	f_depsover_conv_min[i] = 0.0;
	f_depsover_conv_plus[i] = 0.0;

	// 'down' convolutions
	if (n == 1) {

	  for (int j = 0; j < Set.epsilon[0].Npts; ++j) {
	    fmin[j] = depsover1plusepluseps[0][j]
	      * fmin_dlambda[n][i][j];
	    f_depsover_conv_min[i] -= fmin[j]; // Careful ! - sign here
	  }
	} // if (n == 1)

	else { // if (n != 1)

	  for (int j = 0; j < Set.epsilon[n - 1].Npts; ++j) {
	    fmin[j] = (depsover1plusemineps[n-1][j] - depsover1pluseminepsinfty[n-1])
	      * fmin_dlambda[n][i][j];
	    f_depsover_conv_min[i] += fmin[j];
	  }
	} // if (n != 1)


	// 'up' convolutions
	if (n < Set.ntypes - 1) {

	  for (int j = 0; j < Set.epsilon[n+1].Npts; ++j) {
	    fplus[j] = (depsover1plusemineps[n+1][j] - depsover1pluseminepsinfty[n+1])
	      * fplus_dlambda[n][i][j];
	    f_depsover_conv_plus[i] += fplus[j];
	  }

	} // if (n < Set.ntypes - 1...

	// otherwise, we put the function to 1/2 times depsover... of epsilon[n+1] at infinity, minus the same:
	else f_depsover_conv_plus[i] = 0.0;

	// Do some damping:
	DSet.epsilon[n].value[i] = 0.1 * DSet.epsilon[n].prev_value[i]
	  + 0.9 * (DSet.epsilon[n].value_infty + f_depsover_conv_min[i] + f_depsover_conv_plus[i]);
	// Force boundary values to asymptotes:  force boundary 10 points on each side
	if (i < 10)
	  DSet.epsilon[n].value[i] = (1.0 - 0.1 * i) * DSet.epsilon[n].value_infty + 0.1 * i * DSet.epsilon[n].value[i];
	if (i > DSet.epsilon[n].Npts - 11)
	  DSet.epsilon[n].value[i] = (1.0 - 0.1 * (DSet.epsilon[n].Npts-1 - i)) * DSet.epsilon[n].value_infty
	    + 0.1 * (DSet.epsilon[n].Npts-1 - i) * DSet.epsilon[n].value[i];

      } // for (int i = 0...

    } // for (int n = 1...

    // All functions have now been iterated.

    // Now calculate diff:

    DSet.diff = 0.0;

    for (int n = 0; n < DSet.ntypes; ++n) {
      DSet.epsilon[n].diff = 0.0;
      for (int i = 0; i < DSet.epsilon[n].Npts; ++i) {
	DSet.epsilon[n].diff += DSet.epsilon[n].dlambda[i] *
	  fabs((DSet.epsilon[n].value[i] - DSet.epsilon[n].prev_value[i]) * depsover1plusepluseps[n][i]);
      }
      DSet.diff += DSet.epsilon[n].diff;
    }

    return;
  }

  Root_Density_Set Solve_2CBG_deps_dchempot (int option, Root_Density_Set& TBA_Set, DP c_int, DP mu, DP Omega, DP kBT,
					     int Max_Secs, ofstream& LOG_outfile, bool Save_data)
  {
    // This solves the 2CBG deps/dmu (option == 0) or deps/dOmega (option == 1).

    clock_t StartTime, StopTime;

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

    // Set basic precision needed:
    DP running_prec = TBA_Set.diff;

    // We start by converging simplified sets, with fewer points:

    Root_Density_Set TBA_Set_comp = TBA_Set.Return_Compressed_and_Matched_Set(1.0);

    Root_Density_Set DSet = TBA_Set;  // this will be the final function return
    Root_Density_Set DSet_comp = TBA_Set_comp;

    // Set the asymptotics of the TBA_fns:
    Set_2CBG_deps_dchempot_Asymptotics (option, TBA_Set_comp, DSet_comp, mu, Omega, kBT);
    Set_2CBG_deps_dchempot_Asymptotics (option, TBA_Set, DSet, mu, Omega, kBT);

    // Now, start by 'converging' the comp set:

    // Initiate the functions:
    Initiate_2CBG_deps_dchempot_Functions (DSet_comp);

    Vect<Vect<Vect_DP> > a_n_dlambda_comp = Build_a_n_dlambda (TBA_Set_comp, c_int, mu, Omega, kBT);
    Vect<Vect<Vect_DP> > fmin_dlambda_comp = Build_fmin_dlambda (TBA_Set_comp, c_int, mu, Omega, kBT);
    Vect<Vect<Vect_DP> > fplus_dlambda_comp = Build_fplus_dlambda (TBA_Set_comp, c_int, mu, Omega, kBT);

    int CPU_ticks = 0;

    int niter_comp = 0;

    do {

      StartTime = clock();
      Iterate_2CBG_deps_dchempot (option, DSet_comp, TBA_Set_comp,
				  a_n_dlambda_comp, fmin_dlambda_comp, fplus_dlambda_comp, c_int, mu, Omega, kBT);
      niter_comp++;
      StopTime = clock();
      CPU_ticks += StopTime - StartTime;

    } while (CPU_ticks < Max_CPU_ticks/2 && (niter_comp < 100 || DSet_comp.diff > running_prec));
    // use at most half the time, keep rest for later.

    DSet.Match_Densities(DSet_comp);

    Vect<Vect<Vect_DP> > a_n_dlambda = Build_a_n_dlambda (TBA_Set, c_int, mu, Omega, kBT);
    Vect<Vect<Vect_DP> > fmin_dlambda = Build_fmin_dlambda (TBA_Set, c_int, mu, Omega, kBT);
    Vect<Vect<Vect_DP> > fplus_dlambda = Build_fplus_dlambda (TBA_Set, c_int, mu, Omega, kBT);

    int niter = 0;

    do {
      StartTime = clock();
      Iterate_2CBG_deps_dchempot (option, DSet, TBA_Set, a_n_dlambda, fmin_dlambda, fplus_dlambda, c_int, mu, Omega, kBT);
      niter++;
      StopTime = clock();
      CPU_ticks += StopTime - StartTime;
    } while (CPU_ticks < Max_CPU_ticks && DSet_comp.diff > 1.0e-4 * running_prec);


    // We're done !!

    if (Save_data) {
      LOG_outfile << "deps_dchempot, option == " << option << endl;
      LOG_outfile << "c_int " << c_int << "\tmu " << mu << "\tOmega " << Omega << "\tkBT " << kBT
		  << "\trunning_prec " << running_prec << " niter_comp " << niter_comp << " niter_full " << niter
		  << "\tntypes " << DSet.ntypes << "\tdiff " << DSet.diff << endl;
    }

    return(DSet);
  }


  TBA_Data_2CBG Solve_2CBG_TBAE_via_refinements (DP c_int, DP mu, DP Omega, DP kBT, int Max_Secs, bool Save_data)
  {
    // This solves the 2CBG TBAE as best as possible given the time constraint.

    stringstream TBA_stringstream;
    string TBA_string;
    stringstream Dmu_stringstream;
    string Dmu_string;
    stringstream DOmega_stringstream;
    string DOmega_string;
    stringstream LOG_stringstream;
    string LOG_string;
    stringstream GFE_stringstream;
    string GFE_string;

    TBA_stringstream << "EPS_2CBG_c_" << c_int << "_mu_" << mu << "_Omega_" << Omega
		     << "_kBT_" << kBT << "_Secs_" << Max_Secs << ".dat";
    Dmu_stringstream << "EPS_2CBG_c_" << c_int << "_mu_" << mu << "_Omega_" << Omega
		     << "_kBT_" << kBT << "_Secs_" << Max_Secs << ".dmu";
    DOmega_stringstream << "EPS_2CBG_c_" << c_int << "_mu_" << mu << "_Omega_" << Omega
			<< "_kBT_" << kBT << "_Secs_" << Max_Secs << ".dOm";
    LOG_stringstream << "EPS_2CBG_c_" << c_int << "_mu_" << mu << "_Omega_" << Omega
		     << "_kBT_" << kBT << "_Secs_" << Max_Secs << ".log";
    GFE_stringstream << "GFE_2CBG_c_" << c_int << "_mu_" << mu << "_Omega_" << Omega
		     << "_kBT_" << kBT << "_Secs_" << Max_Secs << ".dat";

    TBA_string = TBA_stringstream.str();
    const char* TBA_Cstr = TBA_string.c_str();

    Dmu_string = Dmu_stringstream.str();
    const char* Dmu_Cstr = Dmu_string.c_str();

    DOmega_string = DOmega_stringstream.str();
    const char* DOmega_Cstr = DOmega_string.c_str();

    LOG_string = LOG_stringstream.str();
    const char* LOG_Cstr = LOG_string.c_str();

    GFE_string = GFE_stringstream.str();
    const char* GFE_Cstr = GFE_string.c_str();

    ofstream LOG_outfile;
    ofstream GFE_outfile;

    if (Save_data) {
      LOG_outfile.open(LOG_Cstr);
      LOG_outfile.precision(6);

      GFE_outfile.open(GFE_Cstr);
      GFE_outfile.precision(16);
    }

    Root_Density_Set TBA_Set = Solve_2CBG_TBAE_via_refinements (c_int, mu, Omega, kBT, Max_Secs/3, LOG_outfile, Save_data);

    if (Save_data) {
      // Output the functions:
      TBA_Set.Save(TBA_Cstr);
    }

    Root_Density_Set DSet_dmu =
      Solve_2CBG_deps_dchempot (0, TBA_Set, c_int, mu, Omega, kBT, Max_Secs/3, LOG_outfile, Save_data);
    if (Save_data) DSet_dmu.Save(Dmu_Cstr);

    Root_Density_Set DSet_dOmega =
      Solve_2CBG_deps_dchempot (1, TBA_Set, c_int, mu, Omega, kBT, Max_Secs/3, LOG_outfile, Save_data);
    if (Save_data) DSet_dOmega.Save(DOmega_Cstr);

    DP f = Calculate_Gibbs_Free_Energy (TBA_Set, c_int, mu, Omega, kBT);
    DP dfdmu = Calculate_dGibbs_dchempot (DSet_dmu, TBA_Set, c_int, mu, Omega, kBT);
    DP dfdOm = Calculate_dGibbs_dchempot (DSet_dOmega, TBA_Set, c_int, mu, Omega, kBT);

    if (Save_data)
      GFE_outfile << f << "\t" << TBA_Set.diff
		  << "\t" << dfdmu << "\t" << dfdOm << "\t" << 0.5 * (-dfdmu - dfdOm) << "\t" << 0.5 * (-dfdmu + dfdOm)
		  << endl;

    else cout << setprecision(16) << f << "\t" << 0.5 * (-dfdmu - dfdOm) << "\t" << 0.5 * (-dfdmu + dfdOm);

    if (Save_data) {
      LOG_outfile.close();
      GFE_outfile.close();
    }

    TBA_Data_2CBG answer;
    answer.c_int = c_int;
    answer.mu = mu;
    answer.Omega = Omega;
    answer.kBT = kBT;
    answer.f = f;
    answer.n1 = 0.5 * (-dfdmu - dfdOm);
    answer.n2 = 0.5 * (-dfdmu + dfdOm);

    return(answer);
  }

  void GFE_muscan_2CBG (DP c_int, DP mu_min, DP mu_max, DP Omega, DP kBT, int Npts_mu, int Max_Secs, bool Save_data)
  {
    DP dmu = (mu_max - mu_min)/(Npts_mu - 1);

    stringstream LOG_stringstream;
    string LOG_string;
    stringstream GFE_stringstream;
    string GFE_string;

    LOG_stringstream << "GFE_2CBG_c_" << c_int << "_mu_min_" << mu_min << "_mu_max_" << mu_max
		     << "_Npts_mu_" << Npts_mu << "_Omega_" << Omega << "_kBT_" << kBT << "_Secs_" << Max_Secs << ".log";
    GFE_stringstream << "GFE_2CBG_c_" << c_int << "_mu_min_" << mu_min << "_mu_max_" << mu_max
		     << "_Npts_mu_" << Npts_mu << "_Omega_" << Omega << "_kBT_" << kBT << "_Secs_" << Max_Secs << ".dat";

    LOG_string = LOG_stringstream.str();
    const char* LOG_Cstr = LOG_string.c_str();

    GFE_string = GFE_stringstream.str();
    const char* GFE_Cstr = GFE_string.c_str();

    ofstream LOG_outfile;
    ofstream GFE_outfile;

    LOG_outfile.open(LOG_Cstr);
    LOG_outfile.precision(6);

    GFE_outfile.open(GFE_Cstr);
    GFE_outfile.precision(16);


    Root_Density_Set Scan_Set (10, 10, 10.0);;
    Root_Density_Set Scan_dSet_dmu (10, 10, 10.0);;
    Root_Density_Set Scan_dSet_dOmega (10, 10, 10.0);;
    DP mu;
    int Max_Secs_per_mu_pt = Max_Secs/Npts_mu;

    for (int imu = 0; imu < Npts_mu; ++imu) {
      mu = mu_min + imu * dmu;
      Scan_Set = Solve_2CBG_TBAE_via_refinements (c_int, mu, Omega, kBT, Max_Secs_per_mu_pt/3, LOG_outfile, Save_data);

      Scan_dSet_dmu = Solve_2CBG_deps_dchempot (0, Scan_Set, c_int, mu, Omega, kBT, Max_Secs_per_mu_pt/3,
						LOG_outfile, Save_data);

      Scan_dSet_dOmega = Solve_2CBG_deps_dchempot (1, Scan_Set, c_int, mu, Omega, kBT, Max_Secs_per_mu_pt/3,
						   LOG_outfile, Save_data);

      DP dfdmu = Calculate_dGibbs_dchempot (Scan_dSet_dmu, Scan_Set, c_int, mu, Omega, kBT);
      DP dfdOm = Calculate_dGibbs_dchempot (Scan_dSet_dOmega, Scan_Set, c_int, mu, Omega, kBT);

      GFE_outfile << mu << "\t" << Calculate_Gibbs_Free_Energy (Scan_Set, c_int, mu, Omega, kBT) << "\t" << Scan_Set.diff
		  << "\t" << dfdmu << "\t" << dfdOm << "\t" << 0.5 * (-dfdmu - dfdOm) << "\t" << 0.5 * (-dfdmu + dfdOm)
		  << endl;
    } // for imu

    LOG_outfile.close();
    GFE_outfile.close();

    return;
  }

} // namespace ABACUS