112 lines
4.0 KiB
C++
112 lines
4.0 KiB
C++
/**********************************************************
|
|
|
|
This software is part of J.-S. Caux's ABACUS library.
|
|
|
|
Copyright (c) J.-S. Caux.
|
|
|
|
-----------------------------------------------------------
|
|
|
|
File: LiebLin_Twisted_ln_Overlap.cc
|
|
|
|
Purpose: Calculates Nikita Slavnov's determinant, case RHS is Bethe
|
|
and LHS is twisted Bethe
|
|
|
|
***********************************************************/
|
|
|
|
#include "ABACUS.h"
|
|
|
|
using namespace std;
|
|
using namespace ABACUS;
|
|
|
|
namespace ABACUS {
|
|
|
|
complex<DP> Kernel_Twisted (complex<DP> expbeta, complex<DP> lambdaoc)
|
|
{
|
|
return(1.0/(lambdaoc + II) - expbeta/(lambdaoc - II));
|
|
}
|
|
|
|
complex<DP> Kernel_Twisted (complex<DP> expbeta, DP lambdaoc)
|
|
{
|
|
return(1.0/(lambdaoc + II) - expbeta/(lambdaoc - II));
|
|
}
|
|
|
|
complex<DP> Fn_V (int j, int sign, Vect<complex<DP> >& lstate_lambdaoc, LiebLin_Bethe_State& rstate)
|
|
{
|
|
complex<DP> result_num = 1.0;
|
|
complex<DP> result_den = 1.0;
|
|
|
|
complex<DP> signcx = complex<DP>(sign);
|
|
|
|
for (int m = 0; m < rstate.N; ++m) {
|
|
result_num *= (lstate_lambdaoc[m] - rstate.lambdaoc[j] + signcx * II);
|
|
result_den *= (rstate.lambdaoc[m] - rstate.lambdaoc[j] + signcx * II);
|
|
}
|
|
|
|
return(result_num/result_den);
|
|
}
|
|
|
|
complex<DP> Fn_V_Nikita (int j, int sign, Vect<complex<DP> >& lstate_lambdaoc, LiebLin_Bethe_State& rstate)
|
|
{
|
|
// To match with Nikita's new conventions
|
|
return(1.0/Fn_V (j, -sign, lstate_lambdaoc, rstate));
|
|
}
|
|
|
|
complex<DP> LiebLin_Twisted_ln_Overlap (DP expbeta, Vect<DP> lstate_lambdaoc,
|
|
DP lstate_lnnorm, LiebLin_Bethe_State& rstate)
|
|
{
|
|
Vect<complex<DP> > lstate_lambdaoc_CX(lstate_lambdaoc.size());
|
|
for (int i = 0; i < lstate_lambdaoc.size(); ++i) lstate_lambdaoc_CX[i] = complex<DP>(lstate_lambdaoc[i]);
|
|
return(LiebLin_Twisted_ln_Overlap (complex<DP>(expbeta), lstate_lambdaoc_CX, lstate_lnnorm, rstate));
|
|
}
|
|
|
|
complex<DP> LiebLin_Twisted_ln_Overlap (complex<DP> expbeta, Vect<complex<DP> > lstate_lambdaoc,
|
|
DP lstate_lnnorm, LiebLin_Bethe_State& rstate)
|
|
{
|
|
// Computes the log of the overlap between the left state and the Bethe rstate
|
|
|
|
// If momentum difference is zero but states are different, then form factor is zero:
|
|
|
|
SQMat_CX one_plus_U (0.0, rstate.N);
|
|
|
|
Vect_CX Vplus_Nikita (0.0, rstate.N); // contains V^+_j
|
|
Vect_CX Vminus_Nikita (0.0, rstate.N); // contains V^-_j
|
|
Vect_CX Fn_Prod (0.0, rstate.N); // product_{m\neq j} (\mu_m - \lambda_j)/(\lambda_m - \lambda_j)
|
|
Vect_CX rKern (0.0, rstate.N); // K(lambda_j - lambda_p)
|
|
|
|
int p = 0;
|
|
|
|
complex<DP> Kout = lstate_lambdaoc.sum() - rstate.K;
|
|
|
|
for (int a = 0; a < rstate.N; ++a) {
|
|
Vplus_Nikita[a] = Fn_V_Nikita (a, 1, lstate_lambdaoc, rstate);
|
|
Vminus_Nikita[a] = Fn_V_Nikita (a, -1, lstate_lambdaoc, rstate);
|
|
Fn_Prod[a] = 1.0;
|
|
for (int m = 0; m < rstate.N; ++m)
|
|
if (m != a) Fn_Prod[a] *= (lstate_lambdaoc[m] - rstate.lambdaoc[a])/(rstate.lambdaoc[m] - rstate.lambdaoc[a]);
|
|
rKern[a] = Kernel_Twisted (expbeta, rstate.lambdaoc[p] - rstate.lambdaoc[a]);
|
|
}
|
|
|
|
for (int a = 0; a < rstate.N; ++a)
|
|
for (int b = 0; b < rstate.N; ++b)
|
|
one_plus_U[a][b] = (a == b ? 1.0 : 0.0)
|
|
+ ((lstate_lambdaoc[a] - rstate.lambdaoc[a])/(1.0/Vplus_Nikita[a] - 1.0/Vminus_Nikita[a]))
|
|
* Fn_Prod[a] * (Kernel_Twisted(expbeta, rstate.lambdaoc[a] - rstate.lambdaoc[b]) - rKern[b]);
|
|
|
|
complex<DP> ln_ddalpha_sigma = lndet_LU_CX_dstry(one_plus_U);
|
|
|
|
complex<DP> ln_prod_V = 0.0;
|
|
for (int a = 0; a < rstate.N; ++a) ln_prod_V += log(expbeta * Vplus_Nikita[a]/Vminus_Nikita[a] - 1.0);
|
|
ln_prod_V -= 0.5 * rstate.N * log(expbeta);
|
|
|
|
complex<DP> ln_prod_2 = 0.0;
|
|
for (int a = 0; a < rstate.N; ++a)
|
|
for (int b = 0; b < rstate.N; ++b)
|
|
ln_prod_2 += log((lstate_lambdaoc[a] - rstate.lambdaoc[b] - II)/(rstate.lambdaoc[a] - lstate_lambdaoc[b]));
|
|
|
|
ln_ddalpha_sigma += ln_prod_V + ln_prod_2 - log(Vplus_Nikita[p] - expbeta * Vminus_Nikita[p]);
|
|
|
|
return (log(1.0 - exp(-II * Kout)) + ln_ddalpha_sigma - 0.5 * (lstate_lnnorm + rstate.lnnorm));
|
|
}
|
|
|
|
}
|