Dense FOE for reference purposes

Source/CPlusPlus/CMakeLists.txt

@@ -6,6 +6,7 @@ set(Csrc
     EigenBounds.cc
     EigenSolvers.cc
     ExponentialSolvers.cc
+    FermiOperator.cc
     GeometryOptimization.cc
     HermiteSolvers.cc
     InverseSolvers.cc
@@ -37,6 +38,7 @@ set(Chead
     EigenBounds.h
     EigenSolvers.h
     ExponentialSolvers.h
+    FermiOperator.h
     GeometryOptimization.h
     HermiteSolvers.h
     InverseSolvers.h

Source/Fortran/CMakeLists.txt

@@ -11,6 +11,7 @@ set(Fsrc
     EigenSolversModule.F90
     ErrorModule.F90
     ExponentialSolversModule.F90
+    FermiOperatorModule.F90
     GemmTasksModule.F90
     GeometryOptimizationModule.F90
     HermiteSolversModule.F90

Source/Fortran/DensityMatrixSolversModule.F90

@@ -3,7 +3,7 @@
 MODULE DensityMatrixSolversModule
   USE DataTypesModule, ONLY : NTREAL, MPINTREAL
   USE EigenBoundsModule, ONLY : GershgorinBounds
-  USE EigenSolversModule, ONLY : EigenDecomposition
+  USE FermiOperatorModule, ONLY : ComputeDenseFOE
   USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix
   USE LoggingModule, ONLY : WriteElement, WriteListElement, WriteHeader, &
        & EnterSubLog, ExitSubLog
@@ -14,11 +14,9 @@ MODULE DensityMatrixSolversModule
        & DotMatrix, MatrixTrace, ScaleMatrix
   USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, DestructMatrix, &
        & CopyMatrix, PrintMatrixInformation, FillMatrixIdentity, &
-       & TransposeMatrix, FillMatrixFromTripletList, ConjugateMatrix, &
-       & GetMatrixTripletList
+       & TransposeMatrix
   USE SolverParametersModule, ONLY : SolverParameters_t, PrintParameters, &
        & DestructSolverParameters
-  USE TripletListModule, ONLY : TripletList_r, DestructTripletList
   IMPLICIT NONE
   PRIVATE
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1181,16 +1179,7 @@ SUBROUTINE DenseDensity(H, ISQ, nel, K, &
     TYPE(SolverParameters_t), INTENT(IN), OPTIONAL :: solver_parameters_in
     !! Handling Optional Parameters
     TYPE(SolverParameters_t) :: params
-    !! Local Variables
-    TYPE(Matrix_ps) :: ISQT, WH
-    TYPE(Matrix_ps) :: WD
-    TYPE(Matrix_ps) :: vecs, vecsT, vals, Temp
-    TYPE(MatrixMemoryPool_p) :: pool
-    TYPE(TripletList_r) :: tlist
     REAL(NTREAL) :: chemical_potential, energy_value
-    REAL(NTREAL) :: homo, lumo
-    INTEGER :: II
-    INTEGER :: ierr
 
     !! Optional Parameters
     IF (PRESENT(solver_parameters_in)) THEN
@@ -1199,93 +1188,20 @@ SUBROUTINE DenseDensity(H, ISQ, nel, K, &
        params = SolverParameters_t()
     END IF
 
-    IF (params%be_verbose) THEN
-       CALL WriteHeader("Density Matrix Solver")
-       CALL EnterSubLog
-       CALL WriteElement(key="Method", VALUE="Dense")
-       CALL PrintParameters(params)
-    END IF
-
-    !! Compute the working hamiltonian.
-    CALL TransposeMatrix(ISQ, ISQT)
-    CALL MatrixMultiply(ISQ, H, Temp, &
-         & threshold_in=params%threshold, memory_pool_in=pool)
-    CALL MatrixMultiply(Temp, ISQT, WH, &
-         & threshold_in=params%threshold, memory_pool_in=pool)
-
-    !! Perform the eigendecomposition
-    CALL EigenDecomposition(WH, vecs, vals, params)
-
-    !! Convert to a triplet list, map, and get homo/lumo + energy.
-    CALL GetMatrixTripletList(vals, tlist)
-    homo = 0.0_NTREAL
-    lumo = 0.0_NTREAL
-    energy_value = 0.0_NTREAL
-    DO II = 1, tlist%CurrentSize
-       !! HOMO/LUMO
-       IF (tlist%DATA(II)%index_column .EQ. INT(nel/2)) THEN
-          homo = tlist%DATA(II)%point_value
-       ELSE IF (tlist%DATA(II)%index_column .EQ. INT(nel/2) + 1) THEN
-          lumo = tlist%DATA(II)%point_value
-       END IF
-       !! Mapping
-       IF (tlist%DATA(II)%index_column .LE. INT(nel/2)) THEN
-          energy_value = energy_value + tlist%DATA(II)%point_value
-          tlist%DATA(II)%point_value = 1.0_NTREAL
-       ELSE
-          tlist%DATA(II)%point_value = 0.0_NTREAL
-       ENDIF
-    END DO
-
-    !! Compute MU
-    CALL MPI_ALLREDUCE(MPI_IN_PLACE, homo, 1, MPINTREAL, MPI_SUM, &
-         & H%process_grid%within_slice_comm, ierr)
-    CALL MPI_ALLREDUCE(MPI_IN_PLACE, lumo, 1, MPINTREAL, MPI_SUM, &
-         & H%process_grid%within_slice_comm, ierr)
-    CALL MPI_ALLREDUCE(MPI_IN_PLACE, energy_value, 1, MPINTREAL, MPI_SUM, &
-         & H%process_grid%within_slice_comm, ierr)
-    chemical_potential = homo + 0.5_NTREAL * (lumo - homo)
-
-    !! Fill
-    CALL ConstructEmptyMatrix(vals, H)
-    CALL FillMatrixFromTripletList(vals, tlist, preduplicated_in=.TRUE.)
-
-    !! Multiply Back Together
-    CALL MatrixMultiply(vecs, vals, temp, threshold_in=params%threshold)
-    CALL TransposeMatrix(vecs, vecsT)
-    CALL ConjugateMatrix(vecsT)
-    CALL MatrixMultiply(temp, vecsT, WD, &
-         & threshold_in=params%threshold)
-
-    !! Compute the density matrix in the non-orthogonalized basis
-    CALL MatrixMultiply(ISQT, WD, Temp, &
-         & threshold_in=params%threshold, memory_pool_in=pool)
-    CALL MatrixMultiply(Temp, ISQ, K, &
-         & threshold_in=params%threshold, memory_pool_in=pool)
+    !! Call the unified routine.
+    CALL ComputeDenseFOE(H, ISQ, nel, K, energy_value_out=energy_value, &
+         & chemical_potential_out=chemical_potential, &
+         & solver_parameters_in=params)
 
     !! Optional out variables.
     IF (PRESENT(energy_value_out)) THEN
-       energy_value_out = 2.0_NTREAL * energy_value
+       energy_value_out = energy_value
     END IF
     IF (PRESENT(chemical_potential_out)) THEN
        chemical_potential_out = chemical_potential
     END IF
 
     !! Cleanup
-    CALL DestructMatrix(WH)
-    CALL DestructMatrix(WD)
-    CALL DestructMatrix(ISQT)
-    CALL DestructMatrix(vecs)
-    CALL DestructMatrix(vecst)
-    CALL DestructMatrix(vals)
-    CALL DestructMatrix(temp)
-    CALL DestructTripletList(tlist)
-    CALL DestructMatrixMemoryPool(pool)
-
-    IF (params%be_verbose) THEN
-       CALL ExitSubLog
-    END IF
-
     CALL DestructSolverParameters(params)
   END SUBROUTINE DenseDensity
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Source/Fortran/EigenExaModule.F90

@@ -15,7 +15,7 @@ MODULE EigenExaModule
        & AppendToTripletList, GetTripletAt, ConstructTripletList, &
        & DestructTripletList, RedistributeTripletLists
   USE eigen_libs_mod
-  USE MPI
+  USE NTMPIModule
   IMPLICIT NONE
   PRIVATE
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Source/Fortran/EigenSolversModule.F90

@@ -19,7 +19,7 @@ MODULE EigenSolversModule
        & DestructMatrix
   USE TripletListModule, ONLY : TripletList_r, TripletList_c, &
        & ConstructTripletList, DestructTripletList
-  USE MPI
+  USE NTMPIModule
   IMPLICIT NONE
   PRIVATE
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Source/Fortran/SingularValueSolversModule.F90

@@ -8,7 +8,6 @@ MODULE SingularValueSolversModule
   USE SignSolversModule, ONLY : PolarDecomposition
   USE SolverParametersModule, ONLY : SolverParameters_t, PrintParameters, &
        & DestructSolverParameters
-  USE MPI
   IMPLICIT NONE
   PRIVATE
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Source/Swig/NTPolySwig.i

@@ -15,6 +15,7 @@
 #include "EigenBounds.h"
 #include "EigenSolvers.h"
 #include "ExponentialSolvers.h"
+#include "FermiOperator.h"
 #include "GeometryOptimization.h"
 #include "HermiteSolvers.h"
 #include "InverseSolvers.h"
@@ -54,6 +55,7 @@ using namespace NTPoly;
 %include "EigenBounds.h"
 %include "EigenSolvers.h"
 %include "ExponentialSolvers.h"
+%include "FermiOperator.h"
 %include "GeometryOptimization.h"
 %include "HermiteSolvers.h"
 %include "InverseSolvers.h"

Source/Wrapper/CMakeLists.txt

@@ -6,6 +6,7 @@ set(Wsrc
     EigenBoundsModule_wrp.F90
     EigenSolversModule_wrp.F90
     ExponentialSolversModule_wrp.F90
+    FermiOperatorModule_wrp.F90
     GeometryOptimizationModule_wrp.F90
     HermiteSolversModule_wrp.F90
     InverseSolversModule_wrp.F90

Source/Wrapper/DensityMatrixSolversModule_wrp.F90

@@ -174,8 +174,9 @@ SUBROUTINE DenseDensity_wrp(ih_Hamiltonian, ih_InverseSquareRoot, nel, &
     h_solver_parameters = TRANSFER(ih_solver_parameters, h_solver_parameters)
 
     CALL DenseDensity(h_Hamiltonian%DATA, h_InverseSquareRoot%DATA, INT(nel), &
-         & h_Density%DATA, energy_value_out, chemical_potential_out, &
-         & h_solver_parameters%DATA)
+         & h_Density%DATA, energy_value_out=energy_value_out, &
+         & chemical_potential_out=chemical_potential_out, &
+         & solver_parameters_in=h_solver_parameters%DATA)
   END SUBROUTINE DenseDensity_wrp
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Compute the energy-weighted density matrix.

UnitTests/test_chemistry.py

@@ -185,7 +185,7 @@ def check_cp(self, computed):
         else:
             self.assertTrue(False)
 
-    def basic_solver(self, SRoutine, cpcheck=True):
+    def basic_solver(self, SRoutine, cpcheck=True, temp=None):
         '''Test various kinds of density matrix solvers.'''
         from helpers import result_file
 
@@ -202,10 +202,21 @@ def basic_solver(self, SRoutine, cpcheck=True):
         nt.SquareRootSolvers.InverseSquareRoot(overlap_matrix,
                                                inverse_sqrt_matrix,
                                                self.solver_parameters)
-        energy_value, chemical_potential = SRoutine(fock_matrix,
-                                                    inverse_sqrt_matrix,
-                                                    self.nel, density_matrix,
-                                                    self.solver_parameters)
+        if temp is None:
+            energy_value, chemical_potential = SRoutine(fock_matrix,
+                                                        inverse_sqrt_matrix,
+                                                        self.nel,
+                                                        density_matrix,
+                                                        self.solver_parameters)
+        else:
+            inv_temp = 1/(temp * 3.166811563*10**(-6))
+            print("::: Temperature", temp, inv_temp)
+            energy_value, chemical_potential = SRoutine(fock_matrix,
+                                                        inverse_sqrt_matrix,
+                                                        self.nel,
+                                                        density_matrix,
+                                                        inv_temp,
+                                                        self.solver_parameters)
 
         density_matrix.WriteToMatrixMarket(result_file)
         comm.barrier()
@@ -267,6 +278,10 @@ def test_densedensity(self):
         '''Test routines to compute the density matrix with Dense Method.'''
         self.basic_solver(nt.DensityMatrixSolvers.DenseDensity)
 
+    def test_foe_low(self):
+        '''Test the fermi operator expansion at a low temperature.'''
+        self.basic_solver(nt.FermiOperator.ComputeDenseFOE, temp=50)
+
     def test_energy_density(self):
         '''Test the routines to compute the weighted-energy density matrix.'''
         from helpers import THRESHOLD, result_file