[pull] main from dftbplus:main

CHANGELOG.rst

@@ -11,6 +11,8 @@ Unreleased
 Added
 -----
 
+- MPI-parallelization of Waveplot
+
 - Generalization of mixers to also handle complex density matrices
 
 - General range-separated, long-range corrected CAM hybrid functionals for

app/waveplot/gridcache.F90

@@ -13,8 +13,12 @@
 module waveplot_gridcache
   use dftbp_common_accuracy, only : dp
   use dftbp_common_constants, only : pi
+  use dftbp_common_environment, only : TEnvironment
   use dftbp_common_file, only : TFileDescr, openFile, closeFile
   use dftbp_common_globalenv, only : stdOut
+#:if WITH_MPI
+  use dftbp_common_schedule, only : getStartAndEndIndex
+#:endif
   use dftbp_io_message, only : error
   use waveplot_molorb, only : TMolecularOrbital, getValue
   implicit none
@@ -108,14 +112,17 @@ module waveplot_gridcache
 
   !> Initialises a GridCache instance.
   !! Caveat: Level index is not allowed to contain duplicate entries!
-  subroutine TGridCache_init(sf, levelIndex, nOrb, nAllLevel, nAllKPoint, nAllSpin, nCached,&
-      & nPoints, tVerbose, eigvecbin, gridVec, origin, kPointCoords, tReal, molorb)
+  subroutine TGridCache_init(sf, env, levelIndexAll, nOrb, nAllLevel, nAllKPoint, nAllSpin,&
+      & nCached, nPoints, tVerbose, eigvecBin, gridVec, origin, kPointCoords, tReal, molorb)
 
     !> Structure to initialise
     type(TgridCache), intent(inout) :: sf
 
+    !> Environment settings
+    type(TEnvironment), intent(in) :: env
+
     !> Contains indexes (spin, kpoint, state) of the levels which must be calculated
-    integer, intent(in) :: levelIndex(:,:)
+    integer, intent(in) :: levelIndexAll(:,:)
 
     !> Nr. of orbitals (elements) in an eigenvectors
     integer, intent(in) :: nOrb
