Partial clean up of hybrid requirement checks

* Matrix diagonal function for 3x3 matrix
* Error return propagation fixed in init (bug fix)
* Error message clarification and clean up
* Minor logic clean-up in init
* Cutoff adjustment initialised properly if it had been reduced in
  the restart file (assumes this was done for a reason in the
  previous run).

src/dftbp/dftb/hybridxc.F90

@@ -668,7 +668,7 @@ end subroutine THybridXcFunc_init
 
 
   !> Checks if obtained supercell folding matrix meets current requirements.
-  subroutine checkSupercellFoldingMatrix(supercellFoldingMatrix, errStatus, supercellFoldingDiagOut)
+  subroutine checkSupercellFoldingMatrix(supercellFoldingMatrix, errStatus)
 
     !> Coefficients of the lattice vectors in the linear combination for the super lattice vectors
     !! (should be integer values) and shift of the grid along the three small reciprocal lattice
@@ -678,9 +678,6 @@ subroutine checkSupercellFoldingMatrix(supercellFoldingMatrix, errStatus, superc
     !> Error status
     type(TStatus), intent(inout) :: errStatus
 
-    !> Diagonal elements of supercell folding matrix, if present
-    integer, intent(out), optional :: supercellFoldingDiagOut(:)
-
     !! Supercell folding coefficients and shifts
     real(dp), pointer :: coeffs(:,:), shifts(:)
 
@@ -690,10 +687,6 @@ subroutine checkSupercellFoldingMatrix(supercellFoldingMatrix, errStatus, superc
     !! Auxiliary variables
     integer :: ii, jj
 
-    if (present(supercellFoldingDiagOut)) then
-      @:ASSERT(size(supercellFoldingDiagOut) == 3)
-    end if
-
     coeffs => supercellFoldingMatrix(:, 1:3)
     shifts => supercellFoldingMatrix(:, 4)
 
@@ -717,25 +710,16 @@ subroutine checkSupercellFoldingMatrix(supercellFoldingMatrix, errStatus, superc
       end do
     end do lpOuter
     if (tNotMonkhorstPack) then
-      @:RAISE_ERROR(errStatus, -1, "Range-separated calculations with k-points require a&
-          & Monkhorst-Pack-like sampling, i.e. a uniform extension of the lattice.")
+      @:RAISE_ERROR(errStatus, -1, "Hybrid functionals using integration with k-points requires a&
+          & Monkhorst-Pack-like sampling, i.e., a uniform extension of the lattice.")
     end if
 
     ! Check if shifts are zero
     if (any(abs(shifts) > 1e-06_dp)) then
-      @:RAISE_ERROR(errStatus, -1, "Range-separated calculations with k-points require a&
+      @:RAISE_ERROR(errStatus, -1, "Hybrid functionals using integration with k-points requires a&
           & Monkhorst-Pack-like sampling with zero shift.")
     end if
 
-    ! All checks have passed, continue...
-
-    ! Get diagonal elements as integers, if requested
-    if (present(supercellFoldingDiagOut)) then
-      do ii = 1, 3
-        supercellFoldingDiagOut(ii) = nint(coeffs(ii, ii))
-      end do
-    end if
-
   end subroutine checkSupercellFoldingMatrix
 
 

src/dftbp/dftb/sccinit.F90

@@ -6,6 +6,7 @@
 !--------------------------------------------------------------------------------------------------!
 
 #:include 'common.fypp'
+#:include 'error.fypp'
 
 !> Module for initializing SCC part of the calculation.
 module dftbp_dftb_sccinit
@@ -17,6 +18,7 @@ module dftbp_dftb_sccinit
   use dftbp_dftb_hybridxc, only : checkSupercellFoldingMatrix, hybridXcAlgo
   use dftbp_dftb_periodic, only : getSuperSampling
   use dftbp_io_message, only : error
+  use dftbp_math_simplealgebra, only : diagonal
   use dftbp_type_commontypes, only : TOrbitals
   use dftbp_type_multipole, only : TMultipole
   implicit none
@@ -211,7 +213,7 @@ subroutine initQFromFile(qq, fileName, tReadAscii, orb, qBlock, qiBlock, density
     logical, intent(in) :: tRealHS
 
     !> Error status
-    type(TStatus), intent(inout) :: errStatus
+    type(TStatus), intent(out) :: errStatus
 
     !> Magnetisation checksum for regular spin polarization total magnetic moment
     real(dp), intent(in), optional :: magnetisation
@@ -249,10 +251,10 @@ subroutine initQFromFile(qq, fileName, tReadAscii, orb, qBlock, qiBlock, density
     integer :: fileFormat
     real(dp) :: sumQ
 
-    !! Requested to be re-loaded
+    !! Quantities requested to be re-loaded from the file
     logical :: tBlock, tiBlock, tRho, tKpointInfo, isMultipolar
 
-    !! Present in the file itself
+    !! Are these present in the file itself
     logical :: tBlockPresent, tiBlockPresent, tRhoPresent, tKpointInfoPresent
 
     character(len=120) :: error_string
@@ -555,8 +557,9 @@ subroutine initQFromFile(qq, fileName, tReadAscii, orb, qBlock, qiBlock, density
           else
             read(file%unit, iostat=iErr) coeffsAndShifts
           end if
-          call checkSupercellFoldingMatrix(coeffsAndShifts, errStatus,&
-              & supercellFoldingDiagOut=supercellFoldingDiag)
+          call checkSupercellFoldingMatrix(coeffsAndShifts, errStatus)
+          @:PROPAGATE_ERROR(errStatus)
+          supercellFoldingDiag = nint(diagonal(coeffsAndShifts(:,:3)))
           if (hybridXcAlg == hybridXcAlgo%matrixBased) then
             call getSuperSampling(coeffsAndShifts(:,1:3), modulo(coeffsAndShifts(:,4), 1.0_dp),&
                 & densityMatrix%kPointPrime, densityMatrix%kWeightPrime, reduceByInversion=.true.)

src/dftbp/dftbplus/initprogram.F90

@@ -97,7 +97,7 @@ module dftbp_dftbplus_initprogram
   use dftbp_math_lapackroutines, only : matinv
   use dftbp_math_randomgenpool, only : TRandomGenPool, init
   use dftbp_math_ranlux, only : TRanlux, getRandom
-  use dftbp_math_simplealgebra, only : determinant33
+  use dftbp_math_simplealgebra, only : determinant33, diagonal
   use dftbp_md_andersentherm, only : TAndersenThermostat, init
   use dftbp_md_berendsentherm, only :TBerendsenThermostat, init
   use dftbp_md_dummytherm, only : TDummyThermostat, init
@@ -1396,9 +1396,7 @@ subroutine initProgramVariables(this, input, env)
         & this%boundaryCond, this%coord0, this%species0, this%tCoordsChanged, this%tLatticeChanged,&
         & this%latVec, this%origin, this%recVec, this%invLatVec, this%cellVol, this%recCellVol,&
         & errStatus)
-    if (errStatus%hasError()) then
-      call error(errStatus%message)
-    end if
+    if (errStatus%hasError()) call error(errStatus%message)
 
     ! Get species names and output file
     this%geoOutFile = input%ctrl%outFile
@@ -1689,9 +1687,7 @@ subroutine initProgramVariables(this, input, env)
     call initTransport_(this, env, input, this%electronicSolver, this%nSpin, this%tempElec,&
         & this%tNegf, this%isAContactCalc, this%mu, this%negfInt, this%ginfo, this%transpar,&
         & this%writeTunn, this%tWriteLDOS, this%regionLabelLDOS, errStatus)
-    if (errStatus%hasError()) then
-      call error(errStatus%message)
-    end if
+    if (errStatus%hasError()) call error(errStatus%message)
   #:else
     this%tTunn = .false.
     this%tLocalCurrents = .false.
@@ -2432,9 +2428,7 @@ subroutine initProgramVariables(this, input, env)
       allocate(this%electrostatPot)
       call TElStatPotentials_init(this%electrostatPot, input%ctrl%elStatPotentialsInp,&
           & allocated(this%eField) .or. this%tExtChrg, this%hamiltonianType, errStatus)
-      if (errStatus%hasError()) then
-        call error(errStatus%message)
-      end if
+      if (errStatus%hasError()) call error(errStatus%message)
     end if
 
     if (allocated(input%ctrl%pipekMezeyInp)) then
@@ -2799,31 +2793,37 @@ subroutine initProgramVariables(this, input, env)
       allocate(this%symNeighbourList%iPair(0, this%nAtom))
       if ((.not. this%tReadChrg) .and. this%tPeriodic) then
         this%supercellFoldingMatrix = input%ctrl%supercellFoldingMatrix
+        call checkSupercellFoldingMatrix(this%supercellFoldingMatrix(:,:3), errStatus)
+        if (errStatus%hasError()) call error(errStatus%message)
         this%supercellFoldingDiag = input%ctrl%supercellFoldingDiag
+        @:ASSERT(all(this%supercellFoldingDiag ==&
+            & nint(diagonal(this%supercellFoldingMatrix(:,:3)))))
       end if
       if (this%isHybridXc) then
         call ensureHybridXcReqs(this, input%ctrl%tShellResolved, input%ctrl%hybridXcInp)
-        if (.not. this%tReadChrg) then
-          if (this%tPeriodic .and. this%tRealHS) then
+        if (this%tPeriodic .and. .not. this%tReadChrg) then
+          if (this%tRealHS) then
             ! Periodic system (Gamma-point only), dense Hamiltonian and overlap are real-valued
             call getHybridXcCutOff_gamma(this%cutOff, input%geom%latVecs,&
-                & input%ctrl%hybridXcInp%cutoffRed,&
+                & input%ctrl%hybridXcInp%cutoffRed, errStatus,&
                 & gSummationCutoff=input%ctrl%hybridXcInp%gSummationCutoff,&
                 & gammaCutoff=input%ctrl%hybridXcInp%gammaCutoff)
-          elseif (.not. this%tRealHS) then
+          else
             ! Dense Hamiltonian and overlap are complex-valued (general k-point case)
             call getHybridXcCutOff_kpts(this%cutOff, input%geom%latVecs,&
-                & input%ctrl%hybridXcInp%cutoffRed, this%supercellFoldingDiag,&
+                & input%ctrl%hybridXcInp%cutoffRed, this%supercellFoldingDiag, errStatus,&
                 & gammaCutoff=input%ctrl%hybridXcInp%gammaCutoff,&
                 & wignerSeitzReduction=input%ctrl%hybridXcInp%wignerSeitzReduction,&
                 & gSummationCutoff=input%ctrl%hybridXcInp%gSummationCutoff)
           end if
+          if (errStatus%hasError()) call error(errStatus%message)
         end if
 
         ! Non-periodic system (cluster)
         if (.not. this%tPeriodic) then
           call getHybridXcCutOff_cluster(this%cutOff, input%ctrl%hybridXcInp%cutoffRed)
         end if
+
       end if
 
     #:if WITH_SCALAPACK
@@ -2872,14 +2872,16 @@ subroutine initProgramVariables(this, input, env)
     ! initial run is only known after invoking this%initializeCharges(). Inferring the Coulomb
     ! truncation cutoff, therefore calling getHybridXcCutoff(), needs this information.
     if (this%isHybridXc .and. this%tReadChrg) then
-      ! First, check if supercell folding matrix is identical to previous run, if specified in input
+
+      ! First, check if supercell folding matrix is identical to the previous run, if it is
+      ! specified in the input
       if (allocated(input%ctrl%supercellFoldingMatrix)) then
         if (any(abs(input%ctrl%supercellFoldingMatrix&
             & - this%supercellFoldingMatrix) > 1e-06_dp)) then
           write(tmpStr, "(A,3I5,A,3I5,A,3I5,A,3F10.6)")&
               & 'Error while processing k-point sampling for hybrid run.'&
               & // NEW_LINE('A')&
-              & // '   When restarting, only identical k-point samplings to previous run are'&
+              & // '   When restarting, only identical k-point samplings to the previous run are'&
               & // NEW_LINE('A') // '   supported. In this case this would correspond to the&
               & following supercell' // NEW_LINE('A') // '   folding matrix:'&
               & // NEW_LINE('A'),&
@@ -2890,19 +2892,24 @@ subroutine initProgramVariables(this, input, env)
           call error(trim(tmpStr))
         end if
       end if
+
       if (this%tPeriodic .and. this%tRealHS) then
         ! Periodic system (Gamma-point only), dense Hamiltonian and overlap are real-valued
         call getHybridXcCutOff_gamma(this%cutOff, input%geom%latVecs,&
-            & input%ctrl%hybridXcInp%cutoffRed,&
+            & input%ctrl%hybridXcInp%cutoffRed, errStatus,&
+            & gSummationCutoff=input%ctrl%hybridXcInp%gSummationCutoff,&
             & gammaCutoff=input%ctrl%hybridXcInp%gammaCutoff)
       elseif (.not. this%tRealHS) then
         ! Dense Hamiltonian and overlap are complex-valued (general k-point case)
         call getHybridXcCutOff_kpts(this%cutOff, input%geom%latVecs,&
-            & input%ctrl%hybridXcInp%cutoffRed, this%supercellFoldingDiag,&
+            & input%ctrl%hybridXcInp%cutoffRed, this%supercellFoldingDiag, errStatus,&
             & gammaCutoff=input%ctrl%hybridXcInp%gammaCutoff,&
             & wignerSeitzReduction=input%ctrl%hybridXcInp%wignerSeitzReduction,&
             & gSummationCutoff=input%ctrl%hybridXcInp%gSummationCutoff)
       end if
+
+      if (errStatus%hasError()) call error(errStatus%message)
+
     end if
 
     if (this%isHybridXc) then
@@ -2916,9 +2923,7 @@ subroutine initProgramVariables(this, input, env)
           & gSummationCutoff=this%cutOff%gSummationCutoff,&
           & wignerSeitzReduction=this%cutOff%wignerSeitzReduction,&
           & latVecs=input%geom%latVecs)
-      if (errStatus%hasError()) then
-        call error(errStatus%message)
-      end if
+      if (errStatus%hasError()) call error(errStatus%message)
       ! now all information is present to properly allocate density matrices and associate pointers
       call reallocateHybridXc(this, input%ctrl%hybridXcInp%hybridXcAlg, nLocalRows, nLocalCols,&
           & size(this%parallelKS%localKS, dim=2))
@@ -4012,9 +4017,7 @@ subroutine initProgramVariables(this, input, env)
           & this%tPeriodic, this%parallelKS, this%tRealHS, this%kPoint, this%kWeight,&
           & this%isHybridXc, this%scc, this%tblite, this%eFieldScaling, this%hamiltonianType,&
           & errStatus)
-      if (errStatus%hasError()) then
-        call error(errStatus%message)
-      end if
+      if (errStatus%hasError()) call error(errStatus%message)
 
     end if
 
@@ -4229,7 +4232,7 @@ subroutine initializeCharges(this, errStatus, initialSpins, initialCharges, hybr
     class(TDftbPlusMain), intent(inout) :: this
 
     !> Error status
-    type(TStatus), intent(inout) :: errStatus
+    type(TStatus), intent(out) :: errStatus
 
     !> Initial spins
     real(dp), optional, intent(in) :: initialSpins(:,:)
@@ -4318,6 +4321,7 @@ subroutine initializeCharges(this, errStatus, initialSpins, initialCharges, hybr
 
     ! Charges read from file
     if (this%tReadChrg) then
+
       if (this%tFixEf .or. this%tSkipChrgChecksum) then
         ! do not check charge or magnetisation from file
         call initQFromFile(this%qInput, fCharges, this%tReadChrgAscii, this%orb, this%qBlockIn,&
@@ -4343,23 +4347,20 @@ subroutine initializeCharges(this, errStatus, initialSpins, initialCharges, hybr
         end if
       end if
 
-      ! Check if obtained supercell folding matrix meets current requirements
-      if (this%isHybridXc .and. this%tPeriodic) then
-        allocate(this%supercellFoldingDiag(3))
-        call checkSupercellFoldingMatrix(this%supercellFoldingMatrix, errStatus,&
-            & supercellFoldingDiagOut=this%supercellFoldingDiag)
-        @:PROPAGATE_ERROR(errStatus)
-      end if
-
     #:if WITH_TRANSPORT
       if (this%tUpload) then
         call overrideContactCharges(this%qInput, this%chargeUp, this%transpar, this%qBlockIn,&
             & this%blockUp)
       end if
     #:endif
 
+      if (this%isHybridXc .and. this%tPeriodic) then
+        this%supercellFoldingDiag = nint(diagonal(this%supercellFoldingMatrix(:,:3)))
+      end if
+
     endif
 
+
     if (.not. allocated(this%reks)) then
       !Input charges packed into unique equivalence elements
       #:for NAME in [('this%qDiffRed'),('this%qInpRed'),('this%qOutRed')]
@@ -6081,7 +6082,7 @@ end subroutine getHybridXcCutOff_cluster
 
 
   !> Determine range separated cut-off and also update maximal cutoff
-  subroutine getHybridXcCutOff_gamma(cutOff, latVecs, cutoffRed, gSummationCutoff,&
+  subroutine getHybridXcCutOff_gamma(cutOff, latVecs, cutoffRed, errStatus, gSummationCutoff,&
       & gammaCutoff)
 
     !> Resulting cutoff
@@ -6093,15 +6094,15 @@ subroutine getHybridXcCutOff_gamma(cutOff, latVecs, cutoffRed, gSummationCutoff,
     !> CAM-neighbour list cutoff reduction
     real(dp), intent(in) :: cutoffRed
 
+    !> Operation status, if an error needs to be returned
+    type(TStatus), intent(inout) :: errStatus
+
     !> Cutoff for real-space g-summation
     real(dp), intent(in), optional :: gSummationCutoff
 
     !> Coulomb truncation cutoff for Gamma electrostatics
     real(dp), intent(in), optional :: gammaCutoff
 
-    !! Error status
-    type(TStatus) :: errStatus
-
     if (cutoffRed < 0.0_dp) then
       call error("Cutoff reduction for hybrid xc-functional neighbours should be zero or positive.")
     end if
@@ -6120,7 +6121,6 @@ subroutine getHybridXcCutOff_gamma(cutOff, latVecs, cutoffRed, gSummationCutoff,
       allocate(cutOff%gammaCutoff)
       call getCoulombTruncationCutoff(latVecs, cutOff%gammaCutoff, errStatus)
       @:PROPAGATE_ERROR(errStatus)
-      if (errStatus%hasError()) call error(errStatus%message)
     end if
 
     if (present(gSummationCutoff)) then
@@ -6134,7 +6134,7 @@ end subroutine getHybridXcCutOff_gamma
 
 
   !> Determine range separated cut-off and also update maximal cutoff
-  subroutine getHybridXcCutOff_kpts(cutOff, latVecs, cutoffRed, supercellFoldingDiag,&
+  subroutine getHybridXcCutOff_kpts(cutOff, latVecs, cutoffRed, supercellFoldingDiag, errStatus,&
       & gSummationCutoff, wignerSeitzReduction, gammaCutoff)
 
     !> Resulting cutoff
@@ -6149,6 +6149,9 @@ subroutine getHybridXcCutOff_kpts(cutOff, latVecs, cutoffRed, supercellFoldingDi
     !> Supercell folding coefficients and shifts
     integer, intent(in) :: supercellFoldingDiag(:)
 
+    !> Operation status, if an error needs to be returned
+    type(TStatus), intent(inout) :: errStatus
+
     !> Cutoff for real-space g-summation
     real(dp), intent(in), optional :: gSummationCutoff
 
@@ -6159,9 +6162,6 @@ subroutine getHybridXcCutOff_kpts(cutOff, latVecs, cutoffRed, supercellFoldingDi
     !> Coulomb truncation cutoff for Gamma electrostatics
     real(dp), intent(in), optional :: gammaCutoff
 
-    !! Error status
-    type(TStatus) :: errStatus
-
     if (cutoffRed < 0.0_dp) then
       call error("Cutoff reduction for hybrid xc-functional neighbours should be zero or positive.")
     end if
@@ -6181,7 +6181,6 @@ subroutine getHybridXcCutOff_kpts(cutOff, latVecs, cutoffRed, supercellFoldingDi
       call getCoulombTruncationCutoff(latVecs, cutOff%gammaCutoff, errStatus,&
           & nK=supercellFoldingDiag)
       @:PROPAGATE_ERROR(errStatus)
-      if (errStatus%hasError()) call error(errStatus%message)
     end if
 
     if (present(gSummationCutoff)) then