Skip to content

Commit

Permalink
Move overlap derivative folding (dftbplus#1470)
Browse files Browse the repository at this point in the history 
There are already similar routines in sparse2dense.

Also introduce place-holder citation for k-point based routine:
van der Heide, Unpublished (2024).
Replace this once there is a proper publication/thesis.
  • Loading branch information
@bhourahine
bhourahine authored May 23, 2024
1 parent 292f83e commit 4e4c245
Show file tree
Hide file tree
Showing 3 changed files with 266 additions and 257 deletions.
19 changes: 18 additions & 1 deletion src/dftbp/dftb/dense.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module dftbp_dftb_dense
implicit none

private
public :: buildSquaredAtomIndex
public :: buildSquaredAtomIndex, getDescriptor

contains

Expand Down Expand Up @@ -43,4 +43,21 @@ subroutine buildSquaredAtomIndex(iAtomStart, orb)

end subroutine buildSquaredAtomIndex


!> Finds location of relevant atomic block indices in a dense matrix.
pure function getDescriptor(iAt, iSquare) result(desc)

!> Relevant atom
integer, intent(in) :: iAt

!> Indexing array for start of atom orbitals
integer, intent(in) :: iSquare(:)

!> Resulting location ranges
integer :: desc(3)

desc(:) = [iSquare(iAt), iSquare(iAt + 1) - 1, iSquare(iAt + 1) - iSquare(iAt)]

end function getDescriptor

end module dftbp_dftb_dense
258 changes: 3 additions & 255 deletions src/dftbp/dftb/hybridxc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ module dftbp_dftb_hybridxc
use dftbp_common_environment, only : TEnvironment, globalTimers
use dftbp_common_schedule, only : getStartAndEndIndex
use dftbp_common_status, only : TStatus
use dftbp_dftb_dense, only : getDescriptor
use dftbp_dftb_densitymatrix, only : TDensityMatrix
use dftbp_dftb_nonscc, only : TNonSccDiff, buildS
use dftbp_dftb_periodic, only : TNeighbourList, TSymNeighbourList, getCellTranslations, cart2frac
Expand All @@ -24,7 +25,8 @@ module dftbp_dftb_hybridxc
& getdHfAnalyticalGammaValue_workhorse, getdLrAnalyticalGammaValue_workhorse,&
& getddLrNumericalGammaValue_workhorse
use dftbp_dftb_slakocont, only : TSlakoCont
use dftbp_dftb_sparse2dense, only : unpackHS
use dftbp_dftb_sparse2dense, only : unpackHS, getUnpackedOverlapPrime_real,&
& getUnpackedOverlapPrime_kpts
use dftbp_math_blasroutines, only : gemm, symm, hemm
use dftbp_math_matrixops, only : adjointLowerTriangle
use dftbp_math_simplealgebra, only : determinant33
Expand Down Expand Up @@ -3530,23 +3532,6 @@ subroutine getHybridEnergy_kpts(this, env, localKS, iKiSToiGlobalKS, kWeights, d
end subroutine getHybridEnergy_kpts


!> Finds location of relevant atomic block indices in a dense matrix.
pure function getDescriptor(iAt, iSquare) result(desc)

!> Relevant atom
integer, intent(in) :: iAt

!> Indexing array for start of atom orbitals
integer, intent(in) :: iSquare(:)

!> Resulting location ranges
integer :: desc(3)

desc(:) = [iSquare(iAt), iSquare(iAt + 1) - 1, iSquare(iAt + 1) - iSquare(iAt)]

end function getDescriptor


!> Evaluates energy from Hamiltonian and density matrix (real version).
pure function evaluateEnergy_real(hamiltonian, densityMat) result(energy)

Expand Down Expand Up @@ -3942,95 +3927,6 @@ function getCamGammaGSum(this, iAt1, iAt2, iSp1, iSp2, rCellVecsG) result(gamma)
end function getCamGammaGSum


!> Calculates the derivative of the square, dense, unpacked, dual-space overlap matrix w.r.t.
!! given atom.
subroutine getUnpackedOverlapPrime_kpts(iAtomPrime, skOverCont, orb, derivator, symNeighbourList,&
& nNeighbourCamSym, iSquare, cellVec, rCoords, kPoint, overSqrPrime)

!> Overlap derivative calculated w.r.t. this atom in the central cell
integer, intent(in) :: iAtomPrime

!> Sparse overlap container
type(TSlakoCont), intent(in) :: skOverCont

!> Orbital information.
type(TOrbitals), intent(in) :: orb

!> Differentiation object
class(TNonSccDiff), intent(in) :: derivator

!> List of neighbours for each atom (symmetric version)
type(TSymNeighbourList), intent(in) :: symNeighbourList

!> Nr. of neighbours for each atom.
integer, intent(in) :: nNeighbourCamSym(:)

!> Position of each atom in the rows/columns of the square matrices. Shape: (nAtom)
integer, intent(in) :: iSquare(:)

!> Vectors to neighboring unit cells in relative units
real(dp), intent(in) :: cellVec(:,:)

!> Atomic coordinates in absolute units, potentially including periodic images
real(dp), intent(in) :: rCoords(:,:)

!> Relative coordinates of the k-point where the overlap should be constructed
real(dp), intent(in) :: kPoint(:)

!> Overlap derivative
complex(dp), intent(out) :: overSqrPrime(:,:,:)

!! Temporary overlap derivative storage
real(dp) :: overPrime(orb%mOrb, orb%mOrb, 3)

!! Dense matrix descriptor indices
integer, parameter :: descLen = 3, iStart = 1, iEnd = 2, iNOrb = 3

!! Stores start/end index and number of orbitals of square matrices
integer :: descAt1(descLen), descAt2(descLen)

!! Atom to calculate energy gradient components for
integer :: iAt2, iAt2fold, iNeigh

!! Iterates over coordinates
integer :: iCoord

!! The k-point in relative coordinates, multiplied by 2pi
real(dp) :: kPoint2p(3)

!! Phase factor
complex(dp) :: phase

!! Auxiliary variables
integer :: iVec

overSqrPrime(:,:,:) = (0.0_dp, 0.0_dp)

kPoint2p(:) = 2.0_dp * pi * kPoint

descAt1 = getDescriptor(iAtomPrime, iSquare)
do iNeigh = 0, nNeighbourCamSym(iAtomPrime)
iAt2 = symNeighbourList%neighbourList%iNeighbour(iNeigh, iAtomPrime)
iAt2fold = symNeighbourList%img2CentCell(iAt2)
if (iAtomPrime == iAt2fold) cycle
descAt2 = getDescriptor(iAt2fold, iSquare)
iVec = symNeighbourList%iCellVec(iAt2)
overPrime(:,:,:) = 0.0_dp
call derivator%getFirstDeriv(overPrime, skOverCont, rCoords, symNeighbourList%species,&
& iAtomPrime, iAt2, orb)
phase = exp(imag * dot_product(kPoint2p, cellVec(:, iVec)))
do iCoord = 1, 3
overSqrPrime(iCoord, descAt2(iStart):descAt2(iStart) + descAt2(iNOrb) - 1,&
& descAt1(iStart):descAt1(iStart) + descAt1(iNOrb) - 1)&
& = overSqrPrime(iCoord, descAt2(iStart):descAt2(iStart) + descAt2(iNOrb) - 1,&
& descAt1(iStart):descAt1(iStart) + descAt1(iNOrb) - 1)&
& + overPrime(1:descAt2(iNOrb), 1:descAt1(iNOrb), iCoord) * phase
end do
end do

end subroutine getUnpackedOverlapPrime_kpts


!> Calculates pseudo Fourier transform of square CAM y-matrix with shape [nOrb, nOrb].
subroutine getCamGammaFourierAO(this, iSquare, cellVecsG, rCellVecsG, kPoint, kPointPrime,&
& camGammaSqr)
Expand Down Expand Up @@ -5799,82 +5695,6 @@ subroutine getUnpackedCamGammaAOPrime(env, denseDesc, orb, iAtomPrime, camdGamma

end subroutine getUnpackedCamGammaAOPrime


!> Calculates the derivative of the square, dense, unpacked, Gamma-point overlap matrix w.r.t.
!! given atom.
!!
!! 1st term of Eq.(B5) of Phys. Rev. Materials 7, 063802 (DOI: 10.1103/PhysRevMaterials.7.063802)
subroutine getUnpackedOverlapPrime_real(env, denseDesc, iAtomPrime, skOverCont, orb, derivator,&
& symNeighbourList, nNeighbourCamSym, iSquare, rCoords, overSqrPrime)

!> Environment settings
type(TEnvironment), intent(in) :: env

!> Dense matrix descriptor
type(TDenseDescr), intent(in) :: denseDesc

!> Overlap derivative calculated w.r.t. this atom in the central cell
integer, intent(in) :: iAtomPrime

!> Sparse overlap container
type(TSlakoCont), intent(in) :: skOverCont

!> Orbital information.
type(TOrbitals), intent(in) :: orb

!> Differentiation object
class(TNonSccDiff), intent(in) :: derivator

!> List of neighbours for each atom (symmetric version)
type(TSymNeighbourList), intent(in) :: symNeighbourList

!> Nr. of neighbours for each atom.
integer, intent(in) :: nNeighbourCamSym(:)

!> Position of each atom in the rows/columns of the square matrices. Shape: (nAtom)
integer, intent(in) :: iSquare(:)

!> Atomic coordinates in absolute units, potentially including periodic images
real(dp), intent(in) :: rCoords(:,:)

!> Overlap derivative
real(dp), intent(out) :: overSqrPrime(:,:,:)

!! Temporary overlap derivative storage
real(dp) :: overPrime(orb%mOrb, orb%mOrb, 3)

!! Dense matrix descriptor indices
integer, parameter :: descLen = 3, iStart = 1, iEnd = 2, iNOrb = 3

!! Stores start/end index and number of orbitals of square matrices
integer :: descAt1(descLen), descAt2(descLen)

!! Atom to calculate energy gradient components for
integer :: iAt2, iAt2fold, iNeigh

!! Iterates over coordinates
integer :: iCoord

overSqrPrime(:,:,:) = 0.0_dp

descAt1 = getDescriptor(iAtomPrime, iSquare)
do iNeigh = 0, nNeighbourCamSym(iAtomPrime)
iAt2 = symNeighbourList%neighbourList%iNeighbour(iNeigh, iAtomPrime)
iAt2fold = symNeighbourList%img2CentCell(iAt2)
if (iAtomPrime == iAt2fold) cycle
descAt2 = getDescriptor(iAt2fold, iSquare)
overPrime(:,:,:) = 0.0_dp
call derivator%getFirstDeriv(overPrime, skOverCont, rCoords, symNeighbourList%species,&
& iAtomPrime, iAt2, orb)
do iCoord = 1, 3
call scalafx_addl2g(env%blacs%orbitalGrid, overPrime(1:descAt2(iNOrb), 1:descAt1(iNOrb),&
& iCoord), denseDesc%blacsOrbSqr, descAt2(iStart), descAt1(iStart),&
& overSqrPrime(iCoord, :,:))
end do
end do

end subroutine getUnpackedOverlapPrime_real

#:else

!> Adds CAM gradients due to CAM range-separated contributions, using a matrix-based formulation.
Expand Down Expand Up @@ -6079,78 +5899,6 @@ subroutine getUnpackedCamGammaAOPrime(iAtomPrime, camdGammaEval0, iSquare, camdG

end subroutine getUnpackedCamGammaAOPrime


!> Calculates the derivative of the square, dense, unpacked, Gamma-point overlap matrix w.r.t.
!! given atom.
!!
!! 1st term of Eq.(B5) of Phys. Rev. Materials 7, 063802 (DOI: 10.1103/PhysRevMaterials.7.063802)
subroutine getUnpackedOverlapPrime_real(iAtomPrime, skOverCont, orb, derivator, symNeighbourList,&
& nNeighbourCamSym, iSquare, rCoords, overSqrPrime)

!> Overlap derivative calculated w.r.t. this atom in the central cell
integer, intent(in) :: iAtomPrime

!> Sparse overlap container
type(TSlakoCont), intent(in) :: skOverCont

!> Orbital information.
type(TOrbitals), intent(in) :: orb

!> Differentiation object
class(TNonSccDiff), intent(in) :: derivator

!> List of neighbours for each atom (symmetric version)
type(TSymNeighbourList), intent(in) :: symNeighbourList

!> Nr. of neighbours for each atom.
integer, intent(in) :: nNeighbourCamSym(:)

!> Position of each atom in the rows/columns of the square matrices. Shape: (nAtom)
integer, intent(in) :: iSquare(:)

!> Atomic coordinates in absolute units, potentially including periodic images
real(dp), intent(in) :: rCoords(:,:)

!> Overlap derivative
real(dp), intent(out) :: overSqrPrime(:,:,:)

!! Temporary overlap derivative storage
real(dp) :: overPrime(orb%mOrb, orb%mOrb, 3)

!! Dense matrix descriptor indices
integer, parameter :: descLen = 3, iStart = 1, iEnd = 2, iNOrb = 3

!! Stores start/end index and number of orbitals of square matrices
integer :: descAt1(descLen), descAt2(descLen)

!! Atom to calculate energy gradient components for
integer :: iAt2, iAt2fold, iNeigh

!! Iterates over coordinates
integer :: iCoord

overSqrPrime(:,:,:) = 0.0_dp

descAt1 = getDescriptor(iAtomPrime, iSquare)
do iNeigh = 0, nNeighbourCamSym(iAtomPrime)
iAt2 = symNeighbourList%neighbourList%iNeighbour(iNeigh, iAtomPrime)
iAt2fold = symNeighbourList%img2CentCell(iAt2)
if (iAtomPrime == iAt2fold) cycle
descAt2 = getDescriptor(iAt2fold, iSquare)
overPrime(:,:,:) = 0.0_dp
call derivator%getFirstDeriv(overPrime, skOverCont, rCoords, symNeighbourList%species,&
& iAtomPrime, iAt2, orb)
do iCoord = 1, 3
overSqrPrime(iCoord, descAt2(iStart):descAt2(iStart) + descAt2(iNOrb) - 1,&
& descAt1(iStart):descAt1(iStart) + descAt1(iNOrb) - 1)&
& = overSqrPrime(iCoord, descAt2(iStart):descAt2(iStart) + descAt2(iNOrb) - 1,&
& descAt1(iStart):descAt1(iStart) + descAt1(iNOrb) - 1)&
& + overPrime(1:descAt2(iNOrb), 1:descAt1(iNOrb), iCoord)
end do
end do

end subroutine getUnpackedOverlapPrime_real

#:endif


Expand Down
Loading

0 comments on commit 4e4c245

Please sign in to comment.