[pull] main from dftbplus:main

AUTHORS.rst

@@ -25,6 +25,8 @@ The DFTB+ development is being lead and coordinated by:
 The following people (in alphabetic order by their family names) have
 contributed to DFTB+ :
 
+* Matías Berdakin (Universidad Nacional de Córdoba, Argentina)
+
 * Franco Bonafé (Max Planck Institute for the Structure and Dynamics
   of Matter, Germany)
 

doc/dftb+/manual/dftbp.tex

@@ -5913,7 +5913,7 @@ \section{ElectronDynamics}
 Calculate the real-time propagation of both the electronic state of the system
 and the nuclei within the Ehfenfest approximation in the presence of an optional
 external time-dependent perturbation, as described in Ref.~\cite{realtime1}.
-The real-time propagation is available with DFTB as well as xTB Hamiltonians.
+The real-time propagation is available with DFTB as well as xTB Hamiltonians, and it can be run in parallel using OpenMP as well as MPI. However, calculations for periodic systems, hybrid functionals, DFTB+U, onsite corrections and ion dynamics are not implemented with MPI yet (only OpenMP parallelization can be used on those cases).
 
 \begin{ptable}
   \kw{Steps} &i &  & - & \\

src/dftbp/api/mm/dftbplus.h

@@ -351,7 +351,7 @@ void dftbp_set_coords_lattice_origin(DftbPlus *instance, const double *coords,
  * \param[in] cutoff Cutoff used to compute the neighbour list Unit: Bohr.
  *
  * \param[in] coordNeighbours Coordinates of all neighbour atom in all cells
- *     Shape: [3, nAllAtom]. Unit: Bohr.
+ *     Shape: [nAllAtom, 3]. Unit: Bohr.
  *
  * \param[in] neighbour2CentCell Index of the atom in the central cell a neighbour atom corresponds
  *     to Shape: [nAllAtom].

src/dftbp/api/mm/mmapi.F90

@@ -1243,7 +1243,7 @@ subroutine TDftbPlus_finalizeTimeProp(this)
     !> Instance
     class(TDftbPlus), intent(inout) :: this
 
-    call finalizeTimeProp(this%main)
+    call finalizeTimeProp(this%env, this%main)
 
   end subroutine TDftbPlus_finalizeTimeProp
 

src/dftbp/dftb/populations.F90

@@ -40,7 +40,8 @@ module dftbp_dftb_populations
 #:if WITH_SCALAPACK
   public :: denseMullikenRealBlacs
   public :: denseSubtractDensityOfAtomsRealNonperiodicBlacs,&
-      & denseSubtractDensityOfAtomsRealPeriodicBlacs
+      & denseSubtractDensityOfAtomsRealPeriodicBlacs,&
+      & denseSubtractDensityOfAtomsCmplxNonperiodicBlacs
 #:endif
 
 
@@ -1046,6 +1047,48 @@ subroutine denseSubtractDensityOfAtomsRealNonperiodicBlacs(env, parallelKS, q0,
   end subroutine denseSubtractDensityOfAtomsRealNonperiodicBlacs
 
 
+  !> Subtracts superposition of atomic densities from distributed, dense, square, complex-valued
+  !! density matrix.
+  !!
+  !! For spin-polarized calculations, q0 is distributed equally to alpha and beta density matrices.
+  subroutine denseSubtractDensityOfAtomsCmplxNonperiodicBlacs(env, parallelKS, q0, denseDesc, rho)
+
+    !> Environment settings
+    type(TEnvironment), intent(in) :: env
+
+    !> The k-points and spins to process
+    type(TParallelKS), intent(in) :: parallelKS
+
+    !> Reference atom populations
+    real(dp), intent(in) :: q0(:,:,:)
+
+    !> Dense matrix descriptor
+    type(TDenseDescr), intent(in) :: denseDesc
+
+    !> Spin polarized (lower triangular) matrix
+    complex(dp), intent(inout) :: rho(:,:,:)
+
+    integer :: nAtom, iKS, iAt, iOrbStart, nOrb, iOrb
+    complex(dp) :: tmp(size(q0, dim=1), size(q0, dim=1))
+
+    nAtom = size(q0, dim=2)
+
+    do iKS = 1, parallelKS%nLocalKS
+      do iAt = 1, nAtom
+        iOrbStart = denseDesc%iAtomStart(iAt)
+        nOrb = denseDesc%iAtomStart(iAt + 1) - iOrbStart
+        tmp(:,:) = (0.0_dp, 0.0_dp)
+        do iOrb = 1, nOrb
+          tmp(iOrb, iOrb) = -q0(iOrb, iAt, 1)
+        end do
+        call scalafx_addl2g(env%blacs%orbitalGrid, tmp(1:nOrb, 1:nOrb), denseDesc%blacsOrbSqr,&
+            & iOrbStart, iOrbStart, rho(:,:,iKS))
+      end do
+    end do
+
+  end subroutine denseSubtractDensityOfAtomsCmplxNonperiodicBlacs
+
+
   !> Subtracts superposition of atomic densities from distributed, dense, square, real-valued
   !! density matrix.
   !!

src/dftbp/dftbplus/initprogram.F90

@@ -4074,12 +4074,20 @@ subroutine initProgramVariables(this, input, env)
         call error("Electron dynamics is not compatibile with this spinor Hamiltonian")
       end if
 
-      if (withMpi) then
-        call error("Electron dynamics does not work with MPI yet")
+      if (input%ctrl%elecDynInp%tIons .and. withMPI) then
+        call error("Ion dynamics time propagation does not work with MPI yet")
+      end if
+
+      if (.not. this%tRealHS .and. withMpi) then
+        call error("Electron dynamics of periodic systems does not work with MPI yet")
+      end if
+
+      if (.not. this%tRealHS .and. withMpi) then
+        call error("Electron dynamics of periodic systems does not work with MPI yet")
       end if
 
-      if (this%tFixEf) then
-        call error("Electron dynamics does not work with fixed Fermi levels yet")
+      if ((allocated(this%dftbU) .or. allocated(this%onSiteElements)) .and. withMpi) then
+        call error("Electron dynamics with DFTB+U or onsite corrections not implemented with MPI yet")
       end if
 
       if (this%tSpinSharedEf) then
@@ -4118,7 +4126,7 @@ subroutine initProgramVariables(this, input, env)
           & this%mass, this%nAtom, this%atomEigVal, this%dispersion, this%nonSccDeriv,&
           & this%tPeriodic, this%parallelKS, this%tRealHS, this%kPoint, this%kWeight,&
           & this%isHybridXc, this%scc, this%tblite, this%eFieldScaling, this%hamiltonianType,&
-          & errStatus)
+          & this%denseDesc, errStatus)
       if (errStatus%hasError()) call error(errStatus%message)
 
     end if

