332 lines
9.4 KiB
Fortran
332 lines
9.4 KiB
Fortran
program bcc
|
|
|
|
implicit none
|
|
|
|
integer :: nside,nsite,nelec
|
|
integer :: ix,iy,iz,kx,ky,kz
|
|
integer :: i,k, ielec,kelec,info
|
|
real*8 :: start,finish
|
|
real*8 :: factor
|
|
real*8 :: zpe
|
|
real*8 :: denom1,denom2,numxy,numxz,numyz,fx,fy,fz
|
|
real*8 :: work(8)
|
|
real*8,allocatable :: sin2(:),cos2(:)
|
|
real*8,allocatable :: cmat(:,:,:,:,:)
|
|
real*8,allocatable :: caux(:,:,:,:,:)
|
|
real*8,allocatable :: cmatk(:,:,:,:,:)
|
|
real*8,allocatable :: omega2(:,:)
|
|
|
|
open(unit=10,file='input')
|
|
read(10,*) nside
|
|
close(10)
|
|
|
|
call cpu_time(start)
|
|
|
|
nsite = 2*nside
|
|
nelec = nside**3*2
|
|
|
|
write(6,*) "total number of electrons",nelec
|
|
|
|
allocate (cmat(3,3,nsite,nsite,nsite))
|
|
allocate (caux(3,3,nsite,nsite,nsite))
|
|
allocate (cmatk(3,3,nsite,nsite,nsite))
|
|
allocate (sin2(nsite))
|
|
allocate (cos2(nsite))
|
|
allocate (omega2(3,nelec))
|
|
|
|
factor=acos(-1d0)/nsite
|
|
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
|
|
do iz = 0, nsite-1
|
|
if (ix+iy+iz==0) cycle
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
denom1 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**(1.5d0)
|
|
denom2 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**(2.5d0)
|
|
cmat(1,1,1,1,1) = cmat(1,1,1,1,1) + 3d0*sin2(ix+1)*cos2(ix+1)/denom2 + (sin2(ix+1) - cos2(ix+1))/denom1
|
|
cmat(2,2,1,1,1) = cmat(2,2,1,1,1) + 3d0*sin2(iy+1)*cos2(iy+1)/denom2 + (sin2(iy+1) - cos2(iy+1))/denom1
|
|
cmat(3,3,1,1,1) = cmat(3,3,1,1,1) + 3d0*sin2(iz+1)*cos2(iz+1)/denom2 + (sin2(iz+1) - cos2(iz+1))/denom1
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
!
|
|
! calculate C(i) for i not equal to 0
|
|
!
|
|
|
|
do ix = 0, nsite-1
|
|
do iy = 0, nsite-1
|
|
do iz = 0, nsite-1
|
|
if (ix+iy+iz==0) cycle
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
denom1 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**1.5d0
|
|
denom2 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**2.5d0
|
|
numxy = -3d0*sin(factor*ix)*cos(factor*ix)*sin(factor*iy)*cos(factor*iy)
|
|
numxz = -3d0*sin(factor*ix)*cos(factor*ix)*sin(factor*iz)*cos(factor*iz)
|
|
numyz = -3d0*sin(factor*iy)*cos(factor*iy)*sin(factor*iz)*cos(factor*iz)
|
|
cmat(1,1,ix+1,iy+1,iz+1) = -3d0*sin2(ix+1)*cos2(ix+1)/denom2 + (cos2(ix+1) - sin2(ix+1))/denom1
|
|
cmat(2,2,ix+1,iy+1,iz+1) = -3d0*sin2(iy+1)*cos2(iy+1)/denom2 + (cos2(iy+1) - sin2(iy+1))/denom1
|
|
cmat(3,3,ix+1,iy+1,iz+1) = -3d0*sin2(iz+1)*cos2(iz+1)/denom2 + (cos2(iz+1) - sin2(iz+1))/denom1
|
|
cmat(1,2,ix+1,iy+1,iz+1) = numxy/denom2
|
|
cmat(1,3,ix+1,iy+1,iz+1) = numxz/denom2
|
|
cmat(2,3,ix+1,iy+1,iz+1) = numyz/denom2
|
|
enddo
|
|
enddo
|
|
enddo
|
|
cmat(2,1,:,:,:) = cmat(1,2,:,:,:)
|
|
cmat(3,1,:,:,:) = cmat(1,3,:,:,:)
|
|
cmat(3,2,:,:,:) = cmat(2,3,:,:,:)
|
|
|
|
factor = 3d0*(acos(-1d0))**2/nsite**3
|
|
do ix = 0, nsite-1
|
|
do iy = 0, nsite-1
|
|
do iz = 0, nsite-1
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
cmat(:,:,ix+1,iy+1,iz+1) = cmat(:,:,ix+1,iy+1,iz+1)*factor
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
!
|
|
! 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
|
|
|
|
caux=0d0
|
|
do ix = 0, nsite-1
|
|
fx=factor*ix
|
|
do iy = 0, nsite-1
|
|
do iz = 0, nsite-1
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
do kx = 0, nsite-1
|
|
caux(:,:,kx+1,iy+1,iz+1) = caux(:,:,kx+1,iy+1,iz+1) + cmat(:,:,ix+1,iy+1,iz+1)*cos(kx*fx)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
cmat=0d0
|
|
do iy = 0, nsite-1
|
|
fy=factor*iy
|
|
do iz = 0, nsite-1
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
if (mod(iy+iz,2)==1.or.mod(kx+ky,2)==1) cycle
|
|
cmat(:,:,kx+1,ky+1,iz+1) = cmat(:,:,kx+1,ky+1,iz+1) + caux(:,:,kx+1,iy+1,iz+1)*cos(ky*fy)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
cmatk=0d0
|
|
do iz = 0, nsite-1
|
|
fz=factor*iz
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
do kz = 0, nsite-1
|
|
if (mod(kx+ky,2)==1.or.mod(kx+kz,2)==1.or.mod(ky+kz,2)==1) cycle
|
|
cmatk(:,:,kx+1,ky+1,kz+1) = cmatk(:,:,kx+1,ky+1,kz+1) + cmat(:,:,kx+1,ky+1,iz+1)*cos(kz*fz)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
cmat=0d0
|
|
do iy = 0, nsite-1
|
|
fy=factor*iy
|
|
do iz = 0, nsite-1
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
if (mod(iy+iz,2)==1.or.mod(kx+ky,2)==1) cycle
|
|
cmat(:,:,kx+1,ky+1,iz+1) = cmat(:,:,kx+1,ky+1,iz+1) + caux(:,:,kx+1,iy+1,iz+1)*sin(ky*fy)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
do iz = 0, nsite-1
|
|
fz=factor*iz
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
do kz = 0, nsite-1
|
|
if (mod(kx+ky,2)==1.or.mod(kx+kz,2)==1.or.mod(ky+kz,2)==1) cycle
|
|
cmatk(:,:,kx+1,ky+1,kz+1) = cmatk(:,:,kx+1,ky+1,kz+1) - cmat(:,:,kx+1,ky+1,iz+1)*sin(kz*fz)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
factor=acos(-1d0)/nsite
|
|
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
|
|
do iz = 0, nsite-1
|
|
if (ix+iy+iz==0) cycle
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
denom1 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**(1.5d0)
|
|
denom2 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**(2.5d0)
|
|
cmat(1,1,1,1,1) = cmat(1,1,1,1,1) + 3d0*sin2(ix+1)*cos2(ix+1)/denom2 + (sin2(ix+1) - cos2(ix+1))/denom1
|
|
cmat(2,2,1,1,1) = cmat(2,2,1,1,1) + 3d0*sin2(iy+1)*cos2(iy+1)/denom2 + (sin2(iy+1) - cos2(iy+1))/denom1
|
|
cmat(3,3,1,1,1) = cmat(3,3,1,1,1) + 3d0*sin2(iz+1)*cos2(iz+1)/denom2 + (sin2(iz+1) - cos2(iz+1))/denom1
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
!
|
|
! calculate C(i) for i not equal to 0
|
|
!
|
|
|
|
do ix = 0, nsite-1
|
|
do iy = 0, nsite-1
|
|
do iz = 0, nsite-1
|
|
if (ix+iy+iz==0) cycle
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
denom1 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**1.5d0
|
|
denom2 = (sin2(ix+1)+sin2(iy+1)+sin2(iz+1))**2.5d0
|
|
numxy = -3d0*sin(factor*ix)*cos(factor*ix)*sin(factor*iy)*cos(factor*iy)
|
|
numxz = -3d0*sin(factor*ix)*cos(factor*ix)*sin(factor*iz)*cos(factor*iz)
|
|
numyz = -3d0*sin(factor*iy)*cos(factor*iy)*sin(factor*iz)*cos(factor*iz)
|
|
cmat(1,1,ix+1,iy+1,iz+1) = -3d0*sin2(ix+1)*cos2(ix+1)/denom2 + (cos2(ix+1) - sin2(ix+1))/denom1
|
|
cmat(2,2,ix+1,iy+1,iz+1) = -3d0*sin2(iy+1)*cos2(iy+1)/denom2 + (cos2(iy+1) - sin2(iy+1))/denom1
|
|
cmat(3,3,ix+1,iy+1,iz+1) = -3d0*sin2(iz+1)*cos2(iz+1)/denom2 + (cos2(iz+1) - sin2(iz+1))/denom1
|
|
cmat(1,2,ix+1,iy+1,iz+1) = numxy/denom2
|
|
cmat(1,3,ix+1,iy+1,iz+1) = numxz/denom2
|
|
cmat(2,3,ix+1,iy+1,iz+1) = numyz/denom2
|
|
enddo
|
|
enddo
|
|
enddo
|
|
cmat(2,1,:,:,:) = cmat(1,2,:,:,:)
|
|
cmat(3,1,:,:,:) = cmat(1,3,:,:,:)
|
|
cmat(3,2,:,:,:) = cmat(2,3,:,:,:)
|
|
|
|
factor = 3d0*(acos(-1d0))**2/nsite**3
|
|
do ix = 0, nsite-1
|
|
do iy = 0, nsite-1
|
|
do iz = 0, nsite-1
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
cmat(:,:,ix+1,iy+1,iz+1) = cmat(:,:,ix+1,iy+1,iz+1)*factor
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
factor = 2*acos(-1d0)/nsite
|
|
|
|
caux=0d0
|
|
do ix = 0, nsite-1
|
|
fx=factor*ix
|
|
do iy = 0, nsite-1
|
|
do iz = 0, nsite-1
|
|
if (mod(ix+iy,2)==1.or.mod(ix+iz,2)==1.or.mod(iy+iz,2)==1) cycle
|
|
do kx = 0, nsite-1
|
|
caux(:,:,kx+1,iy+1,iz+1) = caux(:,:,kx+1,iy+1,iz+1) + cmat(:,:,ix+1,iy+1,iz+1)*sin(kx*fx)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
cmat=0d0
|
|
do iy = 0, nsite-1
|
|
fy=factor*iy
|
|
do iz = 0, nsite-1
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
if (mod(iy+iz,2)==1.or.mod(kx+ky,2)==1) cycle
|
|
cmat(:,:,kx+1,ky+1,iz+1) = cmat(:,:,kx+1,ky+1,iz+1) + caux(:,:,kx+1,iy+1,iz+1)*sin(ky*fy)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
do iz = 0, nsite-1
|
|
fz=factor*iz
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
do kz = 0, nsite-1
|
|
if (mod(kx+ky,2)==1.or.mod(kx+kz,2)==1.or.mod(ky+kz,2)==1) cycle
|
|
cmatk(:,:,kx+1,ky+1,kz+1) = cmatk(:,:,kx+1,ky+1,kz+1) - cmat(:,:,kx+1,ky+1,iz+1)*cos(kz*fz)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
cmat=0d0
|
|
do iy = 0, nsite-1
|
|
fy=factor*iy
|
|
do iz = 0, nsite-1
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
if (mod(iy+iz,2)==1.or.mod(kx+ky,2)==1) cycle
|
|
cmat(:,:,kx+1,ky+1,iz+1) = cmat(:,:,kx+1,ky+1,iz+1) + caux(:,:,kx+1,iy+1,iz+1)*cos(ky*fy)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
do iz = 0, nsite-1
|
|
fz=factor*iz
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
do kz = 0, nsite-1
|
|
if (mod(kx+ky,2)==1.or.mod(kx+kz,2)==1.or.mod(ky+kz,2)==1) cycle
|
|
cmatk(:,:,kx+1,ky+1,kz+1) = cmatk(:,:,kx+1,ky+1,kz+1) - cmat(:,:,kx+1,ky+1,iz+1)*sin(kz*fz)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
!
|
|
! diagonalize 3x3 matrices
|
|
!
|
|
|
|
omega2 = 0d0
|
|
kelec=1
|
|
do kx = 0, nsite-1
|
|
do ky = 0, nsite-1
|
|
do kz = 0, nsite-1
|
|
if (mod(kx+ky,2)==1.or.mod(kx+kz,2)==1.or.mod(ky+kz,2)==1) cycle
|
|
call dsyev ("v","u", 3, cmatk(:,:,kx+1,ky+1,kz+1), 3, omega2(:,kelec), work, 8, info)
|
|
kelec = kelec + 1
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
zpe=0d0
|
|
do kelec = 1, nelec
|
|
do i = 1, 3
|
|
! 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,omega2)
|
|
deallocate(cmat,caux,cmatk)
|
|
|
|
call cpu_time(finish)
|
|
print '("Time = ",f6.3," seconds.")',finish-start
|
|
|
|
end program bcc
|