@@ -139,7 +146,7 @@ subroutine TGridCache_init(sf, levelIndex, nOrb, nAllLevel, nAllKPoint, nAllSpin
     logical, intent(in) :: tVerbose
 
     !> Name of the binary eigenvector file
-    character(len=*), intent(in) :: eigvecbin
+    character(len=*), intent(in) :: eigvecBin
 
     !> Grid vectors
     real(dp), intent(in) :: gridVec(:,:)
@@ -156,18 +163,34 @@ subroutine TGridCache_init(sf, levelIndex, nOrb, nAllLevel, nAllKPoint, nAllSpin
     !> Molecular orbital calculator
     type(TMolecularOrbital), pointer, intent(in) :: molorb
 
+    !! Contains indexes (spin, kpoint, state) to be calculated by the current MPI process
+    integer, allocatable :: levelIndex(:,:)
+
+    !! Start and end level index for MPI parallelization, if applicable
+    integer :: iLvlStart, iLvlEnd
+
+    !! Auxiliary variables
     integer ::nAll
     integer :: iSpin, iKPoint, iLevel, ind, ii, iostat
     integer :: curVec(3)
     logical :: tFound
 
+  #:if WITH_MPI
+    call getStartAndEndIndex(env, size(levelIndexAll, dim=2), iLvlStart, iLvlEnd)
+  #:else
+    iLvlStart = 1
+    iLvlEnd = size(levelIndexAll, dim=2)
+  #:endif
+
+    levelIndex = levelIndexAll(:, iLvlStart:iLvlEnd)
+
     @:ASSERT(.not. sf%tInitialised)
     @:ASSERT(size(levelIndex, dim=1) == 3)
     @:ASSERT(size(levelIndex, dim=2) > 0)
     @:ASSERT(minval(levelIndex) > 0)
-    @:ASSERT(maxval(levelIndex(1,:)) <= nAllLevel)
-    @:ASSERT(maxval(levelIndex(2,:)) <= nAllKPoint)
-    @:ASSERT(maxval(levelIndex(3,:)) <= nAllSpin)
+    @:ASSERT(maxval(levelIndex(1, :)) <= nAllLevel)
+    @:ASSERT(maxval(levelIndex(2, :)) <= nAllKPoint)
+    @:ASSERT(maxval(levelIndex(3, :)) <= nAllSpin)
     @:ASSERT(associated(molorb))
     @:ASSERT(all(shape(gridVec) == [3, 3]))
     @:ASSERT(size(origin) == 3)
@@ -204,14 +227,12 @@ subroutine TGridCache_init(sf, levelIndex, nOrb, nAllLevel, nAllKPoint, nAllSpin
           curVec = [iLevel, iKPoint, iSpin]
           tFound = .false.
           lpLevelIndex: do ii = 1, size(levelIndex, dim=2)
-            tFound = all(levelIndex(:,ii) == curVec)
-            if (tFound) then
-              exit lpLevelIndex
-            end if
+            tFound = all(levelIndex(:, ii) == curVec)
+            if (tFound) exit lpLevelIndex
           end do lpLevelIndex
           if (tFound) then
             sf%levelIndex(:, ind) = curVec
-            ind = ind +1
+            ind = ind + 1
           end if
         end do
       end do
@@ -222,7 +243,7 @@ subroutine TGridCache_init(sf, levelIndex, nOrb, nAllLevel, nAllKPoint, nAllSpin
     sf%cachePos = 1
     sf%nReadEigVec = 0
     sf%tFinished = .false.
-    call openFile(sf%fdEigVec, eigvecbin, mode="rb", iostat=iostat)
+    call openFile(sf%fdEigVec, eigvecBin, mode="rb", iostat=iostat)
     if (iostat /= 0) then
       call error("Can't open file '" // trim(eigvecBin) // "'.")
     end if
@@ -249,7 +270,7 @@ subroutine TGridCache_next_real(sf, gridValReal, levelIndex, tFinished)
 
     complex(dp), pointer, save :: gridValCmpl(:,:,:) => null()
 
-    call local_next(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
+    call localNext(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
 
   end subroutine TGridCache_next_real
 
@@ -271,13 +292,13 @@ subroutine TGridCache_next_cmpl(sf, gridValCmpl, levelIndex, tFinished)
 
     real(dp), pointer, save :: gridValReal(:,:,:) => null()
 
-    call local_next(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
+    call localNext(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
 
   end subroutine TGridCache_next_cmpl
 
 
   !> Working subroutine for the TGridCache_next_* subroutines.
-  subroutine local_next(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
+  subroutine localNext(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
 
     !> Gridcache instance
     type(TgridCache), intent(inout), target :: sf
@@ -362,6 +383,6 @@ subroutine local_next(sf, gridValReal, gridValCmpl, levelIndex, tFinished)
     end if
     tFinished = sf%tFinished
 
-  end subroutine local_next
+  end subroutine localNext
 
 end module waveplot_gridcache

app/waveplot/initwaveplot.F90

@@ -10,8 +10,9 @@
 !> Contains the routines for initialising Waveplot.
 module waveplot_initwaveplot
   use dftbp_common_accuracy, only : dp
+  use dftbp_common_environment, only : TEnvironment
   use dftbp_common_file, only : TFileDescr, openFile, closeFile, setDefaultBinaryAccess
-  use dftbp_common_globalenv, only : stdOut
+  use dftbp_common_globalenv, only : stdOut, tIoProc
   use dftbp_common_release, only : releaseYear
   use dftbp_common_status, only : TStatus
   use dftbp_common_unitconversion, only : lengthUnits
@@ -28,10 +29,10 @@ module waveplot_initwaveplot
   use dftbp_io_hsdutils2, only : convertUnitHsd, readHSDAsXML, warnUnprocessedNodes, renameChildren
   use dftbp_io_message, only : warning, error
   use dftbp_io_xmlutils, only : removeChildNodes
+  use dftbp_math_simplealgebra, only : determinant33
   use dftbp_type_linkedlist, only : TListIntR1, TListReal, init, destruct, len, append, asArray
   use dftbp_type_typegeometryhsd, only : TGeometry, readTGeometryGen, readTGeometryHSD,&
       & readTGeometryVasp, readTGeometryXyz, writeTGeometryHSD
-  use dftbp_math_simplealgebra, only : determinant33
   use waveplot_gridcache, only : TGridCache, TGridCache_init
   use waveplot_molorb, only : TMolecularOrbital, TMolecularOrbital_init, TSpeciesBasis
   use waveplot_slater, only : TSlaterOrbital_init
@@ -243,11 +244,14 @@ module waveplot_initwaveplot
 
 
   !> Initialises the program variables.
-  subroutine TProgramVariables_init(this)
+  subroutine TProgramVariables_init(this, env)
 
     !> Container of program variables
     type(TProgramVariables), intent(out), target :: this
 
+    !> Environment settings
+    type(TEnvironment), intent(inout) :: env
+
     !! Pointers to input nodes
     type(fnode), pointer :: root, tmp, detailed, hsdTree
 
@@ -327,12 +331,18 @@ subroutine TProgramVariables_init(this)
     call warnUnprocessedNodes(root, .true.)
 
     ! Finish parsing, dump parsed and processed input
-    call dumpHSD(hsdTree, hsdParsedInput)
+    if (tIoProc) then
+      call dumpHSD(hsdTree, hsdParsedInput)
+    end if
     write(stdout, "(A)") "Processed input written as HSD to '" // hsdParsedInput &
         &//"'"
     write(stdout, "(A,/)") repeat("-", 80)
     call destroyNode(hsdTree)
 
+  #:if WITH_MPI
+    call env%initMpi(1)
+  #:endif
+
     if (this%input%geo%tPeriodic) then
       call TBoundaryConditions_init(this%boundaryCond, boundaryConditions%pbc3d, errStatus)
     else if (this%input%geo%tHelical) then
@@ -368,11 +378,10 @@ subroutine TProgramVariables_init(this)
     call TMolecularOrbital_init(this%loc%molOrb, this%input%geo, this%boundaryCond,&
         & this%basis%basis)
 
-    call TGridCache_init(this%loc%grid, levelIndex=this%loc%levelIndex, nOrb=this%input%nOrb,&
-        & nAllLevel=this%eig%nState, nAllKPoint=nKPoint, nAllSpin=nSpin, nCached=nCached,&
-        & nPoints=this%opt%nPoints, tVerbose=this%opt%tVerbose, eigVecBin=eigVecBin,&
-        & gridVec=this%loc%gridVec, origin=this%opt%gridOrigin,&
-        & kPointCoords=kPointsWeights(1:3, :), tReal=this%input%tRealHam, molorb=this%loc%pMolOrb)
+    call TGridCache_init(this%loc%grid, env, this%loc%levelIndex, this%input%nOrb, this%eig%nState,&
+        & nKPoint, nSpin, nCached, this%opt%nPoints, this%opt%tVerbose, eigVecBin,&
+        & this%loc%gridVec, this%opt%gridOrigin, kPointsWeights(1:3, :), this%input%tRealHam,&
+        & this%loc%pMolOrb)
 
   end subroutine TProgramVariables_init
 
@@ -530,8 +539,8 @@ subroutine readOptions(this, node, nLevel, nKPoint, nSpin, nCached, tShiftGrid)
 
     !! Warning issued, if the detailed.xml id does not match the eigenvector id
     character(len=63) :: warnId(3) = [&
-        & "The external files you are providing differ from those provided", &
-        & "when this input file was generated. The results you obtain with", &
+        & "The external files you are providing differ from those provided",&
+        & "when this input file was generated. The results you obtain with",&
         & "the current files could therefore be different.                "]
 
     !! Auxiliary variables

app/waveplot/molorb.F90

@@ -170,17 +170,17 @@ subroutine TMolecularOrbital_init(this, geometry, boundaryCond, basis)
     ! Get cells to look for when adding STOs from periodic images
     this%tPeriodic = geometry%tPeriodic
     if (this%tPeriodic) then
-      allocate(this%latVecs(3,3))
-      allocate(this%recVecs2p(3,3))
+      allocate(this%latVecs(3, 3))
+      allocate(this%recVecs2p(3, 3))
       this%latVecs(:,:) = geometry%latVecs
       call invert33(this%recVecs2p, this%latVecs)
-      this%recVecs2p = reshape(this%recVecs2p, [3, 3], order=[2, 1])
+      this%recVecs2p(:,:) = reshape(this%recVecs2p, [3, 3], order=[2, 1])
       mCutoff = maxval(this%cutoffs)
       call getCellTranslations(this%cellVec, rCellVec, this%latVecs, this%recVecs2p, mCutoff)
       this%nCell = size(this%cellVec,dim=2)
     else
-      allocate(this%latVecs(3,0))
-      allocate(this%recVecs2p(3,0))
+      allocate(this%latVecs(3, 0))
+      allocate(this%recVecs2p(3, 0))
       allocate(this%cellVec(3, 1))
       this%cellVec(:,:) = 0.0_dp
       allocate(rCellVec(3, 1))
@@ -227,10 +227,10 @@ subroutine TMolecularOrbital_getValue_real(this, origin, gridVecs, eigVecsReal,
     !> Add densities instead of wave functions
     logical, intent(in), optional :: addDensities
 
-    real(dp), save :: kPoints(3,0)
+    real(dp), save :: kPoints(3, 0)
     integer, save :: kIndexes(0)
-    complex(dp), save :: valueCmpl(0,0,0,0)
-    complex(dp), save :: eigVecsCmpl(0,0)
+    complex(dp), save :: valueCmpl(0, 0, 0, 0)
+    complex(dp), save :: eigVecsCmpl(0, 0)
     logical :: tAddDensities
 
     @:ASSERT(this%tInitialised)
@@ -246,7 +246,7 @@ subroutine TMolecularOrbital_getValue_real(this, origin, gridVecs, eigVecsReal,
       tAddDensities = .false.
     end if
 
-    call local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, this%nAtom, this%nOrb,&
+    call localGetValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, this%nAtom, this%nOrb,&
         & this%coords, this%species, this%cutoffs, this%iStos, this%angMoms, this%stos,&
         & this%tPeriodic, .true., this%latVecs, this%recVecs2p, kPoints, kIndexes, this%nCell,&
         & this%cellVec, tAddDensities, valueOnGrid, valueCmpl)
@@ -295,7 +295,7 @@ subroutine TMolecularOrbital_getValue_cmpl(this, origin, gridVecs, eigVecsCmpl,
     @:ASSERT(maxval(kIndexes) <= size(kPoints, dim=2))
     @:ASSERT(minval(kIndexes) > 0)
 
-    call local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, this%nAtom, this%nOrb,&
+    call localGetValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, this%nAtom, this%nOrb,&
         & this%coords, this%species, this%cutoffs, this%iStos, this%angMoms, this%stos,&
         & this%tPeriodic, .false., this%latVecs, this%recVecs2p, kPoints, kIndexes, this%nCell,&
         & this%cellVec, tAddDensities, valueReal, valueOnGrid)
@@ -306,7 +306,7 @@ end subroutine TMolecularOrbital_getValue_cmpl
   !> Returns the values of several molecular orbitals on grids.
   !! Caveat: The flag tPeriodic decides if the complex or the real version is read/written for the
   !! various parameters.
-  subroutine local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, nAtom, nOrb, coords,&
+  subroutine localGetValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, nAtom, nOrb, coords,&
       & species, cutoffs, iStos, angMoms, stos, tPeriodic, tReal, latVecs, recVecs2p, kPoints,&
       & kIndexes, nCell, cellVec, tAddDensities, valueReal, valueCmpl)
 
@@ -409,12 +409,12 @@ subroutine local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, nAtom, nOr
 
     ! Loop over all grid points
     lpI3: do i3 = 1, nPoints(3)
-      curCoords(:, 3) = real(i3 - 1, dp) * gridVecs(:,3)
+      curCoords(:, 3) = real(i3 - 1, dp) * gridVecs(:, 3)
       lpI2: do i2 = 1, nPoints(2)
-        curCoords(:, 2) =  real(i2 - 1, dp) * gridVecs(:,2)
+        curCoords(:, 2) =  real(i2 - 1, dp) * gridVecs(:, 2)
         lpI1: do i1 = 1, nPoints(1)
-          curCoords(:, 1) = real(i1 - 1, dp) * gridVecs(:,1)
-          xyz(:) = sum(curCoords, dim=2) + origin(:)
+          curCoords(:, 1) = real(i1 - 1, dp) * gridVecs(:, 1)
+          xyz(:) = sum(curCoords, dim=2) + origin
           if (tPeriodic) then
             frac(:) = matmul(xyz, recVecs2p)
             xyz(:) = matmul(latVecs, frac - real(floor(frac), dp))
@@ -425,9 +425,9 @@ subroutine local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, nAtom, nOr
             ind = 1
             lpAtom: do iAtom = 1, nAtom
               iSpecies = species(iAtom)
-              diff(:) = xyz - coords(:,iAtom, iCell)
-              xx = sqrt(sum(diff**2))
-              lpOrb: do iOrb = iStos(iSpecies), iStos(iSpecies+1)-1
+              diff(:) = xyz - coords(:, iAtom, iCell)
+              xx = norm2(diff)
+              lpOrb: do iOrb = iStos(iSpecies), iStos(iSpecies + 1) - 1
                 iL = angMoms(iOrb)
                 ! Calculate wave function only if atom is inside the cutoff
                 if (xx <= cutoffs(iOrb)) then
@@ -469,7 +469,7 @@ subroutine local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, nAtom, nOr
           ! the eigenvector)
           if (tReal) then
             if (tAddDensities) then
-              atomAllOrbVal = atomAllOrbVal**2
+              atomAllOrbVal(:,:) = atomAllOrbVal**2
             end if
             atomOrbValReal(:) = sum(atomAllOrbVal, dim=2)
             do iEig = 1, nPoints(4)
@@ -495,6 +495,6 @@ subroutine local_getValue(origin, gridVecs, eigVecsReal, eigVecsCmpl, nAtom, nOr
       end do lpI2
     end do lpI3
 
-  end subroutine local_getValue
+  end subroutine localGetValue
 
 end module waveplot_molorb

app/waveplot/slater.F90

@@ -81,7 +81,7 @@ function realTessY(ll, mm, coord, rrOpt) result(rty)
     if (present(rrOpt)) then
       rr = rrOpt
     else
-      rr = sqrt(sum(coord**2))
+      rr = norm2(coord)
     end if
 
     @:ASSERT(ll >= 0 .and. ll <= 3)