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

627 lines
20 KiB
C++
Raw Permalink Blame History

/**********************************************************
This software is part of J.-S. Caux's ABACUS library.
Copyright (c) J.-S. Caux.
-----------------------------------------------------------
File: XXZ_VOA.cc
Purpose: Defines all class procedures used for the infinite XXZ chain in zero field.
***********************************************************/
#include "ABACUS.h"
using namespace std;
namespace ABACUS {
DP Fermi_velocity_XXZ_h0 (DP Delta)
{
// This uses the J > 0 conventions, with H = J \sum S.S
if (Delta <= -1.0 || Delta >= 1.0) ABACUSerror("Use Delta in AFM gapless regime for Fermi_velocity_XXZ_h0.");
return(0.5 * PI * sqrt(1.0 - Delta * Delta)/acos(Delta));
}
inline DP Integrand_xi_11 (Vect_DP args)
{
// args[0] corresponds to t, args[1] to rho, args[2] to xi
// For increased accuracy: separate in regions 0 < t <= 1 and 1 < t, and 0 < xi t <= 1 and 1 < xi t:
DP answer;
if (args[0] * args[2] <= 1.0) {
if (args[0] <= 1.0) {
answer = sinh(args[0] * (1.0 + args[2])) * sinh(args[0]) * cos(4.0 * args[0] * args[1])
/(args[0] * sinh(args[0] * args[2]) * pow(cosh(args[0]), 2.0));
}
else { // factorize exp t
DP expm2t = exp(-2.0*args[0]);
answer = sinh(args[0] * (1.0 + args[2])) * 2.0 * (1.0 - expm2t) * cos(4.0 * args[0] * args[1])
/(args[0] * sinh(args[0] * args[2]) * (1.0 + expm2t) * (1.0 + expm2t));
}
}
else { // factorize exp (xi t)
if (args[0] <= 1.0) {
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0*args[0]*args[2]);
answer = (1.0 - expm2t * expm2txi) * sinh(args[0]) * cos(4.0 * args[0] * args[1])
/(args[0] * (1.0 - expm2txi) * 0.5 * (1.0 + expm2t) * cosh(args[0]));
}
else { // factorize exp t
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0*args[0]*args[2]);
answer = 2.0 * (1.0 - expm2t * expm2txi) * (1.0 - expm2t) * cos(4.0 * args[0] * args[1])
/(args[0] * (1.0 - expm2txi) * (1.0 + expm2t) * (1.0 + expm2t));
}
}
return(answer);
}
DP I_xi_11_integral (DP xi, DP rho, DP req_prec, int max_nr_pts, DP t1bar)
{
DP rho_used = fabs(rho);
Vect_DP args(3);
args[0] = 0.0;
args[1] = rho_used;
args[2] = xi;
return((Integrate_optimal (Integrand_xi_11, args, 0, 0.0, t1bar, req_prec, req_prec, max_nr_pts)).integ_est);
}
inline DP Integrand_xi_12 (Vect_DP args)
{
// args[0] corresponds to t, args[1] to rho, args[2] to xi
DP expm2t = exp(-2.0*args[0]);
DP expm2txi = exp(-2.0*args[0]*args[2]);
return(2.0 * cos(4.0 * args[0] * args[1]) * (-expm2t * (3.0 + expm2t) + expm2txi * (1.0 + expm2t * (1.0 + 2.0*expm2t)))
/ (args[0] * (1.0 - expm2txi) * (1.0 + expm2t) * (1.0 + expm2t)));
}
DP I_xi_12_integral (DP xi, DP rho, DP req_prec, int max_nr_pts, DP t1bar)
{
DP tinf = 24.0; // such that exp(-2tinf) << machine_eps
DP rho_used = fabs(rho);
Vect_DP args(3);
args[0] = 0.0;
args[1] = rho_used;
args[2] = xi;
return((Integrate_optimal (Integrand_xi_12, args, 0, t1bar, tinf, req_prec, req_prec, max_nr_pts)).integ_est);
}
/*
inline DP Integrand_xi_2 (Vect_DP args)
{
// This version is used for rho <= 1
// args[0] corresponds to t, args[1] to rho, args[2] to xi
DP answer;
if (args[0] * args[2] <= 1.0) {
if (args[0] <= 1.0) {
answer = sinh(args[0] * (args[2] + 1.0)) * pow(sin(2.0 * args[0] * args[1]), 2.0)
/(args[0] * sinh(args[2] * args[0]) * sinh(args[0]) * pow(cosh(args[0]), 2.0));
}
else if (args[0] >= 1.0) {
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0*args[0]*args[2]);
//answer = 8.0 * expm2t * pow(sin(2.0 * args[0] * args[1]), 2.0)
/(args[0] * (1.0 - expm2t) * (1.0 + expm2t) * (1.0 + expm2t));
answer = 8.0 * (1.0 - expm2t*expm2txi) * expm2t *
pow(sin(2.0 * args[0] * args[1]), 2.0)/(args[0] * (1.0 - expm2txi) * (1.0 - expm2t)
* (1.0 + expm2t) * (1.0 + expm2t));
}
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0*args[0]*args[2]);
//answer = 8.0 * expm2t * pow(sin(2.0 * args[0] * args[1]), 2.0)
/(args[0] * (1.0 - expm2t) * (1.0 + expm2t) * (1.0 + expm2t));
answer = 8.0 * (1.0 - expm2t*expm2txi) * expm2t *
pow(sin(2.0 * args[0] * args[1]), 2.0)/(args[0] * (1.0 - expm2txi) * (1.0 - expm2t)
* (1.0 + expm2t) * (1.0 + expm2t));
return(answer);
... NOT USED
}
*/
inline DP Integrand_xi_22 (Vect_DP args)
{
// This version is used when rho > 1.
// args[0] corresponds to t, args[1] to rho, args[2] to xi
DP answer = 0.0;
if (args[0] < 1.0)
answer = pow(sin(2.0 * args[0] * args[1]), 2.0)/(args[0] * tanh(args[0] * args[2]) * pow(cosh(args[0]), 2.0));
else if (args[0] >= 1.0) {
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0*args[0]*args[2]);
answer = 4.0 * expm2t * (1.0 + expm2txi) * pow(sin(2.0 * args[0] * args[1]), 2.0)
/(args[0] * (1.0 - expm2txi) * (1.0 + expm2t) * (1.0 + expm2t));
}
return(answer);
}
DP I_xi_2_integral (DP xi, DP rho, DP req_prec, int max_nr_pts)
{
DP tinf = 24.0; // such that exp(-2tinf) << machine_eps
DP rho_used = fabs(rho);
Vect_DP args(3);
args[0] = 0.0;
args[1] = rho_used;
args[2] = xi;
DP answer = 0.0;
answer = PI * rho + log(0.5 * (1.0 + exp(-2.0 * PI * rho_used))) // This is I^{(21)}
+ (Integrate_optimal (Integrand_xi_22, args, 0, 0.0, tinf, req_prec, req_prec, max_nr_pts)).integ_est;
return(answer);
}
DP I_xi_integral (DP xi, DP rho, DP req_prec, int max_nr_pts)
{
DP t1bar = 1.0; // delimiter between the I_xi_11 and I__xi_12 integrals
DP answer = I_xi_11_integral (xi, rho, req_prec, max_nr_pts, t1bar)
+ I_xi_12_integral (xi, rho, req_prec, max_nr_pts, t1bar)
- 2.0 * Cosine_Integral (4.0 * rho * t1bar)
- I_xi_2_integral (xi, rho, req_prec, max_nr_pts);
return(answer);
}
/********************* TWO SPINONS ********************/
DP Szz_XXZ_h0_2spinons (DP Delta, DP k, DP omega, Integral_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 vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP wu = 4.0 * vF * sin(0.5 * k);
DP wl = 2.0 * vF * fabs(sin(k));
DP rho = acosh(sqrt((wu * wu - wl * wl)/(w * w - wl * wl)))/PI;
DP xi = PI/acos(Delta) - 1.0;
// 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^{+-}
DP expmtwoPIrhooverxi = exp(-twoPI * rho/xi);
return(w < wu && w > wl ? 2.0 * pow(1.0 + 1.0/xi, 2.0) * exp(-Itable.Return_val (rho))
* expmtwoPIrhooverxi/(sqrt(wu * wu - w * w)
* (0.5 * (1.0 + expmtwoPIrhooverxi * expmtwoPIrhooverxi) + expmtwoPIrhooverxi * cos(PI/xi)))
: 0.0);
}
DP Szz_XXZ_h0_2spinons (Vect_DP args, Integral_table Itable)
{
// To be called by integrating functions
// Careful ! This is S(k, omega) !
// This uses args[0] = k, args[1] = omega, args[2] = xi
DP Delta = cos(PI/(args[2] + 1.0));
return(Szz_XXZ_h0_2spinons (Delta, args[0], args[1], Itable));
}
DP Szz_XXZ_h0_2spinons_alt (Vect_DP args, Integral_table Itable)
{
// This returns S(k, omega) \times Jacobian (d omega)/(d alpha) where omega = omega_{lower threshold} cosh \alpha
// so returns S(k, omega) \sqrt{omega^2 - omega_l^2}
// This uses args[0] = k, args[1] = alpha, args[2] = xi
DP Delta = cos(PI/(args[2] + 1.0));
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegaup = 2.0 * vF * sin(0.5 * args[0]);
DP omegalow = vF * fabs(sin(args[0]));
DP omega = omegalow * cosh(args[1]);
if (omega >= omegaup || omega <= omegalow) return(0.0);
DP Jacobian = sqrt(omega * omega - omegalow * omegalow);
return(Jacobian * Szz_XXZ_h0_2spinons (Delta, args[0], omega, Itable));
}
DP Szz_XXZ_h0_2spinons_omega (Vect_DP args, Integral_table Itable)
{
// This uses args[0] = k, args[1] = omega, args[2] = xi
return(args[1] * Szz_XXZ_h0_2spinons (args, Itable));
}
DP Szz_XXZ_h0_2spinons_omega_alt (Vect_DP args, Integral_table Itable)
{
// This uses args[0] = k, args[1] = alpha, args[2] = xi
DP Delta = cos(PI/(args[2] + 1.0));
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegalow = vF * fabs(sin(args[0]));
DP omega = omegalow * cosh(args[1]);
args[1] = omega;
return(Szz_XXZ_h0_2spinons_omega (args, Itable));
}
DP Szz_XXZ_h0_2spinons_intomega (Vect_DP args, Integral_table Itable)
{
// This returns \int_0^2PI domega/2PI S(k, omega)
DP k = args[0];
DP req_prec = args[1];
DP xi = args[2];
int max_nr_pts = int(args[3]);
DP Delta = cos(PI/(args[2] + 1.0));
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegaup = 2.0 * vF * sin(0.5 * args[0]);
DP omegalow = vF * fabs(sin(args[0]));
Vect_DP args_to_SF_2p(4);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be omega
args_to_SF_2p[2] = xi;
args_to_SF_2p[3] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return((Integrate_optimal_using_table (Szz_XXZ_h0_2spinons, args_to_SF_2p, 1, Itable,
omegalow, omegaup, req_prec, req_prec, max_nr_pts)).integ_est/twoPI);
}
DP Szz_XXZ_h0_2spinons_intomega_alt (Vect_DP args, Integral_table Itable)
{
// This returns \int_0^2PI domega/2PI S(k, omega)
DP k = args[0];
DP req_prec = args[1];
DP xi = args[2];
int max_nr_pts = int(args[3]);
DP Delta = cos(PI/(args[2] + 1.0));
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegaup = 2.0 * vF * sin(0.5 * args[0]);
DP omegalow = vF * fabs(sin(args[0]));
Vect_DP args_to_SF_2p(4);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be omega
args_to_SF_2p[2] = xi;
args_to_SF_2p[3] = 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, acos(wl/wu), req_prec, max_rec)/twoPI);
return((Integrate_optimal_using_table (Szz_XXZ_h0_2spinons_alt, args_to_SF_2p, 1, Itable,
0.0, acosh(omegaup/omegalow), req_prec, req_prec, max_nr_pts)).integ_est/twoPI);
}
DP Szz_XXZ_h0_2spinons_check_sumrule (DP Delta, DP req_prec, int max_nr_pts, Integral_table Itable)
{
// Near XX: it's better to use this function.
// Near XXX: it's better to use the ..._alt function below.
Vect_DP args_to_Szz_intomega (4);
args_to_Szz_intomega[0] = 0.0; // this will be k
args_to_Szz_intomega[1] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_Szz_intomega[2] = PI/acos(Delta) - 1.0;
args_to_Szz_intomega[3] = DP(max_nr_pts);
// Factor 4: we expect \int_0^\infty \frac{d\omega}{2\pi} \int_0^{2\pi} \frac{dk}{2\pi} S^zz (k, omega) = 1/4
// Factor 2: int[0, 2PI] = 2 int[0, PI]
return(4.0 * 2.0 * (Integrate_optimal_using_table (Szz_XXZ_h0_2spinons_intomega, args_to_Szz_intomega, 0, Itable,
0.0, PI, req_prec, req_prec, max_nr_pts)).integ_est/twoPI);
}
DP Szz_XXZ_h0_2spinons_check_sumrule_alt (DP Delta, DP req_prec, int max_nr_pts, Integral_table Itable)
{
// This is the preferred version if near XXX.
Vect_DP args_to_Szz_intomega (4);
args_to_Szz_intomega[0] = 0.0; // this will be k
args_to_Szz_intomega[1] = ABACUS::max(1.0e-14, 0.01 * req_prec);
args_to_Szz_intomega[2] = PI/acos(Delta) - 1.0;
args_to_Szz_intomega[3] = DP(max_nr_pts);
// Factor 2: int[0, 2PI] = 2 int[0, PI]
return(4.0 * 2.0 * (Integrate_optimal_using_table (Szz_XXZ_h0_2spinons_intomega_alt, args_to_Szz_intomega, 0, Itable,
0.0, PI, req_prec, req_prec, max_nr_pts)).integ_est/twoPI);
}
DP Fixed_k_sumrule_omega_Szz_XXZ_h0_N (DP Delta, DP k)
{
// This is K_1 (k) = \int domega/2PI omega S(k, omega)
//return(4.0 * (log(2.0) - 0.25) * (1.0 - cos(k))/3.0); // XXX_h0
// Calculate the <S^x S^x> average:
int N = 600;
int M = 300;
DP Xavg = X_avg ('x', Delta, N, M);
return (-2.0 * (1.0 - cos(k)) * Xavg/N);
}
DP Integrand_GSE_XXZ_h0 (Vect_DP args)
{
// args[0] corresponds to t, args[1] to xi
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0 * args[0] * args[1]);
return(2.0 * expm2t * (1.0 - expm2txi)/((1.0 - expm2t * expm2txi) * (1.0 + expm2t)));
}
DP GSE_XXZ_h0 (DP Delta, DP req_prec, int max_nr_pts)
{
DP xi = PI/acos(Delta) - 1.0;
DP tinf = 24.0; // such that exp(-2tinf) << machine_eps
Vect_DP args(2);
args[0] = 0.0;
args[1] = xi;
return(-0.25 * sqrt(1.0 - Delta * Delta) * (2.0/acos(Delta)) * 2.0
* (Integrate_optimal (Integrand_GSE_XXZ_h0, args, 0, 0.0, tinf, req_prec, req_prec, max_nr_pts).integ_est));
}
DP Integrand_2_fSR_XXZ_h0 (Vect_DP args)
{
// args[0] corresponds to t, args[1] to xi
DP expm2t = exp(-2.0 * args[0]);
DP expm2txi = exp(-2.0 * args[0] * args[1]);
return(4.0 * args[0] * expm2t * (1.0 + expm2t * expm2txi)
/((1.0 - expm2t * expm2txi) * (1.0 + expm2t) * (1.0 + expm2t)));
}
DP Fixed_k_sumrule_omega_Szz_XXZ_h0 (DP Delta, DP k, DP req_prec, int max_nr_pts)
{
//
DP xi = PI/acos(Delta) - 1.0;
DP tinf = 24.0; // such that exp(-2tinf) << machine_eps
Vect_DP args(2);
args[0] = 0.0;
args[1] = xi;
DP Xavg = -(0.25/(acos(Delta) * sqrt(1.0 - Delta * Delta))) * 2.0
* (Integrate_optimal (Integrand_GSE_XXZ_h0, args, 0, 0.0, tinf, req_prec, req_prec, max_nr_pts).integ_est)
+ 0.25 * (Delta/pow(acos(Delta), 2.0)) * 2.0
* (Integrate_optimal (Integrand_2_fSR_XXZ_h0, args, 0, 0.0, tinf, req_prec, req_prec, max_nr_pts).integ_est);
return (-2.0 * (1.0 - cos(k)) * Xavg);
}
DP Szz_XXZ_h0_2spinons_check_fixed_k_Szz_sumrule (DP Delta, DP k, DP req_prec, int max_nr_pts, Integral_table Itable)
{
// Near XXX, it's better to use the ..._alt function below.
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegaup = 2.0 * vF * sin(0.5 * k);
DP omegalow = vF * fabs(sin(k));
Vect_DP args_to_SF_2p(4);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be omega
args_to_SF_2p[2] = PI/acos(Delta) - 1.0;
args_to_SF_2p[3] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return(((Integrate_optimal_using_table (Szz_XXZ_h0_2spinons_omega, args_to_SF_2p, 1, Itable,
omegalow, omegaup, req_prec, req_prec, max_nr_pts)).integ_est/twoPI)
/Fixed_k_sumrule_omega_Szz_XXZ_h0(Delta, k, req_prec, max_nr_pts));
}
DP Szz_XXZ_h0_2spinons_check_fixed_k_Szz_sumrule_alt (DP Delta, DP k, DP req_prec, int max_nr_pts, Integral_table Itable)
{
// This is the preferred version if near XXX.
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegaup = 2.0 * vF * sin(0.5 * k);
DP omegalow = vF * fabs(sin(k));
Vect_DP args_to_SF_2p(4);
args_to_SF_2p[0] = k;
args_to_SF_2p[1] = 0.0; // this will be alpha
args_to_SF_2p[2] = PI/acos(Delta) - 1.0;
args_to_SF_2p[3] = ABACUS::max(1.0e-14, 0.01 * req_prec);
return(((Integrate_optimal_using_table (Szz_XXZ_h0_2spinons_omega_alt, args_to_SF_2p, 1, Itable,
0.0, acosh(omegaup/omegalow), req_prec, req_prec, max_nr_pts)).integ_est/twoPI)
/Fixed_k_sumrule_omega_Szz_XXZ_h0(Delta, k, req_prec, max_nr_pts));
}
//******************************** Functions to produce files similar to ABACUS **********************************
void Produce_Szz_XXZ_h0_2spinons_file (DP Delta, int N, int Nomega, DP omegamax, Integral_table Itable)
{
// IMPORTANT NOTE: this produces a file with Szz !!
if (N % 2) ABACUSerror("Please use N even in Produce_Szz_XXZ_h0_2spinons_file.");
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "Szz_XXZ_h0_2spinons_Delta_" << Delta << "_N_" << N << "_Nom_" << Nomega
<< "_ommax_" << 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_XXZ_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_XXZ_2p_dat[dim_K * iw + iK] = Szz_XXZ_h0_2spinons (Delta, K[iK], omega[iw], Itable);
srtot += (iK == N/2 ? 1.0 : 2.0) * SF_XXZ_2p_dat[dim_K * iw + iK];
sr1[iK] += omega[iw] * SF_XXZ_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_XXZ_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_XXZ_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 << "Szz_XXZ_h0_2spinons_Delta_" << Delta << "_N_" << N << "_Nom_" << Nomega
<< "_ommax_" << 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 << "Szz_XXZ_h0_2spinons_Delta_" << Delta << "_N_" << N << "_Nom_" << Nomega
<< "_ommax_" << 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);
// Figure out the f-sumrule factor:
int Nfsr = 600;
int Mfsr = 300;
DP Xavg = 0.0;
if (Delta > 0.0) Xavg = X_avg ('x', Delta, Nfsr, Mfsr);
else Xavg = 0.0;
DP sr1factor = 1.0;
if (Delta > 0.0) sr1factor = -2.0 * Xavg/Nfsr;
else sr1factor = 1.0;
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])) * sr1factor << "\t"
<< sr1[iK] * omegamax/(twoPI * Nomega)/((1.0 - cos(K[iK])) * sr1factor) << endl;
SR1_outfile.close();
return;
}
void Produce_Szz_XXZ_h0_2spinons_fixed_K_file (DP Delta, DP Kover2PI, int Nomega, Integral_table Itable)
{
// IMPORTANT NOTE: this produces a file with Szz !!
stringstream SF_stringstream;
string SF_string;
SF_stringstream << "Szz_XXZ_h0_2spinons_Delta_" << Delta << "_Kover2PI_" << Kover2PI
<< "_Nom_" << Nomega << ".dat";
SF_string = SF_stringstream.str();
const char* SF_Cstr = SF_string.c_str();
DP K = twoPI * Kover2PI;
DP* omega = new DP[Nomega];
DP vF = Fermi_velocity_XXZ_h0 (Delta); // in units of omega
DP omegaup = 2.0 * vF * sin(0.5 * K);
DP omegalow = vF * fabs(sin(K));
for (int iw = 0; iw < Nomega; ++iw)
omega[iw] = 0.99 * omegalow + (1.01 * omegaup - 0.99 * omegalow) * (iw + 0.5)/Nomega;
// factor of 1.01: just to have zeroes on either side of continuum.
DP* SF_XXZ_2p_dat = new DP[Nomega];
DP sr1 = 0.0;
for (int iw = 0; iw < Nomega; ++iw) {
SF_XXZ_2p_dat[iw] = Szz_XXZ_h0_2spinons (Delta, K, omega[iw], Itable);
sr1 += omega[iw] * SF_XXZ_2p_dat[iw];
}
ofstream SF_outfile;
SF_outfile.open(SF_Cstr);
SF_outfile.precision(14);
SF_outfile << omega[0] << "\t" << SF_XXZ_2p_dat[0];
for (int iw = 0; iw < Nomega; ++iw) {
SF_outfile << endl << omega[iw] << "\t" << SF_XXZ_2p_dat[iw];
}
SF_outfile.close();
return;
}
} // namespace ABACUS
Reference in New Issue View Git Blame Copy Permalink