program square

implicit none

integer :: nside,nsite,nelec
integer :: ix,iy,iz,i
real*16  :: factor,bbfac
real*16  :: uee,ubb,eta0,aux
real*16  :: start,finish
real*16,allocatable  :: sin2(:)

open(unit=10,file='input')
   read(10,*) nside
close(10)

call cpu_time(start)

nsite = nside
nelec = nside**2

bbfac = 1.139479116668307748953273191122493968q0

write(6,*) "total number of electrons",nelec

allocate(sin2(nsite))

factor=acos(-1q0)/nsite
do i = 0, nsite-1
   sin2(i+1) = (sin(factor*i))**2
enddo

uee = 0q0
do ix = 0, nsite-1
   do iy = 0, nsite-1
      if (ix+iy==0) cycle
      aux = sqrt(sin2(ix+1)+sin2(iy+1))
      uee = uee + 1q0/aux
   enddo
enddo

factor=sqrt(acos(-1q0))/(nsite*2q0)
uee=uee*factor

ubb = nelec*bbfac/nside

eta0 = uee - ubb

print*,'eta0 (in Ry and Ha) =', 2q0*eta0, eta0
print*, uee,ubb

deallocate(sin2)

call cpu_time(finish)
print '("Time = ",f12.3," seconds.")',finish-start

end program square