Wigner/3D_eta1_bcc.f90
2021-03-11 12:28:00 +01:00

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