Skip to content

Commit

Permalink
Non-collinear hybrids for molecules
Browse files Browse the repository at this point in the history
  • Loading branch information
@bhourahine
bhourahine committed Jun 4, 2024
1 parent 1a31e8a commit 7bd0ab8
Show file tree
Hide file tree
Showing 11 changed files with 686 additions and 107 deletions.
148 changes: 93 additions & 55 deletions src/dftbp/dftb/hybridxc.F90
View file

Large diffs are not rendered by default.

134 changes: 129 additions & 5 deletions src/dftbp/dftb/populations.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,11 +30,12 @@ module dftbp_dftb_populations
implicit none

private
public :: mulliken, skewMulliken, denseMullikenReal
public :: mulliken, skewMulliken, denseMullikenReal, denseMullikenPauli
public :: getChargePerShell, denseBlockMulliken
public :: getOnsitePopulation, getAtomicMultipolePopulation
public :: denseSubtractDensityOfAtomsReal, denseSubtractDensityOfAtomsCmplxNonperiodic,&
& denseSubtractDensityOfAtomsCmplxPeriodic, denseSubtractDensityOfAtomsCmplxPeriodicGlobal,&
public :: denseSubtractDensityOfAtomsReal, denseSubtractDensityOfAtomsPauli,&
& denseSubtractDensityOfAtomsCmplxNonperiodic, denseSubtractDensityOfAtomsCmplxPeriodic,&
& denseSubtractDensityOfAtomsCmplxPeriodicGlobal,&
& denseSubtractDensityOfAtomsNospinRealNonperiodicReks,&
& denseSubtractDensityOfAtomsSpinRealNonperiodicReks
#:if WITH_SCALAPACK
Expand Down Expand Up @@ -558,6 +559,81 @@ subroutine denseMullikenReal(rhoSqr, overSqr, iSquare, qq)
end subroutine denseMullikenReal


!> Mulliken analysis with dense, square, real-space two-component matrices.
subroutine denseMullikenPauli(rhoSqr, overSqr, iSquare, qq)

!> Square (lower triangular) two component density matrix
complex(dp), intent(in) :: rhoSqr(:,:,:)

!> Square (lower triangular) overlap matrix
real(dp), intent(in) :: overSqr(:,:)

!> Atom positions in the row/column of square matrices
integer, intent(in) :: iSquare(:)

!> Mulliken charges on output (mOrb, nAtom, nSpin)
real(dp), intent(out) :: qq(:,:,:)

real(dp) :: tmpSqr(size(qq, dim=1), size(qq, dim=1))
integer :: iAtom1, iAtom2, iStart1, iEnd1, iStart2, iEnd2, nAtom, nSpin, nOrb1, nOrb2, nOrb

nAtom = size(qq, dim=2)
nOrb = size(overSqr, dim=1)
qq(:,:,:) = 0.0_dp

do iAtom1 = 1, nAtom
iStart1 = iSquare(iAtom1)
iEnd1 = iSquare(iAtom1 + 1) - 1
nOrb1 = iEnd1 - iStart1 + 1
do iAtom2 = iAtom1, nAtom
iStart2 = iSquare(iAtom2)
iEnd2 = iSquare(iAtom2 + 1) - 1
nOrb2 = iEnd2 - iStart2 + 1

! 1 matrix
tmpSqr(1:nOrb2, 1:nOrb1) = &
& real(rhoSqr(iStart2:iEnd2, iStart1:iEnd1, 1)&
& + rhoSqr(iStart2+nOrb:iEnd2+nOrb, iStart1+nOrb:iEnd1+nOrb, 1),dp)&
& * overSqr(iStart2:iEnd2, iStart1:iEnd1)
qq(1:nOrb1, iAtom1, 1) = qq(1:nOrb1, iAtom1, 1) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=1)
if (iAtom1 /= iAtom2) then
qq(1:nOrb2, iAtom2, 1) = qq(1:nOrb2, iAtom2, 1) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=2)
end if

