From 3887bffc89d6ff55fd9f8faec22a6512cb12e980 Mon Sep 17 00:00:00 2001 From: William Dawson Date: Mon, 6 Dec 2021 19:03:58 +0900 Subject: [PATCH] Forgot to add files --- Source/C/FermiOperator_c.h | 10 ++ Source/CPlusPlus/FermiOperator.cc | 24 +++ Source/CPlusPlus/FermiOperator.h | 30 ++++ Source/Fortran/FermiOperatorModule.F90 | 177 +++++++++++++++++++++ Source/Wrapper/FermiOperatorModule_wrp.F90 | 44 +++++ 5 files changed, 285 insertions(+) create mode 100644 Source/C/FermiOperator_c.h create mode 100644 Source/CPlusPlus/FermiOperator.cc create mode 100644 Source/CPlusPlus/FermiOperator.h create mode 100644 Source/Fortran/FermiOperatorModule.F90 create mode 100644 Source/Wrapper/FermiOperatorModule_wrp.F90 diff --git a/Source/C/FermiOperator_c.h b/Source/C/FermiOperator_c.h new file mode 100644 index 00000000..d4bbba7a --- /dev/null +++ b/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 diff --git a/Source/CPlusPlus/FermiOperator.cc b/Source/CPlusPlus/FermiOperator.cc new file mode 100644 index 00000000..fd3ddac0 --- /dev/null +++ b/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 diff --git a/Source/CPlusPlus/FermiOperator.h b/Source/CPlusPlus/FermiOperator.h new file mode 100644 index 00000000..2ca5765f --- /dev/null +++ b/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 diff --git a/Source/Fortran/FermiOperatorModule.F90 b/Source/Fortran/FermiOperatorModule.F90 new file mode 100644 index 00000000..61f52cac --- /dev/null +++ b/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 diff --git a/Source/Wrapper/FermiOperatorModule_wrp.F90 b/Source/Wrapper/FermiOperatorModule_wrp.F90 new file mode 100644 index 00000000..43c190bb --- /dev/null +++ b/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