Why does the ZHEEV function give different eigenvalues with the same matrix?

Question

I wrote a Fortran code below to use the ZHEEV function to diagonalize the matrix and obtain the eigenvalue.

SUBROUTINE DOTPRODUCT(v1,v2,v3)

IMPLICIT NONE

INTEGER, PARAMETER                :: dp = SELECTED_REAL_KIND(15,14)
REAL (KIND=dp)                     :: v1(3)
REAL (KIND=dp)                     :: v2(3)
REAL (KIND=dp)                     :: v3(3)

v3(1) = v1(1)*v2(1)
v3(2) = v1(2)*v2(2)
v3(3) = v1(3)*v2(3)

RETURN
END SUBROUTINE DOTPRODUCT

SUBROUTINE CROSSPRODUCT(v1,v2,v3)

IMPLICIT NONE

INTEGER, PARAMETER                 :: dp = SELECTED_REAL_KIND(15,14)
REAL (KIND=dp)                     :: v1(3)
REAL (KIND=dp)                     :: v2(3)
REAL (KIND=dp)                     :: v3(3)

v3(1) = v1(2)*v2(3)-v1(3)*v2(2)
v3(2) = v1(3)*v2(1)-v1(1)*v2(3)
v3(3) = v1(1)*v2(2)-v1(2)*v2(1)
RETURN
END SUBROUTINE CROSSPRODUCT

SUBROUTINE LATTICEINK(lr,lk)

IMPLICIT NONE

INTEGER, PARAMETER                 :: dp = SELECTED_REAL_KIND(15,14)
REAL (KIND=dp)                     :: lr(3,2) !matrix from the main code
REAL (KIND=dp)                     :: lk(3,2) !matrix back to the main code
REAL (KIND=dp)                     :: em(3,3) !meta results
REAL (KIND=dp)                     :: cr(3,3) !meta results
REAL (KIND=dp)                     :: vl(3)   !Volume of the unit cell in the real space
REAL (KIND=dp)                     :: pi      !constant pi

em = 0.0d0
em(:,1:2) = lr
em(3,3) = 1.5d1
pi = DACOS(-1.0d0)
CALL CROSSPRODUCT(em(:,2),em(:,3),cr(:,1))
CALL CROSSPRODUCT(em(:,3),em(:,1),cr(:,2))
CALL CROSSPRODUCT(em(:,1),em(:,2),cr(:,3))
CALL DOTPRODUCT(em(:,1),cr(:,1),vl)
vl(1) = DSQRT(vl(1)**2+vl(2)**2+vl(3)**2)
CALL DOTPRODUCT(em(:,2),cr(:,2),vl)
vl(1) = DSQRT(vl(1)**2+vl(2)**2+vl(3)**2)
CALL DOTPRODUCT(em(:,3),cr(:,3),vl)
vl(1) = DSQRT(vl(1)**2+vl(2)**2+vl(3)**2)
lk(:,1) = 2.0d0*pi*cr(:,1)/vl(1)
lk(:,2) = 2.0d0*pi*cr(:,2)/vl(1)
RETURN
END SUBROUTINE LATTICEINK

PROGRAM GRAPHENE_TIGHTBINDING

USE LAPACK
IMPLICIT NONE

INTEGER, PARAMETER                 :: dp = SELECTED_REAL_KIND(15,14)
INTEGER                            :: i, j
REAL (KIND=dp)                     :: dc      !Distant between carbon atom in pristine graphene unit cell
REAL (KIND=dp)                     :: lr(3,2) !Lattice constant in real space
REAL (KIND=dp)                     :: lk(3,2) !Lattice constant in k space
REAL (KIND=dp)                     :: cr(3,3) !matrix to store the meta data
REAL (KIND=dp)                     :: pi
REAL (KIND=dp)                     :: kh(3,3) !High symmetry k point
INTEGER                            :: km(3)   !k point sampling
REAL (KIND=dp), ALLOCATABLE        :: kc(:,:) !k point coordinate
REAL (KIND=dp)                     :: mk(3,2)
REAL (KIND=dp)                     :: ep(2)   !onsite
REAL (KIND=dp)                     :: hp      !hopping
COMPLEX (KIND=dp)                  :: ha(2,2) !Hamiltonian
REAL (KIND=dp), ALLOCATABLE        :: bd(:,:) !band
INTEGER                            :: nu, lda, lwork, info
REAL (KIND=dp), ALLOCATABLE        :: ei(:), rwork(:)
COMPLEX (KIND=dp), ALLOCATABLE     :: work(:)