! x matrix
tmpSqr(1:nOrb2, 1:nOrb1) = 2.0_dp * &
& real(rhoSqr(iStart2+nOrb:iEnd2+nOrb, iStart1:iEnd1, 1), dp)&
& * overSqr(iStart2:iEnd2, iStart1:iEnd1)
qq(1:nOrb1, iAtom1, 2) = qq(1:nOrb1, iAtom1, 2) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=1)
if (iAtom1 /= iAtom2) then
qq(1:nOrb2, iAtom2, 2) = qq(1:nOrb2, iAtom2, 2) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=2)
end if

! y matrix
tmpSqr(1:nOrb2, 1:nOrb1) = 2.0_dp * &
& aimag(rhoSqr(iStart2+nOrb:iEnd2+nOrb, iStart1:iEnd1, 1))&
& * overSqr(iStart2:iEnd2, iStart1:iEnd1)
qq(1:nOrb1, iAtom1, 3) = qq(1:nOrb1, iAtom1, 3) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=1)
if (iAtom1 /= iAtom2) then
qq(1:nOrb2, iAtom2, 3) = qq(1:nOrb2, iAtom2, 3) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=2)
end if

! z matrix
tmpSqr(1:nOrb2, 1:nOrb1) = &
& real(rhoSqr(iStart2:iEnd2, iStart1:iEnd1, 1)&
& - rhoSqr(iStart2+nOrb:iEnd2+nOrb, iStart1+nOrb:iEnd1+nOrb, 1),dp)&
& * overSqr(iStart2:iEnd2, iStart1:iEnd1)
qq(1:nOrb1, iAtom1, 4) = qq(1:nOrb1, iAtom1, 4) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=1)
if (iAtom1 /= iAtom2) then
qq(1:nOrb2, iAtom2, 4) = qq(1:nOrb2, iAtom2, 4) + sum(tmpSqr(1:nOrb2, 1:nOrb1), dim=2)
end if

end do
end do

end subroutine denseMullikenPauli


!> Returns pre-factor of Mulliken populations, depending on the number of spin-channels.
pure function populationScalingFactor(nSpin) result(scale)

Expand All @@ -567,20 +643,23 @@ pure function populationScalingFactor(nSpin) result(scale)
!> Pre-factor of Mulliken populations
real(dp) :: scale

! distinguish cases: spin-unpolarized, colinear-spin
! distinguish cases: spin-unpolarized, colinear-spin, non-collinear/Pauli
select case(nSpin)
case(1)
scale = 1.0_dp
case(2)
scale = 0.5_dp
case(4)
scale = 0.5_dp
end select

end function populationScalingFactor


!> Subtracts superposition of atomic densities from dense, square, real-space density matrix.
!!
!! For spin-polarized calculations, q0 is distributed equally to alpha and beta density matrices.
!! For spin-polarized calculations, q0 is distributed equally to alpha and beta density matrices
!! for collinear spin.
!! Note: For periodic systems, q0 is normalized by the overlap that includes periodic images.
subroutine denseSubtractDensityOfAtomsReal(q0, iSquare, overSqr, rho)

Expand Down Expand Up @@ -618,6 +697,51 @@ subroutine denseSubtractDensityOfAtomsReal(q0, iSquare, overSqr, rho)
end subroutine denseSubtractDensityOfAtomsReal


!> Subtracts superposition of atomic densities from dense, square, real-space density matrix.
!!
!! For non-collinear spin-polarized calculations, q0 is distributed equally over the identity
!! spinor for Pauli matrices.
!! Note: For periodic systems, q0 is normalized by the overlap that includes periodic images.
subroutine denseSubtractDensityOfAtomsPauli(q0, iSquare, overSqr, rho)

!> Reference atom populations
real(dp), intent(in) :: q0(:,:,:)

