jscaux
/
ABACUS
2
0
Fork 0
Code Issues 1 Pull Requests Releases Wiki Activity
ABACUS/src/TBA/TBA_2CBG.cc

1132 lines
44 KiB
C++
Raw Permalink Blame History

/**********************************************************
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
Reference in New Issue View Git Blame Copy Permalink