Forgot to add files

Source/C/FermiOperator_c.h

@@ -0,0 +1,10 @@
+#ifndef DENSITYMATRIXSOLVERS_ch
+#define DENSITYMATRIXSOLVERS_ch
+
+void ComputeDenseFOE_wrp(const int *ih_Hamiltonian,
+                         const int *ih_InverseSquareRoot, const int *nel,
+                         int *ih_Density, const double *inv_temp_in,
+                         const double *energy_value_out,
+                         const double *chemical_potential_out,
+                         const int *ih_solver_parameters);
+#endif

Source/CPlusPlus/FermiOperator.cc

@@ -0,0 +1,24 @@
+#include "FermiOperator.h"
+#include "PSMatrix.h"
+#include "SolverParameters.h"
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" {
+#include "FermiOperator_c.h"
+}
+
+////////////////////////////////////////////////////////////////////////////////
+namespace NTPoly {
+////////////////////////////////////////////////////////////////////////////////
+void FermiOperator::ComputeDenseFOE(const Matrix_ps &Hamiltonian,
+                                    const Matrix_ps &Overlap, int nel,
+                                    Matrix_ps &Density, double inv_temp,
+                                    double &energy_value_out,
+                                    double &chemical_potential_out,
+                                    const SolverParameters &solver_parameters) {
+  ComputeDenseFOE_wrp(GetIH(Hamiltonian), GetIH(Overlap), &nel, GetIH(Density),
+                      &inv_temp, &energy_value_out, &chemical_potential_out,
+                      GetIH(solver_parameters));
+}
+
+} // namespace NTPoly

Source/CPlusPlus/FermiOperator.h

@@ -0,0 +1,30 @@
+#ifndef FERMIOPERATOREXPANSION_h
+#define FERMIOPERATOREXPANSION_h
+
+#include "SolverBase.h"
+
+////////////////////////////////////////////////////////////////////////////////
+namespace NTPoly {
+class SolverParameters;
+class Matrix_ps;
+//! A Class For Solving Chemistry Systems Using the Fermi Operator Expansion.
+class FermiOperator : public SolverBase {
+public:
+  //! Compute the density matrix using a dense solver.
+  //!\param Hamiltonian the matrix to compute the corresponding density from.
+  //!\param InverseSquareRoot of the overlap matrix.
+  //!\param nel the number of electrons.
+  //!\param Density the density matrix computed by this routine.
+  //!\param inv_temp the inverse temperature.
+  //!\param energy_value_out the energy of the system.
+  //!\param chemical_potential_out the chemical potential calculated.
+  //!\param solver_parameters parameters for the solver
+  static void ComputeDenseFOE(const Matrix_ps &Hamiltonian,
+                              const Matrix_ps &InverseSquareRoot, int nel,
+                              Matrix_ps &Density, double inv_temp,
+                              double &energy_value_out,
+                              double &chemical_potential_out,
+                              const SolverParameters &solver_parameters);
+};
+} // namespace NTPoly
+#endif

Source/Fortran/FermiOperatorModule.F90

@@ -0,0 +1,177 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!> A Module For Computing The Density Matrix Using the Fermi Operator Expansion
+MODULE FermiOperatorModule
+  USE DataTypesModule, ONLY : NTREAL, MPINTREAL
+  USE EigenSolversModule, ONLY : EigenDecomposition
+  USE LoadBalancerModule, ONLY : PermuteMatrix, UndoPermuteMatrix
+  USE LoggingModule, ONLY : WriteElement, WriteHeader, &
+       & EnterSubLog, ExitSubLog
+  USE PSMatrixAlgebraModule, ONLY : MatrixMultiply
+  USE PSMatrixModule, ONLY : Matrix_ps, ConstructEmptyMatrix, &
+       & FillMatrixFromTripletList, GetMatrixTripletList, &
+       & TransposeMatrix, ConjugateMatrix, DestructMatrix
+  USE PMatrixMemoryPoolModule, ONLY : MatrixMemoryPool_p, &
+       & DestructMatrixMemoryPool
+  USE SolverParametersModule, ONLY : SolverParameters_t, &
+       & PrintParameters, DestructSolverParameters
+  USE TripletListModule, ONLY : TripletList_r, DestructTripletList
+  USE NTMPIModule
+  IMPLICIT NONE
+  PRIVATE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  PUBLIC :: ComputeDenseFOE
+CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Compute the density matrix using a dense routine.
+  SUBROUTINE ComputeDenseFOE(H, ISQ, nel, K, inv_temp_in, &
+       & energy_value_out, chemical_potential_out, solver_parameters_in)
+    !> The matrix to compute the corresponding density from.
+    TYPE(Matrix_ps), INTENT(IN) :: H
+    !> The inverse square root of the overlap matrix.
+    TYPE(Matrix_ps), INTENT(IN) :: ISQ
+    !> The number of electrons.
+    INTEGER, INTENT(IN) :: nel
+    !> The density matrix computed by this routine.
+    TYPE(Matrix_ps), INTENT(INOUT) :: K
+    !> The inverse temperature for smearing (optional).
+    REAL(NTREAL), INTENT(IN), OPTIONAL :: inv_temp_in
+    !> The energy of the system (optional).
+    REAL(NTREAL), INTENT(OUT), OPTIONAL :: energy_value_out
+    !> The chemical potential (optional).
+    REAL(NTREAL), INTENT(OUT), OPTIONAL :: chemical_potential_out
+    !> Parameters for the solver (optional).
+    TYPE(SolverParameters_t), INTENT(IN), OPTIONAL :: solver_parameters_in
+    !! Handling Optional Parameters
+    TYPE(SolverParameters_t) :: params
+    REAL(NTREAL) :: inv_temp
+    LOGICAL :: do_smearing
+    !! 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
+    REAL(NTREAL) :: sval, occ
+    INTEGER :: II
+    INTEGER :: ierr
+
+    !! Optional Parameters
+    IF (PRESENT(solver_parameters_in)) THEN
+       params = solver_parameters_in
+    ELSE
+       params = SolverParameters_t()
+    END IF
+    IF (PRESENT(inv_temp_in)) THEN
+       inv_temp = inv_temp_in
+       do_smearing = .TRUE.
+    ELSE
+       do_smearing = .FALSE.
+    END IF
+
+    IF (params%be_verbose) THEN
+       CALL WriteHeader("Density Matrix Solver")
+       CALL EnterSubLog
+       IF (do_smearing) THEN
+          CALL WriteElement(key="Method", VALUE="Dense FOE")
+          CALL WriteElement(key="InverseTemperature", VALUE=inv_temp)
+       ELSE
+          CALL WriteElement(key="Method", VALUE="Dense Step Function")
+       END IF
+       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 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
+       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
+    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)
+    chemical_potential = homo + 0.5_NTREAL * (lumo - homo)
+
+    !! Map
+    DO II = 1, tlist%CurrentSize
+       IF (.NOT. do_smearing) THEN
+          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
+       ELSE
+          sval = tlist%DATA(II)%point_value - chemical_potential
+          occ = 1.0 / (1 + EXP(inv_temp * sval))
+          WRITE(*,*) II, occ
+          energy_value = energy_value + occ * tlist%DATA(II)%point_value
+          tlist%DATA(II)%point_value = occ
+       END IF
+    END DO
+    CALL MPI_ALLREDUCE(MPI_IN_PLACE, energy_value, 1, MPINTREAL, MPI_SUM, &
+         & H%process_grid%within_slice_comm, ierr)
+
+    !! 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)
+
+    !! Optional out variables.
+    IF (PRESENT(energy_value_out)) THEN
+       energy_value_out = 2.0_NTREAL * 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 ComputeDenseFOE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+END MODULE FermiOperatorModule

