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

1434 lines
45 KiB
C++
Raw Blame History

/**********************************************************
This software is part of J.-S. Caux's ABACUS library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: src/HEIS/Heis.cc
Purpose: defines functions in all HEIS classes.
***********************************************************/
#include "ABACUS.h"
using namespace std;
namespace ABACUS {
//***************************************************************************************************
// Function definitions: class Heis_Chain
Heis_Chain::Heis_Chain()
: J(0.0), Delta(0.0), anis(0.0), hz(0.0), Nsites(0), Nstrings(0), Str_L(0), par(0),
si_n_anis_over_2(new DP[1]), co_n_anis_over_2(new DP[1]), ta_n_anis_over_2(new DP[1]), prec(ITER_REQ_PREC) {}
Heis_Chain::Heis_Chain (DP JJ, DP DD, DP hh, int NN)
: J(JJ), Delta(DD), anis(0.0), hz(hh), Nsites(NN), Nstrings(MAXSTRINGS), Str_L(new int[MAXSTRINGS]),
par(new int[MAXSTRINGS]), si_n_anis_over_2(new DP[10*MAXSTRINGS]), co_n_anis_over_2(new DP[10*MAXSTRINGS]),
ta_n_anis_over_2(new DP[10*MAXSTRINGS]), prec(ITER_REQ_PREC)
{
// We restrict to even chains everywhere
if (NN % 2) ABACUSerror("Please use an even-length chain.");
if ((Delta > 0.0) && (Delta < 1.0)) {
// gapless XXZ case
anis = acos(DD);
// Set the Str_L and par vectors:
DP gammaoverpi = acos(DD)/PI;
Vect<int> Nu(MAXSTRINGS);
DP gammaoverpi_eff = gammaoverpi;
DP gammaoverpi_reached = 0.0;
// First, define the series Nu[i], like in TakahashiBOOK eqn 9.35
// CAREFUL: the labelling is different, Nu[0] is \nu_1, etc.
int l = 0;
int ml_temp = 0; // checks the sum of all Nu's
while (fabs(gammaoverpi - gammaoverpi_reached) > 1.0e-8){
Nu[l] = int(1.0/gammaoverpi_eff);
ml_temp += Nu[l];
if (Nu[l] == 0) gammaoverpi_eff = 0.0;
else gammaoverpi_eff = 1.0/gammaoverpi_eff - DP(Nu[l]);
// compute gammaoverpi_reached:
gammaoverpi_reached = Nu[l];
for (int p = 0; p < l; ++p) gammaoverpi_reached = Nu[l - p - 1] + 1.0/gammaoverpi_reached;
gammaoverpi_reached = 1.0/gammaoverpi_reached;
l++;
if (ml_temp > MAXSTRINGS) break; // we defined Str_L and par as arrays of at most MAXSTRINGS elements, so we cut off here...
}
// Check: make sure the last Nu is greater than one: if one, add 1 to previous Nu
if (l > 1) {
if (Nu[l-1] == 1) {
Nu[l-1] = 0;
Nu[l-2] += 1;
l -= 1;
}
}
// Length of continued fraction is l-1, which is denoted l in Takahashi
// Second, define the series y[i] and m[i] as in TakahashiBOOK eqn 9.36
// y_{-1} is treated separately, and here y[i] = y_i, m[i] = m_i, i = 0, ..., l
Vect<int> y(0, l+1);
Vect<int> m(0, l+1);
y[0] = 1;
m[0] = 0;
y[1] = Nu[0];
m[1] = Nu[0];
for (int i = 2; i <= l; ++i) {
y[i] = y[i - 2] + Nu[i-1]*y[i-1];
m[i] = 0;
for (int k = 0; k < i; ++k) m[i] += Nu[k];
}
// Now determine the lengths and parity, following TakahashiBOOK eqn 9.37
// Nstrings = ABACUS::min(m[l] + 1, MAXSTRINGS); // number of different strings that are possible for this chain
// Always set to MAXSTRINGS
Str_L[0] = 1;
par[0] = 1;
Str_L[1] = 1;
par[1] = -1;
int max_j = 0;
for (int i = 0; i < l; ++i) {
for (int j = ABACUS::max(1, m[i]) + 1; j < ABACUS::min(m[i+1] + 1, MAXSTRINGS); ++j) {
if (i == 0) Str_L[j] = (j - m[0])* y[0];
else if (i >= 1) Str_L[j] = y[i-1] + (j - m[i])*y[i];
par[j] = (int(floor((Str_L[j] - 1)*gammaoverpi)) % 2) ? -1 : 1;
max_j = j;
}
}
// Set the rest of Str_L and par vector elements to zero
for (int i = max_j + 1; i < MAXSTRINGS; ++i) {
Str_L[i] = 0;
par[i] = 0;
}
// Set the sin, cos and tan_n_zeta_over_2 vectors:
for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = sin(i * anis/2.0);
for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = cos(i * anis/2.0);
for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = tan(i * anis/2.0);
} // if XXZ gapless
else if (Delta == 1.0) {
// Isotropic antiferromagnet
anis = 0.0;
// Set Str_L and par:
for (int i = 0; i < MAXSTRINGS; ++i) {
Str_L[i] = i + 1;
par[i] = 1;
}
} // if XXX AFM
else if (Delta > 1.0) {
// Gapped antiferromagnet
anis = acosh(DD);
// Set the Str_L and par vectors:
for (int i = 0; i < MAXSTRINGS; ++i) {
Str_L[i] = i + 1;
par[i] = 1;
}
// Set the sinh, cosh and tanh_n_eta_over_2 vectors:
for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = sinh(i * anis/2.0);
for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = cosh(i * anis/2.0);
for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = tanh(i * anis/2.0);
} // if XXZ_gpd
}
Heis_Chain::Heis_Chain (const Heis_Chain& RefChain) // copy constructor
{
J = RefChain.J;
Delta = RefChain.Delta;
anis = RefChain.anis;
hz = RefChain.hz;
Nsites = RefChain.Nsites;
Nstrings = RefChain.Nstrings;
Str_L = new int[RefChain.Nstrings];
for (int i = 0; i < RefChain.Nstrings; ++i) Str_L[i] = RefChain.Str_L[i];
par = new int[RefChain.Nstrings];
for (int i = 0; i < RefChain.Nstrings; ++i) par[i] = RefChain.par[i];
si_n_anis_over_2 = new DP[10*MAXSTRINGS];
for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = RefChain.si_n_anis_over_2[i];
co_n_anis_over_2 = new DP[10*MAXSTRINGS];
for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = RefChain.co_n_anis_over_2[i];
ta_n_anis_over_2 = new DP[10*MAXSTRINGS];
for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = RefChain.ta_n_anis_over_2[i];
prec = RefChain.prec;
}
Heis_Chain& Heis_Chain::operator= (const Heis_Chain& RefChain) // assignment operator
{
if (this != &RefChain) {
J = RefChain.J;
Delta = RefChain.Delta;
anis = RefChain.anis;
hz = RefChain.hz;
Nsites = RefChain.Nsites;
Nstrings = RefChain.Nstrings;
if (Str_L != 0) delete[] Str_L;
Str_L = new int[RefChain.Nstrings];
for (int i = 0; i < RefChain.Nstrings; ++i) Str_L[i] = RefChain.Str_L[i];
if (par != 0) delete[] par;
par = new int[RefChain.Nstrings];
for (int i = 0; i < RefChain.Nstrings; ++i) par[i] = RefChain.par[i];
if (si_n_anis_over_2 != 0) delete[] si_n_anis_over_2;
si_n_anis_over_2 = new DP[10*MAXSTRINGS];
for (int i = 0; i < 10*MAXSTRINGS; ++i) si_n_anis_over_2[i] = RefChain.si_n_anis_over_2[i];
if (co_n_anis_over_2 != 0) delete[] co_n_anis_over_2;
co_n_anis_over_2 = new DP[10*MAXSTRINGS];
for (int i = 0; i < 10*MAXSTRINGS; ++i) co_n_anis_over_2[i] = RefChain.co_n_anis_over_2[i];
if (ta_n_anis_over_2 != 0) delete[] ta_n_anis_over_2;
ta_n_anis_over_2 = new DP[10*MAXSTRINGS];
for (int i = 0; i < 10*MAXSTRINGS; ++i) ta_n_anis_over_2[i] = RefChain.ta_n_anis_over_2[i];
prec = RefChain.prec;
}
return *this;
}
bool Heis_Chain::operator== (const Heis_Chain& RefChain)
{
bool answer = true;
if ((J != RefChain.J) || (Nsites != RefChain.Nsites) || (Delta != RefChain.Delta)) answer = false;
return(answer);
}
bool Heis_Chain::operator!= (const Heis_Chain& RefChain)
{
bool answer = false;
if ((J != RefChain.J) || (Nsites != RefChain.Nsites) || (Delta != RefChain.Delta)) answer = true;
return(answer);
}
Heis_Chain::~Heis_Chain()
{
if (Str_L != 0) delete[] Str_L;
if (par != 0) delete[] par;
if (si_n_anis_over_2 != 0) delete[] si_n_anis_over_2;
if (co_n_anis_over_2 != 0) delete[] co_n_anis_over_2;
if (ta_n_anis_over_2 != 0) delete[] ta_n_anis_over_2;
}
//***************************************************************************************************
// Function definitions: class Heis_Base
Heis_Base::Heis_Base () : Mdown(0), Nrap(Vect<int>()), Nraptot(0), Ix2_infty(Vect<DP>()),
Ix2_min(Vect<int>()), Ix2_max(Vect<int>()), dimH(0.0), baselabel("") {}
Heis_Base::Heis_Base (const Heis_Base& RefBase) // copy constructor
: Mdown(RefBase.Mdown), Nrap(Vect<int>(RefBase.Nrap.size())), Nraptot(RefBase.Nraptot),
Ix2_infty(Vect<DP>(RefBase.Nrap.size())), Ix2_min(Vect<int>(RefBase.Nrap.size())),
Ix2_max(Vect<int>(RefBase.Nrap.size())), baselabel(RefBase.baselabel)
{
for (int i = 0; i < Nrap.size(); ++i) {
Nrap[i] = RefBase.Nrap[i];
Ix2_infty[i] = RefBase.Ix2_infty[i];
Ix2_min[i] = RefBase.Ix2_min[i];
Ix2_max[i] = RefBase.Ix2_max[i];
dimH = RefBase.dimH;
}
}
Heis_Base::Heis_Base (const Heis_Chain& RefChain, int M)
: Mdown(M), Nrap(Vect<int>(RefChain.Nstrings)), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
Ix2_min(Vect<int>(RefChain.Nstrings)), Ix2_max(Vect<int>(RefChain.Nstrings))
{
for (int i = 0; i < RefChain.Nstrings; ++i) Nrap[i] = 0;
Nrap[0] = M;
Nraptot = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];
stringstream M0out;
M0out << M;
baselabel = M0out.str();
// Now compute the Ix2_infty numbers
(*this).Compute_Ix2_limits(RefChain);
// Compute dimensionality of this sub-Hilbert space:
complex<double> ln_dimH_cx = 0.0;
for (int i = 0; i < RefChain.Nstrings; ++i)
if (Nrap[i] > 0) ln_dimH_cx += ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2))
- ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2 - Nrap[i]))
- ln_Gamma(complex<double>(Nrap[i] + 1));
dimH = exp(real(ln_dimH_cx));
}
Heis_Base::Heis_Base (const Heis_Chain& RefChain, const Vect<int>& Nrapidities)
: Mdown(0), Nrap(Nrapidities), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
Ix2_min(Vect<int>(RefChain.Nstrings)), Ix2_max(Vect<int>(RefChain.Nstrings))
{
// Check consistency of Nrapidities vector with RefChain
if (RefChain.Nstrings != Nrapidities.size())
ABACUSerror("Incompatible Nrapidities vector used in Heis_Base constructor.");
int Mcheck = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
Mdown = Mcheck;
Nraptot = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];
stringstream M0out;
M0out << Nrapidities[0];
baselabel = M0out.str();
for (int itype = 1; itype < Nrapidities.size(); ++itype)
if (Nrapidities[itype] > 0) {
baselabel += TYPESEP;
stringstream typeout;
typeout << itype;
baselabel += typeout.str();
baselabel += EXCSEP;
stringstream Mout;
Mout << Nrapidities[itype];
baselabel += Mout.str();
}
// Now compute the Ix2_infty numbers
(*this).Compute_Ix2_limits(RefChain);
// Compute dimensionality of this sub-Hilbert space:
complex<double> ln_dimH_cx = 0.0;
for (int i = 0; i < RefChain.Nstrings; ++i)
if (Nrap[i] > 0) ln_dimH_cx += ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2))
- ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2 - Nrap[i]))
- ln_Gamma(complex<double>(Nrap[i] + 1));
dimH = exp(real(ln_dimH_cx));
}
Heis_Base::Heis_Base (const Heis_Chain& RefChain, string baselabel_ref)
: Mdown(0), Nrap(Vect<int>(0, RefChain.Nstrings)), Nraptot(0), Ix2_infty(Vect<DP>(RefChain.Nstrings)),
Ix2_min(Vect<int>(RefChain.Nstrings)), Ix2_max(Vect<int>(RefChain.Nstrings)), baselabel (baselabel_ref)
{
// Build Nrapidities vector from baselabel_ref.
// This is simply done by using the state label standard reading function after conveniently
// making baselabel into a label (as for a state with no excitations):
string label_ref = baselabel + "_0_";
Vect<Vect<int> > dummyOriginIx2(1);
//State_Label_Data labeldata = Read_State_Label (label_ref, dummyOriginIx2);
State_Label_Data labeldata = Read_Base_Label (label_ref);
// Initialize Nrap:
for (int i = 0; i < labeldata.type.size(); ++i)
Nrap[labeldata.type[i] ] = labeldata.M[i];
int Mcheck = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += RefChain.Str_L[i] * Nrap[i];
Mdown = Mcheck;
Nraptot = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Nraptot += Nrap[i];
// Now compute the Ix2_infty numbers
(*this).Compute_Ix2_limits(RefChain);
// Compute dimensionality of this sub-Hilbert space:
complex<double> ln_dimH_cx = 0.0;
for (int i = 0; i < RefChain.Nstrings; ++i)
if (Nrap[i] > 0) ln_dimH_cx += ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2))
- ln_Gamma(complex<double>((Ix2_max[i] - Ix2_min[i])/2 + 2 - Nrap[i]))
- ln_Gamma(complex<double>(Nrap[i] + 1));
dimH = exp(real(ln_dimH_cx));
}
Heis_Base& Heis_Base::operator= (const Heis_Base& RefBase)
{
if (this != & RefBase) {
Mdown = RefBase.Mdown;
Nrap = RefBase.Nrap;
Nraptot = RefBase.Nraptot;
Ix2_infty = RefBase.Ix2_infty;
Ix2_min = RefBase.Ix2_min;
Ix2_max = RefBase.Ix2_max;
dimH = RefBase.dimH;
baselabel = RefBase.baselabel;
}
return(*this);
}
bool Heis_Base::operator== (const Heis_Base& RefBase)
{
bool answer = (Nrap == RefBase.Nrap);
return (answer);
}
bool Heis_Base::operator!= (const Heis_Base& RefBase)
{
bool answer = (Nrap != RefBase.Nrap);
return (answer);
}
void Heis_Base::Compute_Ix2_limits (const Heis_Chain& RefChain)
{
if ((RefChain.Delta > 0.0) && (RefChain.Delta < 1.0)) {
// Compute the Ix2_infty numbers
DP sum1 = 0.0;
DP sum2 = 0.0;
for (int j = 0; j < RefChain.Nstrings; ++j) {
sum1 = 0.0;
for (int k = 0; k < RefChain.Nstrings; ++k) {
sum2 = 0.0;
sum2 += (RefChain.Str_L[j] == RefChain.Str_L[k]) ? 0.0 :
2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
- 0.5 * fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) * RefChain.anis));
sum2 += 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
- 0.5 * (RefChain.Str_L[j] + RefChain.Str_L[k]) * RefChain.anis));
for (int a = 1; a < ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]); ++a)
sum2 += 2.0 * 2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j] * RefChain.par[k])
- 0.5 * (fabs(RefChain.Str_L[j] - RefChain.Str_L[k]) + 2.0*a) * RefChain.anis));
sum1 += (Nrap[k] - ((j == k) ? 1 : 0)) * sum2;
}
Ix2_infty[j] = (1.0/PI) * fabs(RefChain.Nsites *
2.0 * atan(tan(0.25 * PI * (1.0 + RefChain.par[j])
- 0.5 * RefChain.Str_L[j] * RefChain.anis)) - sum1);
} // The Ix2_infty are now set.
// Now compute the Ix2_max limits
for (int j = 0; j < RefChain.Nstrings; ++j) {
Ix2_max[j] = int(floor(Ix2_infty[j])); // sets basic integer
// Reject formally infinite rapidities (i.e. if Delta is root of unity)
//cout << "Ix2_infty - Ix2_max = " << Ix2_infty[j] - Ix2_max[j] << endl;
//if (Ix2_infty[j] == Ix2_max[j]) {
//Ix2_max[j] -= 2;
//}
// If Nrap is even, Ix2_max must be odd. If odd, then even.
if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;
while (Ix2_max[j] > RefChain.Nsites) {
Ix2_max[j] -= 2;
}
Ix2_min[j] = -Ix2_max[j];
}
} // if XXZ gapless
else if (RefChain.Delta == 1.0) {
// Compute the Ix2_infty numbers
int sum1 = 0;
for (int j = 0; j < RefChain.Nstrings; ++j) {
sum1 = 0;
for (int k = 0; k < RefChain.Nstrings; ++k) {
sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
}
Ix2_infty[j] = (RefChain.Nsites + 1.0 - sum1); // to get counting right...
} // The Ix2_infty are now set.
// Now compute the Ix2_max limits
for (int j = 0; j < RefChain.Nstrings; ++j) {
Ix2_max[j] = int(floor(Ix2_infty[j])); // sets basic integer
// Give the correct parity to Ix2_max
// If Nrap is even, Ix2_max must be odd. If odd, then even.
if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;
// If Ix2_max equals Ix2_infty, we reduce it by 2:
if (Ix2_max[j] == int(Ix2_infty[j])) Ix2_max[j] -= 2;
while (Ix2_max[j] > RefChain.Nsites) {
Ix2_max[j] -= 2;
}
Ix2_min[j] = -Ix2_max[j];
}
} // if XXX AFM
else if (RefChain.Delta > 1.0) {
// Compute the Ix2_infty numbers
int sum1 = 0;
for (int j = 0; j < RefChain.Nstrings; ++j) {
sum1 = 0;
for (int k = 0; k < RefChain.Nstrings; ++k) {
sum1 += Nrap[k] * (2 * ABACUS::min(RefChain.Str_L[j], RefChain.Str_L[k]) - ((j == k) ? 1 : 0));
}
Ix2_infty[j] = (RefChain.Nsites - 1 + 2 * RefChain.Str_L[j] - sum1);
} // The Ix2_infty are now set.
// Now compute the Ix2_max limits
for (int j = 0; j < RefChain.Nstrings; ++j) {
Ix2_max[j] = int(floor(Ix2_infty[j])); // sets basic integer
// Give the correct parity to Ix2_max
// If Nrap is even, Ix2_max must be odd. If odd, then even.
if (!((Nrap[j] + Ix2_max[j]) % 2)) Ix2_max[j] -= 1;
// If Ix2_max equals Ix2_infty, we reduce it by 2:
//if (Ix2_max[j] == Ix2_infty[j]) Ix2_max[j] -= 2;
while (Ix2_max[j] > RefChain.Nsites) {
Ix2_max[j] -= 2;
}
Ix2_min[j] = -Ix2_max[j];
}
} // if XXZ_gpd
}
//***************************************************************************************************
// Function definitions: class Lambda
Lambda::Lambda () : lambda(NULL) {}
Lambda::Lambda (const Heis_Chain& RefChain, int M)
: Nstrings(1), Nrap(Vect<int>(M,1)), Nraptot(M), lambda(new DP*[1]) // single type of string here
{
lambda[0] = new DP[M];
for (int j = 0; j < M; ++j) lambda[0][j] = 0.0;
}
Lambda::Lambda (const Heis_Chain& RefChain, const Heis_Base& base)
: Nstrings(RefChain.Nstrings), Nrap(base.Nrap), Nraptot(base.Nraptot), lambda(new DP*[RefChain.Nstrings])
{
lambda[0] = new DP[base.Nraptot];
for (int i = 1; i < RefChain.Nstrings; ++i) lambda[i] = lambda[i-1] + base[i-1];
for (int i = 0; i < RefChain.Nstrings; ++i) {
for (int j = 0; j < base[i]; ++j) lambda[i][j] = 0.0;
}
}
Lambda& Lambda::operator= (const Lambda& RefLambda)
{
if (this != &RefLambda) {
Nstrings = RefLambda.Nstrings;
Nrap = RefLambda.Nrap;
Nraptot = RefLambda.Nraptot;
if (lambda != 0) {
delete[] (lambda[0]);
delete[] (lambda);
}
lambda = new DP*[Nstrings];
lambda[0] = new DP[Nraptot];
for (int i = 1; i < Nstrings; ++i) lambda[i] = lambda[i-1] + Nrap[i-1];
for (int i = 0; i < Nstrings; ++i)
for (int j = 0; j < Nrap[i]; ++j) lambda[i][j] = RefLambda.lambda[i][j];
}
return(*this);
}
Lambda::~Lambda()
{
if (lambda != 0) {
delete[] (lambda[0]);
delete[] lambda;
}
}
//***************************************************************************************************
// Function definitions: class Heis_Bethe_State
Heis_Bethe_State::Heis_Bethe_State ()
: chain(Heis_Chain()), base(Heis_Base()), Ix2(Vect<Vect<int> > (1)),
lambda(Lambda(chain, 1)), BE(Lambda(chain, 1)), diffsq(0.0), conv(0), dev(1.0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0), label("")
{
};
Heis_Bethe_State::Heis_Bethe_State (const Heis_Bethe_State& RefState) // copy constructor
: chain(RefState.chain), base(RefState.base), Ix2(RefState.Ix2),
lambda(Lambda(RefState.chain, RefState.base)), BE(Lambda(RefState.chain, RefState.base)),
diffsq(RefState.diffsq), conv(RefState.conv), dev(RefState.dev),
iter(RefState.iter), iter_Newton(RefState.iter_Newton),
E(RefState.E), iK(RefState.iK), K(RefState.K), lnnorm(RefState.lnnorm),
label(RefState.label)
{
// copy arrays into new ones
for (int j = 0; j < RefState.chain.Nstrings; ++j) {
for (int alpha = 0; alpha < RefState.base[j]; ++j) {
lambda[j][alpha] = RefState.lambda[j][alpha];
}
}
}
Heis_Bethe_State::Heis_Bethe_State (const Heis_Chain& RefChain, int M)
: chain(RefChain), base(RefChain, M), Ix2(Vect<Vect<int> > (1)),
lambda(Lambda(RefChain, M)),
BE(Lambda(RefChain, M)), diffsq(1.0), conv(0), dev(1.0), iter(0), iter_Newton(0),
E(0.0), iK(0), K(0.0), lnnorm(-100.0)
{
Ix2[0] = Vect<int> (M);
for (int j = 0; j < M; ++j) Ix2[0][j] = -(M - 1) + 2*j;
stringstream M0out;
M0out << M;
label = M0out.str() + "_0_";
}
Heis_Bethe_State::Heis_Bethe_State (const Heis_Chain& RefChain, const Heis_Base& RefBase)
: chain(RefChain), base(RefBase), Ix2 (Vect<Vect<int> > (RefChain.Nstrings)),
lambda(Lambda(RefChain, RefBase)),
BE(Lambda(RefChain, RefBase)), diffsq(1.0), conv(0), dev(1.0),
iter(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0)
{
// Check that the number of rapidities is consistent with Mdown
int Mcheck = 0;
for (int i = 0; i < RefChain.Nstrings; ++i) Mcheck += base[i] * RefChain.Str_L[i];
if (RefBase.Mdown != Mcheck) ABACUSerror("Wrong M from Nrapidities input vector, in Heis_Bethe_State constructor.");
for (int i = 0; i < RefChain.Nstrings; ++i) Ix2[i] = Vect<int> (base[i]);
// And now for the other string types...
for (int i = 0; i < RefChain.Nstrings; ++i) {
for (int j = 0; j < base[i]; ++j) Ix2[i][j] = -(base[i] - 1) + 2*j;
}
label = Return_State_Label (Ix2, Ix2);
}
void Heis_Bethe_State::Set_to_Label (string label_ref, const Vect<Vect<int> >& OriginIx2)
{
State_Label_Data labeldata = Read_State_Label (label_ref, OriginIx2);
label = label_ref;
Vect<Vect<int> > OriginIx2ordered = OriginIx2;
for (int i = 0; i < OriginIx2.size(); ++i) OriginIx2ordered[i].QuickSort();
// Set all Ix2 to OriginIx2:
for (int il = 0; il < OriginIx2.size(); ++il)
for (int i = 0; i < OriginIx2[il].size(); ++i) Ix2[il][i] = OriginIx2ordered[il][i];
// Now set the excitations:
for (int it = 0; it < labeldata.type.size(); ++it)
for (int iexc = 0; iexc < labeldata.nexc[it]; ++iexc)
for (int i = 0; i < labeldata.M[it]; ++i)
if (Ix2[labeldata.type[it] ][i] == labeldata.Ix2old[it][iexc]) {
Ix2[labeldata.type[it] ][i] = labeldata.Ix2exc[it][iexc];
}
// Now reorder the Ix2 to follow convention:
for (int il = 0; il < Ix2.size(); ++il) Ix2[il].QuickSort();
}
void Heis_Bethe_State::Set_Label_from_Ix2 (const Vect<Vect<int> >& OriginIx2)
{
// This function does not assume any ordering of the Ix2.
// ASSUMPTIONS:
// (*this) has a base already identical to base of OriginIx2
// First check that bases are consistent
if ((*this).chain.Nstrings != OriginIx2.size())
ABACUSerror("Inconsistent base sizes in Heis_Bethe_State::Set_Label_from_Ix2.");
// Then check that the filling at each level is equal
for (int il = 0; il < (*this).chain.Nstrings; ++il)
if ((*this).base.Nrap[il] != OriginIx2[il].size())
ABACUSerror("Inconsistent base filling in Heis_Bethe_State::Set_Label_from_Ix2.");
// Determine how many types of particles are present
int ntypes = 1; // level 0 always assumed present
for (int il = 1; il < chain.Nstrings; ++il) if (base.Nrap[il] > 0) ntypes++;
// Set the state label
Vect<int> type_ref(0, ntypes);
Vect<int> M_ref(ntypes);
M_ref[0] = OriginIx2[0].size();
type_ref[0] = 0;
int itype = 1;
for (int il = 1; il < chain.Nstrings; ++il)
if (base.Nrap[il] > 0) {
type_ref[itype] = il;
M_ref[itype++] = OriginIx2[il].size();
}
Vect<int> nexc_ref(0, ntypes);
Vect<Vect<int> > Ix2old_ref(ntypes);
Vect<Vect<int> > Ix2exc_ref(ntypes);
// Count nr of particle-holes at each level:
itype = -1;
for (int il = 0; il < chain.Nstrings; ++il) {
if (il == 0 || base.Nrap[il] > 0) {
itype++;
for (int alpha = 0; alpha < base.Nrap[il]; ++alpha)
if (!OriginIx2[il].includes(Ix2[il][alpha])) nexc_ref[itype] += 1;
Ix2old_ref[itype] = Vect<int>(ABACUS::max(nexc_ref[itype], 1));
Ix2exc_ref[itype] = Vect<int>(ABACUS::max(nexc_ref[itype], 1));
int nexccheck = 0;
for (int alpha = 0; alpha < base.Nrap[il]; ++alpha)
if (!OriginIx2[il].includes(Ix2[il][alpha])) Ix2exc_ref[itype][nexccheck++] = Ix2[il][alpha];
if (nexccheck != nexc_ref[itype]) {
cout << "il = " << il << "\titype = " << itype << "\tnexccheck = " << nexccheck << "\tnexc_ref[itype] = " << nexc_ref[itype] << endl;
cout << OriginIx2[il] << endl << Ix2[il] << endl;
ABACUSerror("Counting excitations wrong (1) in Heis_Bethe_State::Set_Label_from_Ix2.");
}
nexccheck = 0;
for (int alpha = 0; alpha < base.Nrap[il]; ++alpha)
if (!Ix2[il].includes (OriginIx2[il][alpha])) Ix2old_ref[itype][nexccheck++] = OriginIx2[il][alpha];
if (nexccheck != nexc_ref[itype]) {
ABACUSerror("Counting excitations wrong (2) in Heis_Bethe_State::Set_Label_from_Ix2.");
}
// Now order the Ix2old_ref and Ix2exc_ref:
Ix2old_ref[itype].QuickSort();
Ix2exc_ref[itype].QuickSort();
} // if (il == 0 || ...
} // for il
State_Label_Data labeldata(type_ref, M_ref, nexc_ref, Ix2old_ref, Ix2exc_ref);
label = Return_State_Label (labeldata, OriginIx2);
}
bool Heis_Bethe_State::Check_Symmetry ()
{
// Checks whether the I's are symmetrically distributed.
bool symmetric_state = true;
int arg, test1, test3;
for (int j = 0; j < chain.Nstrings; ++j) {
test1 = 0;
test3 = 0;
for (int alpha = 0; alpha < base[j]; ++alpha) {
arg = (Ix2[j][alpha] != chain.Nsites) ? Ix2[j][alpha] : 0; // since Ix2 = N is same as Ix2 = -N by periodicity, this is symmetric.
test1 += arg;
test3 += arg * arg * arg; // to make sure that all I's are symmetrical...
}
if (!(symmetric_state && (test1 == 0) && (test3 == 0))) symmetric_state = false;
}
return symmetric_state;
}
// SOLVING BETHE EQUATIONS:
void Heis_Bethe_State::Compute_diffsq ()
{
DP maxterm = 0.0;
int jmax, alphamax;
diffsq = 0.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
diffsq += pow(BE[j][alpha], 2.0);
if (pow(BE[j][alpha], 2.0)/chain.Nsites > maxterm) {
jmax = j;
alphamax = alpha;
maxterm = pow(BE[j][alpha], 2.0)/chain.Nsites;
}
}
diffsq /= DP(chain.Nsites);
}
void Heis_Bethe_State::Find_Rapidities (bool reset_rapidities)
{
// This function finds the rapidities of the eigenstate
lnnorm = -100.0; // sentinel value, recalculated if Newton method used in the last step of iteration.
Lambda lambda_ref(chain, base);
diffsq = 1.0;
if (reset_rapidities)
(*this).Set_Free_lambdas();
(*this).Compute_BE();
iter = 0;
iter_Newton = 0;
// Start with conventional iterations
DP iter_prec = 1.0e-2/DP(chain.Nsites);
DP iter_factor = 0.99;
DP iter_factor_extrap = 0.99;
clock_t extrap_start;
clock_t extrap_stop;
clock_t Newton_start;
clock_t Newton_stop;
int iter_extrap_start, iter_extrap_stop, iter_Newton_start, iter_Newton_stop;
DP diffsq_start;
DP diffsq_extrap_start, diffsq_extrap_stop, diffsq_Newton_start, diffsq_Newton_stop;
bool straight_improving = true;
bool extrap_improving = true;
bool Newton_improving = true;
DP diffsq_iter_aim = 1.0e-5;
bool info_findrap = false;
if (info_findrap) cout << "Find_Rapidities called for state with label " << (*this).label << endl;
while (diffsq > chain.prec && (iter < 500 && iter_Newton < 100
|| straight_improving || extrap_improving || Newton_improving)) {
// If we haven't reset, first try a few Newton steps...
if (!Newton_improving) diffsq_iter_aim *= 1.0e-5;
// Start with a few straight iterations:
if (!straight_improving) iter_factor *= 2.0/3.0;
else iter_factor = 0.99;
diffsq_start = diffsq;
do {
extrap_start = clock();
iter_extrap_start = iter;
diffsq_extrap_start = diffsq;
if (diffsq > iter_prec) (*this).Solve_BAE_straight_iter (iter_prec, 10, iter_factor);
extrap_stop = clock();
iter_extrap_stop = iter;
diffsq_extrap_stop = diffsq;
iter_prec = diffsq * 0.1;
if (info_findrap) cout << "Straight iter: iter " << iter << "\titer_factor " << iter_factor
<< "\tdiffsq " << diffsq << endl;
} while (diffsq > diffsq_iter_aim && !is_nan(diffsq) && diffsq_extrap_stop/diffsq_extrap_start < 0.01);
straight_improving = (diffsq < diffsq_start);
// Now try to extrapolate to infinite nr of iterations:
if (!extrap_improving) iter_factor_extrap *= 2.0/3.0;
else iter_factor_extrap = 0.99;
diffsq_start = diffsq;
if (diffsq > chain.prec)
do {
extrap_start = clock();
iter_extrap_start = iter;
diffsq_extrap_start = diffsq;
if (diffsq > iter_prec) (*this).Solve_BAE_extrap (iter_prec, 10, iter_factor_extrap);
extrap_stop = clock();
iter_extrap_stop = iter;
diffsq_extrap_stop = diffsq;
iter_prec = diffsq * 0.1;
if (info_findrap) cout << "Extrap: iter " << iter << "\tdiffsq " << diffsq << endl;
} while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_extrap_stop/diffsq_extrap_start < 0.01);
extrap_improving = (diffsq < diffsq_extrap_start);
// Now try Newton
diffsq_Newton_start = diffsq;
if (diffsq > chain.prec)
do {
Newton_start = clock();
iter_Newton_start = iter_Newton;
diffsq_Newton_start = diffsq;
if (diffsq > iter_prec) (*this).Solve_BAE_Newton (chain.prec, 5);
Newton_stop = clock();
iter_Newton_stop = iter_Newton;
diffsq_Newton_stop = diffsq;
if (info_findrap) cout << "Newton: iter_Newton " << iter_Newton << "\tdiffsq " << diffsq << endl;
} while (diffsq > chain.prec && !is_nan(diffsq) && diffsq_Newton_stop/diffsq_Newton_start < 0.01);
Newton_improving = (diffsq < diffsq_Newton_start);
// If none of the methods are improving the result, try the silk gloves...
if (!(straight_improving || extrap_improving || Newton_improving)) {
if (info_findrap) cout << "Before silk gloves: diffsq " << diffsq << endl;
(*this).Solve_BAE_with_silk_gloves (chain.prec, chain.Nsites, 0.9);
if (info_findrap) cout << "After silk gloves: diffsq " << diffsq << endl;
}
if (is_nan(diffsq)) {
(*this).Set_Free_lambdas(); // good start if we've messed up with Newton
diffsq = 1.0;
}
iter_prec *= 1.0e-4;
iter_prec = ABACUS::max(iter_prec, chain.prec);
} // while (diffsq > chain.prec && (iter < 300 && iter_Newton < 50...
// Check convergence:
conv = (diffsq < chain.prec && (*this).Check_Rapidities()) ? 1 : 0;
dev = (*this).String_delta();
return;
}
void Heis_Bethe_State::Solve_BAE_bisect (int j, int alpha, DP req_prec, int itermax)
{
// Finds the root lambda[j][alpha] to precision req_prec using bisection
DP prec_obtained = 1.0;
int niter_here = 0;
DP lambdajalphastart = lambda[j][alpha];
DP lambdamax = (fabs(chain.Delta) <= 1.0 ? pow(log(chain.Nsites), 1.1) : PI - 1.0e-4);
// Find values of lambda such that BE changes sign:
DP lambdaleft = -lambdamax;
DP lambdamid = 0.0;
DP lambdaright = lambdamax;
DP BEleft, BEmid, BEright;
lambda[j][alpha] = lambdaleft;
(*this).Compute_BE (j, alpha);
BEleft= BE[j][alpha];
lambda[j][alpha] = lambdaright;
(*this).Compute_BE (j, alpha);
BEright= BE[j][alpha];
if (BEleft * BEright > 0.0) {
lambda[j][alpha] = lambdajalphastart;
return;
}
while (prec_obtained > req_prec && niter_here < itermax) {
lambdamid = 0.5 * (lambdaleft + lambdaright);
lambda[j][alpha] = lambdamid;
(*this).Compute_BE (j, alpha);
BEmid = BE[j][alpha];
if (BEmid * BEleft < 0.0) { // root is to the left of mid
lambdaright = lambdamid;
BEright = BEmid;
}
else if (BEmid * BEright < 0.0) { // root is to the right of mid
lambdaleft = lambdamid;
BEleft = BEmid;
}
else if (BEmid * BEmid < req_prec) return;
else {
return; // this procedure has failed
}
prec_obtained = BEmid * BEmid;
niter_here++;
}
return;
}
void Heis_Bethe_State::Iterate_BAE (DP iter_factor)
{
for (int j = 0; j < chain.Nstrings; ++j) {
for (int alpha = 0; alpha < base[j]; ++alpha)
{
lambda[j][alpha] += iter_factor * (Iterate_BAE (j, alpha) - lambda[j][alpha]);
}
}
iter++;
(*this).Compute_BE();
(*this).Compute_diffsq();
}
void Heis_Bethe_State::Solve_BAE_straight_iter (DP straight_prec, int max_iter, DP iter_factor)
{
// This function attempts to get convergence diffsq <= straight_prec in at most max_iter steps.
// If this fails, the lambda's are reset to original values, defined as...
Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
diffsq = 1.0;
int iter_done_here = 0;
while ((diffsq > straight_prec) && (max_iter > iter_done_here)) {
(*this).Iterate_BAE(iter_factor);
iter_done_here++;
}
if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
(*this).Compute_BE();
diffsq = diffsq_ref;
}
}
void Heis_Bethe_State::Solve_BAE_extrap (DP extrap_prec, int max_iter_extrap, DP iter_factor)
{
// This function attempts to get convergence diffsq <= extrap_prec in at most max_iter_extrap steps.
// If this fails, the lambda's are reset to original values, defined as...
Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
diffsq = 1.0;
DP diffsq1, diffsq2, diffsq3, diffsq4;
int iter_done_here = 0;
Lambda lambda1(chain, base);
Lambda lambda2(chain, base);
Lambda lambda3(chain, base);
Lambda lambda4(chain, base);
Lambda lambda5(chain, base);
while ((diffsq > extrap_prec) && (max_iter_extrap > iter_done_here)) {
(*this).Iterate_BAE(iter_factor);
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda1[j][alpha] = lambda[j][alpha];
diffsq1 = diffsq;
(*this).Iterate_BAE(iter_factor);
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda2[j][alpha] = lambda[j][alpha];
diffsq2 = diffsq;
(*this).Iterate_BAE(iter_factor);
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda3[j][alpha] = lambda[j][alpha];
diffsq3 = diffsq;
(*this).Iterate_BAE(iter_factor);
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda4[j][alpha] = lambda[j][alpha];
diffsq4 = diffsq;
iter_done_here += 4;
// now extrapolate the result to infinite number of iterations...
if (diffsq4 < diffsq3 && diffsq3 < diffsq2 && diffsq2 < diffsq1 && diffsq1 < diffsq_ref) {
Vect_DP rap(0.0, 4);
Vect_DP oneoverP(0.0, 4);
DP deltalambda = 0.0;
for (int i = 0; i < 4; ++i) oneoverP[i] = 1.0/(1.0 + i*i);
for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) {
rap[0] = lambda1[j][alpha];
rap[1] = lambda2[j][alpha];
rap[2] = lambda3[j][alpha];
rap[3] = lambda4[j][alpha];
polint (oneoverP, rap, 0.0, lambda[j][alpha], deltalambda);
}
// Iterate once to stabilize result
(*this).Iterate_BAE(iter_factor);
}
} // ((diffsq > extrap_prec) && (max_iter_extrap > iter_done_here))
if ((diffsq >= diffsq_ref) || (is_nan(diffsq))) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
(*this).Compute_BE();
diffsq = diffsq_ref;
}
return;
}
void Heis_Bethe_State::Solve_BAE_with_silk_gloves (DP silk_prec, int max_iter_silk, DP iter_factor)
{
// Logic: do iterations, but change only one rapidity at a time.
// This is slow, and called only if straight_iter, extrap and Newton methods don't work.
// If this fails, the lambda's are reset to original values, defined as...
Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
diffsq = 1.0;
int iter_done_here = 0;
while ((diffsq > silk_prec) && (max_iter_silk > iter_done_here)) {
// Find the highest `deviant' rapidity
int jmax = 0;
int alphamax = 0;
DP absBEmax = 0.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
//cout << j << "\t" << alpha << "\t" << BE[j][alpha] << endl;
if (fabs(BE[j][alpha]) > absBEmax) {
jmax = j;
alphamax = alpha;
absBEmax = fabs(BE[j][alpha]);
}
}
Solve_BAE_bisect (jmax, alphamax, silk_prec, max_iter_silk);
iter_done_here++;
// and reset all important arrays.
(*this).Compute_diffsq();
}
if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j) for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
(*this).Compute_BE();
diffsq = diffsq_ref;
}
return;
}
void Heis_Bethe_State::Iterate_BAE_Newton ()
{
// does one step of a Newton method on the rapidities...
// Assumes that BE[j][alpha] have been computed
Vect_CX dlambda (0.0, base.Nraptot); // contains delta lambda computed from Newton's method
SQMat_CX Gaudin (0.0, base.Nraptot);
Vect_INT indx (base.Nraptot);
Lambda lambda_ref(chain, base);
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
(*this).Build_Reduced_Gaudin_Matrix (Gaudin);
int index = 0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
dlambda[index] = - BE[j][alpha] * chain.Nsites;
index++;
}
DP d;
try {
ludcmp_CX (Gaudin, indx, d);
lubksb_CX (Gaudin, indx, dlambda);
}
catch (Divide_by_zero) {
diffsq = log(-1.0); // reset to nan to stop Newton method
return;
}
// Regularize dlambda: max step is +-1.0 to prevent rapidity flying off into outer space.
for (int i = 0; i < base.Nraptot; ++i)
if (fabs(real(dlambda[i])) > 1.0) dlambda[i] = 0.0;//(real(dlambda[i]) > 0) ? 1.0 : -1.0;
index = 0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) {
lambda[j][alpha] = lambda_ref[j][alpha] + real(dlambda[index]);
index++;
}
(*this).Compute_BE();
(*this).Iterate_BAE(1.0);
(*this).Compute_diffsq();
// if we've converged, calculate the norm here, since the work has been done...
if (diffsq < chain.prec) {
lnnorm = 0.0;
for (int j = 0; j < base.Nraptot; j++) lnnorm += log(abs(Gaudin[j][j]));
}
iter_Newton++;
return;
}
void Heis_Bethe_State::Solve_BAE_Newton (DP Newton_prec, int max_iter_Newton)
{
// This function attempts to get convergence diffsq <= Newton_prec in at most max_iter_Newton steps.
// The results are accepted if diffsq has decreased,
// otherwise the lambda's are reset to original values, defined as...
Lambda lambda_ref(chain, base);
DP diffsq_ref = 1.0;
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda_ref[j][alpha] = lambda[j][alpha];
diffsq_ref = diffsq;
// Now begin solving...
int iter_done_here = 0;
while ((diffsq > Newton_prec) && (diffsq < 10.0) && (!is_nan(diffsq)) && (iter_done_here < max_iter_Newton)) {
(*this).Iterate_BAE_Newton();
iter_done_here++;
(*this).Iterate_BAE(1.0);
}
if ((diffsq > diffsq_ref) || (is_nan(diffsq))) {
// This procedure has failed. We reset everything to begin values.
for (int j = 0; j < chain.Nstrings; ++j)
for (int alpha = 0; alpha < base[j]; ++alpha) lambda[j][alpha] = lambda_ref[j][alpha];
(*this).Compute_BE();
diffsq = diffsq_ref;
}
return;
}
void Heis_Bethe_State::Compute_lnnorm ()
{
if (true || lnnorm == -100.0) { // else norm already calculated by Newton method
// Actually, compute anyway to increase accuracy !
SQMat_CX Gaudin_Red(base.Nraptot);
(*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);
complex<DP> lnnorm_check = lndet_LU_CX_dstry(Gaudin_Red);
lnnorm = real(lnnorm_check);
}
return;
}
void Heis_Bethe_State::Compute_Momentum ()
{
int sum_Ix2 = 0;
DP sum_M = 0.0;
for (int j = 0; j < chain.Nstrings; ++j) {
sum_M += 0.5 * (1.0 + chain.par[j]) * base[j];
for (int alpha = 0; alpha < base[j]; ++alpha) {
sum_Ix2 += Ix2[j][alpha];
}
}
iK = (chain.Nsites/2) * int(sum_M + 0.1) - (sum_Ix2/2); // + 0.1: for safety...
while (iK >= chain.Nsites) iK -= chain.Nsites;
while (iK < 0) iK += chain.Nsites;
K = PI * sum_M - PI * sum_Ix2/chain.Nsites;
while (K >= 2.0*PI) K -= 2.0*PI;
while (K < 0.0) K += 2.0*PI;
return;
}
void Heis_Bethe_State::Compute_All (bool reset_rapidities) // solves BAE, computes E, K and lnnorm
{
(*this).Find_Rapidities (reset_rapidities);
if (conv == 1) {
(*this).Compute_Energy ();
(*this).Compute_Momentum ();
(*this).Compute_lnnorm ();
}
return;
}
void Heis_Bethe_State::Set_to_Closest_Matching_Ix2_fixed_Base (const Heis_Bethe_State& StateToMatch)
{
// Given a state with given Ix2 distribution, set the Ix2 to closest match.
// The base of (*this) is fixed, and does not necessarily match that of StateToMatch.
if ((*this).chain != StateToMatch.chain)
ABACUSerror("Heis_Bethe_State::Find_Closest_Matching_Ix2_fixed_Base: "
"trying to match Ix2 for two states with different chains.");
// Check level by level, match quantum numbers from center up.
for (int il = 0; il < chain.Nstrings; ++il) {
// Careful: if parity of Nraps is different, the Ix2 live on different lattices!
int Ix2shift = ((StateToMatch.base.Nrap[il] - (*this).base.Nrap[il]) % 2);
if (StateToMatch.base.Nrap[il] >= (*this).base.Nrap[il]) {
// We match the rapidities from the center towards the sides.
int ashift = (StateToMatch.base.Nrap[il] - (*this).base.Nrap[il])/2;
for (int a = 0; a < (*this).base.Nrap[il]; ++a)
(*this).Ix2[il][a] = StateToMatch.Ix2[il][a + ashift] + Ix2shift;
}
else { // There are less Ix2 in StateToMatch than in (*this)
// We start by filling all the (*this) Ix2 symmetrically from the middle.
// We thereafter start from the left and identify half of the StateToMatch Ix2 with the leftmost (*this)Ix2
// and then do the same thing from the right.
for (int a = 0; a < (*this).base.Nrap[il]; ++a)
(*this).Ix2[il][a] = -(*this).base.Nrap[il] + 1 + 2*a;
int nleft = StateToMatch.base.Nrap[il]/2;
for (int a = 0; a < nleft; ++a)
if (StateToMatch.Ix2[il][a] - Ix2shift < (*this).Ix2[il][a])
(*this).Ix2[il][a] = StateToMatch.Ix2[il][a] - Ix2shift;
for (int a = 0; a < StateToMatch.base.Nrap[il] - 1 - nleft; ++a)
if (StateToMatch.Ix2[il][StateToMatch.base.Nrap[il] - 1 - a] - Ix2shift
> (*this).Ix2[il][(*this).base.Nrap[il] - 1 - a])
(*this).Ix2[il][(*this).base.Nrap[il] - 1 - a] = (StateToMatch.Ix2[il][StateToMatch.base.Nrap[il] - 1 - a]
- Ix2shift);
}
} // for il
}
std::ostream& operator<< (std::ostream& s, const Heis_Bethe_State& state)
{
// sends all the state data to output stream
s << endl << "******** Chain with Delta = " << state.chain.Delta
<< " Nsites = " << state.chain.Nsites << " Mdown = " << state.base.Mdown
<< ": eigenstate with label " << state.label << endl
<< "E = " << state.E << " K = " << state.K << " iK = " << state.iK << " lnnorm = " << state.lnnorm << endl
<< "conv = " << state.conv << " dev = " << state.dev << " iter = " << state.iter
<< " iter_Newton = " << state.iter_Newton << "\tdiffsq " << state.diffsq << endl;
for (int j = 0; j < state.chain.Nstrings; ++j) {
if (state.base.Nrap[j] > 0) {
s << "Type " << j << " Str_L = " << state.chain.Str_L[j] << " par = " << state.chain.par[j]
<< " M_j = " << state.base.Nrap[j]
<< " Ix2_infty = " << state.base.Ix2_infty[j] << " Ix2_min = " << state.base.Ix2_min[j]
<< " Ix2_max = " << state.base.Ix2_max[j] << endl;
Vect_INT qnumbers(state.base.Nrap[j]);
Vect_DP rapidities(state.base.Nrap[j]);
for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) {
qnumbers[alpha] = state.Ix2[j][alpha];
rapidities[alpha] = state.lambda[j][alpha];
}
qnumbers.QuickSort();
rapidities.QuickSort();
s << "Ix2 quantum numbers: " << endl;
for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << qnumbers[alpha] << " ";
s << endl;
s << "Rapidities: " << endl;
for (int alpha = 0; alpha < state.base.Nrap[j]; ++alpha) s << rapidities[alpha] << " ";
s << endl;
}
}
s << endl;
return s;
}
} // namespace ABACUS
Reference in New Issue View Git Blame Copy Permalink