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

1618 lines
51 KiB
C++
Executable File
Raw Permalink Blame History

/****************************************************************
This software is part of J.-S. Caux's C++ library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: XXX_VOA.cc
Defines all class procedures used for the XXX chain in zero field.
******************************************************************/
#include "ABACUS.h"
using namespace std;
namespace ABACUS {
DP I_integral (DP rho, DP req_prec)
{
DP t1 = 1.0; // delimiter between two integrals
DP tinf = 24.0; // such that exp(-2tinf) << machine_eps
DP Euler_Mascheroni = 0.577215664901532860606;
DP rho_used = fabs(rho);
if (rho_used > 10000.0) return(-PI * rho_used - 2.0 * Euler_Mascheroni);
Vect_DP args(2);
args[0] = 0.0;
args[1] = rho_used;
DP answer = -2.0 * Euler_Mascheroni - 2.0 * log(4.0 * rho_used * t1)
+ 2.0 * Integrate_rec (Integrand_11, args, 0, 0.0, t1, req_prec, 12)
- 2.0 * Integrate_rec (Integrand_12, args, 0, t1, tinf, req_prec, 12)
- Integrate_rec (Integrand_2, args, 0, 0.0, tinf, req_prec, 12);
return(answer);
}
/********************* TWO SPINONS ********************/
DP SF_2p (DP k, DP omega, I_table Itable)
{
// Careful ! This is S(k, omega) = S (k, w) |dw/domega| = 2 S(k, w)
DP w = 2.0 * omega;
// Rescale energies by factor 2 because of definitions of H_XXX (omega: S.S; w: 0.5 * sigma.sigma = 2 S.S)
DP wu = twoPI * sin(0.5 * k);
DP wl = PI * fabs(sin(k));
// Factor of 2: return S(k, omega), not S(k, w)
// 0.25 factor: 1/4 * 2 * 1/2, where 1/4 comes from Bougourzi, 2 is the Jacobian |dw/domega|
// and 1/2 is S^{zz} = 1/2 * S^{+-}
return(w < wu && w > wl ? 2.0 * 0.5
* exp(-Itable.Return_val (acosh(sqrt((wu * wu - wl * wl)/(w * w - wl * wl)))/PI))/sqrt(wu * wu - w * w) : 0.0);
}
DP SF_2p (Vect_DP args, I_table Itable)
{
// Careful ! This is S(k, w) !
// This uses args[0] = k, args[1] = w.
DP wu = twoPI * sin(0.5 * args[0]);
DP wl = PI * fabs(sin(args[0]));
// 0.5 factor: 1 from Bougourzi, and 1/2 is S^{zz} = 1/2 * S^{+-}
return(args[1] < wu && args[1] > wl ?
0.5 * exp(-Itable.Return_val (acosh(sqrt((wu * wu - wl * wl)/(args[1] * args[1] - wl * wl)))
/PI))/sqrt(wu * wu - args[1] * args[1]) : 0.0);
}
DP SF_2p_alt (Vect_DP args, I_table Itable)
{
// Careful ! This is S(k, w) !
// This uses args[0] = k, args[1] = alpha.
DP wu = twoPI * sin(0.5 * args[0]);
DP wl = PI * fabs(sin(args[0]));
DP w = wl * cosh(args[1]);
if (w >= wu || w <= wl) return(0.0);
DP factor = sqrt((w * w - wl * wl)/(wu * wu - w * w));
// 0.5 factor: 1 from Bougourzi, and 1/2 is S^{zz} = 1/2 * S^{+-}
return(factor * 0.5 * exp(-Itable.Return_val (acosh(sqrt((wu * wu - wl * wl)/(w * w - wl * wl)))/PI)));
}
DP SF_2p_w (Vect_DP args, I_table Itable)
{
return(args[1] * SF_2p (args, Itable));
}
DP SF_2p_w_alt (Vect_DP args, I_table Itable)
{
DP wu = twoPI * sin(0.5 * args[0]);
DP wl = PI * fabs(sin(args[0]));
DP w = wl * cosh(args[1]);
DP factor = sqrt((w * w - wl * wl)/(wu * wu - w * w));
// 0.5 factor: 1 from Bougourzi, and 1/2 is S^{zz} = 1/2 * S^{+-}
return(w * factor * 0.5 * exp(-Itable.Return_val (acosh(sqrt((wu * wu - wl * wl)/(w * w - wl * wl)))/PI)));
}
DP SF_2p_intw (Vect_DP args, I_table Itable)
{
// This returns \int_0^2PI dw/2PI S(k, w)
DP k = args[0];
DP req_prec = args[1];
int max_rec = int(args[2]);
DP wu = twoPI * sin(0.5 * k);
DP wl = PI * fabs(sin(k));
Vect_DP args_to_SF_2p(2);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be w
args_to_SF_2p[2] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return(Integrate_rec_using_table (SF_2p, args_to_SF_2p, 1, Itable, wl, wu, req_prec, max_rec)/twoPI);
}
DP SF_2p_intw_alt (Vect_DP args, I_table Itable)
{
// This returns \int_0^2PI dw/2PI S(k, w)
DP k = args[0];
DP req_prec = args[1];
int max_rec = int(args[2]);
DP wu = twoPI * sin(0.5 * k);
DP wl = PI * fabs(sin(k));
Vect_DP args_to_SF_2p(2);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be w
args_to_SF_2p[2] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return(Integrate_rec_using_table (SF_2p_alt, args_to_SF_2p, 1, Itable, 0.0, acosh(wu/wl), req_prec, max_rec)/twoPI);
}
DP SF_2p_check_sumrule (DP req_prec, int max_rec, I_table Itable)
{
// It's better to use the ..._alt function below.
Vect_DP args_to_SF_2p_intw (3);
args_to_SF_2p_intw[0] = 0.0; // this will be k
args_to_SF_2p_intw[1] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_SF_2p_intw[2] = DP(max_rec);
// Factor 2: int[0, 2PI] = 2 int[0, PI]
return(4.0 * 2.0 * Integrate_rec_using_table (SF_2p_intw, args_to_SF_2p_intw, 0, Itable,
0.0, PI, req_prec, max_rec)/twoPI);
// 4 : because full integral gives 1/4, return value here is sr fraction obtained.
}
DP SF_2p_check_sumrule_alt (DP req_prec, int max_rec, I_table Itable)
{
// This is the preferred version.
Vect_DP args_to_SF_2p_intw (3);
args_to_SF_2p_intw[0] = 0.0; // this will be k
args_to_SF_2p_intw[1] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_SF_2p_intw[2] = DP(max_rec);
// Factor 2: int[0, 2PI] = 2 int[0, PI]
return(4.0 * 2.0 * Integrate_rec_using_table (SF_2p_intw_alt, args_to_SF_2p_intw, 0, Itable,
0.0, PI, req_prec, max_rec)/twoPI);
// 4 : because full integral gives 1/4, return value here is sr fraction obtained.
}
DP Fixed_k_sumrule_w (DP k)
{
// This is K_1 (k) = \int dw/2PI w S(k, w) = 2 K_1^{KarbachPRB55} (k) = 4 E_G (1-cosk)/3N with E_G = -N(ln2 - 1/4).
return(4.0 * (log(2.0) - 0.25) * (1.0 - cos(k))/3.0);
}
DP Fixed_k_sumrule_omega (DP k)
{
return(0.5 * Fixed_k_sumrule_w(k));
}
DP SF_2p_check_fixed_k_sumrule (DP k, DP req_prec, int max_rec, I_table Itable)
{
// It's better to use the ..._alt function below.
DP wu = twoPI * sin(0.5 * k);
DP wl = PI * fabs(sin(k));
Vect_DP args_to_SF_2p(2);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be w
args_to_SF_2p[2] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return((Integrate_rec_using_table (SF_2p_w, args_to_SF_2p, 1,
Itable, wl, wu, req_prec, max_rec)/twoPI)/Fixed_k_sumrule_w(k));
}
DP SF_2p_check_fixed_k_sumrule_alt (DP k, DP req_prec, int max_rec, I_table Itable)
{
// This is the preferred version.
DP wu = twoPI * sin(0.5 * k);
DP wl = PI * fabs(sin(k));
Vect_DP args_to_SF_2p(2);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be alpha
args_to_SF_2p[2] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return((Integrate_rec_using_table (SF_2p_w_alt, args_to_SF_2p, 1, Itable, 0.0, acosh(wu/wl),
req_prec, max_rec)/twoPI)/Fixed_k_sumrule_w(k));
}
DP SF_2p_check_fixed_k_sumrule_opt (DP k, DP req_prec, int Npts, I_table Itable)
{
// This is the preferred version.
DP wu = twoPI * sin(0.5 * k);
DP wl = PI * fabs(sin(k));
Vect_DP args_to_SF_2p(2);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be alpha
args_to_SF_2p[2] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return(((Integrate_optimal_using_table (SF_2p_w_alt, args_to_SF_2p, 1, Itable, 0.0, acosh(wu/wl),
req_prec, 1.0e-32, Npts)).integ_est/twoPI)/Fixed_k_sumrule_w(k));
}
/********************** FOUR SPINONS **********************/
DP Sum_norm_gl (Vect_DP rho, DP req_prec)
{
complex<DP> g[4];
for (int l = 0; l < 4; ++l) g[l] = 0.0;
complex<DP> Plm[4];
complex<DP> Pm[4];
DP den = 0.0;
for (int j = 0; j < 4; ++j) {
Pm[j] = cosh(twoPI * rho[j]);
den = 1.0;
for (int i = 0; i < 4; ++i) if (i != j) den *= sinh(PI * (rho[j] - rho[i]));
Pm[j] /= den;
}
complex<DP> irhoj0[4];
complex<DP> irhoj1[4];
complex<DP> irhoj2[4];
complex<DP> irhoj3[4];
for (int j = 0; j < 4; ++j) {
irhoj0[j] = II * (rho[j] - rho[0]);
irhoj1[j] = II * (rho[j] - rho[1]);
irhoj2[j] = II * (rho[j] - rho[2]);
irhoj3[j] = II * (rho[j] - rho[3]);
}
// Do m = 0 terms:
for (int j = 0; j < 4; ++j) {
// Calling only Gamma (z) for Re(z) >= 0.5, in view of Lanczos method:
Pm[j] *= exp(ln_Gamma(0.5 + irhoj0[j]) - ln_Gamma(1.0 + irhoj0[j])
+ ln_Gamma(0.5 + irhoj1[j]) - ln_Gamma(1.0 + irhoj1[j])
+ ln_Gamma(0.5 + irhoj2[j]) - ln_Gamma(1.0 + irhoj2[j])
+ ln_Gamma(0.5 + irhoj3[j]) - ln_Gamma(1.0 + irhoj3[j]))
/((-0.5 + irhoj0[j]) * (-0.5 + irhoj1[j]) * (-0.5 + irhoj2[j]) * (-0.5 + irhoj3[j]));
for (int l = 0; l < 4; ++l) {
Plm[j] = 1.0;
for (int i = 0; i < 4; ++i) if (i != l) Plm[j] *= ((l > i ? -0.5 : 0.0) + II * (rho[j] - rho[i]));
if (j <= l) g[l] += Plm[j] * Pm[j]; // otherwise no m = 0 term
}
}
DP sum_norm_gl = norm(g[0]) + norm(g[1]) + norm(g[2]) + norm(g[3]);
DP old_sum_norm_gl = sum_norm_gl;
// Do m = 1, 2, ... terms:
int m = 1;
int m_to_reach = 1;
do {
old_sum_norm_gl = sum_norm_gl;
// We increase m by ten steps before checking sum_norm
m_to_reach = m + 10;
do {
for (int j = 0; j < 4; ++j) {
Pm[j] *= (m - 1.5 + irhoj0[j]) * (m - 1.5 + irhoj1[j]) * (m - 1.5 + irhoj2[j]) * (m - 1.5 + irhoj3[j])
/ ((DP(m) + irhoj0[j]) * (DP(m) + irhoj1[j]) * (DP(m) + irhoj2[j]) * (DP(m) + irhoj3[j]));
// FASTER: unwrap l, i loops
// l = 0:
g[0] += (DP(m) + irhoj1[j]) * (DP(m) + irhoj2[j]) * (DP(m) + irhoj3[j]) * Pm[j];
// l = 1;
g[1] += (m - 0.5 + irhoj0[j]) * (DP(m) + irhoj2[j]) * (DP(m) + irhoj3[j]) * Pm[j];
// l = 2;
g[2] += (m - 0.5 + irhoj0[j]) * (m - 0.5 + irhoj1[j]) * (DP(m) + irhoj3[j]) * Pm[j];
// l = 3;
g[3] += (m - 0.5 + irhoj0[j]) * (m - 0.5 + irhoj1[j]) * (m - 0.5 + irhoj2[j]) * Pm[j];
}
m++;
} while (m < m_to_reach);
sum_norm_gl = norm(g[0]) + norm(g[1]) + norm(g[2]) + norm(g[3]);
} while (m < 10 || sum_norm_gl/old_sum_norm_gl - 1.0 > req_prec && m < 100000);
return(norm(g[0]) + norm(g[1]) + norm(g[2]) + norm(g[3]));
}
DP Compute_C4 (DP req_prec)
{
Vect_DP args(2);
DP answer = exp(-8.0 * real(ln_Gamma (0.25)) - 9.0 * log(2.0)
+ 8.0 * Integrate_rec (Integrand_A, args, 0, 0.0, 50.0, req_prec, 16))/3.0;
return(answer);
}
DP SF_contrib (Vect_DP p, DP req_prec, I_table Itable)
{
Vect_DP rho(4);
DP W, Wu, sum_I_integrals, sum_norm_g;
for (int i = 0; i < 4; ++i) rho[i] = asinh(1.0/tan(p[i]))/twoPI;
Wu = twoPI* fabs(sin(0.5*(p[0] + p[1])));
W = -PI* (sin(p[0]) + sin(p[1]));
sum_I_integrals = 0.0;
for (int i1 = 0; i1 < 3; ++i1) for (int i2 = i1+1; i2 < 4; ++i2) {
sum_I_integrals += Itable.Return_val (rho[i1] - rho[i2]);
}
sum_norm_g = Sum_norm_gl (rho, req_prec);
return(exp(-sum_I_integrals) * sum_norm_g/sqrt(Wu * Wu - W * W));
}
DP J_fn (Vect_DP p, DP req_prec, I_table Itable)
{
Vect_DP rho(4);
DP sum_I_integrals, sum_norm_g;
for (int i = 0; i < 4; ++i) rho[i] = asinh(1.0/tan(p[i]))/twoPI;
sum_I_integrals = 0.0;
for (int i1 = 0; i1 < 3; ++i1) for (int i2 = i1+1; i2 < 4; ++i2) {
if (fabs(rho[i1] - rho[i2]) < 1.0e-10 || fabs(rho[i1] - rho[i2]) >= 1000.0) return(0.0); // safety here
sum_I_integrals += Itable.Return_val (rho[i1] - rho[i2]);
}
sum_norm_g = Sum_norm_gl (rho, req_prec);
return(exp(-sum_I_integrals) * sum_norm_g);
}
bool Set_p_given_kwKW (DP k, DP w, DP K, DP W, Vect_DP& p)
{
// Returns false if any of the p_i are out of -PI, 0
DP argacos1, argacos2;
if (fabs(argacos1 = W/(twoPI * sin(0.5*K))) > 1.0
|| fabs(argacos2 = (w - W)/(twoPI * sin (0.5 * fabs(k - K)))) > 1.0) return(false);
DP acos1 = acos(argacos1);
DP acos2 = acos(argacos2);
p[0] = -0.5 * K + acos1;
p[1] = -0.5 * K - acos1;
if (K <= k) {
p[2] = 0.5 * (K-k) + acos2;
p[3] = 0.5 * (K-k) - acos2;
}
else {
p[2] = 0.5 * (K-k) - PI + acos2;
p[3] = 0.5 * (K-k) - PI - acos2;
}
for (int i = 0; i < 4; ++i) if (p[i] < -PI || p[i] > 0.0) return(false);
return(true);
}
DP G_fn (Vect_DP args_to_G, I_table Itable)
{
Vect_DP p(4);
if (!Set_p_given_kwKW (args_to_G[0], args_to_G[1], args_to_G[2], args_to_G[3], p)) return(0.0);
DP answer = Jacobian_p3p4_KW (args_to_G[0], args_to_G[1], args_to_G[2], args_to_G[3])
* SF_contrib (p, args_to_G[4], Itable);
return(answer);
}
DP G1_fn (Vect_DP args_to_G, I_table Itable)
{
Vect_DP p(4);
DP W = twoPI * sin(0.5*args_to_G[2]) * cos(args_to_G[3]);
if (!Set_p_given_kwKW (args_to_G[0], args_to_G[1], args_to_G[2], W, p)) return(0.0);
return(J_fn (p, args_to_G[4], Itable)/sqrt(pow(twoPI * sin(0.5*(args_to_G[0] - args_to_G[2])), 2.0)
- pow(args_to_G[1] - W, 2.0)));
}
DP G2_fn (Vect_DP args_to_G, I_table Itable)
{
Vect_DP p(4);
DP W = args_to_G[1] - twoPI * fabs(sin(0.5*(args_to_G[0] - args_to_G[2]))) * cos(args_to_G[3]);
if (!Set_p_given_kwKW (args_to_G[0], args_to_G[1], args_to_G[2], W, p)) return(0.0);
return(J_fn (p, args_to_G[4], Itable)/sqrt(pow(twoPI * sin(0.5*args_to_G[2]), 2.0) - W * W));
}
DP G1_fn_mid (Vect_DP args_to_G, I_table Itable)
{
// Called by H_fn_mid.
// For the lower half of W interval
// Translation of arguments to G:
// args_to_G[0] = k;
// args_to_G[1] = w;
// args_to_G[2] = K;
// args_to_G[3] = alpha;
// args_to_G[4] = req_prec;
// args_to_G[5] = max_rec;
// args_to_G[6] = Wmid;
// args_to_G[7] = Wmax_used;
// args_to_G[8] = Wu_sq;
// args_to_G[9] = Wut_sq;
Vect_DP p(4);
DP W = args_to_G[6] * (args_to_G[3] > 1.0e-4 ? 1.0 - cos(args_to_G[3]) : 2.0 * pow(sin(0.5 * args_to_G[3]), 2));
// W = Wmid (1 - cos(alpha)), ensure some precision if alpha small
if (!Set_p_given_kwKW (args_to_G[0], args_to_G[1], args_to_G[2], W, p)) return(0.0);
DP answer = J_fn (p, args_to_G[4], Itable)
* sqrt(W * (2.0 * args_to_G[6] - W)/((args_to_G[8] - W * W) * (args_to_G[9] - pow(args_to_G[1] - W, 2.0))));
if (is_nan(answer)) {
cerr << setprecision(10) << "args_to_G1_fn_mid = " << args_to_G << "G1 = " << answer << "\tPut to zero..." << endl;
answer = 0.0;
}
return(answer);
}
DP G2_fn_mid (Vect_DP args_to_G, I_table Itable)
{
// Called by H_fn_mid.
// For the upper half of W interval
// See above for translation of arguments to G.
Vect_DP p(4);
DP W = args_to_G[7] * cos(args_to_G[3]); // W = Wmax cos(alpha)
if (!Set_p_given_kwKW (args_to_G[0], args_to_G[1], args_to_G[2], W, p)) return(0.0);
DP answer = J_fn (p, args_to_G[4], Itable)
* args_to_G[7] * sin(args_to_G[3]) /sqrt((args_to_G[8] - W * W) * (args_to_G[9] - pow(args_to_G[1] - W, 2.0)));
if (is_nan(answer)) {
cerr << setprecision(10) << "args_to_G2_fn_mid = " << args_to_G << "G2 = " << answer << endl;
cerr << W << "\t" << (args_to_G[7] * args_to_G[7] - W * W) << "\t" << (args_to_G[8] - W * W)
<< "\t" << (args_to_G[9] - pow(args_to_G[1] - W, 2.0)) << endl;
answer = 0.0;
}
return(answer);
}
DP G_fn_alt (Vect_DP args_to_G, I_table Itable)
{
Vect_DP p(4);
DP Wmin = args_to_G[4];
DP W = Wmin * cosh(args_to_G[3]);
if (!Set_p_given_kwKW (args_to_G[0], args_to_G[1], args_to_G[2], W, p)) return(0.0);
DP Wu1 = args_to_G[6];
DP Wu2 = args_to_G[7];
return(J_fn (p, args_to_G[8], Itable)
* sqrt((W * W - Wmin * Wmin)/((Wu1 * Wu1 - W * W) * (Wu2 * Wu2 - (args_to_G[1] - W) * (args_to_G[1] - W)))));
}
DP H_fn (Vect_DP args_to_H, I_table Itable)
{
// Translate arguments to more readable form
DP k = args_to_H[0];
DP w = args_to_H[1];
DP K = args_to_H[2];
DP req_prec = ABACUS::max(1.0e-14, args_to_H[3]);
int max_rec = int(args_to_H[4]);
Vect_DP args_to_G(6);
args_to_G[0] = k;
args_to_G[1] = w;
args_to_G[2] = K;
args_to_G[3] = 0.0;
args_to_G[4] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_G[5] = DP(max_rec);
DP Wmin_used = Wmin (k, w, K);
DP Wmax_used = Wmax (k, w, K);
return(Wmax_used > Wmin_used ?
Integrate_rec_using_table (G_fn, args_to_G, 3, Itable, Wmin_used, Wmax_used, req_prec, max_rec) : 0.0);
}
DP H2_fn (Vect_DP args_to_H, I_table Itable)
{
// Translate arguments to more readable form
DP k = args_to_H[0];
DP w = args_to_H[1];
DP K = args_to_H[2];
DP req_prec = ABACUS::max(1.0e-14, args_to_H[3]);
int max_rec = int(args_to_H[4]);
DP Wmin_used = Wmin (k, w, K);
DP Wmax_used = Wmax (k, w, K);
if (Wmax_used <= Wmin_used) return(0.0);
Vect_DP args_to_G(6);
args_to_G[0] = k;
args_to_G[1] = w;
args_to_G[2] = K;
args_to_G[3] = 0.0; // this will be the alpha parameter
args_to_G[4] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_G[5] = DP(max_rec);
DP Wmid = 0.5 * (Wmin_used + Wmax_used);
DP answer = 0.0;
DP alpha_L1 = acos(ABACUS::min(1.0, Wmax_used/(twoPI * sin(0.5*K)))); // to prevent nan
DP alpha_U1 = acos(Wmid/(twoPI * sin(0.5*K)));
DP alpha_L2 = acos(ABACUS::min(1.0, (w - Wmin_used)/(twoPI * fabs(sin(0.5*(k - K))))));
DP alpha_U2 = acos((w - Wmid)/(twoPI * fabs(sin(0.5*(k - K)))));
answer += Integrate_rec_using_table (G1_fn, args_to_G, 3, Itable, alpha_L1, alpha_U1, req_prec, max_rec);
answer += Integrate_rec_using_table (G2_fn, args_to_G, 3, Itable, alpha_L2, alpha_U2, req_prec, max_rec);
return(answer);
}
DP H_fn_mid (Vect_DP args_to_H, I_table Itable)
{
// For W in [Wmin, Wmid] we use the parametrization W = Wmid sin(alpha), alpha in [0, PI/2]
// such that dW = Wmid cos(alpha) dalpha is approx. 0 around alpha = PI/2 (W = 0).
// For W in [Wmid, Wmax] we use W = Wmax cos(alpha), alpha in [0, acos(Wmid/Wmax)]
// such that dW = -Wmax sin(alpha) dalpha is approx 0 around alpha = 0 (W = Wmax).
// Translation of args_to_H:
// args_to_H[0] = k;
// args_to_H[1] = w;
// args_to_H[2] = K; <------ integrated over, so newly set here
// args_to_H[3] = req_prec;
// args_to_H[4] = max_rec;
DP Wmin_used = Wmin (args_to_H[0], args_to_H[1], args_to_H[2]);
DP Wmax_used = Wmax (args_to_H[0], args_to_H[1], args_to_H[2]);
DP Wu_sq = pow(twoPI * sin(0.5 * args_to_H[2]), 2.0);
DP Wut_sq = pow(twoPI * sin(0.5 * (args_to_H[0] - args_to_H[2])), 2.0);
if (Wmax_used <= Wmin_used) return(0.0);
DP Wmid = 0.5 * (Wmin_used + Wmax_used);
Vect_DP args_to_G(10);
args_to_G[0] = args_to_H[0];
args_to_G[1] = args_to_H[1];
args_to_G[2] = args_to_H[2];
args_to_G[3] = 0.0; // this will be the alpha parameter
args_to_G[4] = ABACUS::max(1.0e-14, 0.01 * args_to_H[3]);
args_to_G[5] = args_to_H[4];
args_to_G[6] = Wmid;
args_to_G[7] = Wmax_used;
args_to_G[8] = Wu_sq;
args_to_G[9] = Wut_sq;
DP answer = 0.0;
DP alpha_L1 = 0.0;
DP alpha_U1 = 0.5 * PI;
DP alpha_L2 = 0.0;
DP alpha_U2 = acos(Wmid/Wmax_used);
answer += Integrate_rec_using_table (G1_fn_mid, args_to_G, 3, Itable, alpha_L1,
alpha_U1, args_to_H[3], int(args_to_H[4]));
answer += Integrate_rec_using_table (G2_fn_mid, args_to_G, 3, Itable, alpha_L2,
alpha_U2, args_to_H[3], int(args_to_H[4]));
return(answer);
}
DP H_fn_mid_opt (Vect_DP args_to_H, I_table Itable)
{
// For W in [Wmin, Wmid] we use the parametrization W = Wmid sin(alpha), alpha in [0, PI/2]
// such that dW = Wmid cos(alpha) dalpha is approx. 0 around alpha = PI/2 (W = 0).
// For W in [Wmid, Wmax] we use W = Wmax cos(alpha), alpha in [0, acos(Wmid/Wmax)]
// such that dW = -Wmax sin(alpha) dalpha is approx 0 around alpha = 0 (W = Wmax).
// Translation of args_to_H:
// args_to_H[0] = k;
// args_to_H[1] = w;
// args_to_H[2] = K; <------ integrated over, so newly set here
// args_to_H[3] = req_prec;
// args_to_H[4] = Npts;
DP Wmin_used = Wmin (args_to_H[0], args_to_H[1], args_to_H[2]);
DP Wmax_used = Wmax (args_to_H[0], args_to_H[1], args_to_H[2]);
DP Wu_sq = pow(twoPI * sin(0.5 * args_to_H[2]), 2.0);
DP Wut_sq = pow(twoPI * sin(0.5 * (args_to_H[0] - args_to_H[2])), 2.0);
if (Wmax_used <= Wmin_used) return(0.0);
DP Wmid = 0.5 * (Wmin_used + Wmax_used);
Vect_DP args_to_G(10);
args_to_G[0] = args_to_H[0];
args_to_G[1] = args_to_H[1];
args_to_G[2] = args_to_H[2];
args_to_G[3] = 0.0; // this will be the alpha parameter
args_to_G[4] = ABACUS::max(1.0e-14, 0.01 * args_to_H[3]);
args_to_G[5] = args_to_H[4];
args_to_G[6] = Wmid;
args_to_G[7] = Wmax_used;
args_to_G[8] = Wu_sq;
args_to_G[9] = Wut_sq;
DP answer = 0.0;
DP alpha_L1 = 0.0;
DP alpha_U1 = 0.5 * PI;
DP alpha_L2 = 0.0;
DP alpha_U2 = acos(Wmid/Wmax_used);
answer += (Integrate_optimal_using_table (G1_fn_mid, args_to_G, 3, Itable, alpha_L1, alpha_U1,
args_to_H[3], 1.0e-32, int(args_to_H[4]))).integ_est;
answer += (Integrate_optimal_using_table (G2_fn_mid, args_to_G, 3, Itable, alpha_L2, alpha_U2,
args_to_H[3], 1.0e-32, int(args_to_H[4]))).integ_est;
return(answer);
}
DP H_fn_alt (Vect_DP args_to_H, I_table Itable)
{
// Translate arguments to more readable form
DP k = args_to_H[0];
DP w = args_to_H[1];
DP K = args_to_H[2];
DP req_prec = ABACUS::max(1.0e-14, args_to_H[3]);
int max_rec = int(args_to_H[4]);
DP Wmin_used = Wmin (k, w, K);
DP Wmax_used = Wmax (k, w, K);
if (Wmax_used <= Wmin_used) return(0.0);
Vect_DP args_to_G(10);
args_to_G[0] = k;
args_to_G[1] = w;
args_to_G[2] = K;
args_to_G[3] = 0.0; // this will become alpha
args_to_G[4] = Wmin_used;
args_to_G[5] = Wmax_used;
args_to_G[6] = twoPI * fabs(sin(0.5*args_to_G[2])); // Wu1
args_to_G[7] = twoPI * fabs(sin(0.5*(args_to_G[0] - args_to_G[2]))); // Wu2
args_to_G[8] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_G[9] = DP(max_rec);
return(Integrate_rec_using_table (G_fn_alt, args_to_G, 3, Itable, 0.0, acosh(Wmax_used/Wmin_used), req_prec, max_rec));
}
DP SF_4p_kwKW (Vect_DP args, I_table Itable)
{
// Translate:
// args[0] = k;
// args[1] = omega;
// args[2] = req_prec;
// args[3] = max_rec;
DP k = args[0];
DP omega = args[1];
DP req_prec = args[2];
int max_rec = int(args[3]);
Vect_DP args_to_H(5);
args_to_H[0] = k; // shift of PI in Bougourzi: because they do FM case.
// We want AFM, so SF_4p (k, omega) is correctly obtained directly from the RHS of their formula.
DP w = 2.0 * omega;
// Rescale energies by factor 2 because of definitions of H_XXX (omega: S.S; w: 0.5 * sigma.sigma = 2 S.S)
args_to_H[1] = w;
args_to_H[2] = 0.0; // this is K
args_to_H[3] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_H[4] = DP(max_rec);
if (w > wmax_4p(k) || w < wmin_4p(k)) {
return(0.0);
}
DP prefactor = 2.0 * 0.5 * 4.0 * Compute_C4 (req_prec); // 4 comes from using p1 > p2 & p3 > p4 instead of whole interval.
// 2 from Jacobian |dw/domega|
// 0.5 from S^{zz} = S^{pm}/2
// Define the K integral domain
Domain<DP> Kdomain;
// First, the inclusions:
if (w <= twoPI * sin(0.5*k)) Kdomain.Include (0.0, twoPI);
else {
if (w < 4.0*PI * sin(0.25 * k)) {
DP K1bm = 0.5 * k - 2.0 * acos (w/(4.0*PI * sin(0.25 * k)));
DP K1bp = 0.5 * k + 2.0 * acos (w/(4.0*PI * sin(0.25 * k)));
Kdomain.Include (K1bm, K1bp);
}
if (w < 4.0*PI * cos(0.25 * k)) {
DP K1am = 0.5 * k + PI - 2.0 * acos (w/(4.0*PI * cos(0.25 * k)));
DP K1ap = 0.5 * k + PI + 2.0 * acos (w/(4.0*PI * cos(0.25 * k)));
Kdomain.Include (K1am, K1ap);
}
}
// Now the exclusions:
if (w < twoPI * sin(0.5*k)) {
DP K2dm = 0.5 * k - acos (w/(twoPI * sin (0.5 * k)));
DP K2dp = 0.5 * k + acos (w/(twoPI * sin (0.5 * k)));
Kdomain.Exclude (K2dm, K2dp);
DP K3cm = K2dm + PI;
DP K3cp = K2dp + PI;
Kdomain.Exclude (K3cm, K3cp);
}
if (w < twoPI * cos(0.5*k)) {
DP K2cm = 0.5 * k + asin(w/(twoPI * cos(0.5*k)));
DP K2cp = 0.5 * k + PI - asin(w/(twoPI * cos(0.5*k)));
Kdomain.Exclude (K2cm, K2cp);
DP K3em = K2cm + PI;
DP K3ep = K2cp + PI;
Kdomain.Exclude (K3em, K3ep);
}
// Use (K,W) -> (k-K, w-W) symmetry to restrict to K in [k/2, k/2+PI]
Kdomain.Exclude (0.0, 0.5 * k);
Kdomain.Exclude (0.5 * k + PI, twoPI);
prefactor *= 2.0;
DP answer = 0.0;
for (int idom = 0; idom < Kdomain.Ndomains(); ++idom)
answer += Integrate_rec_using_table (H_fn_mid, args_to_H, 2, Itable, Kdomain.xmin(idom), Kdomain.xmax(idom),
req_prec, max_rec);
return (prefactor * answer);
}
DP SF_4p_kwKW_opt (Vect_DP args, I_table Itable)
{
// Translate:
// args[0] = k;
// args[1] = omega;
// args[2] = req_prec;
// args[3] = Npts_K;
// args[4] = Npts_W;
DP k = args[0];
DP omega = args[1];
DP req_prec = args[2];
int Npts_K = int(args[3]);
int Npts_W = int(args[4]);
Vect_DP args_to_H(5);
args_to_H[0] = k; // shift of PI in Bougourzi: because they do FM case.
// We want AFM, so SF_4p (k, omega) is correctly obtained directly from the RHS of their formula.
DP w = 2.0 * omega;
// Rescale energies by factor 2 because of definitions of H_XXX (omega: S.S; w: 0.5 * sigma.sigma = 2 S.S)
args_to_H[1] = w;
args_to_H[2] = 0.0; // this is K
args_to_H[3] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_H[4] = DP(Npts_W);
if (w > wmax_4p(k) || w < wmin_4p(k)) {
return(0.0);
}
DP prefactor = 2.0 * 0.5 * 4.0 * Compute_C4 (req_prec);
// 4 comes from using p1 > p2 & p3 > p4 instead of whole interval.
// 2 from Jacobian |dw/domega|
// 0.5 from S^{zz} = S^{pm}/2
// Define the K integral domain
Domain<DP> Kdomain;
// First, the inclusions:
if (w <= twoPI * sin(0.5*k)) Kdomain.Include (0.0, twoPI);
else {
if (w < 4.0*PI * sin(0.25 * k)) {
DP K1bm = 0.5 * k - 2.0 * acos (w/(4.0*PI * sin(0.25 * k)));
DP K1bp = 0.5 * k + 2.0 * acos (w/(4.0*PI * sin(0.25 * k)));
Kdomain.Include (K1bm, K1bp);
}
if (w < 4.0*PI * cos(0.25 * k)) {
DP K1am = 0.5 * k + PI - 2.0 * acos (w/(4.0*PI * cos(0.25 * k)));
DP K1ap = 0.5 * k + PI + 2.0 * acos (w/(4.0*PI * cos(0.25 * k)));
Kdomain.Include (K1am, K1ap);
}
}
// Now the exclusions:
if (w < twoPI * sin(0.5*k)) {
DP K2dm = 0.5 * k - acos (w/(twoPI * sin (0.5 * k)));
DP K2dp = 0.5 * k + acos (w/(twoPI * sin (0.5 * k)));
Kdomain.Exclude (K2dm, K2dp);
DP K3cm = K2dm + PI;
DP K3cp = K2dp + PI;
Kdomain.Exclude (K3cm, K3cp);
}
if (w < twoPI * cos(0.5*k)) {
DP K2cm = 0.5 * k + asin(w/(twoPI * cos(0.5*k)));
DP K2cp = 0.5 * k + PI - asin(w/(twoPI * cos(0.5*k)));
Kdomain.Exclude (K2cm, K2cp);
DP K3em = K2cm + PI;
DP K3ep = K2cp + PI;
Kdomain.Exclude (K3em, K3ep);
}
// Use (K,W) -> (k-K, w-W) symmetry to restrict to K in [k, k+PI]
Kdomain.Exclude (0.0, 0.5 * k);
Kdomain.Exclude (0.5 * k + PI, twoPI);
prefactor *= 2.0;
DP answer = 0.0;
for (int idom = 0; idom < Kdomain.Ndomains(); ++idom)
answer += (Integrate_optimal_using_table (H_fn_mid_opt, args_to_H, 2, Itable, Kdomain.xmin(idom),
Kdomain.xmax(idom), req_prec, 1.0e-32, Npts_K)).integ_est;
return (prefactor * answer);
}
DP SF_4p_kwKW_alpha (Vect_DP args, I_table Itable)
{
// Better version for fixed k sum rule
// Translate:
// args[0] = k;
// args[1] = alpha; <-- integration variable, omega = omegamin + (omegamax - omegamin) (1-cos(alpha))
// args[2] = req_prec;
// args[3] = max_rec;
Vect_DP args_to_SF_4p_kwKW = args;
DP omegamin = 0.5 * wmin_4p (args[0]);
DP omegamax = 0.5 * wmax_4p (args[0]);
args_to_SF_4p_kwKW[1] = omegamin + (omegamax - omegamin) * (1.0 - cos(args[1]));
return((omegamax - omegamin) * sin(args[1]) * SF_4p_kwKW (args_to_SF_4p_kwKW, Itable));
}
DP SF_4p_kwKW_alpha_opt (Vect_DP args, I_table Itable)
{
// Better version for fixed k sum rule
// Translate:
// args[0] = k;
// args[1] = alpha; <-- integration variable, omega = omegamin + (omegamax - omegamin) (1-cos(alpha))
// args[2] = req_prec;
// args[3] = Npts_K;
// args[4] = Npts_W;
Vect_DP args_to_SF_4p_kwKW = args;
DP omegamin = 0.5 * wmin_4p (args[0]);
DP omegamax = 0.5 * wmax_4p (args[0]);
args_to_SF_4p_kwKW[1] = omegamin + (omegamax - omegamin) * (1.0 - cos(args[1]));
return((omegamax - omegamin) * sin(args[1]) * SF_4p_kwKW_opt (args_to_SF_4p_kwKW, Itable));
}
DP SF_4p_kwKW_cosh_alpha_opt (Vect_DP args, I_table Itable)
{
// Better version for fixed k sum rule
// Translate:
// args[0] = k;
// args[1] = alpha; <-- integration variable, omega = omegamin * cosh(alpha)
// args[2] = req_prec;
// args[3] = Npts_K;
// args[4] = Npts_W;
Vect_DP args_to_SF_4p_kwKW = args;
DP omegamin = 0.5 * wmin_4p (args[0]);
args_to_SF_4p_kwKW[1] = omegamin * cosh(args[1]);
return(omegamin * sinh(args[1]) * SF_4p_kwKW_opt (args_to_SF_4p_kwKW, Itable));
}
/******************************************************************************************/
// Interface to used version:
DP SF_4p_rec (DP k, DP omega, DP req_prec, int max_rec, I_table Itable)
{
// CAREFUL !! This is S(k, omega) = 2 S(k, w)
Vect_DP args_to_SF(4);
args_to_SF[0] = k;
args_to_SF[1] = omega;
args_to_SF[2] = req_prec;
args_to_SF[3] = DP(max_rec);
return(2.0 * SF_4p_kwKW (args_to_SF, Itable));
}
DP SF_4p (DP k, DP omega, I_table Itable)
{
// Fixes req_prec and max_rec to default values
return(SF_4p_rec (k, omega, default_req_prec, default_max_rec, Itable));
}
// Interface to used version:
DP SF_4p_opt (DP k, DP omega, DP req_prec, int Npts_K, int Npts_W, I_table Itable)
{
// CAREFUL !! This is S(k, omega) = 2 S(k, w)
Vect_DP args_to_SF(5);
args_to_SF[0] = k;
args_to_SF[1] = omega;
args_to_SF[2] = req_prec;
args_to_SF[3] = DP(Npts_K);
args_to_SF[4] = DP(Npts_W);
return(2.0 * SF_4p_kwKW_opt (args_to_SF, Itable));
}
/******************************************************************************************/
void Translate_raw_4p_data (DP k, int dim_w, const char* SFraw_Cstr, const char* SF_Cstr,
const char* SFsrc_Cstr, I_table Itable)
{
DP omegamin = 0.5 * wmin_4p (k); // Correct for factor of 2 in E between me & Bougourzi
DP omegamax = 0.5 * wmax_4p (k);
DP alpha_in;
DP SF_in;
DP alpha_in_old = -1.0;
DP SF_in_old = -1.0;
DP* alpha = new DP[dim_w];
DP* omega = new DP[dim_w];
DP* SF_4p_dat = new DP[dim_w];
DP* SF_2p_dat = new DP[dim_w];
int* index = new int[dim_w];
ifstream SFraw;
SFraw.open(SFraw_Cstr);
if (SFraw.fail()) ABACUSerror("Couldn't open SFraw file in Translate.");
int i = 0;
for (i = 0; i < dim_w; ++i) {
if (SFraw.peek() == EOF) ABACUSerror("Not enough data points in file...");
index[i] = i;
SFraw >> alpha_in >> SF_in;
alpha[i] = alpha_in;
omega[i] = omegamin + (omegamax - omegamin) * (1.0 - cos(alpha_in));
// CAREFUL !!! SF_in is S (k, w), and we want S (k, omega) = 2 S(k, w)
SF_4p_dat[i] = 2.0 * SF_in/((omegamax - omegamin) * sin(alpha_in));
SF_2p_dat[i] = SF_2p (k, omega[i], Itable); // This already is S(k, omega)
alpha_in_old = alpha_in;
SF_in_old = SF_in;
}
SFraw.close();
if (i != dim_w) {
ABACUSerror("Incorrect number of data points in file.");
}
QuickSort (omega, index, 0, dim_w - 1);
DP fixed_k_sr_2p = 0.0;
DP fixed_k_sr_4p = 0.0;
DP full_sr_2p = 0.0;
DP full_sr_4p = 0.0;
DP Jac_dalpha = 0.0; // This is domega = (omegamax - omegamin) sin alpha dalpha
ofstream SF;
SF.open(SF_Cstr);
for (int j = 0; j < dim_w; ++j)
SF << setprecision(16) << omega[j] << "\t" << SF_4p_dat[index[j] ] << "\t" << SF_2p_dat[index[j] ] << "\t"
<< SF_4p_dat[index[j] ] + SF_2p_dat[index[j] ] << endl;
SF.close();
// Compute first moment sum rule
Jac_dalpha = (omegamax - omegamin) * sin(alpha[index[1] ]) * 0.5 * (alpha[index[2] ] - alpha[index[0] ]);
fixed_k_sr_4p += Jac_dalpha * (omega[0] * SF_4p_dat[index[0] ] + omega[1] * SF_4p_dat[index[1] ]);
fixed_k_sr_2p += Jac_dalpha * (omega[0] * SF_2p_dat[index[0] ] + omega[1] * SF_2p_dat[index[1] ]);
full_sr_4p += Jac_dalpha * (SF_4p_dat[index[0] ] + SF_4p_dat[index[1] ]);
full_sr_2p += Jac_dalpha * (SF_2p_dat[index[0] ] + SF_2p_dat[index[1] ]);
for (int i = 2; i < dim_w - 2; ++i) {
Jac_dalpha = (omegamax - omegamin) * sin(alpha[index[i] ]) * 0.5 * (alpha[index[i + 1] ] - alpha[index[i - 1] ]);
fixed_k_sr_4p += Jac_dalpha * omega[i] * SF_4p_dat[index[i] ];
fixed_k_sr_2p += Jac_dalpha * omega[i] * SF_2p_dat[index[i] ];
full_sr_4p += Jac_dalpha * SF_4p_dat[index[i] ];
full_sr_2p += Jac_dalpha * SF_2p_dat[index[i] ];
}
Jac_dalpha = (omegamax - omegamin) * sin(alpha[index[dim_w - 2] ])
* 0.5 * (alpha[index[dim_w - 1] ] - alpha[index[dim_w - 3] ]);
fixed_k_sr_4p += Jac_dalpha * (omega[dim_w - 2] * SF_4p_dat[index[dim_w - 2] ]
+ omega[dim_w - 1] * SF_4p_dat[index[dim_w - 1] ]);
fixed_k_sr_2p += Jac_dalpha * (omega[dim_w - 2] * SF_2p_dat[index[dim_w - 2] ]
+ omega[dim_w - 1] * SF_2p_dat[index[dim_w - 1] ]);
full_sr_4p += Jac_dalpha * (SF_4p_dat[index[dim_w - 2] ] + SF_4p_dat[index[dim_w - 1] ]);
full_sr_2p += Jac_dalpha * (SF_2p_dat[index[dim_w - 2] ] + SF_2p_dat[index[dim_w - 1] ]);
ofstream SFsrc;
SFsrc.open(SFsrc_Cstr);
// Reintroduce 1/PI since \int domega/2PI
full_sr_4p /= twoPI;
full_sr_2p /= twoPI;
fixed_k_sr_4p /= twoPI;
fixed_k_sr_2p /= twoPI;
SFsrc << setprecision(16) << fixed_k_sr_4p << "\t"
<< fixed_k_sr_2p << "\t"
<< fixed_k_sr_4p + fixed_k_sr_2p << "\t"
<< fixed_k_sr_4p/Fixed_k_sumrule_omega(k) << "\t"
<< fixed_k_sr_2p/Fixed_k_sumrule_omega(k) << "\t"
<< (fixed_k_sr_4p + fixed_k_sr_2p)/Fixed_k_sumrule_omega(k) << endl;
SFsrc << setprecision(16) << full_sr_4p << "\t" << full_sr_2p << "\t" << full_sr_4p + full_sr_2p << "\t"
<< 0.25 * full_sr_4p << "\t" << 0.25 * full_sr_2p << "\t" << 0.25 * (full_sr_4p + full_sr_2p) << endl;
SFsrc.close();
delete[] omega;
delete[] SF_4p_dat;
delete[] SF_2p_dat;
delete[] index;
return;
}
void Translate_raw_4p_data_cosh (DP k, int dim_w, const char* SFraw_Cstr, const char* SF_Cstr,
const char* SFsrc_Cstr, I_table Itable)
{
// Here, omega = omegamin * cosh(alpha)
DP omegamin = 0.5 * wmin_4p (k); // Correct for factor of 2 in E between me & Bougourzi
DP alpha_in;
DP SF_in;
DP alpha_in_old = -1.0;
DP SF_in_old = -1.0;
DP* alpha = new DP[dim_w];
DP* omega = new DP[dim_w];
DP* SF_4p_dat = new DP[dim_w];
DP* SF_2p_dat = new DP[dim_w];
int* index = new int[dim_w];
ifstream SFraw;
SFraw.open(SFraw_Cstr);
if (SFraw.fail()) ABACUSerror("Couldn't open SFraw file in Translate.");
int i = 0;
for (i = 0; i < dim_w; ++i) {
if (SFraw.peek() == EOF) ABACUSerror("Not enough data points in file...");
index[i] = i;
SFraw >> alpha_in >> SF_in;
alpha[i] = alpha_in;
omega[i] = omegamin * cosh(alpha_in);
// CAREFUL !!! SF_in is S (k, w), and we want S (k, omega) = 2 S(k, w)
SF_4p_dat[i] = 2.0 * SF_in/(omegamin * sinh(alpha_in));
SF_2p_dat[i] = SF_2p (k, omega[i], Itable); // This already is S(k, omega)
alpha_in_old = alpha_in;
SF_in_old = SF_in;
}
SFraw.close();
if (i != dim_w) {
ABACUSerror("Incorrect number of data points in file.");
}
QuickSort (omega, index, 0, dim_w - 1);
DP fixed_k_sr_2p = 0.0;
DP fixed_k_sr_4p = 0.0;
DP full_sr_2p = 0.0;
DP full_sr_4p = 0.0;
DP Jac_dalpha = 0.0; // This is domega = omegamin sinh alpha dalpha
ofstream SF;
SF.open(SF_Cstr);
for (int j = 0; j < dim_w; ++j)
SF << setprecision(16) << omega[j] << "\t" << SF_4p_dat[index[j] ] << "\t" << SF_2p_dat[index[j] ] << "\t"
<< SF_4p_dat[index[j] ] + SF_2p_dat[index[j] ] << endl;
SF.close();
// Compute first moment sum rule
Jac_dalpha = omegamin * sinh(alpha[index[1] ]) * 0.5 * (alpha[index[2] ] - alpha[index[0] ]);
fixed_k_sr_4p += Jac_dalpha * (omega[0] * SF_4p_dat[index[0] ] + omega[1] * SF_4p_dat[index[1] ]);
fixed_k_sr_2p += Jac_dalpha * (omega[0] * SF_2p_dat[index[0] ] + omega[1] * SF_2p_dat[index[1] ]);
full_sr_4p += Jac_dalpha * (SF_4p_dat[index[0] ] + SF_4p_dat[index[1] ]);
full_sr_2p += Jac_dalpha * (SF_2p_dat[index[0] ] + SF_2p_dat[index[1] ]);
for (int i = 2; i < dim_w - 2; ++i) {
Jac_dalpha = omegamin * sinh(alpha[index[i] ]) * 0.5 * (alpha[index[i + 1] ] - alpha[index[i - 1] ]);
fixed_k_sr_4p += Jac_dalpha * omega[i] * SF_4p_dat[index[i] ];
fixed_k_sr_2p += Jac_dalpha * omega[i] * SF_2p_dat[index[i] ];
full_sr_4p += Jac_dalpha * SF_4p_dat[index[i] ];
full_sr_2p += Jac_dalpha * SF_2p_dat[index[i] ];
}
Jac_dalpha = omegamin * sinh(alpha[index[dim_w - 2] ]) * 0.5 * (alpha[index[dim_w - 1] ] - alpha[index[dim_w - 3] ]);
fixed_k_sr_4p += Jac_dalpha * (omega[dim_w - 2] * SF_4p_dat[index[dim_w - 2] ]
+ omega[dim_w - 1] * SF_4p_dat[index[dim_w - 1] ]);
fixed_k_sr_2p += Jac_dalpha * (omega[dim_w - 2] * SF_2p_dat[index[dim_w - 2] ]
+ omega[dim_w - 1] * SF_2p_dat[index[dim_w - 1] ]);
full_sr_4p += Jac_dalpha * (SF_4p_dat[index[dim_w - 2] ] + SF_4p_dat[index[dim_w - 1] ]);
full_sr_2p += Jac_dalpha * (SF_2p_dat[index[dim_w - 2] ] + SF_2p_dat[index[dim_w - 1] ]);
ofstream SFsrc;
SFsrc.open(SFsrc_Cstr);
// Reintroduce 1/PI since \int domega/2PI
full_sr_4p /= twoPI;
full_sr_2p /= twoPI;
fixed_k_sr_4p /= twoPI;
fixed_k_sr_2p /= twoPI;
SFsrc << setprecision(16) << fixed_k_sr_4p << "\t"
<< fixed_k_sr_2p << "\t"
<< fixed_k_sr_4p + fixed_k_sr_2p << "\t"
<< fixed_k_sr_4p/Fixed_k_sumrule_omega(k) << "\t"
<< fixed_k_sr_2p/Fixed_k_sumrule_omega(k) << "\t"
<< (fixed_k_sr_4p + fixed_k_sr_2p)/Fixed_k_sumrule_omega(k) << endl;
SFsrc << setprecision(16) << full_sr_4p << "\t" << full_sr_2p << "\t" << full_sr_4p + full_sr_2p << "\t"
<< 0.25 * full_sr_4p << "\t" << 0.25 * full_sr_2p << "\t" << 0.25 * (full_sr_4p + full_sr_2p) << endl;
SFsrc.close();
delete[] omega;
delete[] SF_4p_dat;
delete[] SF_2p_dat;
delete[] index;
return;
}
// Function producing a fixed k scan, with data file
DP SF_4p_rec (DP k, DP req_prec, int max_rec_w, int max_rec, I_table Itable)
{
stringstream SFraw_stringstream;
string SFraw_string;
SFraw_stringstream << "SF_4p_k_" << k << "_prec_" << req_prec << "_max_rec_w_" << max_rec_w
<< "_max_rec_" << max_rec << ".raw";
SFraw_string = SFraw_stringstream.str();
const char* SFraw_Cstr = SFraw_string.c_str();
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "SF_4p_k_" << k << "_prec_" << req_prec << "_max_rec_w_" << max_rec_w
<< "_max_rec_" << max_rec << ".dat";
SF_string = SF_stringstream.str();
const char* SF_Cstr = SF_string.c_str();
stringstream SFsrc_stringstream;
string SFsrc_string;
SFsrc_stringstream << "SF_4p_k_" << k << "_prec_" << req_prec << "_max_rec_w_" << max_rec_w
<< "_max_rec_" << max_rec << ".src";
SFsrc_string = SFsrc_stringstream.str();
const char* SFsrc_Cstr = SFsrc_string.c_str();
ofstream SFraw_outfile;
SFraw_outfile.open(SFraw_Cstr);
ofstream SFsrc_outfile;
SFsrc_outfile.open(SFsrc_Cstr, ofstream::app);
Vect_DP args_to_SF(4);
args_to_SF[0] = k;
args_to_SF[1] = 0.0; // integration variable
args_to_SF[2] = req_prec;
args_to_SF[3] = DP(max_rec);
// Version using omega = omegamin + (omegamax - omegamin) * (1-cos(alpha))
DP answer = Integrate_rec_using_table (SF_4p_kwKW_alpha, args_to_SF, 1, Itable, 0.0, 0.5*PI, req_prec, max_rec_w, SFraw_outfile)/twoPI;
SFraw_outfile.close();
SFsrc_outfile << answer << endl;
SFsrc_outfile.close();
// Translate raw data into SF_4p (k,omega) data
Translate_raw_4p_data (k, int(pow(3.0, max_rec_w + 2)), SFraw_Cstr, SF_Cstr, SFsrc_Cstr, Itable);
return(answer);
}
Integral_result SF_4p_opt (DP k, DP req_prec, int Npts_w, int Npts_K, int Npts_W, I_table Itable)
{
stringstream SFraw_stringstream;
string SFraw_string;
SFraw_stringstream << "SF_4p_k_" << k << "_prec_" << req_prec << "_Npts_" << Npts_w << "_"
<< Npts_K << "_" << Npts_W << ".raw";
SFraw_string = SFraw_stringstream.str();
const char* SFraw_Cstr = SFraw_string.c_str();
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "SF_4p_k_" << k << "_prec_" << req_prec << "_Npts_" << Npts_w << "_"
<< Npts_K << "_" << Npts_W << ".dat";
SF_string = SF_stringstream.str();
const char* SF_Cstr = SF_string.c_str();
stringstream SFsrc_stringstream;
string SFsrc_string;
SFsrc_stringstream << "SF_4p_k_" << k << "_prec_" << req_prec << "_Npts_" << Npts_w << "_"
<< Npts_K << "_" << Npts_W << ".src";
SFsrc_string = SFsrc_stringstream.str();
const char* SFsrc_Cstr = SFsrc_string.c_str();
ofstream SFraw_outfile;
SFraw_outfile.open(SFraw_Cstr);
ofstream SFsrc_outfile;
SFsrc_outfile.open(SFsrc_Cstr);
Vect_DP args_to_SF(5);
args_to_SF[0] = k;
args_to_SF[1] = 0.0; // integration variable
args_to_SF[2] = req_prec;
args_to_SF[3] = DP(Npts_K);
args_to_SF[4] = DP(Npts_W);
// Version using omega = omegamin + (omegamax - omegamin) * (1-cos(alpha))
Integral_result answer = Integrate_optimal_using_table (SF_4p_kwKW_alpha_opt, args_to_SF, 1,
Itable, 0.0, 0.5*PI, req_prec, 1.0e-32, Npts_w, SFraw_outfile);
// Version using omega = omegamin * cosh(alpha)
//Integral_result answer = Integrate_optimal_using_table (SF_4p_kwKW_cosh_alpha_opt, args_to_SF, 1, Itable, 0.0,
//acosh(wmax_4p(k)/wmin_4p(k)), req_prec, 1.0e-32, Npts_w, SFraw_outfile);
answer.integ_est /= twoPI;
answer.abs_prec /= twoPI;
SFraw_outfile.close();
SFsrc_outfile << answer << endl;
SFsrc_outfile.close();
// Translate raw data into SF_4p (k,omega) data
Translate_raw_4p_data (k, answer.n_vals, SFraw_Cstr, SF_Cstr, SFsrc_Cstr, Itable);
return(answer);
}
Integral_result SF_4p_opt (DP k, DP req_prec, int Npts_w, int Npts_KW, I_table Itable)
{
return(SF_4p_opt (k, req_prec, Npts_w, Npts_KW, Npts_KW, Itable));
}
Integral_result SF_4p_opt (DP k, DP req_prec, int Npts, I_table Itable)
{
return(SF_4p_opt (k, req_prec, Npts, Npts, Npts, Itable));
}
//******************************** Functions to produce files similar to ABACUS **********************************
void Produce_SF_2p_file (int N, int Nomega, DP omegamax, I_table Itable)
{
// IMPORTANT NOTE: this produces a file with Szz !!
DP ABACUS_factor = 1.0;
if (N % 2) ABACUSerror("Please use N even in Produce_SF_2p_file.");
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "SF_2p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << ".dat";
SF_string = SF_stringstream.str();
const char* SF_Cstr = SF_string.c_str();
Write_K_File (N, 0, N);
Write_Omega_File (Nomega, 0.0, omegamax);
int dim_K = N/2 + 1;
DP* K = new DP[dim_K];
for (int iK = 0; iK < dim_K; ++iK) K[iK] = (twoPI * iK)/N;
DP* omega = new DP[Nomega];
for (int iw = 0; iw < Nomega; ++iw) omega[iw] = omegamax * (iw + 0.5)/Nomega;
DP* SF_2p_dat = new DP[dim_K * Nomega];
DP srtot = 0.0;
Vect_DP sr1(0.0, dim_K);
for (int iK = 0; iK < dim_K; ++iK)
for (int iw = 0; iw < Nomega; ++iw) {
SF_2p_dat[dim_K * iw + iK] = ABACUS_factor * SF_2p (K[iK], omega[iw], Itable);
srtot += (iK == N/2 ? 1.0 : 2.0) * SF_2p_dat[dim_K * iw + iK];
sr1[iK] += omega[iw] * SF_2p_dat[dim_K * iw + iK];
}
ofstream SF_outfile;
SF_outfile.open(SF_Cstr);
SF_outfile.precision(14);
for (int iw = 0; iw < Nomega; ++iw) {
for (int iK = 0; iK < dim_K; ++iK)
SF_outfile << SF_2p_dat[dim_K * iw + iK] << "\t";
for (int iKt = dim_K - 2; iKt >= 0; --iKt) // put K in [PI, 2PI] back in
SF_outfile << SF_2p_dat[dim_K * iw + iKt] << "\t";
SF_outfile << endl;
}
SF_outfile.close();
// Do sum rule files:
stringstream SRC_stringstream;
string SRC_string;
SRC_stringstream << "SF_2p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << ".src";
SRC_string = SRC_stringstream.str();
const char* SRC_Cstr = SRC_string.c_str();
ofstream SRC_outfile;
SRC_outfile.open(SRC_Cstr);
SRC_outfile.precision(14);
SRC_outfile << srtot * omegamax/(twoPI * Nomega * N) << "\t" << srtot * 4.0 * omegamax/(twoPI * Nomega * N) << endl;
SRC_outfile.close();
stringstream SR1_stringstream;
string SR1_string;
SR1_stringstream << "SF_2p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << ".sr1";
SR1_string = SR1_stringstream.str();
const char* SR1_Cstr = SR1_string.c_str();
ofstream SR1_outfile;
SR1_outfile.open(SR1_Cstr);
SR1_outfile.precision(14);
for (int iK = 1; iK < dim_K; ++iK)
SR1_outfile << iK << "\t" << K[iK] << "\t" << sr1[iK] * omegamax/(twoPI * Nomega)
<< "\t" << -((1.0 - cos(K[iK])) * 2.0 * (0.25 - log(2.0))/3.0) << "\t"
<< -sr1[iK] * omegamax/(twoPI * Nomega)/((1.0 - cos(K[iK])) * 2.0 * (0.25 - log(2.0))/3.0) << endl;
SR1_outfile.close();
return;
}
void Produce_SF_4p_file (int N, int Nomega, DP omegamax, DP req_prec, int max_rec, I_table Itable)
{
// IMPORTANT NOTE: this produces a file with the same normalization as Smp from ABACUS,
// so we use a factor of 2 (for Szz -> Smp) and 2PI.
DP ABACUS_factor = 1.0;
if (N % 2) ABACUSerror("Please use N even in Produce_SF_2p_file.");
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "SF_4p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << "_prec_" << req_prec
<< "_max_rec_" << max_rec << ".dat";
SF_string = SF_stringstream.str();
const char* SF_Cstr = SF_string.c_str();
Write_K_File (N, 0, N);
Write_Omega_File (Nomega, 0.0, omegamax);
int dim_K = N/2 + 1;
DP* K = new DP[dim_K];
for (int iK = 0; iK < dim_K; ++iK) K[iK] = (twoPI * iK)/N;
DP* omega = new DP[Nomega];
for (int iw = 0; iw < Nomega; ++iw) omega[iw] = omegamax * (iw + 0.5)/Nomega;
DP* SF_4p_dat = new DP[dim_K * Nomega];
DP srtot = 0.0;
Vect_DP sr1(0.0, dim_K);
for (int iK = 0; iK < dim_K; ++iK)
for (int iw = 0; iw < Nomega; ++iw) {
SF_4p_dat[dim_K * iw + iK] = ABACUS_factor * SF_4p_rec (K[iK], omega[iw], req_prec, max_rec, Itable);
srtot += (iK == N/2 ? 1.0 : 2.0) * SF_4p_dat[dim_K * iw + iK];
sr1[iK] += omega[iw] * SF_4p_dat[dim_K * iw + iK];
}
ofstream SF_outfile;
SF_outfile.open(SF_Cstr);
SF_outfile.precision(14);
for (int iw = 0; iw < Nomega; ++iw) {
for (int iK = 0; iK < dim_K; ++iK)
SF_outfile << SF_4p_dat[dim_K * iw + iK] << "\t";
for (int iKt = dim_K - 2; iKt >= 0; --iKt) // put K in [PI, 2PI] back in
SF_outfile << SF_4p_dat[dim_K * iw + iKt] << "\t";
SF_outfile << endl;
}
SF_outfile.close();
// Do sum rule files:
stringstream SRC_stringstream;
string SRC_string;
SRC_stringstream << "SF_4p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << "_prec_" << req_prec
<< "_max_rec_" << max_rec << ".src";
SRC_string = SRC_stringstream.str();
const char* SRC_Cstr = SRC_string.c_str();
ofstream SRC_outfile;
SRC_outfile.open(SRC_Cstr);
SRC_outfile.precision(14);
SRC_outfile << srtot * omegamax/(twoPI * Nomega * N) << "\t" << srtot * 4.0 * omegamax/(twoPI * Nomega * N) << endl;
SRC_outfile.close();
stringstream SR1_stringstream;
string SR1_string;
SR1_stringstream << "SF_4p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << "_prec_" << req_prec
<< "_max_rec_" << max_rec << ".sr1";
SR1_string = SR1_stringstream.str();
const char* SR1_Cstr = SR1_string.c_str();
ofstream SR1_outfile;
SR1_outfile.open(SR1_Cstr);
SR1_outfile.precision(14);
for (int iK = 1; iK < dim_K; ++iK)
SR1_outfile << iK << "\t" << K[iK] << "\t" << sr1[iK] * omegamax/(twoPI * Nomega)
<< "\t" << -((1.0 - cos(K[iK])) * 2.0 * (0.25 - log(2.0))/3.0) << "\t"
<< -sr1[iK] * omegamax/(twoPI * Nomega)/((1.0 - cos(K[iK])) * 2.0 * (0.25 - log(2.0))/3.0) << endl;
SR1_outfile.close();
return;
}
void Produce_SF_4p_file (int N, int Nomega, DP omegamax, DP req_prec, int Npts_K, int Npts_W, I_table Itable)
{
// IMPORTANT NOTE: this produces a file with the same normalization as Smp from ABACUS,
// so we use a factor of 2 (for Szz -> Smp) and 2PI.
DP ABACUS_factor = 1.0;
if (N % 2) ABACUSerror("Please use N even in Produce_SF_2p_file.");
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "SF_4p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << "_prec_" << req_prec
<< "_Npts_K_" << Npts_K << "_Npts_W_" << Npts_W << ".dat";
SF_string = SF_stringstream.str();
const char* SF_Cstr = SF_string.c_str();
Write_K_File (N, 0, N);
Write_Omega_File (Nomega, 0.0, omegamax);
int dim_K = N/2 + 1;
DP* K = new DP[dim_K];
for (int iK = 0; iK < dim_K; ++iK) K[iK] = (twoPI * iK)/N;
DP* omega = new DP[Nomega];
for (int iw = 0; iw < Nomega; ++iw) omega[iw] = omegamax * (iw + 0.5)/Nomega;
DP* SF_4p_dat = new DP[dim_K * Nomega];
DP srtot = 0.0;
Vect_DP sr1(0.0, dim_K);
for (int iK = 0; iK < dim_K; ++iK)
for (int iw = 0; iw < Nomega; ++iw) {
SF_4p_dat[dim_K * iw + iK] = ABACUS_factor * SF_4p_opt (K[iK], omega[iw], req_prec, Npts_K, Npts_W, Itable);
srtot += (iK == N/2 ? 1.0 : 2.0) * SF_4p_dat[dim_K * iw + iK];
sr1[iK] += omega[iw] * SF_4p_dat[dim_K * iw + iK];
}
// Output SF:
ofstream SF_outfile;
SF_outfile.open(SF_Cstr);
SF_outfile.precision(14);
for (int iw = 0; iw < Nomega; ++iw) {
for (int iK = 0; iK < dim_K; ++iK)
SF_outfile << SF_4p_dat[dim_K * iw + iK] << "\t";
for (int iKt = dim_K - 2; iKt >= 0; --iKt) // put K in [PI, 2PI] back in
SF_outfile << SF_4p_dat[dim_K * iw + iKt] << "\t";
SF_outfile << endl;
}
SF_outfile.close();
// Do sum rule files:
stringstream SRC_stringstream;
string SRC_string;
SRC_stringstream << "SF_4p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << "_prec_" << req_prec
<< "_Npts_K_" << Npts_K << "_Npts_W_" << Npts_W << ".src";
SRC_string = SRC_stringstream.str();
const char* SRC_Cstr = SRC_string.c_str();
ofstream SRC_outfile;
SRC_outfile.open(SRC_Cstr);
SRC_outfile.precision(14);
SRC_outfile << srtot * omegamax/(twoPI * Nomega * N) << "\t" << srtot * 4.0 * omegamax/(twoPI * Nomega * N) << endl;
SRC_outfile.close();
stringstream SR1_stringstream;
string SR1_string;
SR1_stringstream << "SF_4p_N_" << N << "_Nw_" << Nomega << "_wmax_" << omegamax << "_prec_" << req_prec
<< "_Npts_K_" << Npts_K << "_Npts_W_" << Npts_W << ".sr1";
SR1_string = SR1_stringstream.str();
const char* SR1_Cstr = SR1_string.c_str();
ofstream SR1_outfile;
SR1_outfile.open(SR1_Cstr);
SR1_outfile.precision(14);
for (int iK = 1; iK < dim_K; ++iK)
SR1_outfile << iK << "\t" << K[iK] << "\t" << sr1[iK] * omegamax/(twoPI * Nomega)
<< "\t" << -((1.0 - cos(K[iK])) * 2.0 * (0.25 - log(2.0))/3.0) << "\t"
<< -sr1[iK] * omegamax/(twoPI * Nomega)/((1.0 - cos(K[iK])) * 2.0 * (0.25 - log(2.0))/3.0) << endl;
SR1_outfile.close();
return;
}
} // namespace ABACUS
Reference in New Issue View Git Blame Copy Permalink