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

569 lines
21 KiB
C++
Raw Permalink Blame History

/**********************************************************
This software is part of J.-S. Caux's ABACUS library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: ln_Sz_ME_XXZ.cc
Purpose: compute the S^z matrix elemment for XXZ
***********************************************************/
#include "ABACUS.h"
using namespace std;
using namespace ABACUS;
namespace ABACUS {
inline complex<DP> ln_Fn_F (XXZ_Bethe_State& B, int k, int beta, int b)
{
complex<DP> ans = 0.0;
complex<DP> prod_temp = 1.0;
int counter = 0;
int arg = 0;
int absarg = 0;
int par_comb_1, par_comb_2;
for (int j = 0; j < B.chain.Nstrings; ++j) {
par_comb_1 = B.chain.par[j] == B.chain.par[k] ? 1 : 0;
par_comb_2 = B.chain.par[k] == B.chain.par[j] ? 0 : B.chain.par[k];
for (int alpha = 0; alpha < B.base.Nrap[j]; ++alpha) {
for (int a = 1; a <= B.chain.Str_L[j]; ++a) {
if (!((j == k) && (alpha == beta) && (a == b))) {
arg = B.chain.Str_L[j] - B.chain.Str_L[k] - 2 * (a - b);
absarg = abs(arg);
prod_temp *= ((B.sinhlambda[j][alpha] * B.coshlambda[k][beta]
- B.coshlambda[j][alpha] * B.sinhlambda[k][beta])
* (B.chain.co_n_anis_over_2[absarg] * par_comb_1
- sgn_int(arg) * B.chain.si_n_anis_over_2[absarg] * par_comb_2)
+ II * (B.coshlambda[j][alpha] * B.coshlambda[k][beta]
- B.sinhlambda[j][alpha] * B.sinhlambda[k][beta])
* (sgn_int(arg) * B.chain.si_n_anis_over_2[absarg] * par_comb_1
+ B.chain.co_n_anis_over_2[absarg] * par_comb_2));
}
if (counter++ > 100) { // we do at most 100 products before taking a log
ans += log(prod_temp);
prod_temp = 1.0;
counter = 0;
}
}}}
return(ans + log(prod_temp));
}
inline complex<DP> ln_Fn_G (XXZ_Bethe_State& A, XXZ_Bethe_State& B, int k, int beta, int b)
{
complex<DP> ans = 0.0;
complex<DP> prod_temp = 1.0;
int counter = 0;
int arg = 0;
int absarg = 0;
int par_comb_1, par_comb_2;
for (int j = 0; j < A.chain.Nstrings; ++j) {
par_comb_1 = A.chain.par[j] == B.chain.par[k] ? 1 : 0;
par_comb_2 = B.chain.par[k] == A.chain.par[j] ? 0 : B.chain.par[k];
for (int alpha = 0; alpha < A.base.Nrap[j]; ++alpha) {
for (int a = 1; a <= A.chain.Str_L[j]; ++a) {
arg = A.chain.Str_L[j] - B.chain.Str_L[k] - 2 * (a - b);
absarg = abs(arg);
prod_temp *= ((A.sinhlambda[j][alpha] * B.coshlambda[k][beta]
- A.coshlambda[j][alpha] * B.sinhlambda[k][beta])
* (A.chain.co_n_anis_over_2[absarg] * par_comb_1
- sgn_int(arg) * A.chain.si_n_anis_over_2[absarg] * par_comb_2)
+ II * (A.coshlambda[j][alpha] * B.coshlambda[k][beta]
- A.sinhlambda[j][alpha] * B.sinhlambda[k][beta])
* (sgn_int(arg) * A.chain.si_n_anis_over_2[absarg] * par_comb_1
+ A.chain.co_n_anis_over_2[absarg] * par_comb_2));
if (counter++ > 100) { // we do at most 100 products before taking a log
ans += log(prod_temp);
prod_temp = 1.0;
counter = 0;
}
}}}
return(ans + log(prod_temp));
}
inline complex<DP> Fn_K (XXZ_Bethe_State& A, int j, int alpha, int a, XXZ_Bethe_State& B, int k, int beta, int b)
{
int arg1 = A.chain.Str_L[j] - B.chain.Str_L[k] - 2 * (a - b);
int absarg1 = abs(arg1);
int arg2 = arg1 + 2;
int absarg2 = abs(arg2);
return(4.0/(
((A.sinhlambda[j][alpha] * B.coshlambda[k][beta] - A.coshlambda[j][alpha] * B.sinhlambda[k][beta])
* (A.chain.co_n_anis_over_2[absarg1] * (1.0 + A.chain.par[j] * B.chain.par[k])
- sgn_int(arg1) * A.chain.si_n_anis_over_2[absarg1] * (B.chain.par[k] - A.chain.par[j]))
+ II * (A.coshlambda[j][alpha] * B.coshlambda[k][beta] - A.sinhlambda[j][alpha] * B.sinhlambda[k][beta])
* (sgn_int(arg1) * A.chain.si_n_anis_over_2[absarg1] * (1.0 + A.chain.par[j] * B.chain.par[k])
+ A.chain.co_n_anis_over_2[absarg1] * (B.chain.par[k] - A.chain.par[j])) )
*
((A.sinhlambda[j][alpha] * B.coshlambda[k][beta] - A.coshlambda[j][alpha] * B.sinhlambda[k][beta])
* (A.chain.co_n_anis_over_2[absarg2] * (1.0 + A.chain.par[j] * B.chain.par[k])
- sgn_int(arg2) * A.chain.si_n_anis_over_2[absarg2] * (B.chain.par[k] - A.chain.par[j]))
+ II * (A.coshlambda[j][alpha] * B.coshlambda[k][beta] - A.sinhlambda[j][alpha] * B.sinhlambda[k][beta])
* (sgn_int(arg2) * A.chain.si_n_anis_over_2[absarg2] * (1.0 + A.chain.par[j] * B.chain.par[k])
+ A.chain.co_n_anis_over_2[absarg2] * (B.chain.par[k] - A.chain.par[j])) )
));
}
inline complex<DP> Fn_L (XXZ_Bethe_State& A, int j, int alpha, int a, XXZ_Bethe_State& B, int k, int beta, int b)
{
return (sinh(2.0 * (A.lambda[j][alpha] - B.lambda[k][beta]
+ 0.5 * II * B.chain.anis * (A.chain.Str_L[j] - B.chain.Str_L[k] - 2.0 * (a - b - 0.5))
+ 0.25 * II * PI * complex<DP>(-A.chain.par[j] + B.chain.par[k])))
* pow(Fn_K (A, j, alpha, a, B, k, beta, b), 2.0));
}
// Version without phases:
complex<DP> ln_Sz_ME (XXZ_Bethe_State& A, XXZ_Bethe_State& B)
{
// This function returns the natural log of the S^z operator matrix element.
// The A and B states can contain strings.
// Check that the two states refer to the same XXZ_Chain
if (A.chain != B.chain) {
ABACUSerror("Incompatible XXZ_Chains in Sz matrix element.");
}
// Check that A and B are compatible: same Mdown
if (A.base.Mdown != B.base.Mdown) {
cout << "Bra state: " << endl << A << endl;
cout << "Ket state: " << endl << B << endl;
cout << A.base.Mdown << "\t" << B.base.Mdown << endl;
ABACUSerror("Incompatible Mdown between the two states in Sz matrix element!");
}
if (A.iK == B.iK && (A.label != B.label))
return(-300.0); // matrix element identically vanishes
// Compute the sinh and cosh of rapidities
A.Compute_sinhlambda();
A.Compute_coshlambda();
B.Compute_sinhlambda();
B.Compute_coshlambda();
// Some convenient arrays
Lambda re_ln_Fn_F_B_0(B.chain, B.base);
Lambda im_ln_Fn_F_B_0(B.chain, B.base);
Lambda re_ln_Fn_G_0(B.chain, B.base);
Lambda im_ln_Fn_G_0(B.chain, B.base);
Lambda re_ln_Fn_G_2(B.chain, B.base);
Lambda im_ln_Fn_G_2(B.chain, B.base);
complex<DP> ln_prod1 = 0.0;
complex<DP> ln_prod2 = 0.0;
complex<DP> ln_prod3 = 0.0;
complex<DP> ln_prod4 = 0.0;
for (int i = 0; i < A.chain.Nstrings; ++i)
for (int alpha = 0; alpha < A.base.Nrap[i]; ++alpha)
for (int a = 1; a <= A.chain.Str_L[i]; ++a)
ln_prod1 += log(norm(sinh(A.lambda[i][alpha] + 0.5 * II * A.chain.anis * (A.chain.Str_L[i] + 1.0 - 2.0 * a - 1.0)
+ 0.25 * II * PI * (1.0 - A.chain.par[i]))));
for (int i = 0; i < B.chain.Nstrings; ++i)
for (int alpha = 0; alpha < B.base.Nrap[i]; ++alpha)
for (int a = 1; a <= B.chain.Str_L[i]; ++a)
if (norm(sinh(B.lambda[i][alpha] + 0.5 * II * B.chain.anis * (B.chain.Str_L[i] + 1.0 - 2.0 * a - 1.0)
+ 0.25 * II * PI * (1.0 - B.chain.par[i]))) > 100.0 * MACHINE_EPS_SQ)
ln_prod2 += log(norm(sinh(B.lambda[i][alpha] + 0.5 * II * B.chain.anis * (B.chain.Str_L[i] + 1.0 - 2.0 * a - 1.0)
+ 0.25 * II * PI * (1.0 - B.chain.par[i]))));
// Define the F ones earlier...
complex<DP> ln_Fn_F_B_0;
for (int j = 0; j < B.chain.Nstrings; ++j) {
for (int alpha = 0; alpha < B.base.Nrap[j]; ++alpha) {
ln_Fn_F_B_0 = ln_Fn_F(B, j, alpha, 0);
re_ln_Fn_F_B_0[j][alpha] = real(ln_Fn_F_B_0);
im_ln_Fn_F_B_0[j][alpha] = imag(ln_Fn_F_B_0);
re_ln_Fn_G_0[j][alpha] = real(ln_Fn_G(A, B, j, alpha, 0));
im_ln_Fn_G_0[j][alpha] = imag(ln_Fn_G(A, B, j, alpha, 0));
re_ln_Fn_G_2[j][alpha] = real(ln_Fn_G(A, B, j, alpha, 2));
im_ln_Fn_G_2[j][alpha] = imag(ln_Fn_G(A, B, j, alpha, 2));
}
}
DP logabssinzeta = log(abs(sin(A.chain.anis)));
// Define regularized products in prefactors
for (int j = 0; j < A.chain.Nstrings; ++j)
for (int alpha = 0; alpha < A.base.Nrap[j]; ++alpha)
for (int a = 1; a <= A.chain.Str_L[j]; ++a) {
ln_prod3 += ln_Fn_F (A, j, alpha, a - 1);
}
ln_prod3 -= A.base.Mdown * log(abs(sin(A.chain.anis)));
for (int k = 0; k < B.chain.Nstrings; ++k)
for (int beta = 0; beta < B.base.Nrap[k]; ++beta)
for (int b = 1; b <= B.chain.Str_L[k]; ++b) {
if (b == 1) ln_prod4 += re_ln_Fn_F_B_0[k][beta];
else if (b > 1) ln_prod4 += ln_Fn_F(B, k, beta, b - 1);
}
ln_prod4 -= B.base.Mdown * log(abs(sin(B.chain.anis)));
// Now proceed to build the Hm2P matrix
SQMat_CX Hm2P(0.0, A.base.Mdown);
int index_a = 0;
int index_b = 0;
complex<DP> sum1 = 0.0;
complex<DP> sum2 = 0.0;
complex<DP> prod_num = 0.0;
complex<DP> Fn_K_0_G_0 = 0.0;
complex<DP> Prod_powerN = 0.0;
complex<DP> Fn_K_1_G_2 = 0.0;
complex<DP> two_over_A_sinhlambda_sq_plus_sinzetaover2sq;
for (int j = 0; j < A.chain.Nstrings; ++j) {
for (int alpha = 0; alpha < A.base.Nrap[j]; ++alpha) {
for (int a = 1; a <= A.chain.Str_L[j]; ++a) {
index_b = 0;
two_over_A_sinhlambda_sq_plus_sinzetaover2sq
= 2.0/((sinh(A.lambda[j][alpha] + 0.5 * II * A.chain.anis * (A.chain.Str_L[j] + 1.0 - 2.0 * a)
+ 0.25 * II * PI * (1.0 - A.chain.par[j])))
* (sinh(A.lambda[j][alpha] + 0.5 * II * A.chain.anis * (A.chain.Str_L[j] + 1.0 - 2.0 * a)
+ 0.25 * II * PI * (1.0 - A.chain.par[j])))
+ pow(sin(0.5*A.chain.anis), 2.0));
for (int k = 0; k < B.chain.Nstrings; ++k) {
for (int beta = 0; beta < B.base.Nrap[k]; ++beta) {
for (int b = 1; b <= B.chain.Str_L[k]; ++b) {
if (B.chain.Str_L[k] == 1) {
// use simplified code for one-string here: original form of Hm2P matrix
Fn_K_0_G_0 = Fn_K (A, j, alpha, a, B, k, beta, 0) *
exp(re_ln_Fn_G_0[k][beta] + II * im_ln_Fn_G_0[k][beta] - re_ln_Fn_F_B_0[k][beta] + logabssinzeta);
Fn_K_1_G_2 = Fn_K (A, j, alpha, a, B, k, beta, 1) *
exp(re_ln_Fn_G_2[k][beta] + II * im_ln_Fn_G_2[k][beta] - re_ln_Fn_F_B_0[k][beta] + logabssinzeta);
Prod_powerN = pow( B.chain.par[k] == 1 ?
(B.sinhlambda[k][beta] * B.chain.co_n_anis_over_2[1]
+ II * B.coshlambda[k][beta] * B.chain.si_n_anis_over_2[1])
/(B.sinhlambda[k][beta] * B.chain.co_n_anis_over_2[1]
- II * B.coshlambda[k][beta] * B.chain.si_n_anis_over_2[1])
:
(B.coshlambda[k][beta] * B.chain.co_n_anis_over_2[1]
+ II * B.sinhlambda[k][beta] * B.chain.si_n_anis_over_2[1])
/(B.coshlambda[k][beta] * B.chain.co_n_anis_over_2[1]
- II * B.sinhlambda[k][beta] * B.chain.si_n_anis_over_2[1])
, complex<DP> (B.chain.Nsites));
Hm2P[index_a][index_b] = Fn_K_0_G_0 - Prod_powerN * Fn_K_1_G_2
- two_over_A_sinhlambda_sq_plus_sinzetaover2sq
* exp(II*im_ln_Fn_F_B_0[k][beta] + logabssinzeta);
}
else {
if (b <= B.chain.Str_L[k] - 1) Hm2P[index_a][index_b] = Fn_K(A, j, alpha, a, B, k, beta, b);
else if (b == B.chain.Str_L[k]) {
Vect_CX ln_FunctionF(B.chain.Str_L[k] + 2);
for (int i = 0; i < B.chain.Str_L[k] + 2; ++i) ln_FunctionF[i] = ln_Fn_F (B, k, beta, i);
Vect_CX ln_FunctionG(B.chain.Str_L[k] + 2);
for (int i = 0; i < B.chain.Str_L[k] + 2; ++i) ln_FunctionG[i] = ln_Fn_G (A, B, k, beta, i);
sum1 = 0.0;
sum1 += Fn_K (A, j, alpha, a, B, k, beta, 0)
* exp(ln_FunctionG[0] + ln_FunctionG[1] - ln_FunctionF[0] - ln_FunctionF[1]);
sum1 += Fn_K (A, j, alpha, a, B, k, beta, B.chain.Str_L[k])
* exp(ln_FunctionG[B.chain.Str_L[k]] + ln_FunctionG[B.chain.Str_L[k] + 1]
- ln_FunctionF[B.chain.Str_L[k]] - ln_FunctionF[B.chain.Str_L[k] + 1]);
for (int jsum = 1; jsum < B.chain.Str_L[k]; ++jsum)
sum1 -= Fn_L (A, j, alpha, a, B, k, beta, jsum) *
exp(ln_FunctionG[jsum] + ln_FunctionG[jsum + 1] - ln_FunctionF[jsum] - ln_FunctionF[jsum + 1]);
sum2 = 0.0;
for (int jsum = 1; jsum <= B.chain.Str_L[k]; ++jsum)
sum2 += exp(ln_FunctionG[jsum] - ln_FunctionF[jsum]);
prod_num = exp(II * im_ln_Fn_F_B_0[k][beta] + ln_FunctionF[1]
- ln_FunctionG[B.chain.Str_L[k]] + logabssinzeta);
for (int jsum = 2; jsum <= B.chain.Str_L[k]; ++jsum)
prod_num *= exp(ln_FunctionG[jsum] - real(ln_FunctionF[jsum - 1]) + logabssinzeta);
// include all string contributions F_B_0 in this term
Hm2P[index_a][index_b] = prod_num * (sum1 - sum2 * two_over_A_sinhlambda_sq_plus_sinzetaover2sq);
} // else if (b == B.chain.Str_L[k])
} // else
index_b++;
}}} // sums over k, beta, b
index_a++;
}}} // sums over j, alpha, a
DP re_ln_det = real(lndet_LU_CX_dstry(Hm2P));
complex<DP> ln_ME_sq = log(0.25 * A.chain.Nsites) + real(ln_prod1 - ln_prod2) - real(ln_prod3) + real(ln_prod4)
+ 2.0 * re_ln_det - A.lnnorm - B.lnnorm;
return(0.5 * ln_ME_sq); // Return ME, not MEsq
}
// Version with phases:
complex<DP> ln_Sz_ME_with_phase (XXZ_Bethe_State& A, XXZ_Bethe_State& B)
{
// This function returns the natural log of the S^z operator matrix element.
// The A and B states can contain strings.
// Check that the two states refer to the same XXZ_Chain
if (A.chain != B.chain) ABACUSerror("Incompatible XXZ_Chains in Sz matrix element.");
// Check that A and B are compatible: same Mdown
if (A.base.Mdown != B.base.Mdown) {
cout << "Bra state: " << endl << A << endl;
cout << "Ket state: " << endl << B << endl;
cout << A.base.Mdown << "\t" << B.base.Mdown << endl;
ABACUSerror("Incompatible Mdown between the two states in Sz matrix element!");
}
if (A.iK == B.iK && (A.label != B.label))
return(-300.0); // matrix element identically vanishes
// Compute the sinh and cosh of rapidities
A.Compute_sinhlambda();
A.Compute_coshlambda();
B.Compute_sinhlambda();
B.Compute_coshlambda();
// Some convenient arrays
Lambda re_ln_Fn_F_B_0(B.chain, B.base);
Lambda im_ln_Fn_F_B_0(B.chain, B.base);
Lambda re_ln_Fn_G_0(B.chain, B.base);
Lambda im_ln_Fn_G_0(B.chain, B.base);
Lambda re_ln_Fn_G_2(B.chain, B.base);
Lambda im_ln_Fn_G_2(B.chain, B.base);
complex<DP> ln_prod1 = 0.0;
complex<DP> ln_prod2 = 0.0;
complex<DP> ln_prod3 = 0.0;
complex<DP> ln_prod4 = 0.0;
for (int i = 0; i < A.chain.Nstrings; ++i)
for (int alpha = 0; alpha < A.base.Nrap[i]; ++alpha)
for (int a = 1; a <= A.chain.Str_L[i]; ++a)
ln_prod1 += log(sinh(A.lambda[i][alpha] + 0.5 * II * A.chain.anis * (A.chain.Str_L[i] + 1.0 - 2.0 * a - 1.0)
+ 0.25 * II * PI * (1.0 - A.chain.par[i])));
complex<DP> shB;
for (int i = 0; i < B.chain.Nstrings; ++i)
for (int alpha = 0; alpha < B.base.Nrap[i]; ++alpha)
for (int a = 1; a <= B.chain.Str_L[i]; ++a)
if (norm(shB = sinh(B.lambda[i][alpha] + 0.5 * II * B.chain.anis * (B.chain.Str_L[i] + 1.0 - 2.0 * a - 1.0)
+ 0.25 * II * PI * (1.0 - B.chain.par[i]))) > 100.0 * MACHINE_EPS_SQ)
//ln_prod2 += log(norm(sinh(B.lambda[i][alpha] + 0.5 * II * B.chain.anis * (B.chain.Str_L[i] + 1.0 - 2.0 * a - 1.0)
// + 0.25 * II * PI * (1.0 - B.chain.par[i]))));
ln_prod2 += log(shB);
// Define the F ones earlier...
complex<DP> ln_Fn_F_B_0, ln_Fn_G_0, ln_Fn_G_2;
for (int j = 0; j < B.chain.Nstrings; ++j) {
for (int alpha = 0; alpha < B.base.Nrap[j]; ++alpha) {
//re_ln_Fn_F_B_0[j][alpha] = real(ln_Fn_F(B, j, alpha, 0));
//im_ln_Fn_F_B_0[j][alpha] = imag(ln_Fn_F(B, j, alpha, 0));
ln_Fn_F_B_0 = ln_Fn_F(B, j, alpha, 0);
re_ln_Fn_F_B_0[j][alpha] = real(ln_Fn_F_B_0);
im_ln_Fn_F_B_0[j][alpha] = imag(ln_Fn_F_B_0);
//re_ln_Fn_G_0[j][alpha] = real(ln_Fn_G(A, B, j, alpha, 0));
//im_ln_Fn_G_0[j][alpha] = imag(ln_Fn_G(A, B, j, alpha, 0));
ln_Fn_G_0 = ln_Fn_G(A, B, j, alpha, 0);
re_ln_Fn_G_0[j][alpha] = real(ln_Fn_G_0);
im_ln_Fn_G_0[j][alpha] = imag(ln_Fn_G_0);
//re_ln_Fn_G_2[j][alpha] = real(ln_Fn_G(A, B, j, alpha, 2));
//im_ln_Fn_G_2[j][alpha] = imag(ln_Fn_G(A, B, j, alpha, 2));
ln_Fn_G_0 = ln_Fn_G(A, B, j, alpha, 2);
re_ln_Fn_G_2[j][alpha] = real(ln_Fn_G_2);
im_ln_Fn_G_2[j][alpha] = imag(ln_Fn_G_2);
}
}
DP logsinzeta = log(sin(A.chain.anis));
// Define regularized products in prefactors
for (int j = 0; j < A.chain.Nstrings; ++j)
for (int alpha = 0; alpha < A.base.Nrap[j]; ++alpha)
for (int a = 1; a <= A.chain.Str_L[j]; ++a) {
ln_prod3 += ln_Fn_F (A, j, alpha, a - 1);
}
//ln_prod3 -= A.base.Mdown * log(abs(sin(A.chain.anis)));
ln_prod3 -= A.base.Mdown * logsinzeta;
for (int k = 0; k < B.chain.Nstrings; ++k)
for (int beta = 0; beta < B.base.Nrap[k]; ++beta)
for (int b = 1; b <= B.chain.Str_L[k]; ++b) {
if (b == 1) ln_prod4 += re_ln_Fn_F_B_0[k][beta];
else if (b > 1) ln_prod4 += ln_Fn_F(B, k, beta, b - 1);
}
//ln_prod4 -= B.base.Mdown * log(abs(sin(B.chain.anis)));
ln_prod4 -= B.base.Mdown * logsinzeta;
// Now proceed to build the Hm2P matrix
SQMat_CX Hm2P(0.0, A.base.Mdown);
int index_a = 0;
int index_b = 0;
complex<DP> sum1 = 0.0;
complex<DP> sum2 = 0.0;
complex<DP> prod_num = 0.0;
complex<DP> Fn_K_0_G_0 = 0.0;
complex<DP> Prod_powerN = 0.0;
complex<DP> Fn_K_1_G_2 = 0.0;
complex<DP> two_over_A_sinhlambda_sq_plus_sinzetaover2sq;
for (int j = 0; j < A.chain.Nstrings; ++j) {
for (int alpha = 0; alpha < A.base.Nrap[j]; ++alpha) {
for (int a = 1; a <= A.chain.Str_L[j]; ++a) {
index_b = 0;
two_over_A_sinhlambda_sq_plus_sinzetaover2sq = 2.0/((sinh(A.lambda[j][alpha] + 0.5 * II * A.chain.anis * (A.chain.Str_L[j] + 1.0 - 2.0 * a)
+ 0.25 * II * PI * (1.0 - A.chain.par[j])))
* (sinh(A.lambda[j][alpha] + 0.5 * II * A.chain.anis * (A.chain.Str_L[j] + 1.0 - 2.0 * a)
+ 0.25 * II * PI * (1.0 - A.chain.par[j])))
+ pow(sin(0.5*A.chain.anis), 2.0));
for (int k = 0; k < B.chain.Nstrings; ++k) {
for (int beta = 0; beta < B.base.Nrap[k]; ++beta) {
for (int b = 1; b <= B.chain.Str_L[k]; ++b) {
if (B.chain.Str_L[k] == 1) {
// use simplified code for one-string here: original form of Hm2P matrix
Fn_K_0_G_0 = Fn_K (A, j, alpha, a, B, k, beta, 0) *
exp(re_ln_Fn_G_0[k][beta] + II * im_ln_Fn_G_0[k][beta] - re_ln_Fn_F_B_0[k][beta] + logsinzeta);
Fn_K_1_G_2 = Fn_K (A, j, alpha, a, B, k, beta, 1) *
exp(re_ln_Fn_G_2[k][beta] + II * im_ln_Fn_G_2[k][beta] - re_ln_Fn_F_B_0[k][beta] + logsinzeta);
Prod_powerN = pow( B.chain.par[k] == 1 ?
(B.sinhlambda[k][beta] * B.chain.co_n_anis_over_2[1] + II * B.coshlambda[k][beta] * B.chain.si_n_anis_over_2[1])
/(B.sinhlambda[k][beta] * B.chain.co_n_anis_over_2[1] - II * B.coshlambda[k][beta] * B.chain.si_n_anis_over_2[1])
:
(B.coshlambda[k][beta] * B.chain.co_n_anis_over_2[1] + II * B.sinhlambda[k][beta] * B.chain.si_n_anis_over_2[1])
/(B.coshlambda[k][beta] * B.chain.co_n_anis_over_2[1] - II * B.sinhlambda[k][beta] * B.chain.si_n_anis_over_2[1])
, complex<DP> (B.chain.Nsites));
//cout << "Prod_powerN = " << Prod_powerN << "\t" << abs(Prod_powerN) << endl;
Hm2P[index_a][index_b] = Fn_K_0_G_0 - Prod_powerN * Fn_K_1_G_2 - two_over_A_sinhlambda_sq_plus_sinzetaover2sq
* exp(II*im_ln_Fn_F_B_0[k][beta] + logsinzeta);
}
else {
if (b <= B.chain.Str_L[k] - 1) Hm2P[index_a][index_b] = Fn_K(A, j, alpha, a, B, k, beta, b);
else if (b == B.chain.Str_L[k]) {
Vect_CX ln_FunctionF(B.chain.Str_L[k] + 2);
for (int i = 0; i < B.chain.Str_L[k] + 2; ++i) ln_FunctionF[i] = ln_Fn_F (B, k, beta, i);
Vect_CX ln_FunctionG(B.chain.Str_L[k] + 2);
for (int i = 0; i < B.chain.Str_L[k] + 2; ++i) ln_FunctionG[i] = ln_Fn_G (A, B, k, beta, i);
sum1 = 0.0;
sum1 += Fn_K (A, j, alpha, a, B, k, beta, 0) * exp(ln_FunctionG[0] + ln_FunctionG[1] - ln_FunctionF[0] - ln_FunctionF[1]);
sum1 += Fn_K (A, j, alpha, a, B, k, beta, B.chain.Str_L[k])
* exp(ln_FunctionG[B.chain.Str_L[k]] + ln_FunctionG[B.chain.Str_L[k] + 1]
- ln_FunctionF[B.chain.Str_L[k]] - ln_FunctionF[B.chain.Str_L[k] + 1]);
for (int jsum = 1; jsum < B.chain.Str_L[k]; ++jsum)
sum1 -= Fn_L (A, j, alpha, a, B, k, beta, jsum) *
exp(ln_FunctionG[jsum] + ln_FunctionG[jsum + 1] - ln_FunctionF[jsum] - ln_FunctionF[jsum + 1]);
sum2 = 0.0;
for (int jsum = 1; jsum <= B.chain.Str_L[k]; ++jsum) sum2 += exp(ln_FunctionG[jsum] - ln_FunctionF[jsum]);
prod_num = exp(II * im_ln_Fn_F_B_0[k][beta] + ln_FunctionF[1] - ln_FunctionG[B.chain.Str_L[k]] + logsinzeta);
for (int jsum = 2; jsum <= B.chain.Str_L[k]; ++jsum)
prod_num *= exp(ln_FunctionG[jsum] - real(ln_FunctionF[jsum - 1]) + logsinzeta);
// include all string contributions F_B_0 in this term
Hm2P[index_a][index_b] = prod_num * (sum1 - sum2 * two_over_A_sinhlambda_sq_plus_sinzetaover2sq);
} // else if (b == B.chain.Str_L[k])
} // else
index_b++;
}}} // sums over k, beta, b
index_a++;
}}} // sums over j, alpha, a
complex<DP> ln_det = lndet_LU_CX_dstry(Hm2P);
complex<DP> ln_ME = log(0.5 * sqrt(A.chain.Nsites)) + ln_prod1 - ln_prod2 - ln_prod3 + ln_prod4 + 2.0 * ln_det - A.lnnorm - B.lnnorm;
//cout << endl << ln_prod1 << "\t" << ln_prod2 << "\t" << ln_prod3 << "\t" << ln_prod4 << "\t" << A.lnnorm << "\t" << B.lnnorm
//<< "\t" << re_ln_det << "\t" << ln_form_factor_sq << endl;
//return(ln_ME_sq);
return(ln_ME); // Return ME, not MEsq
}
} // namespace ABACUS
Reference in New Issue View Git Blame Copy Permalink