!> Atom positions in the row/column of square matrices
integer, intent(in) :: iSquare(:)

!> Square (lower triangular) overlap matrix, note that for Pauli matrices this factor is
!> multiplied by the identity spinor, so is half the size of rho
real(dp), intent(in) :: overSqr(:,:)

!> Spin polarized (lower triangular) density matrix
complex(dp), intent(inout) :: rho(:,:,:)

integer :: nAtom, iAtom, nSpin, iStart, iEnd, iOrb, iBlock, offset
real(dp) :: scale

nAtom = size(iSquare) - 1
nSpin = size(q0, dim=3)
@:ASSERT(nSpin == 4)

scale = populationScalingFactor(nSpin)

do iBlock = 0, 1
offset = iBlock * (iSquare(nAtom + 1) - 1)
do iAtom = 1, nAtom
iStart = iSquare(iAtom)
iEnd = iSquare(iAtom + 1) - 1
do iOrb = 1, iEnd - iStart + 1
rho(iStart + iOrb - 1 + offset, iStart + iOrb -1 + offset, 1) =&
& rho(iStart + iOrb - 1 + offset, iStart + iOrb - 1 + offset, 1)&
& - scale * q0(iOrb, iAtom, 1) / overSqr(iStart+iOrb-1, iStart+iOrb-1)
end do
end do
end do

end subroutine denseSubtractDensityOfAtomsPauli


!> Subtracts superposition of atomic densities from dense, square, complex-valued density matrix.
!!
!! Note: q0 is equally distributed to alpha and beta density matrices.
Expand Down
43 changes: 33 additions & 10 deletions src/dftbp/dftbplus/initprogram.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1839,7 +1839,8 @@ subroutine initProgramVariables(this, input, env)
iMixer = input%ctrl%iMixSwitch
nGeneration = input%ctrl%iGenerations
mixParam = input%ctrl%almix
tCmplxMixer = (.not. this%tRealHS) .and. (this%hybridXcAlg == hybridXcAlgo%matrixBased)
tCmplxMixer = (.not. this%tRealHS) .and. (this%hybridXcAlg == hybridXcAlgo%matrixBased)&
& .or. (this%t2Component .and. this%isHybridXc)
if (tCmplxMixer) then
allocate(this%pChrgMixerCmplx)
select case (iMixer)
Expand Down Expand Up @@ -2908,7 +2909,7 @@ subroutine initProgramVariables(this, input, env)
call THybridXcFunc_init(this%hybridXc, this%nAtom, this%species0, hubbU(1, :),&
& input%ctrl%hybridXcInp%screeningThreshold, input%ctrl%hybridXcInp%omega,&
& input%ctrl%hybridXcInp%camAlpha, input%ctrl%hybridXcInp%camBeta,&
& this%tSpin, allocated(this%reks), input%ctrl%hybridXcInp%hybridXcAlg,&
& this%tSpin, this%nSpin, allocated(this%reks), input%ctrl%hybridXcInp%hybridXcAlg,&
& input%ctrl%hybridXcInp%hybridXcType, input%ctrl%hybridXcInp%gammaType, this%tPeriodic,&
& this%tRealHS, errStatus, coeffsDiag=this%supercellFoldingDiag,&
& gammaCutoff=this%cutOff%gammaCutoff,&
Expand All @@ -2923,7 +2924,11 @@ subroutine initProgramVariables(this, input, env)
& size(this%parallelKS%localKS, dim=2))
! reset number of mixer elements, so that there is enough space for density matrices
if (this%tRealHS) then
this%nMixElements = size(this%densityMatrix%deltaRhoIn)
if (this%t2Component) then
this%nMixElements = size(this%densityMatrix%deltaRhoInCplx)
else
this%nMixElements = size(this%densityMatrix%deltaRhoIn)
end if
else
if (input%ctrl%hybridXcInp%hybridXcAlg == hybridXcAlgo%matrixBased) then
this%nMixElements = size(this%densityMatrix%deltaRhoInCplx)
Expand Down Expand Up @@ -5853,17 +5858,17 @@ subroutine ensureHybridXcReqs(this, tShellResolved, hybridXcInp)
& moment.")
end if

