Non-collinear hybrids for molecules

H atom test cases

src/dftbp/dftb/populations.F90

@@ -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
@@ -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) = &
+            & 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) = &
+            & 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)
 
@@ -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)
 
@@ -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.

src/dftbp/dftbplus/initprogram.F90

@@ -2897,7 +2897,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,&
@@ -2912,7 +2912,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)
@@ -5842,17 +5846,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
@@ -5869,6 +5873,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
@@ -6195,10 +6208,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