!Parameters
dc = 1.42d0
pi = DACOS(-1.0d0)
lr = 0.0d0
lr(1,1) = DCOS(pi/6.0d0)*dc*2.0d0
lr(1,2) = DCOS(pi*5.0d0/6.0d0)*dc
lr(2,2) = DSIN(pi/6.0d0)*dc+dc
kh = 0.0d0
kh(1,2) = 1.0d0/3.0d0
kh(2,2) = 1.0d0/3.0d0
kh(1,3) = 5.0d-1
km(1) = 100
km(2) = 100
km(3) = 100
ep = 0.0d0
hp = 2.8d0
nu = 2
lda = nu
lwork = -1
!
OPEN (UNIT=3, FILE='bandstructure.dat', STATUS='UNKNOWN')
!Computing the lattice constant in k space
CALL LATTICEINK(lr,lk)
!
!Computing the k point coordinate
ALLOCATE (kc(4,km(1)+km(2)+km(3)))
kc = 0.0d0
mk(1,1) = kh(1,1)*lk(1,1)+kh(2,1)*lk(1,2)
mk(2,1) = kh(1,1)*lk(2,1)+kh(2,1)*lk(2,2)
mk(3,1) = kh(1,1)*lk(3,1)+kh(2,1)*lk(3,2)
mk(1,2) = kh(1,2)*lk(1,1)+kh(2,2)*lk(1,2)
mk(2,2) = kh(1,2)*lk(2,1)+kh(2,2)*lk(2,2)
mk(3,2) = kh(1,2)*lk(3,1)+kh(2,2)*lk(3,2)
DO i = 1, km(1), 1
   kc(1,i) = mk(1,1)+DBLE(i-1)*(mk(1,2)-mk(1,1))/DBLE(km(1))
   kc(2,i) = mk(2,1)+DBLE(i-1)*(mk(2,2)-mk(2,1))/DBLE(km(1))
   kc(3,i) = mk(3,1)+DBLE(i-1)*(mk(3,2)-mk(3,1))/DBLE(km(1))
END DO
mk(1,1) = kh(1,2)*lk(1,1)+kh(2,2)*lk(1,2)
mk(2,1) = kh(1,2)*lk(2,1)+kh(2,2)*lk(2,2)
mk(3,1) = kh(1,2)*lk(3,1)+kh(2,2)*lk(3,2)
mk(1,2) = kh(1,3)*lk(1,1)+kh(2,3)*lk(1,2)
mk(2,2) = kh(1,3)*lk(2,1)+kh(2,3)*lk(2,2)
mk(3,2) = kh(1,3)*lk(3,1)+kh(2,3)*lk(3,2)
DO i = km(1)+1, km(1)+km(2), 1
   kc(1,i) = mk(1,1)+DBLE(i-km(1)-1)*(mk(1,2)-mk(1,1))/DBLE(km(2))
   kc(2,i) = mk(2,1)+DBLE(i-km(1)-1)*(mk(2,2)-mk(2,1))/DBLE(km(2))
   kc(3,i) = mk(3,1)+DBLE(i-km(1)-1)*(mk(3,2)-mk(3,1))/DBLE(km(2))
