Renormalized Susceptibility Code
D. Stanford
<![CDATA[
#include <iostream>
#include <math.h>
#include <fstream>
#include <unistd.h>
//VARIABLES
//number of ky momentum windows for the integral
int N=37;
//number of energy windows
const int M=1000;
float pi=3.14159265,w=1,qx=pi-.02,qy=pi-.02,eps=0.01,Emax=10,Emin=-10,dk=pi/float(N),dw=(Emax-Emin)/float(M);
float t=1;
float tprime=1;
float U=0;
float Delta=0;
int counter=0;
//STRUCTS
struct douglas{
float w;
float c;
float r;
}output[M+1];
//FUNCTIONS
float E(float kx, float ky);
float MatrixElement(float kx, float ky, float qx, float qy);
float Grad(float kx, float ky);
float L(float Edr,float Edl, float Eur, float Eul);
void ComputeBareSusceptibility();
void MakeRetarded();
void LadderSum();
void HilbertTransform();
void Outfile();
float SumRule();
float Grad2(float kx,float ky);
void InfileSusceptibility();
using namespace std;
int main(){
// InfileSusceptibility();
ComputeBareSusceptibility();
MakeRetarded();
HilbertTransform();
LadderSum();
Outfile();
return 0;
}
float E(float kx, float ky){
return(sqrt(4*t*t*(pow(cos(kx),2)+pow(cos(ky),2))+pow(Delta+4*tprime*sin(kx)*sin(ky),2)));
}
float MatrixElement(float kx, float ky, float qx, float qy){
return((E(kx+qx,ky+qy)*E(kx,ky)-4*t*t*(cos(kx+qx)*cos(kx)+cos(ky+qy)*cos(ky))-16*tprime*tprime*sin(kx+qx)*sin(ky+qy)*sin(kx)*sin(ky))/(2*E(kx,ky)*E(kx+qx,ky+qy)));
}
float L(float Edr, float Edl, float Eur, float Eul){
float t1=Edr/(Edr-Edl),t2=Edl/(Edl-Eul),t3=Eul/(Eul-Eur),t4=Eur/(Eur-Edr);
if((0<t1&&t1<1)&&(0<t2&&t2<1)){return(dk*pow(pow(1-t1,2)+pow(t2,2),.5));}
else if((0<t1&&t1<1)&&(0<t3&&t3<1)){return(dk*pow(pow(1-t1-t3,2)+1,.5));}
else if((0<t1&&t1<1)&&(0<t4&&t4<1)){return(dk*pow(pow(1-t4,2)+t1*t1,.5));}
else if((0<t2&&t2<1)&&(0<t3&&t3<1)){return(dk*pow(pow(1-t2,2)+pow(t3,2),.5));}
else if((0<t2&&t2<1)&&(0<t4&&t4<1)){return(dk*pow(pow(1-t4-t2,2)+1,.5));}
else if((0<t3&&t3<1)&&(0<t4&&t4<1)){return(dk*pow(pow(1-t3,2)+pow(t4,2),.5));}
else { return(0);}
}
float Grad(float kx, float ky){
return(.5*pow((pow(8*t*t*cos(ky)*sin(ky)-8*tprime*cos(ky)*sin(kx)*(Delta+4*tprime*sin(kx)*sin(ky)),2)+pow(8*t*t*cos(kx)*sin(kx)-8*tprime*cos(kx)*sin(ky)*(Delta+4*tprime*sin(kx)*sin(ky)),2))/(4*t*t*(pow(cos(kx),2)+pow(cos(ky),2))+pow(Delta+4*tprime*sin(kx)*sin(ky),2)),.5));
}
void ComputeBareSusceptibility(){
//loop over output momenta
for(float w=Emin;w<Emax;w+=dw){
counter+=1;
output[counter].w=w;
//loop over B.Z.
for(float kx=-pi; kx<pi; kx+=dk){
for(float ky=-pi; ky<pi; ky+=dk){
float en1=0, en2=0;
float Edr=E(kx+dk+qx,ky+qy)+E(kx+dk,ky)-w;
float Edl=E(kx+qx,ky+qy)+E(kx,ky)-w;
float Eur=E(kx+dk+qx,ky+dk+qy)+E(kx+dk,ky+dk)-w;
float Eul=E(kx+qx,ky+dk+qy)+E(kx,ky+dk)-w;
en1=-L(Edr,Edl,Eur,Eul);
Edr=Edr+2*w;
Edl=Edl+2*w;
Eur=Eur+2*w;
Eul=Eul+2*w;
en2=-L(Edr,Edl,Eur,Eul);
output[counter].c+=((en1+en2)*2/(8*pi*Grad2(kx,ky)));
}
}
}
}
void MakeRetarded(){
//Reverse sign of negative-frequency imaginary part
for(int i=1; i<=counter; i++){
if(output[i].w<0){
output[i].c=-output[i].c;
}
}
}
void HilbertTransform(){
//Hilbert Transform
for(int i=1; i<=counter; i++){
float sum=0;
for(int j=1; j<=counter; j++){
if(i!=j){
sum+=output[j].c*(1/(pi*(j-i)));
}
}
output[i].r=sum;
}
}
void LadderSum(){
//ladder sum
for(int i=1; i<=M; i++){
float r=output[i].r,im=output[i].c,r2=(1+.5*U*r)/(pow(1+.5*U*r,2)+pow(.5*U*im,2));
float im2=-.5*U*im/(pow(1+.5*U*r,2)+pow(.5*U*im,2));
output[i].r=r*r2-im*im2;
output[i].c=r*im2+r2*im;
}
}
void Outfile(){
ofstream outfile;
outfile.open("output.txt");
for(int i=1; i<=M; i++){
outfile<<output[i].w<<'\t'<<output[i].c<<'\t'<<output[i].r<<endl;
}
outfile.close();
}
float SumRule(){
float sum=0;
for(int i=1; i<=counter; i++){
if(output[i].c<40){
sum+=output[i].c*dw;}
}
return(sum);
}
float Grad2(float kx,float ky){
return(pow(pow((E(kx+eps/2,ky)+E(kx+qx+eps/2,ky+qy)-E(kx-eps/2,ky)-E(kx+qx-eps/2,ky+qy))/eps,2)+pow((E(kx,ky+eps/2)+E(kx+qx,ky+qy+eps/2)-E(kx+qx,ky+qy-eps/2)-E(kx,ky-eps/2))/eps,2),.5));
}
void InfileSusceptibility(){
ifstream infile;
infile.open("bare_suscep.txt");
for(int i=1; i<=M; i++){
infile>>output[i].w;
infile>>output[i].c;
infile>>output[i].r;
}
}
]]>