455 lines
14 KiB
C++
455 lines
14 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c)
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: src/HEIS/XXX_Bethe_State.cc
|
|
|
|
Purpose: Defines all functions for XXX_Bethe_State
|
|
|
|
******************************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
|
|
namespace ABACUS {
|
|
|
|
// Function prototypes
|
|
|
|
DP Theta_XXX (DP lambda, int nj, int nk);
|
|
DP ddlambda_Theta_XXX (DP lambda, int nj, int nk);
|
|
|
|
//***************************************************************************************************
|
|
|
|
// Function definitions: class XXX_Bethe_State
|
|
|
|
XXX_Bethe_State::XXX_Bethe_State ()
|
|
: Heis_Bethe_State()
|
|
{};
|
|
|
|
XXX_Bethe_State::XXX_Bethe_State (const XXX_Bethe_State& RefState) // copy constructor
|
|
: Heis_Bethe_State(RefState)
|
|
{
|
|
}
|
|
|
|
XXX_Bethe_State::XXX_Bethe_State (const Heis_Chain& RefChain, int M)
|
|
: Heis_Bethe_State(RefChain, M)
|
|
{
|
|
if (RefChain.Delta != 1.0) {
|
|
cout << setprecision(16) << RefChain.Delta << endl;
|
|
ABACUSerror("Delta != 1.0 in XXX_Bethe_State constructor");
|
|
}
|
|
}
|
|
|
|
XXX_Bethe_State::XXX_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& RefBase)
|
|
: Heis_Bethe_State(RefChain, RefBase)
|
|
{
|
|
if (RefChain.Delta != 1.0) {
|
|
cout << setprecision(16) << RefChain.Delta << endl;
|
|
ABACUSerror("Delta != 1.0 in XXX_Bethe_State constructor");
|
|
}
|
|
}
|
|
|
|
XXX_Bethe_State& XXX_Bethe_State::operator= (const XXX_Bethe_State& RefState)
|
|
{
|
|
if (this != &RefState) {
|
|
chain = RefState.chain;
|
|
base = RefState.base;
|
|
Ix2 = RefState.Ix2;
|
|
lambda = RefState.lambda;
|
|
BE = RefState.BE;
|
|
diffsq = RefState.diffsq;
|
|
conv = RefState.conv;
|
|
iter = RefState.iter;
|
|
iter_Newton = RefState.iter_Newton;
|
|
E = RefState.E;
|
|
iK = RefState.iK;
|
|
K = RefState.K;
|
|
lnnorm = RefState.lnnorm;
|
|
label = RefState.label;
|
|
}
|
|
return(*this);
|
|
}
|
|
|
|
// Member functions
|
|
|
|
void XXX_Bethe_State::Set_Free_lambdas()
|
|
{
|
|
// Sets all the rapidities to the solutions of the BAEs without scattering terms
|
|
|
|
for (int i = 0; i < chain.Nstrings; ++i) {
|
|
|
|
for (int alpha = 0; alpha < base[i]; ++alpha) {
|
|
|
|
lambda[i][alpha] = chain.Str_L[i] * 0.5 * tan(PI * 0.5 * Ix2[i][alpha]/chain.Nsites);
|
|
|
|
}
|
|
}
|
|
|
|
return;
|
|
}
|
|
|
|
bool XXX_Bethe_State::Check_Admissibility(char option)
|
|
{
|
|
// This function checks the admissibility of the Ix2's of a state:
|
|
// returns false if there are higher strings with Ix2 = 0, a totally symmetric distribution of I's at each level,
|
|
// and strings of equal length modulo 2 and parity with Ix2 = 0, meaning at least two equal roots in BAE.
|
|
|
|
bool answer = true;
|
|
Vect<bool> Zero_at_level(false, chain.Nstrings); // whether there exists an Ix2 == 0 at a given level
|
|
|
|
bool higher_string_on_zero = false;
|
|
|
|
for (int j = 0; j < chain.Nstrings; ++j) {
|
|
// The following line puts answer to true if there is at least one higher string with zero Ix2
|
|
for (int alpha = 0; alpha < base[j]; ++alpha)
|
|
if ((Ix2[j][alpha] == 0) && (chain.Str_L[j] > 2) && !(chain.Str_L[j] % 2))
|
|
higher_string_on_zero = true;
|
|
for (int alpha = 0; alpha < base[j]; ++alpha) if (Ix2[j][alpha] == 0) Zero_at_level[j] = true;
|
|
// NOTE: if base[j] == 0, Zero_at_level[j] remains false.
|
|
}
|
|
|
|
bool symmetric_state = (*this).Check_Symmetry();
|
|
|
|
bool string_coincidence = false;
|
|
// Checks that we have strings of equal length modulo 2 with Ix2 == 0, so equal rapidities, and inadmissibility
|
|
for (int j1 = 0; j1 < chain.Nstrings; ++j1) {
|
|
for (int j2 = j1 + 1; j2 < chain.Nstrings; ++j2)
|
|
if (Zero_at_level[j1] && Zero_at_level[j2] && (!((chain.Str_L[j1] + chain.Str_L[j2])%2)))
|
|
string_coincidence = true;
|
|
}
|
|
answer = !(symmetric_state && (higher_string_on_zero || string_coincidence /*|| onep_onem_on_zero*/));
|
|
|
|
// Now check that no Ix2 is equal to +N (since we take -N into account, and I + N == I by periodicity of exp)
|
|
|
|
for (int j = 0; j < chain.Nstrings; ++j)
|
|
for (int alpha = 0; alpha < base[j]; ++alpha)
|
|
if ((Ix2[j][alpha] < -chain.Nsites) || (Ix2[j][alpha] >= chain.Nsites)) answer = false;
|
|
|
|
if (!answer) {
|
|
E = 0.0;
|
|
K = 0.0;
|
|
conv = 0;
|
|
iter = 0;
|
|
iter_Newton = 0;
|
|
lnnorm = -100.0;
|
|
}
|
|
|
|
return(answer); // answer == true: nothing wrong with this Ix2_config
|
|
}
|
|
|
|
void XXX_Bethe_State::Compute_BE (int j, int alpha)
|
|
{
|
|
// Fills in the BE members with the value of the Bethe equations.
|
|
|
|
DP sumtheta = 0.0;
|
|
|
|
sumtheta = 0.0;
|
|
for (int k = 0; k < chain.Nstrings; ++k) {
|
|
for (int beta = 0; beta < base[k]; ++beta)
|
|
|
|
if ((chain.Str_L[j] == 1) && (chain.Str_L[k] == 1))
|
|
sumtheta += atan(lambda[j][alpha] - lambda[k][beta]);
|
|
|
|
else sumtheta += 0.5 * Theta_XXX((lambda[j][alpha] - lambda[k][beta]), chain.Str_L[j], chain.Str_L[k]);
|
|
}
|
|
sumtheta *= 2.0;
|
|
|
|
BE[j][alpha] = 2.0 * atan(2.0 * lambda[j][alpha]/chain.Str_L[j]) - (sumtheta + PI*Ix2[j][alpha])/chain.Nsites;
|
|
}
|
|
|
|
void XXX_Bethe_State::Compute_BE ()
|
|
{
|
|
// Fills in the BE members with the value of the Bethe equations.
|
|
|
|
DP sumtheta = 0.0;
|
|
|
|
for (int j = 0; j < chain.Nstrings; ++j) {
|
|
for (int alpha = 0; alpha < base[j]; ++alpha) {
|
|
|
|
sumtheta = 0.0;
|
|
for (int k = 0; k < chain.Nstrings; ++k) {
|
|
for (int beta = 0; beta < base[k]; ++beta)
|
|
|
|
if ((chain.Str_L[j] == 1) && (chain.Str_L[k] == 1))
|
|
sumtheta += atan(lambda[j][alpha] - lambda[k][beta]);
|
|
|
|
else sumtheta += 0.5 * Theta_XXX((lambda[j][alpha] - lambda[k][beta]), chain.Str_L[j], chain.Str_L[k]);
|
|
}
|
|
sumtheta *= 2.0;
|
|
|
|
BE[j][alpha] = 2.0 * atan(2.0 * lambda[j][alpha]/chain.Str_L[j]) - (sumtheta + PI*Ix2[j][alpha])/chain.Nsites;
|
|
}
|
|
}
|
|
}
|
|
|
|
DP XXX_Bethe_State::Iterate_BAE (int j, int alpha)
|
|
{
|
|
// Returns a new iteration value for lambda[j][alpha] given BE[j][alpha]
|
|
|
|
return(0.5 * chain.Str_L[j] * tan(0.5 * (2.0 * atan(2.0 * lambda[j][alpha]/chain.Str_L[j]) - BE[j][alpha])));
|
|
}
|
|
|
|
bool XXX_Bethe_State::Check_Rapidities()
|
|
{
|
|
bool nonan = true;
|
|
|
|
for (int j = 0; j < chain.Nstrings; ++j)
|
|
for (int alpha = 0; alpha < base[j]; ++alpha) nonan *= !is_nan(lambda[j][alpha]);
|
|
|
|
return nonan;
|
|
}
|
|
|
|
DP XXX_Bethe_State::String_delta ()
|
|
{
|
|
// Computes the sum of absolute value of \delta^{a, a+1} in string hypothesis, for a given Bethe eigenstate
|
|
|
|
DP delta = 0.0;
|
|
|
|
int occupied_strings = 0;
|
|
for (int i = 0; i < (*this).chain.Nstrings; ++i)
|
|
if ((*this).chain.Str_L[i] > 1) occupied_strings += (*this).base.Nrap[i];
|
|
|
|
if (occupied_strings == 0) delta = 0.0;
|
|
|
|
else {
|
|
|
|
Vect_DP ln_deltadiff(0.0, 1000); // contains ln |delta^{a, a+1}|
|
|
Vect_DP deltadiff(0.0, 1000); // contains |delta^{a, a+1}|
|
|
|
|
complex<DP> log_BAE_reg = 0.0;
|
|
|
|
for (int j = 0; j < (*this).chain.Nstrings; ++j) {
|
|
for (int alpha = 0; alpha < (*this).base[j]; ++alpha) {
|
|
|
|
ln_deltadiff = 0.0;
|
|
|
|
for (int a = 1; a <= (*this).chain.Str_L[j]; ++a) {
|
|
|
|
if ((*this).chain.Str_L[j] > 1) { // else the BAE are already 1
|
|
|
|
log_BAE_reg = DP((*this).chain.Nsites)
|
|
* log(((*this).lambda[j][alpha] + 0.5 * II * ((*this).chain.Str_L[j] + 1.0 - 2.0 * a + 1.0))
|
|
/((*this).lambda[j][alpha] + 0.5 * II * ((*this).chain.Str_L[j] + 1.0 - 2.0 * a - 1.0)));
|
|
|
|
for (int k = 0; k < (*this).chain.Nstrings; ++k)
|
|
for (int beta = 0; beta < (*this).base[k]; ++beta)
|
|
for (int b = 1; b <= (*this).chain.Str_L[k]; ++b) {
|
|
if ((j != k) || (alpha != beta) || (a != b - 1))
|
|
|
|
log_BAE_reg += log((*this).lambda[j][alpha] + 0.5 * II * ((*this).chain.Str_L[j] + 1.0 - 2.0 * a )
|
|
- ((*this).lambda[k][beta] + 0.5 * II * ((*this).chain.Str_L[k] + 1.0 - 2.0 * b )
|
|
) - II );
|
|
|
|
if ((j != k) || (alpha != beta) || (a != b + 1))
|
|
|
|
log_BAE_reg -= log(((*this).lambda[j][alpha] + 0.5 * II * ((*this).chain.Str_L[j] + 1.0 - 2.0 * a))
|
|
- ((*this).lambda[k][beta] + 0.5 * II * ((*this).chain.Str_L[k] + 1.0 - 2.0 * b))
|
|
+ II );
|
|
}
|
|
|
|
// The regular LHS of BAE is now defined. Now sum up the deltas...
|
|
|
|
if (a == 1) ln_deltadiff[0] = - real(log_BAE_reg);
|
|
|
|
else if (a < (*this).chain.Str_L[j]) ln_deltadiff[a - 1] = ln_deltadiff[a-2] - real(log_BAE_reg);
|
|
|
|
else if (a == (*this).chain.Str_L[j]) ln_deltadiff[a-1] = real(log_BAE_reg);
|
|
|
|
} // if ((*this).chain.Str_L[j] > 1)
|
|
|
|
} // for (int a = 1; ...
|
|
|
|
for (int a = 0; a < (*this).chain.Str_L[j]; ++a) {
|
|
deltadiff[a] = ln_deltadiff[a] != 0.0 ? exp(ln_deltadiff[a]) : 0.0;
|
|
delta += fabs(deltadiff[a]);
|
|
}
|
|
|
|
} // alpha sum
|
|
} // j sum
|
|
|
|
if (is_nan(delta)) delta = 1.0; // sentinel
|
|
|
|
} // else
|
|
|
|
return delta;
|
|
}
|
|
|
|
|
|
void XXX_Bethe_State::Compute_Energy ()
|
|
{
|
|
DP sum = 0.0;
|
|
|
|
for (int j = 0; j < chain.Nstrings; ++j) {
|
|
for (int alpha = 0; alpha < base[j]; ++alpha) {
|
|
sum += chain.Str_L[j] / ( 4.0 * lambda[j][alpha] * lambda[j][alpha] + chain.Str_L[j] * chain.Str_L[j]);
|
|
}
|
|
}
|
|
|
|
sum *= - chain.J * 2.0;
|
|
|
|
E = sum;
|
|
|
|
return;
|
|
}
|
|
|
|
void XXX_Bethe_State::Build_Reduced_Gaudin_Matrix (SQMat<complex<DP> >& Gaudin_Red)
|
|
{
|
|
|
|
if (Gaudin_Red.size() != base.Nraptot)
|
|
ABACUSerror("Passing matrix of wrong size in Build_Reduced_Gaudin_Matrix.");
|
|
|
|
int index_jalpha;
|
|
int index_kbeta;
|
|
|
|
DP sum_hbar_XXX = 0.0;
|
|
|
|
index_jalpha = 0;
|
|
for (int j = 0; j < chain.Nstrings; ++j) {
|
|
for (int alpha = 0; alpha < base[j]; ++alpha) {
|
|
index_kbeta = 0;
|
|
for (int k = 0; k < chain.Nstrings; ++k) {
|
|
for (int beta = 0; beta < base[k]; ++beta) {
|
|
|
|
if ((j == k) && (alpha == beta)) {
|
|
|
|
sum_hbar_XXX = 0.0;
|
|
|
|
for (int kp = 0; kp < chain.Nstrings; ++kp) {
|
|
for (int betap = 0; betap < base[kp]; ++betap) {
|
|
if (!((j == kp) && (alpha == betap)))
|
|
sum_hbar_XXX
|
|
+= ddlambda_Theta_XXX (lambda[j][alpha] - lambda[kp][betap], chain.Str_L[j], chain.Str_L[kp]);
|
|
}
|
|
}
|
|
|
|
Gaudin_Red[index_jalpha][index_kbeta]
|
|
= complex<DP> ( chain.Nsites * chain.Str_L[j]/(lambda[j][alpha] * lambda[j][alpha]
|
|
+ 0.25 * chain.Str_L[j] * chain.Str_L[j]) - sum_hbar_XXX);
|
|
}
|
|
|
|
else {
|
|
if ((chain.Str_L[j] == 1) && (chain.Str_L[k] == 1))
|
|
Gaudin_Red[index_jalpha][index_kbeta] =
|
|
complex<DP> ( 2.0/(pow(lambda[j][alpha] - lambda[k][beta], 2.0) + 1.0));
|
|
|
|
else
|
|
Gaudin_Red[index_jalpha][index_kbeta] =
|
|
complex<DP> (ddlambda_Theta_XXX (lambda[j][alpha] - lambda[k][beta], chain.Str_L[j], chain.Str_L[k]));
|
|
}
|
|
index_kbeta++;
|
|
}
|
|
}
|
|
index_jalpha++;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
bool XXX_Bethe_State::Check_Finite_rap ()
|
|
{
|
|
bool answer = true;
|
|
|
|
for (int j = 0; j < chain.Nstrings; ++j) {
|
|
for (int alpha = 0; alpha < base[j]; ++alpha) {
|
|
if (fabs(lambda[j][alpha]) > 1.0e6) answer = false;
|
|
}
|
|
}
|
|
|
|
return(answer);
|
|
}
|
|
|
|
// ****************************************************************************************************
|
|
|
|
// non-member functions
|
|
|
|
DP Theta_XXX (DP lambda, int nj, int nk)
|
|
{
|
|
DP result;
|
|
|
|
if ((nj == 1) && (nk == 1)) result = 2.0 * atan(lambda);
|
|
|
|
else {
|
|
|
|
result = (nj == nk) ? 0.0 : 2.0 * atan(2.0 * lambda/fabs(nj - nk));
|
|
|
|
for (int a = 1; a < ABACUS::min(nj, nk); ++a) result += 4.0 * atan(2.0 * lambda/(fabs(nj - nk) + 2*a));
|
|
|
|
result += 2.0 * atan(2.0 * lambda/(nj + nk));
|
|
}
|
|
|
|
return (result);
|
|
}
|
|
|
|
DP ddlambda_Theta_XXX (DP lambda, int nj, int nk)
|
|
{
|
|
int n = abs(nj - nk);
|
|
|
|
DP result = (nj == nk) ? 0.0 : DP(n)/(lambda * lambda + 0.25 * n * n);
|
|
|
|
for (int a = 1; a < ABACUS::min(nj, nk); ++a) result += 2.0 * (n + 2.0*a)
|
|
/ (lambda * lambda + 0.25 * (n + 2.0*a) * (n + 2.0*a));
|
|
|
|
result += DP(nj + nk)/(lambda * lambda + 0.25 * (nj + nk) * (nj + nk));
|
|
|
|
return (result);
|
|
}
|
|
|
|
XXX_Bethe_State Add_Particle_at_Center (const XXX_Bethe_State& RefState)
|
|
{
|
|
if (2*RefState.base.Mdown == RefState.chain.Nsites)
|
|
ABACUSerror("Trying to add a down spin to a zero-magnetized chain in Add_Particle_at_Center.");
|
|
|
|
Vect<int> newM = RefState.base.Nrap;
|
|
newM[0] = newM[0] + 1;
|
|
|
|
Heis_Base newBase (RefState.chain, newM);
|
|
|
|
XXX_Bethe_State ReturnState (RefState.chain, newBase);
|
|
|
|
for (int il = 1; il < RefState.chain.Nstrings; ++il)
|
|
for (int alpha = 0; alpha < RefState.base.Nrap[il]; ++alpha)
|
|
ReturnState.Ix2[il][alpha] = RefState.Ix2[il][alpha];
|
|
|
|
// Add a quantum number in middle (explicitly: to right of index M[0]/2)
|
|
// and shift quantum numbers by half-integer away from added one:
|
|
ReturnState.Ix2[0][RefState.base.Nrap[0]/2] = RefState.Ix2[0][RefState.base.Nrap[0]/2] - 1;
|
|
for (int i = 0; i < RefState.base.Nrap[0] + 1; ++i)
|
|
ReturnState.Ix2[0][i + (i >= RefState.base.Nrap[0]/2)]
|
|
= RefState.Ix2[0][i] - 1 + 2*(i >= RefState.base.Nrap[0]/2);
|
|
|
|
return(ReturnState);
|
|
}
|
|
|
|
|
|
XXX_Bethe_State Remove_Particle_at_Center (const XXX_Bethe_State& RefState)
|
|
{
|
|
if (RefState.base.Nrap[0] == 0)
|
|
ABACUSerror("Trying to remove a down spin in an empty Nrap[0] state.");
|
|
|
|
Vect<int> newM = RefState.base.Nrap;
|
|
newM[0] = newM[0] - 1;
|
|
|
|
Heis_Base newBase (RefState.chain, newM);
|
|
|
|
XXX_Bethe_State ReturnState (RefState.chain, newBase);
|
|
|
|
for (int il = 1; il < RefState.chain.Nstrings; ++il)
|
|
for (int alpha = 0; alpha < RefState.base.Nrap[il]; ++alpha)
|
|
ReturnState.Ix2[il][alpha] = RefState.Ix2[il][alpha];
|
|
|
|
// Remove midmost and shift quantum numbers by half-integer towards removed one:
|
|
for (int i = 0; i < RefState.base.Nrap[0]-1; ++i)
|
|
ReturnState.Ix2[0][i]
|
|
= RefState.Ix2[0][i + (i >= RefState.base.Nrap[0]/2)] + 1 - 2*(i >= RefState.base.Nrap[0]/2);
|
|
|
|
return(ReturnState);
|
|
}
|
|
|
|
|
|
} // namespace ABACUS
|