END DO
mk(1,1) = kh(1,3)*lk(1,1)+kh(2,3)*lk(1,2)
mk(2,1) = kh(1,3)*lk(2,1)+kh(2,3)*lk(2,2)
mk(3,1) = kh(1,3)*lk(3,1)+kh(2,3)*lk(3,2)
mk(1,2) = kh(1,1)*lk(1,1)+kh(2,1)*lk(1,2)
mk(2,2) = kh(1,1)*lk(2,1)+kh(2,1)*lk(2,2)
mk(3,2) = kh(1,1)*lk(3,1)+kh(2,1)*lk(3,2)
DO i = km(1)+km(2)+1, km(1)+km(2)+km(3), 1
   kc(1,i) = mk(1,1)+DBLE(i-km(1)-km(2)-1)*(mk(1,2)-mk(1,1))/DBLE(km(3))
   kc(2,i) = mk(2,1)+DBLE(i-km(1)-km(2)-1)*(mk(2,2)-mk(2,1))/DBLE(km(3))
   kc(3,i) = mk(3,1)+DBLE(i-km(1)-km(2)-1)*(mk(3,2)-mk(3,1))/DBLE(km(3))
END DO
!
!Building Hamiltonian and computing the band structure
ALLOCATE (bd(nu,km(1)+km(2)+km(3)))
ALLOCATE (ei(nu))
ALLOCATE (work(lwork))
ALLOCATE (rwork(3*nu-2))
ha = DCMPLX(0.0d0, 0.0d0)
ha(1,1) = DCMPLX(ep(1), 0.0d0)
ha(1,2) = hp*(1.0d0+EXP(DCMPLX(0.0d0,-1.0d0)*(kc(1,km(1)+1)*lr(1,1)&
          +kc(2,km(1)+1)*lr(2,1)+kc(3,km(1)+1)*lr(3,1)))+EXP(DCMPLX(0.0d0,-1.0d0)&
          *(kc(1,km(1)+1)*lr(1,2)+kc(2,km(1)+1)*lr(2,2)+kc(3,km(1)+1)*lr(3,2))))
ha(2,2) = DCMPLX(ep(2), 0.0d0)
ha(2,1) = DCONJG(ha(1,2))
print*,'outsideloop'
print*,'ha(1,1)',ha(1,1)
print*,'ha(1,2)',ha(1,2)
print*,'ha(2,1)',ha(2,1)
print*,'ha(2,2)',ha(2,2)
print*,'nu',nu
print*,'lda',lda
print*,'ei',ei
print*,'work',work
print*,'lwork',lwork
print*,'rwork',rwork
print*,'info',info
CALL ZHEEV('V','U',nu,ha,lda,ei,work,lwork,rwork,info)
print*,'ei,kc',ei,kc(:,km(1)+1)
lwork = INT(REAL(work(1)))
DEALLOCATE (work)
ALLOCATE (work(lwork))
DO i = 1, km(1)+km(2)+km(3), 1
   rwork = 0.0d0
   ei = 0.0d0
   ha = DCMPLX(0.0d0, 0.0d0)
   ha(1,1) = DCMPLX(ep(1), 0.0d0)
   ha(1,2) = hp*(1.0d0+EXP(DCMPLX(0.0d0,-1.0d0)*(kc(1,i)*lr(1,1)&
             +kc(2,i)*lr(2,1)+kc(3,i)*lr(3,1)))+EXP(DCMPLX(0.0d0,-1.0d0)&
             *(kc(1,i)*lr(1,2)+kc(2,i)*lr(2,2)+kc(3,i)*lr(3,2))))
   ha(2,2) = DCMPLX(ep(2), 0.0d0)
   ha(2,1) = DCONJG(ha(1,2))
if (i == km(1)+1) then
print*,'within the loop at K'
print*,'ha(1,1)',ha(1,1)
print*,'ha(1,2)',ha(1,2)
print*,'ha(2,1)',ha(2,1)
print*,'ha(2,2)',ha(2,2)
print*,'nu',nu
print*,'lda',lda
print*,'ei',ei
print*,'work',work
print*,'lwork',lwork
print*,'rwork',rwork
print*,'info',info
end if
   CALL ZHEEV('V','U',nu,ha,lda,ei,work,lwork,rwork,info)
