Improve estimate of the gap by using power method

william-dawson · Sep 8, 2023 · 6cdef57 · 6cdef57
1 parent 612fd8c
commit 6cdef57
Show file tree

Hide file tree

Showing 3 changed files with 51 additions and 41 deletions.
diff --git a/Source/Fortran/EigenBoundsModule.F90 b/Source/Fortran/EigenBoundsModule.F90
@@ -97,13 +97,13 @@ SUBROUTINE PowerBounds(this, max_value, solver_parameters_in)
 
     !! Guess Vector
     CALL ConstructTripletList(temp_list)
-    IF (this%process_grid%global_rank .EQ. 0) THEN
+    DO II = this%start_column, this%end_column - 1
        temp_triplet%index_row = 1
-       temp_triplet%index_column = 1
-       temp_triplet%point_value = 1.0_NTREAL
+       temp_triplet%index_column = II
+       temp_triplet%point_value = 1.0_NTREAL / vector%actual_matrix_dimension
        CALL AppendToTripletList(temp_list, temp_triplet)
-    END IF
-    CALL FillMatrixFromTripletList(vector, temp_list)
+    END DO
+    CALL FillMatrixFromTripletList(vector, temp_list, prepartitioned_in=.TRUE.)
 
     !! Iterate
     IF (params%be_verbose) THEN

diff --git a/Source/Fortran/EigenSolversModule.F90 b/Source/Fortran/EigenSolversModule.F90
@@ -8,7 +8,8 @@ MODULE EigenSolversModule
 #if EIGENEXA
   USE EigenExaModule, ONLY : EigenExa_s
 #endif
-  USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteHeader, WriteElement
+  USE LoggingModule, ONLY : EnterSubLog, ExitSubLog, WriteHeader, &
+       & WriteComment, WriteElement
   USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, &
        & DestructMatrixMemoryPool
   USE PSMatrixAlgebraModule, ONLY : IncrementMatrix, MatrixMultiply, &
@@ -176,10 +177,19 @@ SUBROUTINE EstimateGap(H, K, chemical_potential, gap,  solver_parameters_in)
        CALL PrintParameters(params)
     END IF
 
+    !! Use the power method to get the smallest eigenvalue. This works
+    !! because we project out the rydberg states, but if not we fall
+    !! back to the Gershgorin Bounds
+    CALL MatrixMultiply(K, H, KH, &
+         & threshold_in = params%threshold, memory_pool_in = pool)
+    CALL PowerBounds(KH, e_min, params)
+    IF (e_min .GT. 0_NTREAL) THEN
+       CALL WriteComment("Switching to GershgorinBounds")
+       CALL GershgorinBounds(H, e_min, e_max)
+    END IF
+    CALL WriteElement("Estimated e_min", VALUE = e_min)
+
     !! Shift the matrix so the HOMO is the largest magnitude eigenvalue
-    CALL GershgorinBounds(H, e_min, e_max)
-    CALL WriteElement("Gershgorin e_min", VALUE = e_min)
-    CALL WriteElement("Gershgorin e_max", VALUE = e_max)
     CALL ConstructEmptyMatrix(ShiftH, H)
     CALL FillMatrixIdentity(ShiftH)
     CALL ScaleMatrix(ShiftH, -e_min)

diff --git a/UnitTests/test_solvers.py b/UnitTests/test_solvers.py
@@ -1055,38 +1055,38 @@ def test_svd(self):
         comm.barrier()
         self.check_result()
 
-    # def test_estimategap(self):
-    #     '''Test routines to estimate homo-lumo gap.'''
-    #     from scipy.linalg import eigh
-    #     # Starting Matrix
-    #     H = self.create_matrix(add_gap=True)
-    #     self.write_matrix(H, self.input_file)
-
-    #     # Check Matrix
-    #     vals, _ = eigh(H.todense())
-    #     nel = int(self.mat_dim/2)
-    #     gap_gold = vals[nel] - vals[nel - 1]
-    #     cp = vals[nel - 1] + 0.5 * gap_gold
-
-    #     # Compute
-    #     Hmat = nt.Matrix_ps(self.input_file, False)
-    #     Kmat = nt.Matrix_ps(Hmat.GetActualDimension())
-
-    #     # Density Part
-    #     ISQ = nt.Matrix_ps(Hmat.GetActualDimension())
-    #     ISQ.FillIdentity()
-    #     SRoutine = nt.DensityMatrixSolvers.TRS4
-    #     _, _ = SRoutine(Hmat, ISQ, nel, Kmat, self.isp)
-
-    #     # Estimate Driver
-    #     gap = nt.EigenSolvers.EstimateGap(Hmat, Kmat, cp, self.isp)
-
-    #     # Check the value. Threshold has to be lose because of degenerate
-    #     # eigenvalues near the gap.
-    #     threshold = 0.5
-    #     relative_error = abs(gap_gold - gap)
-    #     global_error = comm.bcast(relative_error, root=0)
-    #     self.assertLessEqual(global_error, threshold)
+    def test_estimategap(self):
+        '''Test routines to estimate homo-lumo gap.'''
+        from scipy.linalg import eigh
+        # Starting Matrix
+        H = self.create_matrix(add_gap=True)
+        self.write_matrix(H, self.input_file)
+
+        # Check Matrix
+        vals, _ = eigh(H.todense())
+        nel = int(self.mat_dim/2)
+        gap_gold = vals[nel] - vals[nel - 1]
+        cp = vals[nel - 1] + 0.5 * gap_gold
+
+        # Compute
+        Hmat = nt.Matrix_ps(self.input_file, False)
+        Kmat = nt.Matrix_ps(Hmat.GetActualDimension())
+
+        # Density Part
+        ISQ = nt.Matrix_ps(Hmat.GetActualDimension())
+        ISQ.FillIdentity()
+        SRoutine = nt.DensityMatrixSolvers.TRS4
+        _, _ = SRoutine(Hmat, ISQ, nel, Kmat, self.isp)
+
+        # Estimate Driver
+        gap = nt.EigenSolvers.EstimateGap(Hmat, Kmat, cp, self.isp)
+
+        # Check the value. Threshold has to be lose because of degenerate
+        # eigenvalues near the gap.
+        threshold = 0.5
+        relative_error = abs(gap_gold - gap)
+        global_error = comm.bcast(relative_error, root=0)
+        self.assertLessEqual(global_error, threshold)
 
 
 class TestSolvers_r(TestSolvers):