jscaux/ABACUS-v1
3
0
Fork 0
Code Issues 1 Pull Requests Releases Wiki Activity

Include Usage_Example files; work on LiebLin (code, documentation)

Browse Source
This commit is contained in:
J.-S. Caux
2018-02-15 10:08:36 +01:00
parent f654f45171
commit 2677b62120
10 changed files with 549 additions and 153 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/LIEBLIN/LiebLin_Bethe_State.cc
+131 -41
View File
@@ -24,18 +24,18 @@ namespace ABACUS {
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)
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)),
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)
@@ -93,7 +93,8 @@ namespace ABACUS {
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().");
ABACUSerror("Trying to set an incorrect label on LiebLin_Bethe_State: "
"N != OriginStateIx2.size().");
}
label = label_ref;
@@ -106,7 +107,8 @@ namespace ABACUS {
// 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];
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();
@@ -143,10 +145,13 @@ namespace ABACUS {
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");
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];
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;
@@ -167,7 +172,8 @@ namespace ABACUS {
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.");
if (N != OriginStateIx2.size())
ABACUSerror("N != OriginStateIx2.size() in Set_Label_Internals_from_Ix2.");
Vect<int> OriginStateIx2ordered = OriginStateIx2;
OriginStateIx2ordered.QuickSort();
@@ -207,7 +213,8 @@ namespace ABACUS {
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;
// 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;
}
@@ -218,10 +225,12 @@ namespace ABACUS {
// 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]
// 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
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...");
}
@@ -233,7 +242,8 @@ namespace ABACUS {
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
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...");
}
@@ -245,7 +255,8 @@ namespace ABACUS {
bool LiebLin_Bethe_State::Check_Admissibility (char whichDSF)
{
//if (Ix2.min() < -13 || Ix2.max() > 13) return(false); // For testing with restricted Hilbert space
// For testing with restricted Hilbert space:
//if (Ix2.min() < -13 || Ix2.max() > 13) return(false);
return(true);
}
@@ -253,7 +264,8 @@ namespace ABACUS {
{
// This function finds the rapidities of the eigenstate
lnnorm = -100.0; // sentinel value, recalculated if Newton method used in the last step of iteration.
// sentinel value, recalculated if Newton method used in the last step of iteration.
lnnorm = -100.0;
diffsq = 1.0;
@@ -289,15 +301,20 @@ namespace ABACUS {
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); // no strings (thus no deviations) in replusive LiebLin
return(0.0);
}
/**
Checks whether the quantum numbers are symmetrically distributed: \f$\{ I \} = \{ -I \}\f$.
*/
bool LiebLin_Bethe_State::Check_Symmetry ()
{
// Checks whether the I's are symmetrically distributed.
bool symmetric_state = true;
Vect<int> Ix2check = Ix2;
@@ -309,6 +326,18 @@ namespace ABACUS {
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
@@ -317,7 +346,13 @@ namespace ABACUS {
(*this).Build_Reduced_Gaudin_Matrix(Gaudin_Red);
lnnorm = real(lndet_LU_dstry(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
@@ -329,7 +364,15 @@ namespace ABACUS {
return;
}
void LiebLin_Bethe_State::Compute_All (bool reset_rapidities) // solves BAE, computes E, K and lnnorm
/**
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) {
@@ -345,7 +388,8 @@ namespace ABACUS {
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.
// 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];
@@ -403,7 +447,8 @@ namespace ABACUS {
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];
dlambdaoc[j] = (PI*Ix2[j]/cxL - S[j] + lambdaoc[j] * dSdlambdaoc[j])
/(1.0 + dSdlambdaoc[j]) - lambdaoc[j];
}
@@ -421,16 +466,18 @@ namespace ABACUS {
return;
}
/**
Performs one step of the matrix Newton method on the rapidities.
*/
void LiebLin_Bethe_State::Iterate_BAE_Newton (DP damping)
{
// does one step of a Newton method on the rapidities...
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);
DP sumtheta = 0.0;
int atanintshift = 0; // for large |lambda|, use atan (lambda) = sgn(lambda) pi/2 - atan(1/lambda)
// 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:
@@ -463,11 +510,13 @@ namespace ABACUS {
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;
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.
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 {
@@ -552,16 +601,38 @@ namespace ABACUS {
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)
{
@@ -575,7 +646,8 @@ namespace ABACUS {
if (j == k) {
sum_Kernel = 0.0;
for (int kp = 0; kp < N; ++kp) if (j != kp) sum_Kernel += Kernel (lambdaoc[j] - lambdaoc[kp]);
for (int kp = 0; kp < N; ++kp)
if (j != kp) sum_Kernel += Kernel (lambdaoc[j] - lambdaoc[kp]);
Gaudin_Red[j][k] = cxL + sum_Kernel;
}
@@ -586,6 +658,11 @@ namespace ABACUS {
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
@@ -612,7 +689,8 @@ namespace ABACUS {
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]);
sum_Kernel += Kernel (lambdaoc[j+N/2] - lambdaoc[kp])
+ Kernel (lambdaoc[j+N/2] + lambdaoc[kp]);
Gaudin_Red[j][k] = cxL + sum_Kernel;
}
@@ -625,10 +703,13 @@ namespace ABACUS {
}
void LiebLin_Bethe_State::Annihilate_ph_pair (int ipart, int ihole, const Vect<int>& OriginStateIx2)
/**
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)
{
// 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).
State_Label_Data currentdata = Read_State_Label ((*this).label, OriginStateIx2);
@@ -645,8 +726,10 @@ namespace ABACUS {
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));
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) {
@@ -661,9 +744,12 @@ namespace ABACUS {
(*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 ()
{
// For simplicity, we don't redo base_id, type_id, id.
Vect_INT Ix2buff = Ix2;
Vect_DP lambdaocbuff = lambdaoc;
for (int i = 0; i < N; ++i) Ix2[i] = -Ix2buff[N - 1 - i];
@@ -672,16 +758,20 @@ namespace ABACUS {
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 << 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;
s << "E = " << state.E << " iK = " << state.iK << " K = " << state.K
<< " lnnorm = " << state.lnnorm << endl;
return(s);
}
src/LIEBLIN/ln_Density_ME.cc
+26 -9
View File
@@ -19,35 +19,51 @@ using namespace ABACUS;
namespace ABACUS {
/**
Given a index \f$j\f$ together with a bra and ket respectively parametrized
by sets of rapidities \f$\{ \mu \}\f$ and \f$\{ \lambda \}\f$,
this function returns the product (see e.g.
[J.-S. Caux et al, JSTAT P01008 (2007)](https://doi.org/10.1088/1742-5468/2007/01/P01008))
\f[
V_j^+ \equiv \prod_l
\frac{\tilde{\mu}_l - \tilde{\lambda}_j + i}{\tilde{\lambda}_l - \tilde{\lambda}_j + i}.
\f]
*/
complex<DP> Fn_V (int j, LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate)
{
complex<DP> result = 1.0;
for (int m = 0; m < lstate.N; ++m) {
result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II)/(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II);
result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II)
/(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II);
}
return(result);
}
/**
Computes the log of the density operator \f$\hat{\rho}(x = 0)\f$ matrix element
between lstate (bra) and rstate (ket). If lstate == rstate, density matrix element = N/L;
if momentum difference is zero but states are different, then matrix element is zero
(but this function then returns the log artificially set to -200.0).
*/
complex<DP> ln_Density_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate)
{
// Computes the log of the density operator \rho(x = 0) matrix element between lstate and rstate.
// If we have lstate == rstate, density matrix element = N/L:
if (lstate.Ix2 == rstate.Ix2) return(log(lstate.N/lstate.L));
// If momentum difference is zero but states are different, then matrix element is zero:
else if (lstate.iK == rstate.iK) return(-200.0); // so exp(.) is zero
SQMat_DP one_plus_U (0.0, lstate.N);
Vect_CX Vplus (0.0, lstate.N); // contains V^+_j; V^-_j is the conjugate
Vect_CX ln_Fn_Prod (0.0, lstate.N); // product_{m\neq j} (\mu_m - \lambdaoc_j)/(\lambdaoc_m - \lambdaoc_j)
// prod_{m\neq j} (\mu_m - \lambdaoc_j)/(\lambdaoc_m - \lambdaoc_j):
Vect_CX ln_Fn_Prod (0.0, lstate.N);
Vect_DP rKern (0.0, lstate.N); // K(lambdaoc_j - lambdaoc_(p == arbitrary))
// "Phantom" rapidity in ME expression. Choice doesn't matter, see 1990_Slavnov_TMP_82
// after (3.8). Choose rapidity around the middle.
//int p = 0;
int p = rstate.N/2-1; // choice doesn't matter, see 1990_Slavnov_TMP_82 after (3.8). Choose rapidity around the middle.
int p = rstate.N/2-1;
DP Kout = lstate.K - rstate.K;
@@ -62,8 +78,9 @@ namespace ABACUS {
for (int a = 0; a < lstate.N; ++a)
for (int b = 0; b < lstate.N; ++b)
one_plus_U[a][b] = (a == b ? 1.0 : 0.0) + 0.5 * ((lstate.lambdaoc[a] - rstate.lambdaoc[a])/imag(Vplus[a]))
* real(exp(ln_Fn_Prod[a])) * (rstate.Kernel(a,b) - rKern[b]); // BUGRISK: why real here?
one_plus_U[a][b] = (a == b ? 1.0 : 0.0)
+ 0.5 * ((lstate.lambdaoc[a] - rstate.lambdaoc[a])/imag(Vplus[a]))
* real(exp(ln_Fn_Prod[a])) * (rstate.Kernel(a,b) - rKern[b]);
complex<DP> ln_ddalpha_sigma = lndet_LU_dstry(one_plus_U);
src/LIEBLIN/ln_Psi_ME.cc
+33 -17
View File
@@ -19,62 +19,78 @@ using namespace ABACUS;
namespace ABACUS {
/**
Given a index \f$j\f$ together with a bra and ket respectively parametrized
by sets of rapidities \f$\{ \mu \}\f$ and \f$\{ \lambda \}\f$,
this function returns the product (see e.g.
[J.-S. Caux et al, JSTAT P01008 (2007)](https://doi.org/10.1088/1742-5468/2007/01/P01008))
\f[
V_j^+ \equiv \prod_l
\frac{\tilde{\mu}_l - \tilde{\lambda}_j + i}{\tilde{\lambda}_l - \tilde{\lambda}_j + i}.
\f]
*/
complex<DP> Fn_V_Psi (int j, LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate)
{
// lstate has one particle less than rstate
complex<DP> result_num = 1.0;
complex<DP> result_den = 1.0;
complex<DP> result = 1.0;
for (int m = 0; m < lstate.N; ++m)
result_num *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II);
result *= (lstate.lambdaoc[m] - rstate.lambdaoc[j] + II)
/(rstate.lambdaoc[m] - rstate.lambdaoc[j] + II);
result /= (rstate.lambdaoc[rstate.N - 1] - rstate.lambdaoc[j] + II);
for (int m = 0; m < rstate.N; ++m)
result_den *= (rstate.lambdaoc[m] - rstate.lambdaoc[j] + II);
return(result_num/result_den);
return(result);
}
/**
Computes the log of the annihilation operator matrix element between lstate and rstate.
*/
complex<DP> ln_Psi_ME (LiebLin_Bethe_State& lstate, LiebLin_Bethe_State& rstate)
{
// Computes the log of the annihilation operator matrix element between lstate and rstate.
if (lstate.N + 1 != rstate.N) ABACUSerror("Wrong particle numbers in left and right states for Psi FF.");
if (lstate.N + 1 != rstate.N)
ABACUSerror("Wrong particle numbers in left and right states for Psi FF.");
SQMat_DP U_Psi (0.0, lstate.N);
Vect_CX Vplus (0.0, lstate.N); // contains V^+_j; V^-_j is the conjugate
Vect_DP Fn_Prod (0.0, lstate.N); // product_{m} (\mu_m - \lambdaoc_j)/product_{m neq j} (\lambdaoc_m - \lambdaoc_j)
// product_{m} (\mu_m - \lambdaoc_j)/product_{m neq j} (\lambdaoc_m - \lambdaoc_j)
Vect_DP Fn_Prod (0.0, lstate.N);
Vect_DP rKern (0.0, lstate.N); // K(lambdaoc_j - lambdaoc_(p == N))
int p = rstate.N - 1;
for (int a = 0; a < lstate.N; ++a) {
Vplus[a] = Fn_V_Psi (a, lstate, rstate);
Fn_Prod[a] = (lstate.lambdaoc[a] - rstate.lambdaoc[a])/(rstate.lambdaoc[rstate.N - 1] - rstate.lambdaoc[a]);
Fn_Prod[a] = (lstate.lambdaoc[a] - rstate.lambdaoc[a])
/(rstate.lambdaoc[rstate.N - 1] - rstate.lambdaoc[a]);
for (int m = 0; m < lstate.N; ++m)
if (m != a) Fn_Prod[a] *= (lstate.lambdaoc[m] - rstate.lambdaoc[a])/(rstate.lambdaoc[m] - rstate.lambdaoc[a]);
if (m != a) Fn_Prod[a] *= (lstate.lambdaoc[m] - rstate.lambdaoc[a])
/(rstate.lambdaoc[m] - rstate.lambdaoc[a]);
rKern[a] = rstate.Kernel (a, p);
}
for (int a = 0; a < lstate.N; ++a)
for (int b = 0; b < lstate.N; ++b)
U_Psi[a][b] = (a == b ? 2.0 * imag(Vplus[a]) : 0.0) + Fn_Prod[a] * (rstate.Kernel(a,b) - rKern[b]);
U_Psi[a][b] = (a == b ? 2.0 * imag(Vplus[a]) : 0.0)
+ Fn_Prod[a] * (rstate.Kernel(a,b) - rKern[b]);
complex<DP> ln_det_U_Psi = lndet_LU_dstry(U_Psi);
complex<DP> ln_prod_lambdaocsq_plus_one = 0.0;
for (int a = 0; a < rstate.N - 1; ++a)
for (int b = a; b < rstate.N; ++b)
ln_prod_lambdaocsq_plus_one += log(complex<DP>((rstate.lambdaoc[a] - rstate.lambdaoc[b])
* (rstate.lambdaoc[a] - rstate.lambdaoc[b])) + 1.0);
ln_prod_lambdaocsq_plus_one +=
log(complex<DP>((rstate.lambdaoc[a] - rstate.lambdaoc[b])
* (rstate.lambdaoc[a] - rstate.lambdaoc[b])) + 1.0);
complex<DP> ln_prod_lambdaoca_min_mub = 0.0;
for (int a = 0; a < rstate.N; ++a)
for (int b = 0; b < lstate.N; ++b)
ln_prod_lambdaoca_min_mub += log(complex<DP>(rstate.lambdaoc[a] - lstate.lambdaoc[b]));
return (ln_det_U_Psi + 0.5 * log(lstate.c_int) + ln_prod_lambdaocsq_plus_one - ln_prod_lambdaoca_min_mub
return (ln_det_U_Psi + 0.5 * log(lstate.c_int)
+ ln_prod_lambdaocsq_plus_one - ln_prod_lambdaoca_min_mub
- 0.5 * (lstate.lnnorm + rstate.lnnorm));
}