if (i == km(1)+1) then
print*,'within the loop at K',ei,kc(:,i)
end if
   bd(:,i) = ei
END DO
!
!Plotting band structure
kc(4,1) = 0.0d0
DO i = 2, km(1)+km(2)+km(3), 1
   kc(4,i) = kc(4,i-1)+DSQRT((kc(1,i)-kc(1,i-1))**2+(kc(2,i)-kc(2,i-1))**2+(kc(3,i)-kc(3,i-1))**2)
END DO
DO i = 1, nu, 1
   DO j = 1, km(1)+km(2)+km(3), 1
      WRITE (UNIT=3, FMT='(F15.10,3X,F15.10)') kc(4,j), bd(i,j)
   END DO
   WRITE (UNIT=3, FMT=*)
END DO
!
DEALLOCATE (bd)
DEALLOCATE (ei)
DEALLOCATE (kc)
DEALLOCATE (rwork)
DEALLOCATE (work)
CLOSE (UNIT=3)
END PROGRAM GRAPHENE_TIGHTBINDING

The 'ei' and 'ha' variables store the eigenvalue and matrix respectively. I printed the eigenvalues obtained by the ZHEEV function with the same matrix both outside and inside the loop; however, the output eigenvalues are not same.

outsideloop
ha(1,1)               (0.0000000000000000,0.0000000000000000)
ha(1,2)              (0.0000000000000000,-4.8497422611928558)
ha(2,1)               (0.0000000000000000,4.8497422611928558)
ha(2,2)               (0.0000000000000000,0.0000000000000000)
nu           2
lda           2
ei   0.0000000000000000        0.0000000000000000     
work
lwork          -1
rwork   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
info           0
ei,kc   0.0000000000000000        0.0000000000000000       0.85154899729305999        1.4749261284459125        0.0000000000000000        0.0000000000000000     
within the loop at K
ha(1,1)               (0.0000000000000000,0.0000000000000000)
ha(1,2)              (0.0000000000000000,-4.8497422611928558)
ha(2,1)               (0.0000000000000000,4.8497422611928558)
ha(2,2)               (0.0000000000000000,0.0000000000000000)
nu           2
lda           2
ei   0.0000000000000000        0.0000000000000000     
work               (66.000000000000000,0.0000000000000000)               (0.0000000000000000,0.0000000000000000)               (32.000000000000000,0.0000000000000000)               (0.0000000000000000,0.0000000000000000)               (0.0000000000000000,0.0000000000000000)               (0.0000000000000000,0.0000000000000000)              
......
lwork          66
rwork   0.0000000000000000        0.0000000000000000        0.0000000000000000        0.0000000000000000     
info           0
within the loop at K  -4.8497422611928558        4.8497422611928558       0.85154899729305999        1.4749261284459125        0.0000000000000000        0.0000000000000000

It can be seen that the eigenvalues are not same with the same matrix diagonalized by the ZHEEV function. Before the loop the eigenvalues are '0.0000000000000000 0.0000000000000000'; while, inside the loop, the eigenvalues are '-4.8497422611928558 4.8497422611928558'. The matrix is same as follows:

ha(1,1)               (0.0000000000000000,0.0000000000000000)
ha(1,2)              (0.0000000000000000,-4.8497422611928558)
ha(2,1)               (0.0000000000000000,4.8497422611928558)
ha(2,2)               (0.0000000000000000,0.0000000000000000)

Would anyone please tell me what is wrong with my code and how to solve it?

Thank you in advance.

Answer 1

Ian is right in his second comment, the first call is just enquiery, it allows you to get the right sizes of your working arrays. This call does not do any calculation. You can set the values of the output to some specific values and see if they are retained.

Only the subsequent calls to zheev do any calculation and hence you see some specific values returned in the output arrays.

Since the values of ei are most likely undefined at the point you are printing them, you should be able to diagnose this problem with valgrind or with sanitizations (-fsanitize=undefined in GCC).