189 lines
5.0 KiB
Fortran
189 lines
5.0 KiB
Fortran
program triangle
|
|
|
|
implicit none
|
|
|
|
integer :: nside,nsite,nelec
|
|
integer :: ix,iy,kx,ky
|
|
integer :: i,k, ielec,kelec,info
|
|
real*8 :: start,finish
|
|
real*8 :: factor,height,factor2
|
|
real*8 :: zpe
|
|
real*8 :: denom1,denom2,numxy,fx,fy
|
|
real*8 :: work(8)
|
|
real*8,allocatable :: sin2(:),cos2(:)
|
|
real*8,allocatable :: cmat(:,:,:,:),cmatk(:,:,:,:)
|
|
real*8,allocatable :: caux1(:,:,:,:),caux2(:,:,:,:)
|
|
real*8,allocatable :: caux3(:,:,:,:),caux4(:,:,:,:)
|
|
real*8,allocatable :: omega2(:,:)
|
|
|
|
open(unit=10,file='input')
|
|
read(10,*) nside
|
|
close(10)
|
|
|
|
call cpu_time(start)
|
|
|
|
nsite = 2*nside
|
|
nelec = nside**2*2
|
|
|
|
write(6,*) "total number of electrons",nelec
|
|
|
|
allocate (cmat(2,2,nsite,nsite))
|
|
allocate (cmatk(2,2,nsite,nsite))
|
|
allocate (caux1(2,2,nsite,nsite))
|
|
allocate (caux2(2,2,nsite,nsite))
|
|
allocate (caux3(2,2,nsite,nsite))
|
|
allocate (caux4(2,2,nsite,nsite))
|
|
allocate (sin2(nsite))
|
|
allocate (cos2(nsite))
|
|
allocate (omega2(2,nelec))
|
|
|
|
factor=acos(-1d0)/nsite
|
|
height=2d0*sqrt(3d0)
|
|
|
|
do i = 0, nsite-1
|
|
sin2(i+1) = (sin(factor*i))**2
|
|
cos2(i+1) = (cos(factor*i))**2
|
|
enddo
|
|
|
|
!
|
|
! calculate C(0)
|
|
!
|
|
|
|
cmat = 0d0
|
|
do ix = 0, nsite-1
|
|
do iy = 0, nsite-1
|
|
if (ix+iy==0) cycle
|
|
if (mod(ix+iy,2)==1) cycle
|
|
denom1 = (sin2(ix+1)+0.25d0*height**2*sin2(iy+1))**1.5d0
|
|
denom2 = (sin2(ix+1)+0.25d0*height**2*sin2(iy+1))**2.5d0
|
|
numxy = 1.5d0*height*sin(factor*ix)*cos(factor*ix)*sin(factor*iy)*cos(factor*iy)
|
|
cmat(1,1,1,1) = cmat(1,1,1,1) + 3d0*sin2(ix+1)*cos2(ix+1)/denom2 + &
|
|
(sin2(ix+1) - cos2(ix+1))/denom1
|
|
cmat(2,2,1,1) = cmat(2,2,1,1) + 0.75d0*height**2*sin2(iy+1)*cos2(iy+1)/denom2 + &
|
|
(sin2(iy+1) - cos2(iy+1))/denom1
|
|
cmat(1,2,1,1) = cmat(1,2,1,1) + numxy/denom2
|
|
enddo
|
|
enddo
|
|
|
|
!
|
|
! calculate C(i) for i not equal to 0
|
|
!
|
|
|
|
do ix = 0, nsite-1
|
|
do iy = 0, nsite-1
|
|
if (ix+iy==0) cycle
|
|
if (mod(ix+iy,2)==1) cycle
|
|
denom1 = (sin2(ix+1)+0.25d0*height**2*sin2(iy+1))**1.5d0
|
|
denom2 = (sin2(ix+1)+0.25d0*height**2*sin2(iy+1))**2.5d0
|
|
numxy = -1.5d0*height*sin(factor*ix)*cos(factor*ix)*sin(factor*iy)*cos(factor*iy)
|
|
cmat(1,1,ix+1,iy+1) = -3d0*sin2(ix+1)*cos2(ix+1)/denom2 + &
|
|
(cos2(ix+1) - sin2(ix+1))/denom1
|
|
cmat(2,2,ix+1,iy+1) = -0.75d0*height**2*sin2(iy+1)*cos2(iy+1)/denom2 + &
|
|
(cos2(iy+1) - sin2(iy+1))/denom1
|
|
cmat(1,2,ix+1,iy+1) = numxy/denom2
|
|
enddo
|
|
enddo
|
|
cmat(2,1,:,:) = cmat(1,2,:,:)
|
|
|
|
cmat=cmat*(acos(-1d0)*height/nsite**2)**1.5d0
|
|
|
|
print*,'testje',cmat(1,1,1,1),cmat(2,2,1,1)
|
|
|
|
!
|
|
! calculate sum_{i} C_ab(i) cos(2pi k.i/2n) using index decoupling,
|
|
! i.e., cos (x + y) = cos(x)cos(y)-sin(x)sin(y)
|
|
!
|
|
|
|
factor = 2*acos(-1d0)/nsite
|
|
factor2 = 2*acos(-1d0)/(nside*sqrt(3d0))
|
|
|
|
caux1=0d0
|
|
caux2=0d0
|
|
do ix = 0, nsite-1
|
|
fx=factor*ix
|
|
do iy = 0, nsite-1
|
|
if (mod(ix+iy,2)==1) cycle
|
|
do kx = 0, nsite-1
|
|
caux1(:,:,kx+1,iy+1) = caux1(:,:,kx+1,iy+1) + cmat(:,:,ix+1,iy+1)*cos(kx*fx)
|
|
caux2(:,:,kx+1,iy+1) = caux2(:,:,kx+1,iy+1) + cmat(:,:,ix+1,iy+1)*sin(kx*fx)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
caux3=0d0
|
|
caux4=0d0
|
|
do iy = 0, nsite-1
|
|
fy=factor*iy
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
if (mod(kx+ky,2)==1) cycle
|
|
caux3(:,:,kx+1,ky+1) = caux3(:,:,kx+1,ky+1) + caux1(:,:,kx+1,iy+1)*cos(ky*fy)
|
|
caux4(:,:,kx+1,ky+1) = caux4(:,:,kx+1,ky+1) + caux2(:,:,kx+1,iy+1)*sin(ky*fy)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
cmatk = caux3-caux4
|
|
|
|
!!
|
|
!!debug : implementation without index decoupling
|
|
!!
|
|
!cmatk=0d0
|
|
!do ix = 0, nsite-1
|
|
! do iy = 0, nsite-1
|
|
! if (mod(ix+iy,2)==1) cycle
|
|
! do kx = 0, nsite-1
|
|
! do ky = 0, nsite-1
|
|
! if (mod(kx+ky,2)==1) cycle
|
|
!! print*,'testje2',kx,ky,cmat(1,1,ix+1,iy+1),cos(factor*(ix*kx+iy*ky))
|
|
! cmatk(:,:,kx+1,ky+1) = cmatk(:,:,kx+1,ky+1) + cmat(:,:,ix+1,iy+1)*cos(factor*(ix*kx+iy*ky))
|
|
!! cmatk(:,:,kx+1,ky+1) = cmatk(:,:,kx+1,ky+1) + cmat(:,:,ix+1,iy+1)*cos(factor*ix*kx+factor2*iy*(ky-kx))
|
|
!! cmatk(1,1,kx+1,ky+1) = cmatk(1,1,kx+1,ky+1) + cmat(1,1,ix+1,iy+1)*cos(factor*(ix*kx+iy*ky))
|
|
!! print*,'testje3',kx,ky,cmatk(1,1,kx+1,ky+1)
|
|
! enddo
|
|
! enddo
|
|
! enddo
|
|
!enddo
|
|
|
|
!do kx = 0, nsite-1
|
|
! do ky = 0, nsite-1
|
|
! if (mod(kx+ky,2)==1) cycle
|
|
! print*,'testje',kx,ky,cmatk(1,1,kx+1,ky+1),cmatk(2,2,kx+1,ky+1),cmatk(1,2,kx+1,ky+1),cmatk(2,1,kx+1,ky+1)
|
|
! enddo
|
|
!enddo
|
|
|
|
!
|
|
! diagonalize 2x2 matrices
|
|
!
|
|
|
|
omega2 = 0d0
|
|
kelec=1
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
if (mod(kx+ky,2)==1) cycle
|
|
call dsyev ("v","u", 2, cmatk(:,:,kx+1,ky+1), 2, omega2(:,kelec), work, 8, info)
|
|
kelec = kelec + 1
|
|
enddo
|
|
enddo
|
|
|
|
zpe=0d0
|
|
do kelec = 1, nelec
|
|
do i = 1, 2
|
|
! if (omega2(i,kelec)<0d0) print*,'warning',kelec,i,omega2(i,kelec)
|
|
! print*,omega2(i,kelec)
|
|
if (omega2(i,kelec)<0d0) cycle
|
|
zpe = zpe + sqrt(omega2(i,kelec))
|
|
enddo
|
|
enddo
|
|
zpe=zpe/nelec
|
|
|
|
print*,'zero-point energy (in Ry and Ha) =',zpe,zpe/2
|
|
|
|
deallocate(sin2,cos2,cmat,cmatk,omega2)
|
|
deallocate(caux1,caux2,caux3,caux4)
|
|
|
|
call cpu_time(finish)
|
|
print '("Time = ",f6.3," seconds.")',finish-start
|
|
|
|
end program triangle
|