if (this%nSpin > 2) then
call error("Hybrid calculations not implemented for non-colinear calculations.")
if (this%nSpin > 2 .and. .not. this%tRealHS) then
call error("Hybrid calculations not implemented for non-colinear calculations in this case.")
end if

if ((.not. this%tRealHS) .and. this%nSpin == 2&
if ((.not. this%tRealHS) .and. this%nSpin > 1&
& .and. hybridXcInp%gammaType /= hybridXcGammaTypes%truncated) then
call error("Hybrid functionality does not yet support spin-polarized calculations of periodic&
& systems beyond the Gamma-point for CoulombMatrix settings other than 'Truncated'.")
end if

if ((.not. this%tRealHS) .and. this%nSpin == 2 .and. this%tForces) then
if ((.not. this%tRealHS) .and. this%nSpin > 1 .and. this%tForces) then
call error("Hybrid functionality currently does not yet support spin-polarized gradient&
& evaluation for periodic systems beyond the Gamma-point.")
end if
Expand All @@ -5880,6 +5885,15 @@ subroutine ensureHybridXcReqs(this, tShellResolved, hybridXcInp)
call error("Hybrid calculations not currently implemented for 3rd-order DFTB.")
end if

if (this%nSpin == 4 .and. hybridXcInp%hybridXcAlg /= hybridXcAlgo%matrixBased) then
call error("Non-collinear spin only available for matrix alogorithm hybrids at present.")
end if

if (this%nSpin == 4 .and. this%boundaryCond%iBoundaryCondition /= boundaryConditions%cluster)&
& then
call error("Non-collinear spin only available for hybrids with molecular systems at present.")
end if

if (this%isHybLinResp .and. hybridXcInp%hybridXcType == hybridXcFunc%cam) then
call error("General CAM functionals not currently implemented for linear response.")
end if
Expand Down Expand Up @@ -6206,10 +6220,19 @@ subroutine reallocateHybridXc(this, hybridXcAlg, nLocalRows, nLocalCols, nLocalK

if (this%tRealHS) then
! Prevent for deleting charges read in from file
if (.not. allocated(this%densityMatrix%deltaRhoIn)) then
allocate(this%densityMatrix%deltaRhoIn(nLocalRows, nLocalCols, nLocalKS), source=0.0_dp)
if (this%t2Component) then
if (.not. allocated(this%densityMatrix%deltaRhoInCplx)) then
allocate(this%densityMatrix%deltaRhoInCplx(nLocalRows, nLocalCols, nLocalKS),&
& source=(0.0_dp,0.0_dp))
end if
allocate(this%densityMatrix%deltaRhoOutCplx(nLocalRows, nLocalCols, nLocalKS),&
& source=(0.0_dp,0.0_dp))
else
if (.not. allocated(this%densityMatrix%deltaRhoIn)) then
allocate(this%densityMatrix%deltaRhoIn(nLocalRows, nLocalCols, nLocalKS), source=0.0_dp)
end if
allocate(this%densityMatrix%deltaRhoOut(nLocalRows, nLocalCols, nLocalKS), source=0.0_dp)
end if
allocate(this%densityMatrix%deltaRhoOut(nLocalRows, nLocalCols, nLocalKS), source=0.0_dp)
elseif (this%tReadChrg .and. (.not. allocated(this%supercellFoldingDiag))) then
! in case of k-points and restart from file, we have to wait until charges.bin was read
if (hybridXcAlg == hybridXcAlgo%matrixBased) then
Expand Down
Loading

0 comments on commit 7bd0ab8

Please sign in to comment.