Merge pull request dftbplus#1433 from bhourahine/hybridNeighFix

Fixes array sizing/unallocated access

src/dftbp/dftb/hybridxc.F90

@@ -992,7 +992,9 @@ subroutine updateCoords_kpts(this, env, symNeighbourList, nNeighbourCamSym, skOv
 
     if (this%tScreeningInited) then
       this%hPrevCplxHS(:,:,:) = 0.0_dp
-      this%dRhoPrevCplxHS(:,:,:,:,:,:) = 0.0_dp
+      ! Note, this is only allocated if passing through addCamHamiltonianNeighbour_kpts_mic or
+      ! addCamHamiltonianNeighbour_kpts_ct :
+      if (allocated(this%dRhoPrevCplxHS)) this%dRhoPrevCplxHS(:,:,:,:,:,:) = 0.0_dp
       this%camEnergy = 0.0_dp
     end if
 
@@ -2961,7 +2963,6 @@ subroutine addCamHamiltonianNeighbour_kpts_mic(this, env, deltaRhoSqr, symNeighb
     if (.not. this%tScreeningInited) then
       allocate(this%hprevCplxHS(squareSize, squareSize, nKS))
       this%hprevCplxHS(:,:,:) = (0.0_dp, 0.0_dp)
-      this%dRhoPrevCplxHS = deltaRhoSqr
       ! there is no previous delta density matrix, therefore just copy over
       deltaDeltaRhoSqr = deltaRhoSqr
       this%tScreeningInited = .true.
@@ -3271,7 +3272,6 @@ subroutine addCamHamiltonianNeighbour_kpts_ct(this, env, deltaRhoSqr, symNeighbo
     if (.not. this%tScreeningInited) then
       allocate(this%hprevCplxHS(squareSize, squareSize, nKS))
       this%hprevCplxHS(:,:,:) = (0.0_dp, 0.0_dp)
-      this%dRhoPrevCplxHS = deltaRhoSqr
       ! there is no previous delta density matrix, therefore just copy over
       deltaDeltaRhoSqr = deltaRhoSqr
       this%tScreeningInited = .true.

src/dftbp/dftb/sparse2dense.F90

@@ -28,7 +28,7 @@ module dftbp_dftb_sparse2dense
   public :: unpackHS, packHS, iPackHS, packErho, unpackDQ
   public :: packHSPauli, packHSPauliImag, unpackHPauli, unpackSPauli
   public :: unpackHelicalHS, packHelicalHS
-  public :: getSparseDescriptor, getSparseSize
+  public :: getSparseDescriptor
 
 #:if WITH_SCALAPACK
   public :: unpackHSRealBlacs, unpackHSCplxBlacs, unpackHPauliBlacs, unpackSPauliBlacs
@@ -3227,41 +3227,4 @@ subroutine getSparseDescriptor(iNeighbour, nNeighbour, img2CentCell, orb, iPair,
 
   end subroutine getSparseDescriptor
 
-
-  !> Calculates number of elements in sparse arrays like the real space overlap.
-  subroutine getSparseSize(iNeighbour, nNeighbour, img2CentCell, orb, sparseSize)
-
-    !> Neighbours of each atom
-    integer, intent(in) :: iNeighbour(0:,:)
-
-    !> Number of neighbours of each atom
-    integer, intent(in) :: nNeighbour(:)
-
-    !> Indexing for mapping image atoms to central cell
-    integer, intent(in) :: img2CentCell(:)
-
-    !> Atomic orbital information
-    type(TOrbitals), intent(in) :: orb
-
-    !> Total number of elements in a sparse structure (ignoring extra indices like spin)
-    integer, intent(out) :: sparseSize
-
-    integer :: nAtom, mNeighbour
-    integer :: ind, iAt1, nOrb1, iNeigh1, nOrb2
-
-    nAtom = size(iNeighbour, dim=2)
-    mNeighbour = size(iNeighbour, dim=1)
-
-    ind = 0
-    do iAt1 = 1, nAtom
-      nOrb1 = orb%nOrbAtom(iAt1)
-      do iNeigh1 = 0, nNeighbour(iAt1)
-        nOrb2 = orb%nOrbAtom(img2CentCell(iNeighbour(iNeigh1, iAt1)))
-        ind = ind + nOrb1 * nOrb2
-      end do
-    end do
-    sparseSize = ind
-
-  end subroutine getSparseSize
-
 end module dftbp_dftb_sparse2dense

src/dftbp/dftbplus/main.F90

@@ -53,7 +53,7 @@ module dftbp_dftbplus_main
   use dftbp_dftb_slakocont, only : TSlakoCont
   use dftbp_dftb_sparse2dense, only : unpackHPauli, unpackHS, packHS, packHS, unpackHelicalHS,&
       & packerho, packHSPauli, packHelicalHS, packHSPauliImag, iPackHS, unpackSPauli,&
-      & getSparseDescriptor, getSparseSize
+      & getSparseDescriptor
   use dftbp_dftb_spin, only : ud2qm, qm2ud
   use dftbp_dftb_spinorbit, only : addOnsiteSpinOrbitHam, getOnsiteSpinOrbitEnergy
   use dftbp_dftb_stress, only : getkineticstress, getBlockStress, getBlockiStress, getNonSCCStress
@@ -4760,10 +4760,6 @@ subroutine getNextInputDensityCplx(env, ints, neighbourList, nNeighbourSK, dense
     !! Number of spins and spin index
     integer :: nSpin
 
-    !! Total size of orbitals in the sparse data structures, where the decay of the overlap sets the
-    !! sparsity pattern
-    integer :: sparseSize
-
     !! Sparse density matrix storage
     real(dp), allocatable :: rhoPrim(:,:)
 
@@ -4817,9 +4813,7 @@ subroutine getNextInputDensityCplx(env, ints, neighbourList, nNeighbourSK, dense
         #:endif
 
           ! Construct sparse density matrix for later Mulliken analysis
-          call getSparseSize(neighbourList%iNeighbour, nNeighbourSK, img2CentCell, orb,&
-              & sparseSize)
-          allocate(rhoPrim(sparseSize, nSpin), source=0.0_dp)
+          allocate(rhoPrim(size(ints%overlap), nSpin), source=0.0_dp)
 
           do iKS = 1, parallelKS%nLocalKS
             iK = parallelKS%localKS(1, iKS)