729 lines
25 KiB
C++
729 lines
25 KiB
C++
/**********************************************************
|
|
|
|
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
|