Source/Wrapper/FermiOperatorModule_wrp.F90

@@ -0,0 +1,44 @@
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!> Wraps the density matrix solvers module for calling from other languages.
+MODULE FermiOperatorModule_wrp
+  USE DataTypesModule, ONLY : NTREAL
+  USE FermiOperatorModule, ONLY : ComputeDenseFOE
+  USE PSMatrixModule_wrp, ONLY : Matrix_ps_wrp
+  USE SolverParametersModule_wrp, ONLY : SolverParameters_wrp
+  USE WrapperModule, ONLY : SIZE_wrp
+  USE ISO_C_BINDING, ONLY : c_int
+  IMPLICIT NONE
+  PRIVATE
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  PUBLIC :: ComputeDenseFOE_wrp
+CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Compute the density matrix from a Hamiltonian using the PM method.
+  SUBROUTINE ComputeDenseFOE_wrp(ih_Hamiltonian, ih_InverseSquareRoot, nel, &
+       & ih_Density, inv_temp_in, energy_value_out, chemical_potential_out, &
+       & ih_solver_parameters) BIND(c,name="ComputeDenseFOE_wrp")
+    INTEGER(kind=c_int), INTENT(IN) :: ih_Hamiltonian(SIZE_wrp)
+    INTEGER(kind=c_int), INTENT(IN) :: ih_InverseSquareRoot(SIZE_wrp)
+    INTEGER(kind=c_int), INTENT(IN) :: nel
+    INTEGER(kind=c_int), INTENT(INOUT) :: ih_Density(SIZE_wrp)
+    REAL(NTREAL), INTENT(IN) :: inv_temp_in
+    REAL(NTREAL), INTENT(OUT) :: energy_value_out
+    REAL(NTREAL), INTENT(OUT) :: chemical_potential_out
+    INTEGER(kind=c_int), INTENT(IN) :: ih_solver_parameters(SIZE_wrp)
+    TYPE(Matrix_ps_wrp) :: h_Hamiltonian
+    TYPE(Matrix_ps_wrp) :: h_InverseSquareRoot
+    TYPE(Matrix_ps_wrp) :: h_Density
+    TYPE(SolverParameters_wrp) :: h_solver_parameters
+
+    h_Hamiltonian = TRANSFER(ih_Hamiltonian,h_Hamiltonian)
+    h_InverseSquareRoot = TRANSFER(ih_InverseSquareRoot,h_InverseSquareRoot)
+    h_Density = TRANSFER(ih_Density,h_Density)
+    h_solver_parameters = TRANSFER(ih_solver_parameters, h_solver_parameters)
+
+    CALL ComputeDenseFOE(h_Hamiltonian%DATA, h_InverseSquareRoot%DATA, &
+         & INT(nel), h_Density%DATA, inv_temp_in=inv_temp_in, &
+         & energy_value_out=energy_value_out, &
+         & chemical_potential_out=chemical_potential_out, &
+         & solver_parameters_in=h_solver_parameters%DATA)
+  END SUBROUTINE ComputeDenseFOE_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+END MODULE FermiOperatorModule_wrp