From 2f9e75b894fd8b26be181cf444b043bd5c7287d2 Mon Sep 17 00:00:00 2001
From: Alessandro Pecchia <alessandro.pecchia@ismn.cnr.it>
Date: Fri, 10 May 2024 14:15:49 +0200
Subject: [PATCH] Fix for DM with k-points and add test labels

---
 src/dftbp/dftbplus/mainio.F90   |  4 ++--
 src/dftbp/transport/negfint.F90 | 42 +++++++++++++++++++++------------
 test/app/dftb+/CMakeLists.txt   |  1 +
 3 files changed, 30 insertions(+), 17 deletions(-)

diff --git a/src/dftbp/dftbplus/mainio.F90 b/src/dftbp/dftbplus/mainio.F90
index e7d3b89f56..8d207b4de5 100644
--- a/src/dftbp/dftbplus/mainio.F90
+++ b/src/dftbp/dftbplus/mainio.F90
@@ -3689,9 +3689,9 @@ subroutine writeDetailedOut8(fd, neFermi)
     if (allocated(neFermi)) then
       write(fd,"(A)", advance='no')'Density of states at the Fermi energy (a.u.): '
       if (size(neFermi)==2) then
-        write(fd,"(E12.6,A,E12.6,A)")neFermi(1), ' (up) ', neFermi(2), ' (down)'
+        write(fd,"(E13.6,A,E13.6,A)") neFermi(1), ' (up) ', neFermi(2), ' (down)'
       else
-        write(fd,"(E12.6)")neFermi
+        write(fd,"(E13.6)") neFermi
       end if
     end if
 
diff --git a/src/dftbp/transport/negfint.F90 b/src/dftbp/transport/negfint.F90
index 4bb619e813..47e89beefe 100644
--- a/src/dftbp/transport/negfint.F90
+++ b/src/dftbp/transport/negfint.F90
@@ -931,7 +931,7 @@ subroutine negf_density(negf,spin,nkpoint,HH,SS,mu,DensMat,EnMat)
     if(present(EnMat)) then
        params%DorE = 'E'
        call set_params(negf,params)
-       call pass_DM(negf,rhoE=EnMat)
+       call pass_DM(negf,rho=EnMat)
     endif
     if (present(DensMat).and.present(EnMat)) then
        params%DorE  = 'B'
@@ -976,7 +976,7 @@ subroutine  negf_density_inel(negf, DensMat, EnMat)
     if(present(EnMat)) then
        params%DorE = 'E'
        call set_params(negf,params)
-       call pass_DM(negf,rhoE=EnMat)
+       call pass_DM(negf,rho=EnMat)
     endif
 
     call compute_density_dft(negf)
@@ -1374,8 +1374,9 @@ subroutine calc_density_ela()
       end do
 
      #:if WITH_MPI
+      ! In place reduce of the density matrix along energy (groupComm)
+      ! and k-points (interGroupComm)
       do iS = 1, nSpin
-        ! In place all-reduce of the density matrix
         call mpifx_allreduceip(env%mpi%groupComm, rho(:,iS), MPI_SUM)
       end do
       call mpifx_allreduceip(env%mpi%interGroupComm, rho, MPI_SUM)
@@ -1388,8 +1389,9 @@ end subroutine calc_density_ela
     ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     subroutine calc_density_inel()
 
-      type(z_CSR), target :: csrDens
-      type(z_CSR), pointer :: pcsrDens
+      !type(z_CSR), target :: csrDens
+      !type(z_CSR), pointer :: pcsrDens
+      type(TMatrixArray) :: csrDens(nKS)
 
       call get_params(this%negf, params)
       if (params%Np_real == 0) then
@@ -1400,9 +1402,9 @@ subroutine calc_density_inel()
         call error('Collinear spin not supported with inelastic scattering yet')
       end if
 
-      if (nKS > 1) then
-        call error('Inelastic DM currently works with just one k-point')
-      end if
+      !if (nKS > 1) then
+      !  call error('Inelastic DM currently works with just one k-point')
+      !end if
 
       write(stdOut, *)
       write(stdOut, '(80("="))')
@@ -1429,10 +1431,13 @@ subroutine calc_density_inel()
               & img2CentCell, iCellVec, cellVec, orb)
 
         call copy_HS(this%negf, this%csrHam, this%csrOver, iKS)
+        
+        allocate(csrDens(iKS)%Mat)
+        call pass_DM(this%negf, rho=csrDens(iKS)%Mat, iKS=iKS)
       end do
 
-      pcsrDens => csrDens
-      call pass_DM(this%negf,rho=pcsrDens)
+      !pcsrDens => csrDens
+      !call pass_DM(this%negf,rho=pcsrDens)
 
       call compute_density_dft(this%negf)
 
@@ -1440,13 +1445,20 @@ subroutine calc_density_inel()
         iK = groupKS(1, iKS)
         iS = groupKS(2, iKS)
 
-        !call get_DM(this%negf, csrDens)
-
-        call unfoldFromCSR(rho(:,1), csrDens, kPoints(:,iK), kWeights(iK), iAtomStart, &
+        call unfoldFromCSR(rho(:,1), csrDens(iKS)%Mat, kPoints(:,iK), kWeights(iK), iAtomStart, &
                      & iPair, iNeighbor, nNeighbor, img2CentCell, iCellVec, cellVec, orb)
+     
+        call destruct(csrDens(iKS)%Mat)
       end do
-
-      call destruct(csrDens)
+     
+     #:if WITH_MPI
+      ! In place reduce of the density matrix along energy (groupComm)
+      ! and k-points (interGroupComm)
+      do iS = 1, nSpin
+        call mpifx_allreduceip(env%mpi%groupComm, rho(:,iS), MPI_SUM)
+      end do
+      call mpifx_allreduceip(env%mpi%interGroupComm, rho, MPI_SUM)
+     #:endif
 
       write(stdOut,'(80("="))')
       write(stdOut,*)
diff --git a/test/app/dftb+/CMakeLists.txt b/test/app/dftb+/CMakeLists.txt
index da9dac2a30..b3412b3bad 100644
--- a/test/app/dftb+/CMakeLists.txt
+++ b/test/app/dftb+/CMakeLists.txt
@@ -39,4 +39,5 @@ foreach(test IN LISTS tests)
     dftb+_${test}
     PROPERTIES
     ENVIRONMENT "DFTBPLUS_PARAM_DIR=${PROJECT_SOURCE_DIR}/external")
+  set_tests_properties(dftb+_${test} PROPERTIES LABELS "dftb+/${test}") 
 endforeach()