src/dftbp/dftbplus/mainapi.F90

@@ -753,13 +753,16 @@ end subroutine initializeTimeProp
 
 
   !> Finalizes the dynamics (releases memory, closes eventual open files)
-  subroutine finalizeTimeProp(main)
+  subroutine finalizeTimeProp(env, main)
+
+    !> dftb+ environment
+    type(TEnvironment), intent(inout) :: env
 
     !> Instance
     type(TDftbPlusMain), intent(inout) :: main
 
     if (allocated(main%electronDynamics)) then
-      call finalizeDynamics(main%electronDynamics)
+      call finalizeDynamics(main%electronDynamics, env)
     end if
 
   end subroutine finalizeTimeProp

src/dftbp/dftbplus/parser.F90

@@ -5807,13 +5807,6 @@ subroutine readElecDynamics(node, input, geom, masses)
     real (dp) :: defPpRange(2)
     logical :: defaultWrite
 
-  #:if WITH_MPI
-    if (associated(node)) then
-      call detailedError(node, 'This DFTB+ binary has been compiled with MPI settings and &
-          & electron dynamics are not currently available for distributed parallel calculations.')
-    end if
-  #:endif
-
     call getChildValue(node, "Steps", input%steps)
     call getChildValue(node, "TimeStep", input%dt, modifier=modifier, child=child)
     call convertUnitHsd(char(modifier), timeUnits, child, input%dt)

src/dftbp/extlibs/scalapackfx.F90

@@ -12,9 +12,9 @@ module dftbp_extlibs_scalapackfx
 #:if WITH_SCALAPACK
   use libscalapackfx_module, only : blacsfx_gemr2d, blacsfx_gsum, blacsgrid, blocklist, CSRC_,&
       & DLEN_, linecomm, M_, MB_, N_, NB_, pblasfx_pgemm, pblasfx_phemm, pblasfx_psymm,&
-      & pblasfx_psymv, pblasfx_ptran, pblasfx_ptranc, RSRC_, scalafx_addg2l, scalafx_addl2g,&
+      & pblasfx_psymv, pblasfx_ptran, pblasfx_ptranc, pblasfx_ptranu, RSRC_, scalafx_addg2l, scalafx_addl2g,&
       & scalafx_cpg2l, scalafx_cpl2g, scalafx_getdescriptor, scalafx_getlocalshape,&
-      & scalafx_indxl2g, scalafx_infog2l, scalafx_islocal, scalafx_pgetrf, scalafx_phegv,&
+      & scalafx_indxl2g, scalafx_infog2l, scalafx_islocal, scalafx_pgetri, scalafx_pgetrf, scalafx_phegv,&
       & scalafx_phegvd, scalafx_phegvr, scalafx_pposv, scalafx_ppotrf, scalafx_ppotri,&
       & scalafx_psygv, scalafx_psygvd, scalafx_psygvr, size
 #:endif