Skip to content

Commit

Permalink
Fix block hermitian code and clean up (dftbplus#1337)
Browse files Browse the repository at this point in the history 
* Clean up matrix sym/hermit routines

Add unit testing and remove the atomic block version of routines.
Also add distributed matrix case.

* Module re-name

* Routine rename
  • Loading branch information
@bhourahine
bhourahine authored Feb 6, 2024
1 parent 78e76e9 commit 7530d37
Show file tree
Hide file tree
Showing 16 changed files with 345 additions and 241 deletions.
36 changes: 18 additions & 18 deletions src/dftbp/dftb/hybridxc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,12 +17,12 @@ module dftbp_dftb_hybridxc
use dftbp_common_status, only : TStatus
use dftbp_dftb_nonscc, only : TNonSccDiff
use dftbp_dftb_slakocont, only : TSlakoCont
use dftbp_dftb_sparse2dense, only : blockSymmetrizeHS, symmetrizeHS, hermitianSquareMatrix
use dftbp_math_blasroutines, only : gemm, symm
use dftbp_math_sorting, only : index_heap_sort
use dftbp_type_commontypes, only : TOrbitals, TParallelKS
use dftbp_dftb_periodic, only : TSymNeighbourList, getCellTranslations, cart2frac
use dftbp_dftb_nonscc, only : buildS
use dftbp_math_matrixops, only : adjointLowerTriangle
use dftbp_math_simplealgebra, only : determinant33
use dftbp_common_parallel, only : getStartAndEndIndex
use dftbp_math_wignerseitz, only : generateWignerSeitzGrid
Expand Down Expand Up @@ -1557,14 +1557,14 @@ subroutine allocateAndInit()

allocate(tmpOvr(matrixSize, matrixSize))
tmpOvr(:,:) = overlap
call blockSymmetrizeHS(tmpOvr, iSquare)
call adjointLowerTriangle(tmpOvr)

allocate(tmpDHam(matrixSize, matrixSize))
tmpDHam(:,:) = 0.0_dp

allocate(tmpDRho(matrixSize, matrixSize))
tmpDRho(:,:) = deltaRho
call symmetrizeHS(tmpDRho)
call adjointLowerTriangle(tmpDRho)
call checkAndInitScreening(this, matrixSize, tmpDRho)

allocate(tmpDDRho(matrixSize, matrixSize))
Expand Down Expand Up @@ -1743,7 +1743,7 @@ subroutine addCamHamiltonianNeighbour_real(this, deltaRhoSqr, overSparse, iNeigh
tmpHH(:,:) = 0.125_dp * tmpHH
end if

call symmetrizeHS(tmpHH)
call adjointLowerTriangle(tmpHH)

HSqrReal(:,:) = HSqrReal + tmpHH
this%camEnergy = this%camEnergy + evaluateEnergy_real(tmpHH, tmpDRho)
Expand All @@ -1763,7 +1763,7 @@ subroutine allocateAndInit(tmpHH, tmpDRho)
tmpHH(:,:) = 0.0_dp
allocate(tmpDRho(size(deltaRhoSqr, dim=1), size(deltaRhoSqr, dim=1)))
tmpDRho(:,:) = deltaRhoSqr
call symmetrizeHS(tmpDRho)
call adjointLowerTriangle(tmpDRho)

end subroutine allocateAndInit

Expand Down Expand Up @@ -2162,11 +2162,11 @@ subroutine allocateAndInit(this, iSquare, overlap, densSqr, HH, Smat, Dmat, camG
nAtom = size(this%camGammaEval0, dim=1)

! Symmetrize Hamiltonian, overlap, density matrices
call symmetrizeHS(HH)
call adjointLowerTriangle(HH)
Smat(:,:) = overlap
call symmetrizeHS(Smat)
call adjointLowerTriangle(Smat)
Dmat(:,:) = densSqr
call symmetrizeHS(Dmat)
call adjointLowerTriangle(Dmat)

! Get CAM gamma variable
camGammaAO(:,:) = 0.0_dp
Expand Down Expand Up @@ -2335,11 +2335,11 @@ subroutine allocateAndInit(this, iSquare, overlap, densSqr, HH, Smat, Dmat, camG
nAtom = size(this%camGammaEval0, dim=1)

!! Symmetrize Hamiltonian, overlap, density matrices
call hermitianSquareMatrix(HH)
call adjointLowerTriangle(HH)
Smat(:,:) = overlap
call hermitianSquareMatrix(Smat)
call adjointLowerTriangle(Smat)
Dmat(:,:) = densSqr
call hermitianSquareMatrix(Dmat)
call adjointLowerTriangle(Dmat)

! Get CAM gamma variable
camGammaAO(:,:) = 0.0_dp
Expand Down Expand Up @@ -4374,10 +4374,10 @@ subroutine allocateAndInit(tmpOvr, tmpRho, tmpderiv)
tmpOvr(:,:) = SSqrReal
tmpRho(:,:,:) = deltaRhoSqr

call symmetrizeHS(tmpOvr)
call adjointLowerTriangle(tmpOvr)

do iSpin = 1, size(deltaRhoSqr, dim=3)
call symmetrizeHS(tmpRho(:,:, iSpin))
call adjointLowerTriangle(tmpRho(:,:, iSpin))
end do

end subroutine allocateAndInit
Expand Down Expand Up @@ -4497,7 +4497,7 @@ subroutine addCamGradientsNeighbour_gamma(this, deltaRhoSqr, skOverCont, symNeig

tmpDeltaRhoSqr = deltaRhoSqr
do iSpin = 1, nSpin
call symmetrizeHS(tmpDeltaRhoSqr(:,:, iSpin))
call adjointLowerTriangle(tmpDeltaRhoSqr(:,:, iSpin))
end do

loopK1: do iAtK = 1, nAtom0
Expand Down Expand Up @@ -5265,10 +5265,10 @@ subroutine addCamGradientsMatrix_real(this, deltaRhoSqr, overlap, skOverCont,&

! symmetrize square overlap and density matrix
overlapSym = overlap
call symmetrizeHS(overlapSym)
call adjointLowerTriangle(overlapSym)
deltaRhoSqrSym = deltaRhoSqr
do iSpin = 1, nSpin
call symmetrizeHS(deltaRhoSqrSym(:,:, iSpin))
call adjointLowerTriangle(deltaRhoSqrSym(:,:, iSpin))
end do

! get CAM \tilde{gamma} super-matrix
Expand Down Expand Up @@ -6117,8 +6117,8 @@ function evaluateLrEnergyDirect_cluster(this, env, deltaRho, overlap, iSquare) r
tmpOvr(:,:) = overlap
tmpDRho(:,:) = deltaRho

call symmetrizeHS(tmpOvr)
call symmetrizeHS(tmpDRho)
call adjointLowerTriangle(tmpOvr)
call adjointLowerTriangle(tmpDRho)

energy = 0.0_dp
do iAt1 = 1, nAtom0
Expand Down
9 changes: 4 additions & 5 deletions src/dftbp/dftb/populations.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,15 @@
module dftbp_dftb_populations
use dftbp_common_accuracy, only : dp
use dftbp_common_constants, only : pi
use dftbp_type_commontypes, only : TOrbitals, TParallelKS
use dftbp_dftb_hybridxc, only : THybridXcFunc
use dftbp_dftb_sparse2dense, only : symmetrizeHS
use dftbp_common_environment, only : TEnvironment, globalTimers
use dftbp_common_schedule, only : distributeRangeWithWorkload, assembleChunks
use dftbp_type_commontypes, only : TOrbitals
use dftbp_type_integral, only : TIntegral
use dftbp_dftb_hybridxc, only : THybridXcFunc
use dftbp_dftb_periodic, only : TNeighbourList
use dftbp_dftb_sparse2dense, only : unpackHS
use dftbp_math_matrixops, only : adjointLowerTriangle
use dftbp_type_commontypes, only : TOrbitals, TParallelKS
use dftbp_type_densedescr, only : TDenseDescr
use dftbp_type_integral, only : TIntegral
#:if WITH_SCALAPACK
use dftbp_extlibs_mpifx, only : MPI_SUM, mpifx_allreduceip
use dftbp_extlibs_scalapackfx, only : DLEN_, NB_, MB_, CSRC_, RSRC_, scalafx_indxl2g,&
Expand Down
164 changes: 8 additions & 156 deletions src/dftbp/dftb/sparse2dense.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ module dftbp_dftb_sparse2dense
use dftbp_common_constants, only : pi, imag
use dftbp_dftb_periodic, only : TNeighbourList
use dftbp_math_angmomentum, only : rotateZ
use dftbp_math_matrixops, only : adjointLowerTriangle
use dftbp_type_commontypes, only : TOrbitals
use dftbp_type_densedescr, only : TDenseDescr
#:if WITH_SCALAPACK
Expand All @@ -25,7 +26,6 @@ module dftbp_dftb_sparse2dense

private
public :: unpackHS, packHS, iPackHS, packErho, unpackDQ
public :: blockSymmetrizeHS, blockHermitianHS, symmetrizeHS, hermitianSquareMatrix
public :: packHSPauli, packHSPauliImag, unpackHPauli, unpackSPauli
public :: unpackHelicalHS, packHelicalHS
public :: getSparseDescriptor
Expand Down Expand Up @@ -90,31 +90,10 @@ module dftbp_dftb_sparse2dense
end interface packErho


!> Symmetrize the square matrix except the on-site blocks
interface blockSymmetrizeHS
module procedure blockSymmetrizeHS_real
module procedure blockSymmetrizeHS_cmplx
end interface blockSymmetrizeHS


!> Hermitian the square matrix except the on-site blocks
interface blockHermitianHS
module procedure blockSymmetrizeHS_real
module procedure blockHermitianHS_cmplx
end interface blockHermitianHS


!> Symmetrize the square matrix including the on-site blocks
interface symmetrizeHS
module procedure symmetrizeHS_real
end interface symmetrizeHS


contains

!> Unpacks sparse matrix to square form (complex version) Note the non on-site blocks are only
!> filled in the lower triangle part of the matrix. To fill the matrix completely, apply the
!> blockSymmetrizeHS subroutine.
!> filled in the lower triangle part of the matrix.
subroutine unpackHS_cmplx_kpts(square, orig, kPoint, iNeighbour, nNeighbourSK, iCellVec, cellVec,&
& iAtomStart, iSparseStart, img2CentCell)

Expand Down Expand Up @@ -194,9 +173,8 @@ end subroutine unpackHS_cmplx_kpts


!> Unpacks sparse matrix to square form (real version for Gamma point)
!>
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix. To fill
!> the matrix completely, apply the blockSymmetrizeHS subroutine.
!!
!! Note: The non on-site blocks are only filled in the lower triangle part of the matrix.
subroutine unpackHS_real(square, orig, iNeighbour, nNeighbourSK, iAtomStart, iSparseStart,&
& img2CentCell)

Expand Down Expand Up @@ -256,8 +234,7 @@ end subroutine unpackHS_real


!> Unpacks sparse matrix to square form (complex version) for helical geometries. Note the non
!> on-site blocks are only filled in the lower triangle part of the matrix. To fill the matrix
!> completely, apply the blockSymmetrizeHS subroutine.
!> on-site blocks are only filled in the lower triangle part of the matrix.
subroutine unpackHSHelical_cmplx(square, orig, kPoint, iNeighbour, nNeighbourSK, iCellVec,&
& cellVec, iAtomStart, iSparseStart, img2CentCell, orb, species, coord)

Expand Down Expand Up @@ -346,8 +323,7 @@ end subroutine unpackHSHelical_cmplx

!> Unpacks sparse matrix to square form (real version for Gamma point) for helical geometry
!>
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix. To fill
!> the matrix completely, apply the blockSymmetrizeHS subroutine.
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix.
subroutine unpackHSHelical_real(square, orig, iNeighbour, nNeighbourSK, iAtomStart, iSparseStart,&
& img2CentCell, orb, species, coord)

Expand Down Expand Up @@ -576,8 +552,7 @@ end subroutine unpackDQ_real

!> Unpacks sparse matrices to square form (2 component version for k-points)
!>
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix. To fill
!> the matrix completely, apply the blockSymmetrizeHS subroutine.
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix.
subroutine unpackHPauli(ham, kPoint, iNeighbour, nNeighbourSK, iSparseStart, iAtomStart,&
& img2CentCell, iCellVec, cellVec, HSqrCplx, iHam)

Expand Down Expand Up @@ -705,8 +680,7 @@ end subroutine unpackHPauli

!> Unpacks sparse overlap matrices to square form (2 component version for k-points)
!>
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix. To fill
!> the matrix completely, apply the blockSymmetrizeHS subroutine.
!> Note: The non on-site blocks are only filled in the lower triangle part of the matrix.
subroutine unpackSPauli(over, kPoint, iNeighbour, nNeighbourSK, iAtomStart, iSparseStart,&
& img2CentCell, iCellVec, cellVec, SSqrCplx)

Expand Down Expand Up @@ -1899,128 +1873,6 @@ subroutine packHSPauliERho_kpts(primitive, square, kPoint, kWeight, iNeighbour,
end subroutine packHSPauliERho_kpts


!> Symmetrize a square matrix leaving the on-site atomic blocks alone. (Complex version)
subroutine blockSymmetrizeHS_cmplx(square, iAtomStart)

!> Square form matrix.
complex(dp), intent(inout) :: square(:, :)

!> Returns the offset array for each atom.
integer, intent(in) :: iAtomStart(:)

integer :: nAtom, iAtom, iStart, iEnd, mOrb

nAtom = size(iAtomStart, dim=1) - 1
mOrb = iAtomStart(nAtom+1) - 1

@:ASSERT(nAtom > 0)
@:ASSERT(size(square, dim=1) == size(square, dim=2))
@:ASSERT((size(square, dim=1) == 2*mOrb) .or. (size(square, dim=1) == mOrb))
@:ASSERT(size(square, dim=1) == mOrb)

do iAtom = 1, nAtom
iStart = iAtomStart(iAtom)
iEnd = iAtomStart(iAtom+1) - 1
square(iStart:iEnd, iEnd+1:mOrb) = transpose(square(iEnd+1:mOrb, iStart:iEnd))
end do

end subroutine blockSymmetrizeHS_cmplx


!> Symmetrize a square matrix leaving the on-site atomic blocks alone. (Complex version)
subroutine blockHermitianHS_cmplx(square, iAtomStart)

!> Square form matrix.
complex(dp), intent(inout) :: square(:, :)

!> Returns the offset array for each atom.
integer, intent(in) :: iAtomStart(:)

integer :: nAtom, iAtom, iStart, iEnd, mOrb

nAtom = size(iAtomStart, dim=1) - 1
mOrb = iAtomStart(nAtom+1) - 1

@:ASSERT(nAtom > 0)
@:ASSERT(size(square, dim=1) == size(square, dim=2))
@:ASSERT((size(square, dim=1) == mOrb) .or. (size(square, dim=1) == 2*mOrb))
@:ASSERT(size(square, dim=1) == mOrb .or. size(square, dim=1) == 2*mOrb)

do iAtom = 1, nAtom
iStart = iAtomStart(iAtom)
iEnd = iAtomStart(iAtom+1) - 1
square(iStart:iEnd, iEnd+1:mOrb) = transpose(conjg(square(iEnd+1:mOrb, iStart:iEnd)))
end do

if (size(square, dim=1) == 2*mOrb) then
! 2 component matrix
do iAtom = 1, nAtom
iStart = iAtomStart(iAtom) + mOrb
iEnd = iAtomStart(iAtom+1) + mOrb - 1
square(iStart:iEnd, iEnd+1:) = transpose(conjg(square(iEnd+1:, iStart:iEnd)))
end do
end if

end subroutine blockHermitianHS_cmplx


!> Symmetrize a square matrix leaving the on-site atomic blocks alone. (Real version)
subroutine blockSymmetrizeHS_real(square, iAtomStart)

!> Square form matrix.
real(dp), intent(inout) :: square(:,:)

!> Returns the offset array for each atom.
integer, intent(in) :: iAtomStart(:)

integer :: nAtom, iAtom, iStart, iEnd, mOrb

nAtom = size(iAtomStart) - 1
mOrb = iAtomStart(nAtom+1) - 1

@:ASSERT(nAtom > 0)
@:ASSERT(size(square, dim=1) == size(square, dim=2))
@:ASSERT(size(square, dim=1) == mOrb)

do iAtom = 1, nAtom
iStart = iAtomStart(iAtom)
iEnd = iAtomStart(iAtom+1) - 1
square(iStart:iEnd, iEnd+1:mOrb) = transpose(square(iEnd+1:mOrb, iStart:iEnd))
end do

end subroutine blockSymmetrizeHS_real


!> Copy lower triangle to upper for a square matrix. (Real version)
subroutine symmetrizeHS_real(matrix)

!> matrix to symmetrize
real(dp), intent(inout) :: matrix(:,:)
integer :: ii, matSize

matSize = size(matrix, dim=1)
do ii = 1, matSize - 1
matrix(ii, ii + 1 : matSize) = matrix(ii + 1 : matSize, ii)
end do

end subroutine symmetrizeHS_real


!> Copy lower triangle to upper for a square matrix.
subroutine hermitianSquareMatrix(matrix)

!> matrix to symmetrize
complex(dp), intent(inout) :: matrix(:,:)
integer :: ii, matSize

matSize = size(matrix, dim = 1)
do ii = 1, matSize - 1
matrix(ii, ii + 1 : matSize) = conjg(matrix(ii + 1 : matSize, ii))
end do

end subroutine hermitianSquareMatrix


#:if WITH_SCALAPACK
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Scalapack routines
Expand Down
Loading

0 comments on commit 7530d37

Please sign in to comment.