From 4af47e1074b822e9b0ebf65022ada1221c2eb4eb Mon Sep 17 00:00:00 2001
From: Alessandro Pecchia <alessandro.pecchia@ismn.cnr.it>
Date: Wed, 6 Mar 2024 19:46:02 +0100
Subject: [PATCH] Add local current to tagged out

---
 src/dftbp/dftbplus/main.F90   |  3 ++-
 src/dftbp/dftbplus/mainio.F90 | 19 +++++++++++++------
 src/dftbp/io/taggedoutput.F90 |  6 ++++++
 3 files changed, 21 insertions(+), 7 deletions(-)

diff --git a/src/dftbp/dftbplus/main.F90 b/src/dftbp/dftbplus/main.F90
index d0543113ff..ae50c14078 100644
--- a/src/dftbp/dftbplus/main.F90
+++ b/src/dftbp/dftbplus/main.F90
@@ -475,7 +475,8 @@ subroutine runDftbPlus(this, env)
           & this%tStress, this%totalStress, this%pDynMatrix,&
           & this%dftbEnergy(this%deltaDftb%iFinal), this%extPressure, this%coord0, this%tLocalise,&
           & localisation, this%electrostatPot, this%taggedWriter, this%tunneling, this%ldos,&
-          & this%lCurrArray, this%polarisability, this%dEidE, this%dipoleMoment, this%eFieldScaling)
+          & this%lCurrArray, this%leadCurrents, this%polarisability, this%dEidE, this%dipoleMoment,&
+          & this%eFieldScaling)
     end if
     if (this%tWriteResultsTag) then
       call writeResultsTag(resultsTag, this%dftbEnergy(this%deltaDftb%iFinal), this%derivs,&
diff --git a/src/dftbp/dftbplus/mainio.F90 b/src/dftbp/dftbplus/mainio.F90
index 3b5dec2e41..e7d3b89f56 100644
--- a/src/dftbp/dftbplus/mainio.F90
+++ b/src/dftbp/dftbplus/mainio.F90
@@ -1881,7 +1881,7 @@ end subroutine initOutputFile
   subroutine writeAutotestTag(fileName, electronicSolver, tPeriodic, cellVol, tMulliken, qOutput,&
       & derivs, chrgForces, excitedDerivs, tStress, totalStress, pDynMatrix, energy, pressure,&
       & endCoords, tLocalise, localisation, esp, taggedWriter, tunneling, ldos, lCurrArray,&
-      & polarisability, dEidE, dipoleMoment, eFieldScaling)
+      & currents, polarisability, dEidE, dipoleMoment, eFieldScaling)
 
     !> Name of output file
     character(*), intent(in) :: fileName
@@ -1950,6 +1950,9 @@ subroutine writeAutotestTag(fileName, electronicSolver, tPeriodic, cellVol, tMul
     !> This array is for testing only since it misses info
     real(dp), allocatable, intent(in) :: lCurrArray(:,:)
 
+    !> contact or layer currents
+    real(dp), allocatable, intent(in) :: currents(:)
+
     !> Static electric polarisability
     real(dp), intent(in), allocatable :: polarisability(:,:,:)
 
@@ -2025,6 +2028,10 @@ subroutine writeAutotestTag(fileName, electronicSolver, tPeriodic, cellVol, tMul
       end if
     end if
 
+    if (allocated(currents)) then
+      call taggedWriter%write(fd, tagLabels%contactCurrents, currents)
+    end if
+
     if (allocated(lCurrArray)) then
       call taggedWriter%write(fd, tagLabels%localCurrents, lCurrArray)
     end if
@@ -2498,12 +2505,12 @@ subroutine writeHessianOut(fileName, pDynMatrix, indMovedAtoms)
 
     integer :: ii, fd
     character(10) :: suffix1, suffix2
-    logical :: tPartialHessian = .false. 
+    logical :: tPartialHessian = .false.
 
     ! Sanity check in case some bug is introduced
     if (size(pDynMatrix, dim=2) /= 3*size(indMovedAtoms)) then
-      call error('Internal error: incorrect number of rows of dynamical Matrix')    
-    end if       
+      call error('Internal error: incorrect number of rows of dynamical Matrix')
+    end if
     ! It is a partial Hessian Calculation if DynMatrix is not squared
     if (size(pDynMatrix, dim=1) > size(pDynMatrix, dim=2)) then
       tPartialHessian = .true.
@@ -2511,10 +2518,10 @@ subroutine writeHessianOut(fileName, pDynMatrix, indMovedAtoms)
 
     if (tPartialHessian) then
       write(suffix1,'(I10)') indMovedAtoms(1)
-      write(suffix2,'(I10)') indMovedAtoms(size(indMovedAtoms))     
+      write(suffix2,'(I10)') indMovedAtoms(size(indMovedAtoms))
       open(newunit=fd, file=fileName//"."//trim(adjustl(suffix1))//"-"//trim(adjustl(suffix2)), &
             & action="write", status="replace")
-    else 
+    else
       open(newunit=fd, file=fileName, action="write", status="replace")
     end if
 
diff --git a/src/dftbp/io/taggedoutput.F90 b/src/dftbp/io/taggedoutput.F90
index fc12f1b571..0cf96783f7 100644
--- a/src/dftbp/io/taggedoutput.F90
+++ b/src/dftbp/io/taggedoutput.F90
@@ -152,6 +152,12 @@ module dftbp_io_taggedoutput
     !> total projected DOS vector
     character(lenLabel) :: ldos = 'total_localdos'
 
+    !> contact currents
+    character(lenLabel) :: contactCurrents = 'contact_currents'
+
+    !> layer currents
+    character(lenLabel) :: layerCurrents = 'layer_currents'
+
     !> total bond currents
     character(lenLabel) :: localCurrents = 'local_currents'