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