c calculates spin-spin correlation function near pipi, in ladder approximation,
c for positive frequencies, for one value of momentum.
c sets the accuracy of integral and resolution of graphs.
parameter(nq=40)
parameter(ne=2000)
dimension d1(ne),d2(ne),d3(ne),ha2(ne)
complex g1(ne),g2(ne),g3(ne),h(ne),h22(ne),h12(ne)
complex zz
c establishes frequency domain.
emin=0.0
emax=20.0
de=(emax-emin)/ne
c sets the hopping parameter.
c sets will-be-half-the-band-gap, delta.
c sets a momentum.
write(*,88)
88 format('Enter values for hopping, delta, qqx, qqy')
read(*,*) t,delta,qqx,qqy
c establishes a value for pi.
pi=3.14159265
c sets the region of integration.
dq=pi/nq
c NOW solves SDW self-consistency equation for 2-d Hubbard model.
c opens output file.
open(20,file='self-consistency.dat',status='new')
c zeroes out accumulator.
s=0.0
do 1 j=1,nq
qx=(j-0.5)*dq
do 1 k=1,nq
qy=(k-0.5)*dq
e=-2.0*t*(cos(qx)+cos(qy))
s=s+dq**2/((pi**2)*sqrt(e**2+delta**2))
1 continue
c determines the value of u.
u=2.0/s
c outputs results.
write(20,21) u,delta
21 format(2e18.7)
c NOW calculates imaginary part of HF spin-spin correlation function, divided by -pi.
c zeroes out HFimaginary.
do 2 i=1,ne
d1(i)=0.0
d2(i)=0.0
d3(i)=0.0
2 continue
c zeroes out HFreal.
a1=0.0
a2=0.0
a3=0.0
c establishes regulator for gradient components to avoid subtle cases.
eps=0.00001
c begins looping over momenta.
do 5 j=1,nq
qx=(j-0.5)*dq
do 5 k=1,nq
qy=(k-0.5)*dq
c calculates first energy and gradient.
e1=-2.0*t*(cos(qx)+cos(qy))
ee1=sqrt(e1**2+delta**2)
g1x=2.0*t*sin(qx)*e1/ee1
g1y=2.0*t*sin(qy)*e1/ee1
c calculates second energy and gradient.
e2=-2.0*t*(cos(qx+qqx)+cos(qy+qqy))
ee2=sqrt(e2**2+delta**2)
g2x=2.0*t*sin(qx+qqx)*e2/ee2
g2y=2.0*t*sin(qy+qqy)*e2/ee2
c sums energies.
e=ee1+ee2
c sums gradients and moves into upper right quadrant.
x=abs(g1x+g2x)+eps
y=abs(g1y+g2y)+eps
c makes sure x > y.
if(x.gt.y) go to 3
z=x
x=y
y=z
3 continue
c computes energy window.
im=(e-emin-dq/2.0*(x+y))/de-1
ip=(e-emin+dq/2.0*(x+y))/de+1
if(im.lt.1) im=1
if(ip.gt.ne) ip=ne
c loops over energy bins.
do 4 i=im,ip
c computes two renormalized energies.
em=2.0/dq*((i-0.5)*de+emin-e)
ep=2.0/dq*((i+0.5)*de+emin-e)
z=a(em,x,y)-a(ep,x,y)
d1(i)=d1(i)+z*(ee1*ee2-e1*e2+delta**2)/4.0/ee1/ee2
d2(i)=d2(i)+z*(ee1*ee2+e1*e2+delta**2)/4.0/ee1/ee2
d3(i)=d3(i)+z*(ee1+ee2)*delta/4.0/ee1/ee2
4 continue
5 continue
c corrects normalization of integrations.
z=(dq/pi)**2/de
s1=0.0
s2=0.0
s3=0.0
do 6 i=1,ne
d1(i)=d1(i)*z
d2(i)=d2(i)*z
d3(i)=d3(i)*z
6 continue
c NOW hilbert transforms imaginary parts to get real parts, for positive frequencies.
do 9 i=1,ne
e=i*de+emin
do 10 ii=1,ne
ee=ii*de+emin
if(i.eq.ii) go to 10
a1=a1+d1(ii)*(1.0/(e-ee)-1.0/(e+ee))
a2=a2+d2(ii)*(1.0/(e-ee)-1.0/(e+ee))
a3=a3+d3(ii)*(1.0/(e-ee)+1.0/(e+ee))
10 continue
a1=de*(a1-d1(i)/2.0/e)
a2=de*(a2-d2(i)/2.0/e)
a3=de*(a3+d3(i)/2.0/e)
c gets full HF correlation function.
g1(i)=cmplx(a1,-pi*d1(i))
g2(i)=cmplx(a2,-pi*d2(i))
g3(i)=cmplx(a3,-pi*d3(i))
c NOW gets spin-spin correlation function in Ladder approximation, near pipi.
h(i)=(1+u*g1(i))*(1+u*g2(i))-(u*g3(i))**2
h22(i)=(g2(i)+u*g1(i)*g2(i)-(u*(g3(i)**2)))/h(i)
h12(i)=(g3(i))/h(i)
ha2(i)=abs(h22(i))+abs(h12(i))
9 continue
c outputs results for HFfull.
open(24,file='HFfull.dat',status='new')
do 11 i=1,ne
e=i*de+emin
write(24,25) e,g1(i),g2(i),g3(i)
25 format(f10.4,6e18.7)
11 continue
c outputs results for Ladder.
open(26,file='Ladder.dat',status='new')
do 12 i=1,ne
e=i*de+emin
write(26,27) e,ha2(i)
27 format(f10.4,6e18.7)
12 continue
stop
end
c calculates area of little box above energy e.
function a(e,x,y)
z=x*y*8.0
c zeroes out area.
a=0.0
c did we get the whole square?
if(e.gt.-x-y) go to 3
a=1.0
return
3 continue
c no? did we get any of the square?
if(e.gt.x+y) return
c yes? then computes upper triangle area.
a=(e-x-y)**2/z
c did we overshoot the bottom of the box?
if(e.gt.x-y) return
c yes? then clips off lower right triangle.
a=a-(e-x+y)**2/z
c did we overshoot the left side of the box?
if(e.gt.-x+y) return
c yes? then clips off upper left triangle.
a=a-(e+x-y)**2/z
return
end