80 lines
1.8 KiB
Fortran
80 lines
1.8 KiB
Fortran
program hcp
|
|
|
|
implicit none
|
|
|
|
integer(kind=8) :: nside,nsitexz,nsitey,nelec
|
|
integer(kind=8) :: ix,iy,iz,i
|
|
real*16 :: factorxz,factory,factor,bbfac,pi
|
|
real*16 :: uee,ubb,eta0,aux
|
|
real*16 :: start,finish
|
|
real*16,allocatable :: sin2xz(:),sin2y(:)
|
|
|
|
open(unit=10,file='input')
|
|
read(10,*) nside
|
|
close(10)
|
|
|
|
call cpu_time(start)
|
|
|
|
nsitexz = 2*nside
|
|
nsitey = 6*nside
|
|
nelec = nside**3*4
|
|
|
|
bbfac = 0.54063377545280790285q0
|
|
pi = acos(-1q0)
|
|
|
|
write(6,*) "total number of electrons",nelec
|
|
|
|
allocate(sin2xz(nsitexz))
|
|
allocate(sin2y(nsitey))
|
|
|
|
factorxz=acos(-1q0)/nsitexz
|
|
factory=acos(-1q0)/nsitey
|
|
do i = 0, nsitexz-1
|
|
sin2xz(i+1) = (sin(factorxz*i))**2
|
|
enddo
|
|
do i = 0, nsitey-1
|
|
sin2y(i+1) = (sin(factory*i))**2
|
|
enddo
|
|
|
|
uee = 0q0
|
|
do ix = 0, nsitexz-1
|
|
do iy = 0, nsitey-1
|
|
do iz = 0, nsitexz-1
|
|
if (ix+iy+iz==0) cycle
|
|
aux = sqrt(sin2xz(ix+1)+3q0*sin2y(iy+1)+8q0/3q0*sin2xz(iz+1))
|
|
if (mod(ix,2)==0.and.mod(iz,2)==0.and.mod(iy,6)==0) then
|
|
uee = uee + 1q0/aux
|
|
! print*,'testje1',ix,iy,iz,aux
|
|
endif
|
|
if (mod(ix,2)==1.and.mod(iz,2)==0.and.mod(iy+3,6)==0) then
|
|
uee = uee + 1q0/aux
|
|
! print*,'testje2',ix,iy,iz,aux
|
|
endif
|
|
if (mod(ix,2)==0.and.mod(iz,2)==1.and.mod(iy+4,6)==0) then
|
|
uee = uee + 1q0/aux
|
|
! print*,'testje3',ix,iy,iz,aux
|
|
endif
|
|
if (mod(ix,2)==1.and.mod(iz,2)==1.and.mod(iy+1,6)==0) then
|
|
uee = uee + 1q0/aux
|
|
! print*,'testje4',ix,iy,iz,aux
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
factor=(3q0*pi**(2q0)*sqrt(2q0))**(1q0/3q0)/(2q0*nsitexz)
|
|
uee=uee*factor
|
|
|
|
ubb = nelec*bbfac/nside
|
|
|
|
eta0 = uee - ubb
|
|
|
|
print*,'eta0 (in Ry and Ha) =', 2q0*eta0, eta0
|
|
|
|
deallocate(sin2xz,sin2y)
|
|
|
|
call cpu_time(finish)
|
|
print '("Time = ",f12.3," seconds.")',finish-start
|
|
|
|
end program hcp
|