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

783 lines
25 KiB
C++
Raw Permalink Blame History

/**********************************************************
This software is part of J.-S. Caux's ABACUS library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: LiebLin_Bethe_State.cc
Purpose: Definitions for LiebLin_Bethe_State class.
***********************************************************/
#include "ABACUS.h"
using namespace std;
namespace ABACUS {
//***************************************************************************************************
// Function definitions: class LiebLin_Bethe_State
LiebLin_Bethe_State::LiebLin_Bethe_State ()
: c_int (0.0), L(0.0), cxL(0.0), N(0),
Ix2_available(Vect<int>(0, 1)), index_first_hole_to_right (Vect<int>(0,1)),
displacement (Vect<int>(0,1)), Ix2(Vect<int>(0, 1)), lambdaoc(Vect<DP>(0.0, 1)),
S(Vect<DP>(0.0, 1)), dSdlambdaoc(Vect<DP>(0.0, 1)), diffsq(0.0), prec(ITER_REQ_PREC_LIEBLIN),
conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0)
{
stringstream Nout; Nout << N; label = Nout.str() + LABELSEP + ABACUScoding[0] + LABELSEP;
}
LiebLin_Bethe_State::LiebLin_Bethe_State (DP c_int_ref, DP L_ref, int N_ref)
: c_int(c_int_ref), L(L_ref), cxL(c_int_ref * L_ref), N(N_ref),
Ix2_available(Vect<int>(0, 2)), index_first_hole_to_right (Vect<int>(0,N)),
displacement (Vect<int>(0,N)), Ix2(Vect<int>(0, N)), lambdaoc(Vect<DP>(0.0, N)),
S(Vect<DP>(0.0, N)), dSdlambdaoc(Vect<DP>(0.0, N)),
diffsq(0.0), prec(ABACUS::max(1.0, 1.0/(c_int * c_int)) * ITER_REQ_PREC_LIEBLIN),
conv(0), iter_Newton(0), E(0.0), iK(0), K(0.0), lnnorm(-100.0)
{
if (c_int < 0.0) ABACUSerror("You must use a positive interaction parameter !");
if (N < 0) ABACUSerror("Particle number must be strictly positive.");
stringstream Nout; Nout << N; label = Nout.str() + LABELSEP + ABACUScoding[0] + LABELSEP;
// Set quantum numbers to ground-state configuration:
for (int i = 0; i < N; ++i) Ix2[i] = -(N-1) + 2*i;
Vect<int> OriginIx2 = Ix2;
(*this).Set_Label_from_Ix2 (OriginIx2);
}
LiebLin_Bethe_State& LiebLin_Bethe_State::operator= (const LiebLin_Bethe_State& RefState)
{
if (this != &RefState) {
c_int = RefState.c_int;
L = RefState.L;
cxL = RefState.cxL;
N = RefState.N;
label = RefState.label;
Ix2_available = RefState.Ix2_available;
index_first_hole_to_right = RefState.index_first_hole_to_right;
displacement = RefState.displacement;
Ix2 = RefState.Ix2;
lambdaoc = RefState.lambdaoc;
S = RefState.S;
dSdlambdaoc = RefState.dSdlambdaoc;
diffsq = RefState.diffsq;
prec = RefState.prec;
conv = RefState.conv;
iter_Newton = RefState.iter_Newton;
E = RefState.E;
iK = RefState.iK;
K = RefState.K;
lnnorm = RefState.lnnorm;
}
return(*this);
}
void LiebLin_Bethe_State::Set_to_Label (string label_ref, const Vect<int>& OriginStateIx2)
{
State_Label_Data labeldata = Read_State_Label (label_ref, OriginStateIx2);
if (N != labeldata.M[0]) {
cout << label_ref << endl;
cout << labeldata.M << endl;
ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: N != M[0].");
}
if (N != OriginStateIx2.size()) {
cout << label_ref << endl;
cout << labeldata.M << endl;
ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: "
"N != OriginStateIx2.size().");
}
label = label_ref;
Vect<int> OriginStateIx2ordered = OriginStateIx2;
OriginStateIx2ordered.QuickSort();
// Set all Ix2 to OriginState's Ix2:
for (int i = 0; i < N; ++i) Ix2[i] = OriginStateIx2ordered[i];
// Now set the excitations:
for (int iexc = 0; iexc < labeldata.nexc[0]; ++iexc)
for (int i = 0; i < N; ++i)
if (Ix2[i] == labeldata.Ix2old[0][iexc]) Ix2[i] = labeldata.Ix2exc[0][iexc];
// Now reorder the Ix2 to follow convention:
Ix2.QuickSort();
(*this).Set_Label_from_Ix2 (OriginStateIx2ordered);
}
void LiebLin_Bethe_State::Set_to_Label (string label_ref)
{
// This function assumes that OriginState is the ground state.
Vect<int> OriginStateIx2(N);
for (int i = 0; i < N; ++i) OriginStateIx2[i] = -N + 1 + 2*i;
(*this).Set_to_Label(label_ref, OriginStateIx2);
}
void LiebLin_Bethe_State::Set_Label_from_Ix2 (const Vect<int>& OriginStateIx2)
{
// This function does not assume any ordering of the Ix2.
if (N != OriginStateIx2.size()) ABACUSerror("N != OriginStateIx2.size() in Set_Label_from_Ix2.");
// Set the state label:
Vect<int> type_ref(0,1);
Vect<int> M_ref(N, 1);
Vect<int> nexc_ref(0, 1);
// Count nr of particle-holes:
for (int i = 0; i < N; ++i) if (!OriginStateIx2.includes(Ix2[i])) nexc_ref[0] += 1;
Vect<Vect<int> > Ix2old_ref(1);
Vect<Vect<int> > Ix2exc_ref(1);
Ix2old_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
Ix2exc_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
int nexccheck = 0;
for (int i = 0; i < N; ++i)
if (!OriginStateIx2.includes(Ix2[i])) Ix2exc_ref[0][nexccheck++] = Ix2[i];
if (nexccheck != nexc_ref[0])
ABACUSerror("Counting excitations wrong (1) in LiebLin_Bethe_State::Set_Label_from_Ix2");
nexccheck = 0;
for (int i = 0; i < N; ++i)
if (!Ix2.includes (OriginStateIx2[i])) Ix2old_ref[0][nexccheck++] = OriginStateIx2[i];
if (nexccheck != nexc_ref[0]) {
cout << "nexc_ref[0] = " << nexc_ref[0] << "\tnexccheck = " << nexccheck << endl;
cout << OriginStateIx2 << endl;
cout << Ix2 << endl;
cout << nexc_ref[0] << endl;
cout << Ix2old_ref[0] << endl;
cout << Ix2exc_ref[0] << endl;
ABACUSerror("Counting excitations wrong (2) in LiebLin_Bethe_State::Set_Label_from_Ix2");
}
// Now order the Ix2old_ref and Ix2exc_ref:
Ix2old_ref[0].QuickSort();
Ix2exc_ref[0].QuickSort();
State_Label_Data labeldata(type_ref, M_ref, nexc_ref, Ix2old_ref, Ix2exc_ref);
label = Return_State_Label (labeldata, OriginStateIx2);
}
void LiebLin_Bethe_State::Set_Label_Internals_from_Ix2 (const Vect<int>& OriginStateIx2)
{
if (N != OriginStateIx2.size())
ABACUSerror("N != OriginStateIx2.size() in Set_Label_Internals_from_Ix2.");
Vect<int> OriginStateIx2ordered = OriginStateIx2;
OriginStateIx2ordered.QuickSort();
// Set the state label:
Vect<int> type_ref(0,1);
Vect<int> M_ref(N, 1);
Vect<int> nexc_ref(0, 1);
// Count nr of particle-holes:
for (int i = 0; i < N; ++i) if (!OriginStateIx2.includes(Ix2[i])) nexc_ref[0] += 1;
Vect<Vect<int> > Ix2old_ref(1);
Vect<Vect<int> > Ix2exc_ref(1);
Ix2old_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
Ix2exc_ref[0] = Vect<int>(ABACUS::max(nexc_ref[0],1));
int nexccheck = 0;
for (int i = 0; i < N; ++i)
if (Ix2[i] != OriginStateIx2ordered[i]) {
Ix2old_ref[0][nexccheck] = OriginStateIx2ordered[i];
Ix2exc_ref[0][nexccheck++] = Ix2[i];
}
State_Label_Data labeldata(type_ref, M_ref, nexc_ref, Ix2old_ref, Ix2exc_ref);
label = Return_State_Label (labeldata, OriginStateIx2);
// Construct the Ix2_available vector: we give one more quantum number on left and right:
int navailable = 2 + (ABACUS::max(Ix2.max(), OriginStateIx2.max())
- ABACUS::min(Ix2.min(), OriginStateIx2.min()))/2 - N + 1;
Ix2_available = Vect<int>(navailable);
index_first_hole_to_right = Vect<int>(N);
// First set Ix2_available to all holes from left
for (int i = 0; i < Ix2_available.size(); ++i)
Ix2_available[i] = ABACUS::min(Ix2.min(), OriginStateIx2.min()) - 2 + 2*i;
// Now shift according to Ix2 of OriginState:
for (int j = 0; j < N; ++j) {
int i = 0;
while (Ix2_available[i] < OriginStateIx2ordered[j]) i++;
// We now have Ix2_available[i] == OriginStateIx2[j].
// Shift all Ix2_available to the right of this by 2;
for (int i1 = i; i1 < navailable; ++i1) Ix2_available[i1] += 2;
index_first_hole_to_right[j] = i;
}
// Ix2_available and index_first_hole_to_right are now fully defined.
// Now set displacement vector:
displacement = Vect<int>(0, N);
// Set displacement vector from the Ix2:
for (int j = 0; j < N; ++j) {
if (Ix2[j] < OriginStateIx2ordered[j]) {
// Ix2[j] must be equal to some OriginState_Ix2_available[i]
// for i < OriginState_index_first_hole_to_right[j]
while (Ix2[j] != Ix2_available[index_first_hole_to_right[j] + displacement[j] ]) {
if (index_first_hole_to_right[j] + displacement[j] == 0) {
cout << label << endl << j << endl << OriginStateIx2 << endl
<< Ix2 << endl << Ix2_available
<< endl << index_first_hole_to_right << endl << displacement << endl;
ABACUSerror("Going down too far in Set_Label_Internals...");
}
displacement[j]--;
}
}
if (Ix2[j] > OriginStateIx2ordered[j]) {
// Ix2[j] must be equal to some Ix2_available[i] for i >= index_first_hole_to_right[j]
displacement[j] = 1; // start with this value to prevent segfault
while (Ix2[j] != Ix2_available[index_first_hole_to_right[j] - 1 + displacement[j] ]) {
if (index_first_hole_to_right[j] + displacement[j] == Ix2_available.size() - 1) {
cout << label << endl << j << endl << OriginStateIx2 << endl
<< Ix2 << endl << Ix2_available
<< endl << index_first_hole_to_right << endl << displacement << endl;
ABACUSerror("Going up too far in Set_Label_Internals...");
}
displacement[j]++;
}
}
}
}
bool LiebLin_Bethe_State::Check_Admissibility (char whichDSF)
{
// For testing with restricted Hilbert space:
//if (Ix2.min() < -13 || Ix2.max() > 13) return(false);
return(true);
}
void LiebLin_Bethe_State::Find_Rapidities (bool reset_rapidities)
{
// This function finds the rapidities of the eigenstate
// sentinel value, recalculated if Newton method used in the last step of iteration.
lnnorm = -100.0;
diffsq = 1.0;
if (reset_rapidities) (*this).Set_Free_lambdaocs();
iter_Newton = 0;
DP damping = 1.0;
DP diffsq_prev = 1.0e+6;
Vect_DP RHSBAE (0.0, N); // contains RHS of BAEs
Vect_DP dlambdaoc (0.0, N); // contains delta lambdaoc computed from Newton's method
SQMat_DP Gaudin (0.0, N);
Vect_INT indx (N);
while (diffsq > prec && !is_nan(diffsq) && iter_Newton < 100) {
(*this).Iterate_BAE_Newton(damping, RHSBAE, dlambdaoc, Gaudin, indx);
if (diffsq > diffsq_prev && damping > 0.5) damping /= 2.0;
else if (diffsq < diffsq_prev) damping = 1.0;
diffsq_prev = diffsq;
}
conv = ((diffsq < prec) && (*this).Check_Rapidities()) ? 1 : 0;
if (!conv) {
cout << "Alert! State " << label << " did not converge... diffsq " << diffsq
<< "\titer_Newton " << iter_Newton << (*this) << endl;
}
return;
}
bool LiebLin_Bethe_State::Check_Rapidities()
{
bool nonan = true;
for (int j = 0; j < N; ++j) nonan *= !is_nan(lambdaoc[j]);
return nonan;
}
/**
For the repulsive Lieb-Liniger model, all eigenstates have real rapidities
(there are no strings). All string deviations are thus set to zero.
*/
DP LiebLin_Bethe_State::String_delta()
{
return(0.0);
}
/**
Checks whether the quantum numbers are symmetrically distributed: \f$\{ I \} = \{ -I \}\f$.
*/
bool LiebLin_Bethe_State::Check_Symmetry ()
{
bool symmetric_state = true;
Vect<int> Ix2check = Ix2;
Ix2check.QuickSort();
for (int alpha = 0; alpha <= N/2; ++alpha)
symmetric_state = symmetric_state && (Ix2check[alpha] == -Ix2check[N - 1 - alpha]);
return(symmetric_state);
}
/**
This function computes the log of the norm of a `LiebLin_Bethe_State` instance.
This is obtained from the Gaudin-Korepin norm formula
\f{eqnarray*}{
\langle 0 | \prod_{j=0}^{N-1} C(\lambda_j) \prod_{j=0}^{N-1} B(\lambda_j) | 0 \rangle
= \prod_{j=0}^{N-2}\prod_{k=j+1}^{N-1}
\frac{(\lambda_j - \lambda_k )^2 + c^2}{(\lambda_j - \lambda_k)^2}
\times ~\mbox{Det}_N G
\f}
in which \f$G\f$ is the Gaudin matrix computed by
LiebLin_Bethe_State::Build_Reduced_Gaudin_Matrix.
*/
void LiebLin_Bethe_State::Compute_lnnorm ()
{
if (lnnorm == -100.0) { // else Gaudin part already calculated by Newton method
SQMat_DP Gaudin_Red(N);
(*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);
complex<DP> lnnorm_CX = lndet_LU_dstry(Gaudin_Red);
if (abs(imag(lnnorm_CX)) > 1.0e-6) {
cout << "lnnorm_CX = " << lnnorm_CX << endl;
ABACUSerror("Gaudin norm of LiebLin_Bethe_State should be real and positive.");
}
lnnorm = real(lnnorm_CX);
// Add the pieces outside of Gaudin determinant
for (int j = 0; j < N - 1; ++j) for (int k = j+1; k < N; ++k)
lnnorm += log(1.0 + 1.0/pow(lambdaoc[j] - lambdaoc[k], 2.0));
}
return;
}
/**
This function solves the Bethe equations, computes the energy, momentum and state norm.
It calls LiebLin_Bethe_State::Find_Rapidities, LiebLin_Bethe_State::Compute_Energy,
LiebLin_Bethe_State::Compute_Momentum and LiebLin_Bethe_State::Compute_lnnorm.
Assumptions:
- the set of quantum numbers is admissible.
*/
void LiebLin_Bethe_State::Compute_All (bool reset_rapidities)
{
(*this).Find_Rapidities (reset_rapidities);
if (conv == 1) {
(*this).Compute_Energy ();
(*this).Compute_Momentum ();
(*this).Compute_lnnorm ();
}
return;
}
void LiebLin_Bethe_State::Set_Free_lambdaocs()
{
if (cxL >= 1.0)
for (int a = 0; a < N; ++a) lambdaoc[a] = PI * Ix2[a]/cxL;
// For small values of c, use better approximation using approximate
// zeroes of Hermite polynomials: see Gaudin eqn 4.71.
if (cxL < 1.0) {
DP oneoversqrtcLN = 1.0/pow(cxL * N, 0.5);
for (int a = 0; a < N; ++a) lambdaoc[a] = oneoversqrtcLN * PI * Ix2[a];
}
return;
}
void LiebLin_Bethe_State::Iterate_BAE (DP damping)
{
// does one step of simple iterations
DP sumtheta = 0.0;
Vect_DP dlambdaoc (0.0, N);
for (int j = 0; j < N; ++j) {
sumtheta = 0.0;
for (int k = 0; k < N; ++k) sumtheta += atan((lambdaoc[j] - lambdaoc[k]));
sumtheta *= 2.0;
dlambdaoc[j] = damping * ((PI*Ix2[j] - sumtheta)/cxL - lambdaoc[j]);
}
diffsq = 0.0;
for (int i = 0; i < N; ++i) {
lambdaoc[i] += dlambdaoc[i];
// Normalize the diffsq by the typical value of the rapidities:
if (cxL > 1.0) diffsq += dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
else diffsq += cxL * cxL * dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
}
diffsq /= DP(N);
return;
}
void LiebLin_Bethe_State::Iterate_BAE_S (DP damping)
{
// This is essentially Newton's method but only in one variable.
// The logic is that the derivative of the LHS of the BE_j w/r to lambdaoc_j is much larger
// than with respect to lambdaoc_l with l != j.
Vect_DP dlambdaoc (0.0, N);
// Start by calculating S and dSdlambdaoc:
for (int j = 0; j < N; ++j) {
S[j] = 0.0;
for (int k = 0; k < N; ++k) S[j] += atan((lambdaoc[j] - lambdaoc[k]));
S[j] *= 2.0/cxL;
dSdlambdaoc[j] = 0.0;
for (int k = 0; k < N; ++k)
dSdlambdaoc[j] += 1.0/((lambdaoc[j] - lambdaoc[k]) * (lambdaoc[j] - lambdaoc[k]) + 1.0);
dSdlambdaoc[j] *= 2.0/(PI * cxL);
dlambdaoc[j] = (PI*Ix2[j]/cxL - S[j] + lambdaoc[j] * dSdlambdaoc[j])
/(1.0 + dSdlambdaoc[j]) - lambdaoc[j];
}
diffsq = 0.0;
for (int i = 0; i < N; ++i) {
lambdaoc[i] += damping * dlambdaoc[i];
// Normalize the diffsq by the typical value of the rapidities:
if (cxL > 1.0) diffsq += dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
else diffsq += cxL * cxL * dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
}
diffsq /= DP(N);
iter_Newton++;
return;
}
/**
Performs one step of the matrix Newton method on the rapidities.
*/
void LiebLin_Bethe_State::Iterate_BAE_Newton (DP damping, Vect_DP& RHSBAE, Vect_DP& dlambdaoc, SQMat_DP& Gaudin, Vect_INT& indx)
{
DP sumtheta = 0.0;
// for large |lambda|, use atan (lambda) = sgn(lambda) pi/2 - atan(1/lambda)
int atanintshift = 0;
DP lambdahere = 0.0;
// Compute the RHS of the BAEs:
for (int j = 0; j < N; ++j) {
sumtheta = 0.0;
atanintshift = 0;
for (int k = 0; k < N; ++k)
if (j != k) { // otherwise 0
if (fabs(lambdahere = lambdaoc[j] - lambdaoc[k]) < 1.0) { // use straight atan
sumtheta += atan(lambdahere);
}
else { // for large rapidities, use dual form of atan, extracting pi/2 factors
atanintshift += sgn_DP(lambdahere);
sumtheta -= atan(1.0/lambdahere);
}
}
sumtheta *= 2.0;
RHSBAE[j] = cxL * lambdaoc[j] + sumtheta - PI*(Ix2[j] - atanintshift);
}
(*this).Build_Reduced_Gaudin_Matrix (Gaudin);
for (int j = 0; j < N; ++j) dlambdaoc[j] = - RHSBAE[j];
DP d;
ludcmp (Gaudin, indx, d);
lubksb (Gaudin, indx, dlambdaoc);
bool ordering_changed = false;
for (int j = 0; j < N-1; ++j)
if (lambdaoc[j] + dlambdaoc[j] > lambdaoc[j+1] + dlambdaoc[j+1]) ordering_changed = true;
// To prevent Newton from diverging, we limit the size of the rapidity changes.
// The leftmost and rightmost rapidities can grow by one order of magnitude per iteration step.
if (ordering_changed) {
// We explicitly ensure that the ordering remains correct after the iteration step.
bool ordering_still_changed = false;
DP maxdlambdaoc = 0.0;
do {
ordering_still_changed = false;
if (dlambdaoc[0] < 0.0 && fabs(dlambdaoc[0])
> (maxdlambdaoc = 10.0*ABACUS::max(fabs(lambdaoc[0]), fabs(lambdaoc[N-1]))))
dlambdaoc[0] = -maxdlambdaoc;
if (lambdaoc[0] + dlambdaoc[0] > lambdaoc[1] + dlambdaoc[1]) {
dlambdaoc[0] = 0.25 * (lambdaoc[1] + dlambdaoc[1] - lambdaoc[0] ); // max quarter distance
ordering_still_changed = true;
}
if (dlambdaoc[N-1] > 0.0 && fabs(dlambdaoc[N-1])
> (maxdlambdaoc = 10.0*ABACUS::max(fabs(lambdaoc[0]), fabs(lambdaoc[N-1]))))
dlambdaoc[N-1] = maxdlambdaoc;
if (lambdaoc[N-1] + dlambdaoc[N-1] < lambdaoc[N-2] + dlambdaoc[N-2]) {
dlambdaoc[N-1] = 0.25 * (lambdaoc[N-2] + dlambdaoc[N-2] - lambdaoc[N-1]);
ordering_still_changed = true;
}
for (int j = 1; j < N-1; ++j) {
if (lambdaoc[j] + dlambdaoc[j] > lambdaoc[j+1] + dlambdaoc[j+1]) {
dlambdaoc[j] = 0.25 * (lambdaoc[j+1] + dlambdaoc[j+1] - lambdaoc[j]);
ordering_still_changed = true;
}
if (lambdaoc[j] + dlambdaoc[j] < lambdaoc[j-1] + dlambdaoc[j-1]) {
dlambdaoc[j] = 0.25 * (lambdaoc[j-1] + dlambdaoc[j-1] - lambdaoc[j]);
ordering_still_changed = true;
}
}
} while (ordering_still_changed);
}
diffsq = 0.0;
for (int i = 0; i < N; ++i) {
// Normalize the diffsq by the typical value of the rapidities:
if (cxL > 1.0) diffsq += dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
else diffsq += cxL * cxL * dlambdaoc[i] * dlambdaoc[i]/(lambdaoc[i] * lambdaoc[i] + 1.0e-6);
}
diffsq /= DP(N);
if (ordering_changed) diffsq = 1.0; // reset if Newton wanted to change ordering
for (int j = 0; j < N; ++j) lambdaoc[j] += damping * dlambdaoc[j];
iter_Newton++;
// if we've converged, calculate the norm here, since the work has been done...
if (diffsq < prec && !ordering_changed) {
lnnorm = 0.0;
for (int j = 0; j < N; j++) lnnorm += log(fabs(Gaudin[j][j]));
// Add the pieces outside of Gaudin determinant
for (int j = 0; j < N - 1; ++j)
for (int k = j+1; k < N; ++k)
lnnorm += log(1.0 + 1.0/pow(lambdaoc[j] - lambdaoc[k], 2.0));
}
return;
}
void LiebLin_Bethe_State::Compute_Energy()
{
E = 0.0;
for (int j = 0; j < N; ++j) E += lambdaoc[j] * lambdaoc[j];
E *= c_int * c_int;
}
void LiebLin_Bethe_State::Compute_Momentum()
{
iK = 0;
for (int j = 0; j < N; ++j) {
iK += Ix2[j];
}
if (iK % 2) {
cout << Ix2 << endl;
cout << iK << "\t" << iK % 2 << endl;
ABACUSerror("Sum of Ix2 is not even: inconsistency.");
}
iK /= 2; // sum of Ix2 is guaranteed even.
K = 2.0 * iK * PI/L;
}
/**
The fundamental scattering kernel for Lieb-Liniger,
\f[
K (\tilde{\lambda}) = \frac{ 2 }{\tilde{\lambda}^2 + 1}
\f]
given two indices for the rapidities as arguments.
*/
DP LiebLin_Bethe_State::Kernel (int a, int b)
{
return(2.0/(pow(lambdaoc[a] - lambdaoc[b], 2.0) + 1.0));
}
/**
The fundamental scattering kernel for Lieb-Liniger,
\f[
K (\tilde{\lambda}) = \frac{ 2 }{\tilde{\lambda}^2 + 1}
\f]
given the rapidity difference as argument.
*/
DP LiebLin_Bethe_State::Kernel (DP lambdaoc_ref)
{
return(2.0/(lambdaoc_ref * lambdaoc_ref + 1.0));
}
/**
This function constructs the Gaudin matrix \f$G\f$
which is defined as the Hessian of the Yang-Yang action \f$S^{YY}\f$,
\f[
G_{jk} = \delta_{jk} \left\{ cL + \sum_{l} K (\tilde{\lambda}_j - \tilde{\lambda}_l) \right\}
- K (\tilde{\lambda}_j - \tilde{\lambda}_k).
\f]
*/
void LiebLin_Bethe_State::Build_Reduced_Gaudin_Matrix (SQMat<DP>& Gaudin_Red)
{
if (Gaudin_Red.size() != N)
ABACUSerror("Passing matrix of wrong size in Build_Reduced_Gaudin_Matrix.");
DP sum_Kernel = 0.0;
for (int j = 0; j < N; ++j)
for (int k = 0; k < N; ++k) {
if (j == k) {
sum_Kernel = 0.0;
for (int kp = 0; kp < N; ++kp)
if (j != kp) sum_Kernel += Kernel (lambdaoc[j] - lambdaoc[kp]);
Gaudin_Red[j][k] = cxL + sum_Kernel;
}
else Gaudin_Red[j][k] = - Kernel (lambdaoc[j] - lambdaoc[k]);
}
return;
}
/**
For a summetric `LiebLin_Bethe_State`, this function computes the Gaudin matrix
as an \f$N/2 \times N/2\f$ matrix to accelerate computations.
*/
void LiebLin_Bethe_State::Build_Reduced_BEC_Quench_Gaudin_Matrix (SQMat<DP>& Gaudin_Red)
{
// Passing a matrix of dimension N/2
if (N % 2 != 0) ABACUSerror("Choose a state with even numer of particles please");
// Check Parity invariant
bool ck = true;
for (int j = 0; j < N/2; ++j){ if(Ix2[j] != - Ix2[N-j-1]) ck = false;}
if (!ck) ABACUSerror("Choose a parity invariant state please");
if (Gaudin_Red.size() != N/2)
ABACUSerror("Passing matrix of wrong size in Build_Reduced_Gaudin_Matrix.");
DP sum_Kernel = 0.0;
for (int j = 0; j < N/2; ++j)
for (int k = 0; k < N/2; ++k) {
if (j == k) {
sum_Kernel = 0.0;
for (int kp = N/2; kp < N; ++kp)
if (j + N/2 != kp)
sum_Kernel += Kernel (lambdaoc[j+N/2] - lambdaoc[kp])
+ Kernel (lambdaoc[j+N/2] + lambdaoc[kp]);
Gaudin_Red[j][k] = cxL + sum_Kernel;
}
else Gaudin_Red[j][k] = - (Kernel (lambdaoc[j+ N/2] - lambdaoc[k+ N/2])
+ Kernel (lambdaoc[j+ N/2] + lambdaoc[k+ N/2]) );
}
return;
}
/**
This function changes the Ix2 of a given state by annihilating a particle and hole
pair specified by ipart and ihole (counting from the left, starting with index 0).
*/
void LiebLin_Bethe_State::Annihilate_ph_pair (int ipart, int ihole,
const Vect<int>& OriginStateIx2)
{
State_Label_Data currentdata = Read_State_Label ((*this).label, OriginStateIx2);
if (ipart >= currentdata.nexc[0])
ABACUSerror("Particle label too large in LiebLin_Bethe_State::Annihilate_ph_pair.");
if (ihole >= currentdata.nexc[0])
ABACUSerror("Hole label too large in LiebLin_Bethe_State::Annihilate_ph_pair.");
// Simply remove the given pair:
Vect<int> type_new = currentdata.type;
Vect<int> M_new = currentdata.M;
Vect<int> nexc_new = currentdata.nexc;
nexc_new[0] -= 1; // we drill one more particle-hole pair at level 0
int ntypespresent = 1; // only one type for LiebLin
Vect<Vect<int> > Ix2old_new(ntypespresent);
Vect<Vect<int> > Ix2exc_new(ntypespresent);
for (int it = 0; it < ntypespresent; ++it)
Ix2old_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
for (int it = 0; it < ntypespresent; ++it)
Ix2exc_new[it] = Vect<int>(ABACUS::max(nexc_new[it],1));
// Copy earlier data in, leaving out ipart and ihole:
for (int it = 0; it < ntypespresent; ++it) {
for (int i = 0; i < nexc_new[it]; ++i) {
Ix2old_new[it][i] = currentdata.Ix2old[it][i + (i >= ihole)];
Ix2exc_new[it][i] = currentdata.Ix2exc[it][i + (i >= ipart)];
}
}
State_Label_Data newdata (type_new, M_new, nexc_new, Ix2old_new, Ix2exc_new);
(*this).Set_to_Label (Return_State_Label(newdata, OriginStateIx2));
}
/**
Performs a spatial inversion of a `LiebLin_Bethe_State`, mapping all quantum numbers
and rapidities to minus themselves.
*/
void LiebLin_Bethe_State::Parity_Flip ()
{
Vect_INT Ix2buff = Ix2;
Vect_DP lambdaocbuff = lambdaoc;
for (int i = 0; i < N; ++i) Ix2[i] = -Ix2buff[N - 1 - i];
for (int i = 0; i < N; ++i) lambdaoc[i] = -lambdaocbuff[N - 1 - i];
iK = -iK;
K = -K;
}
/**
Send information about a `LiebLin_Bethe_State` to the output stream.
*/
std::ostream& operator<< (std::ostream& s, const LiebLin_Bethe_State& state)
{
s << endl << "******** State for c = " << state.c_int << " L = " << state.L
<< " N = " << state.N << " with label " << state.label << " ********" << endl;
s << "Ix2:" << endl;
for (int j = 0; j < state.N; ++j) s << state.Ix2[j] << " ";
s << endl << "lambdaocs:" << endl;
for (int j = 0; j < state.N; ++j) s << state.lambdaoc[j] << " ";
s << endl << "conv = " << state.conv << " iter_Newton = " << state.iter_Newton << endl;
s << "E = " << state.E << " iK = " << state.iK << " K = " << state.K
<< " lnnorm = " << state.lnnorm << endl;
return(s);
}
} // namespace ABACUS
Reference in New Issue View Git Blame Copy Permalink