diff --git a/src/dirac.f90 b/src/dirac.f90
index 0fc99c0..957aeb7 100644
--- a/src/dirac.f90
+++ b/src/dirac.f90
@@ -234,7 +234,7 @@ subroutine Ffunc(x, y, eng)
     real(dp) :: E_dirac_shift
     integer :: idx
     logical :: accurate_eigensolver
-    accurate_eigensolver = .true.
+    accurate_eigensolver = .false.
     iter = iter + 1
     print *, "SCF iteration:", iter
     Vin = reshape(x, shape(Vin))
diff --git a/src/lapack.f90 b/src/lapack.f90
index 75fee7d..e2f64cc 100644
--- a/src/lapack.f90
+++ b/src/lapack.f90
@@ -77,6 +77,18 @@ SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, &
     REAL(dp)           A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
     END SUBROUTINE
 
+    SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, &
+                       VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, &
+                       LWORK, IWORK, IFAIL, INFO )
+    import :: dp
+    CHARACTER          JOBZ, RANGE, UPLO
+    INTEGER            IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
+    REAL(dp)           ABSTOL, VL, VU
+    INTEGER            IFAIL( * ), IWORK( * )
+    REAL(dp)           A( LDA, * ), B( LDB, * ), W( * ), WORK( * ), &
+                       Z( LDZ, * )
+    END SUBROUTINE
+
     REAL(dp) FUNCTION DLAMCH( CMACH )
     import :: dp
     CHARACTER          CMACH
diff --git a/src/linalg.f90 b/src/linalg.f90
index a12ee3b..6ffa550 100644
--- a/src/linalg.f90
+++ b/src/linalg.f90
@@ -1,6 +1,6 @@
 module linalg
   use types, only: dp
-  use lapack, only: dsyevd, dsygvd, dgesv
+  use lapack, only: dsyevd, dsygvd, dgesv, dsygvx
   implicit none
   private
   public eigh, solve
@@ -24,8 +24,10 @@ subroutine deigh_generalized(Am, Bm, lam, c)
     integer :: n
     ! lapack variables
     integer :: lwork, liwork, info
-    integer, allocatable :: iwork(:)
-    real(dp), allocatable :: Bmt(:,:), work(:)
+    integer, allocatable :: iwork(:), ifail(:)
+    real(dp), allocatable :: Bmt(:,:), work(:), Z(:,:), Amt(:,:)
+    integer :: il, iu, M
+    real(dp) :: abstol
 
     ! solve
     n = size(Am,1)
@@ -34,9 +36,22 @@ subroutine deigh_generalized(Am, Bm, lam, c)
     call assert_shape(c, [n, n], "eigh", "c")
     lwork = 1 + 6*n + 2*n**2
     liwork = 3 + 5*n
-    allocate(Bmt(n,n), work(lwork), iwork(liwork))
-    c = Am; Bmt = Bm  ! Bmt temporaries overwritten by dsygvd
-    call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info)
+    allocate(work(lwork), iwork(liwork))
+    allocate(ifail(n))
+    !call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info)
+    il = 1
+    iu = 7
+    M = iu-il+1
+    allocate(z(n,M))
+    abstol = 1e-4_dp
+    call dsygvx(1,'V','I','L',n,Am,n,Bm,n, &
+        0._dp, 0._dp, 1, 7, abstol, M, lam, c, n, work, &
+        lwork, iwork, ifail, info)
+    !SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, &
+    !                   LWORK, IWORK, LIWORK, INFO )
+    !SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, &
+    !                   VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, &
+    !                   LWORK, IWORK, IFAIL, INFO )
     if (info /= 0) then
        print *, "dsygvd returned info =", info
        if (info < 0) then