152 lines
4.9 KiB
C++
152 lines
4.9 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: LiebLin_Tgt0.cc
|
|
|
|
Purpose: Finite temperature correlations for Lieb-Liniger
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
|
|
|
|
DP Entropy_Fixed_Delta (LiebLin_Bethe_State& RefState, int Delta)
|
|
{
|
|
// This function calculates the discrete entropy of a finite Lieb-Liniger state,
|
|
// counting the possible permutations within windows of fixed width.
|
|
|
|
// We assume that the quantum numbers are ordered.
|
|
|
|
// Fill in vector of occupancies:
|
|
// assume Ix2 are ordered, leave space of Delta on both sides
|
|
int nrIs = (RefState.Ix2[RefState.N-1] - RefState.Ix2[0])/2 + 1 + 2*Delta;
|
|
Vect<int> occupancy(0, nrIs); // leave space of Delta on both sides
|
|
for (int i = 0; i < RefState.N; ++i) occupancy[Delta + (RefState.Ix2[i] -RefState.Ix2[0])/2] = 1;
|
|
// Check:
|
|
int ncheck = 0;
|
|
for (int i = 0; i < nrIs; ++i) if(occupancy[i] == 1) ncheck++;
|
|
if (ncheck != RefState.N) {
|
|
cout << ncheck << "\t" << RefState.N << endl;
|
|
ABACUSerror("Counting q numbers incorrectly in Entropy.");
|
|
}
|
|
|
|
// Define some useful numbers:
|
|
Vect_DP oneoverDeltalnchoose(Delta + 1);
|
|
for (int i = 0; i <= Delta; ++i) oneoverDeltalnchoose[i] = ln_choose(Delta, i)/Delta;
|
|
|
|
// Compute entropy:
|
|
DP entropy = 0.0;
|
|
int np;
|
|
for (int i = 0; i < nrIs - Delta; ++i) {
|
|
np = 0;
|
|
for (int iD = 0; iD < Delta; ++iD) if (occupancy[i + iD] == 1) np++;
|
|
|
|
entropy += oneoverDeltalnchoose[np];
|
|
}
|
|
|
|
return(entropy);
|
|
}
|
|
|
|
DP Entropy (LiebLin_Bethe_State& RefState, int Delta)
|
|
{
|
|
// Perform an average of entropies for regulators from Delta to 2Delta:
|
|
DP entropysum = 0.0;
|
|
for (int ie = 0; ie < Delta; ++ie) entropysum += Entropy_Fixed_Delta (RefState, Delta + ie);
|
|
return(entropysum/Delta);
|
|
}
|
|
|
|
DP Entropy (LiebLin_Bethe_State& RefState)
|
|
{
|
|
int Delta = int(log(RefState.L));
|
|
return(Entropy (RefState, Delta));
|
|
}
|
|
|
|
DP Canonical_Free_Energy (LiebLin_Bethe_State& RefState, DP kBT)
|
|
{
|
|
return(RefState.E - kBT * Entropy (RefState));
|
|
}
|
|
|
|
DP rho_of_lambdaoc_1 (LiebLin_Bethe_State& RefState, DP lambdaoc, DP delta)
|
|
{
|
|
DP answer = 0.0;
|
|
for (int i = 0; i < RefState.N; ++i)
|
|
answer += atan((lambdaoc - RefState.lambdaoc[i])/delta + 0.5)
|
|
- atan((lambdaoc - RefState.lambdaoc[i])/delta - 0.5);
|
|
answer *= 1.0/(PI * delta * RefState.L);
|
|
|
|
return(answer);
|
|
}
|
|
|
|
DP rho_of_lambdaoc_2 (LiebLin_Bethe_State& RefState, DP lambdaoc, DP delta)
|
|
{
|
|
DP answer = 0.0;
|
|
for (int i = 0; i < RefState.N; ++i)
|
|
answer += 1.0/(pow(lambdaoc - RefState.lambdaoc[i], 2.0) + delta*delta);
|
|
answer *= delta/(PI * RefState.L);
|
|
|
|
return(answer);
|
|
}
|
|
|
|
// Better implementation: making use of rediscretized TBA state.
|
|
LiebLin_Bethe_State Canonical_Saddle_Point_State (DP c_int, DP L, int N, DP kBT)
|
|
{
|
|
// This function returns the discretized state minimizing the canonical free energy
|
|
// F = E - T S.
|
|
|
|
// This is obtained by rediscretizing the solution coming from TBA.
|
|
|
|
// ASSUMPTIONS:
|
|
// Periodic boundary conditions (the state which is output is forced to be symmetric Ix2 == -Ix2).
|
|
|
|
// For zero temperature, return the ground state:
|
|
if (fabs(kBT) < 1.0e-4) return(LiebLin_Bethe_State(c_int, L, N));
|
|
|
|
// Otherwise, return the discretized TBA saddle-point state:
|
|
LiebLin_TBA_Solution TBAsol = LiebLin_TBA_Solution_fixed_nbar (c_int, N/L, kBT, 1.0e-4, ABACUS::max(N/10, 10));
|
|
|
|
LiebLin_Bethe_State spstate = Discretized_LiebLin_Bethe_State (c_int, L, N, TBAsol.rho);
|
|
|
|
// Explicitly symmetrize:
|
|
for (int i = 0; i < N/2; ++i) spstate.Ix2[i] = -spstate.Ix2[N-1-i];
|
|
spstate.Compute_All(false);
|
|
|
|
return(spstate);
|
|
}
|
|
|
|
LiebLin_Bethe_State Add_Particle_at_Center (const LiebLin_Bethe_State& RefState)
|
|
{
|
|
LiebLin_Bethe_State ReturnState (RefState.c_int, RefState.L, RefState.N + 1);
|
|
|
|
// Add a quantum number in middle (explicitly: to right of index N/2)
|
|
// and shift quantum numbers by half-integer away from added one:
|
|
ReturnState.Ix2[RefState.N/2] = RefState.Ix2[RefState.N/2] - 1;
|
|
for (int i = 0; i < RefState.N+1; ++i)
|
|
ReturnState.Ix2[i + (i >= RefState.N/2)] = RefState.Ix2[i] - 1 + 2*(i >= RefState.N/2);
|
|
|
|
return(ReturnState);
|
|
}
|
|
|
|
LiebLin_Bethe_State Remove_Particle_at_Center (const LiebLin_Bethe_State& RefState)
|
|
{
|
|
LiebLin_Bethe_State ReturnState (RefState.c_int, RefState.L, RefState.N - 1);
|
|
|
|
// Remove midmost and shift quantum numbers by half-integer towards removed one:
|
|
for (int i = 0; i < RefState.N-1; ++i)
|
|
ReturnState.Ix2[i] = RefState.Ix2[i + (i >= RefState.N/2)] + 1 - 2*(i >= RefState.N/2);
|
|
|
|
return(ReturnState);
|
|
}
|
|
|
|
} // namespace ABACUS
|