From 89b6b4b9f545bc002ae207ed0338cb824f72b764 Mon Sep 17 00:00:00 2001 From: Ben Hourahine Date: Wed, 1 Jan 2025 11:05:44 +0000 Subject: [PATCH 1/2] Missed language compatability in parser --- src/dftbp/dftbplus/parser.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dftbp/dftbplus/parser.F90 b/src/dftbp/dftbplus/parser.F90 index 954403ba64..8311987cb4 100644 --- a/src/dftbp/dftbplus/parser.F90 +++ b/src/dftbp/dftbplus/parser.F90 @@ -5235,6 +5235,7 @@ subroutine readAnalysis(node, ctrl, geo) call getChildValue(node, "WriteBandOut", ctrl%tWriteBandDat, tWriteBandDatDefault) + call renameChildren(node, "Polarizability", "Polarisability") call getChild(node, "Polarisability", child=child, requested=.false.) call getChild(node, "ResponseKernel", child=child2, requested=.false.) if (associated(child) .or. associated(child2)) then From 0422e7504c5a988549ba42cd475cd1da108737bc Mon Sep 17 00:00:00 2001 From: Tammo van der Heide Date: Fri, 10 Jan 2025 10:20:32 +0100 Subject: [PATCH 2/2] Parallelize (MPI) the Waveplot code (#1587) --- CHANGELOG.rst | 2 + app/waveplot/gridcache.F90 | 55 +++++-- app/waveplot/initwaveplot.F90 | 31 ++-- app/waveplot/molorb.F90 | 40 ++--- app/waveplot/slater.F90 | 2 +- app/waveplot/waveplot.F90 | 301 +++++++++++++++++++--------------- 6 files changed, 254 insertions(+), 177 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index bafb54e331..57830f19c0 100644 --- a/CHANGELOG.rst +++ b/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 diff --git a/app/waveplot/gridcache.F90 b/app/waveplot/gridcache.F90 index 2c96345e7c..381e1b7e58 100644 --- a/app/waveplot/gridcache.F90 +++ b/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 diff --git a/app/waveplot/initwaveplot.F90 b/app/waveplot/initwaveplot.F90 index b227042ccd..b037de0b5a 100644 --- a/app/waveplot/initwaveplot.F90 +++ b/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 diff --git a/app/waveplot/molorb.F90 b/app/waveplot/molorb.F90 index 24b886f237..4b9766ef1f 100644 --- a/app/waveplot/molorb.F90 +++ b/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 diff --git a/app/waveplot/slater.F90 b/app/waveplot/slater.F90 index a4d4abb41e..dd63fe0576 100644 --- a/app/waveplot/slater.F90 +++ b/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) diff --git a/app/waveplot/waveplot.F90 b/app/waveplot/waveplot.F90 index 22eac9bb99..90a49f4ad0 100644 --- a/app/waveplot/waveplot.F90 +++ b/app/waveplot/waveplot.F90 @@ -10,10 +10,12 @@ !> Program for plotting molecular orbitals as cube files. program waveplot use dftbp_common_accuracy, only : dp + use dftbp_common_environment, only : TEnvironment, TEnvironment_init use dftbp_common_file, only : TFileDescr, openFile, closeFile - use dftbp_common_globalenv, only : stdOut + use dftbp_common_globalenv, only : stdOut, initGlobalEnv, synchronizeAll, destructGlobalEnv use dftbp_dftb_periodic, only : getCellTranslations use dftbp_io_charmanip, only : i2c + use dftbp_io_message, only : warning, error use dftbp_math_simplealgebra, only : invert33 use dftbp_type_linkedlist, only : TListInt, TListRealR1, len, init, append, asArray use dftbp_type_typegeometry, only : TGeometry @@ -21,9 +23,7 @@ program waveplot use waveplot_initwaveplot, only : TProgramVariables, TProgramVariables_init use waveplot_molorb, only : getValue #:if WITH_MPI - use mpi, only : MPI_THREAD_FUNNELED - use dftbp_common_mpienv, only : TMpiEnv, TMpiEnv_init - use dftbp_extlibs_mpifx, only : mpifx_init_thread, mpifx_finalize + use dftbp_extlibs_mpifx, only : MPI_SUM, MPI_LOR, mpifx_allreduceip, mpifx_bcast #:endif implicit none @@ -31,6 +31,9 @@ program waveplot !> Container of program variables type(TProgramVariables), target :: wp + !> Environment settings + type(TEnvironment) :: env + !> Pointer to real-valued grid real(dp), pointer :: gridValReal(:,:,:) @@ -62,51 +65,43 @@ program waveplot type(TListInt) :: speciesList !> Auxiliary variables - integer :: i1, i2, i3 - integer :: iCell, iLevel, iKPoint, iSpin, iSpinOld, iAtom, iSpecies, iAng, mAng, ind, nBox - logical :: tFinished, tPlotLevel + integer :: i1, i2, i3, ioStat + integer :: iCell, iLevel, iKPoint, iSpin, iAtom, iSpecies, iAng, mAng, ind, nBox + logical :: tFinished, tPlotLevel, hasIoError character(len=80) :: comments(2), fileName real(dp) :: mDist, dist real(dp) :: cellMiddle(3), boxMiddle(3), frac(3), cubeCorner(3), coord(3), shift(3) real(dp) :: invBoxVecs(3,3), recVecs2p(3,3) real(dp), allocatable :: cellVec(:,:), rCellVec(:,:) -#:if WITH_MPI - !> MPI environment - type(TMpiEnv) :: mpiEnv - - ! As this is serial code, trap for run time execution on more than 1 processor with MPI enabled - call mpifx_init_thread(requiredThreading=MPI_THREAD_FUNNELED) - call TMpiEnv_init(mpiEnv) - call mpiEnv%mpiSerialEnv() -#:endif + call initGlobalEnv() + call TEnvironment_init(env) ! Allocate resources - call TProgramVariables_init(wp) - write(stdout, "(/,A,/)") "Starting main program" + call TProgramVariables_init(wp, env) + write(stdOut, "(/,A,/)") "Starting main program" ! Allocating buffer for general grid, total charge and spin up allocate(buffer(wp%opt%nPoints(1), wp%opt%nPoints(2), wp%opt%nPoints(3))) if (wp%opt%tCalcTotChrg) then - allocate(totChrg(wp%opt%nPoints(1), wp%opt%nPoints(2), wp%opt%nPoints(3))) - totChrg = 0.0_dp + allocate(totChrg(wp%opt%nPoints(1), wp%opt%nPoints(2), wp%opt%nPoints(3)), source=0.0_dp) if (wp%opt%tPlotTotSpin) then - allocate(spinUp(wp%opt%nPoints(1), wp%opt%nPoints(2), wp%opt%nPoints(3))) + allocate(spinUp(wp%opt%nPoints(1), wp%opt%nPoints(2), wp%opt%nPoints(3)), source=0.0_dp) end if end if - write (comments(1), "(A)") "Cube file generated by WAVEPLOT from data created by DFTB+" + write(comments(1), "(A)") "Cube file generated by WAVEPLOT from data created by DFTB+" ! Repeat boxes if necessary nBox = product(wp%opt%repeatBox) if (nBox > 1) then - ! If tFillBox is off, coordinates must be repeated here. Otherwise the part for filling with - ! atoms will do that. + ! If tFillBox is off, coordinates must be repeated here. + ! Otherwise the part for filling with atoms will do that. if (.not. wp%opt%tFillBox) then allocate(coords(3, size(wp%input%geo%coords))) allocate(species(size(wp%input%geo%species))) - coords(:,:) = wp%input%geo%coords(:,:) - species(:) = wp%input%geo%species(:) + coords(:,:) = wp%input%geo%coords + species(:) = wp%input%geo%species deallocate(wp%input%geo%coords) deallocate(wp%input%geo%species) allocate(wp%input%geo%coords(3, nBox * wp%input%geo%nAtom)) @@ -117,9 +112,9 @@ program waveplot do i3 = 0, wp%opt%repeatBox(3) - 1 shift(:) = matmul(wp%opt%boxVecs, real([i1, i2, i3], dp)) do iAtom = 1, wp%input%geo%nAtom - wp%input%geo%coords(:,ind+iAtom) = coords(:,iAtom) + shift(:) + wp%input%geo%coords(:,ind+iAtom) = coords(:,iAtom) + shift end do - wp%input%geo%species(ind+1:ind+wp%input%geo%nAtom) = species(:) + wp%input%geo%species(ind+1:ind+wp%input%geo%nAtom) = species ind = ind + wp%input%geo%nAtom end do end do @@ -131,51 +126,61 @@ program waveplot end do end if - write(stdout, "(A)") "Origin" - write(stdout, "(2X,3(F0.5,1X))") wp%opt%origin(:) - write(stdout, "(A)") "Box" + write(stdOut, "(A)") "Origin" + write(stdOut, "(2X,3(F0.5,1X))") wp%opt%origin + write(stdOut, "(A)") "Box" do i1 = 1, 3 - write(stdout, "(2X,3(F0.5,1X))") wp%opt%boxVecs(:, i1) + write(stdOut, "(2X,3(F0.5,1X))") wp%opt%boxVecs(:, i1) end do - write(stdout, "(A)") "Spatial resolution [1/Bohr]:" - write(stdout, "(2X,3(F0.5,1X))") 1.0_dp / sqrt(sum(wp%loc%gridVec**2, dim=1)) - write(stdout, *) + write(stdOut, "(A)") "Spatial resolution [1/Bohr]:" + write(stdOut, "(2X,3(F0.5,1X))") 1.0_dp / norm2(wp%loc%gridVec, dim=1) + write(stdOut, *) ! Create density superposition of the atomic orbitals. Occupation is distributed equally on ! orbitals with the same angular momentum. if (wp%opt%tCalcAtomDens) then allocate(atomicChrg(wp%opt%nPoints(1), wp%opt%nPoints(2), wp%opt%nPoints(3), 1)) allocate(orbitalOcc(wp%input%nOrb, 1)) - ind = 1 - do iAtom = 1, wp%input%geo%nAtom - iSpecies = wp%input%geo%species(iAtom) - do iAng = 1, size(wp%basis%basis(iSpecies)%angMoms) - mAng = 2 * wp%basis%basis(iSpecies)%angMoms(iAng) + 1 - orbitalOcc(ind:ind + mAng - 1,1) = wp%basis%basis(iSpecies)%occupations(iAng)& - & / real(mAng, dp) - ind = ind + mAng + if (env%tGlobalLead) then + ind = 1 + do iAtom = 1, wp%input%geo%nAtom + iSpecies = wp%input%geo%species(iAtom) + do iAng = 1, size(wp%basis%basis(iSpecies)%angMoms) + mAng = 2 * wp%basis%basis(iSpecies)%angMoms(iAng) + 1 + orbitalOcc(ind:ind + mAng - 1,1) = wp%basis%basis(iSpecies)%occupations(iAng)& + & / real(mAng, dp) + ind = ind + mAng + end do end do - end do - call getValue(wp%loc%molorb, wp%opt%gridOrigin, wp%loc%gridVec, orbitalOcc, atomicChrg,& - & addDensities=.true.) - sumAtomicChrg = sum(atomicChrg) * wp%loc%gridVol - buffer(:,:,:) = atomicChrg(:,:,:, 1) + call getValue(wp%loc%molorb, wp%opt%gridOrigin, wp%loc%gridVec, orbitalOcc, atomicChrg,& + & addDensities=.true.) + sumAtomicChrg = sum(atomicChrg) * wp%loc%gridVol + buffer(:,:,:) = atomicChrg(:,:,:, 1) - if (wp%opt%tVerbose) then - write(stdout, "('Total charge of atomic densities:',F12.6,/)") sumAtomicChrg - end if - if (wp%opt%tPlotAtomDens) then - write (comments(2), 9989) wp%input%identity -9989 format('Calc-Id:',I11,', atomdens') - fileName = "wp-atomdens.cube" - call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + if (wp%opt%tVerbose) then + write(stdOut, "('Total charge of atomic densities:',F12.6,/)") sumAtomicChrg + end if + if (wp%opt%tPlotAtomDens) then + write(comments(2), 9989) wp%input%identity +9989 format('Calc-Id:',I11,', atomdens') + fileName = "wp-atomdens.cube" + call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + if (ioStat /= 0) then + call error("Error while writing file '" // trim(fileName) // "'.") + else + write(stdOut, "(A)") "File '" // trim(fileName) // "' written" + end if + end if end if + #:if WITH_MPI + call mpifx_bcast(env%mpi%globalComm, atomicChrg) + call mpifx_bcast(env%mpi%globalComm, sumAtomicChrg) + #:endif end if if (wp%opt%tVerbose) then - write(stdout, "(/,A5,' ',A6,' ',A6,' ',A7,' ',A11,' ',A11)") "Spin", "KPoint", "State",& + write(stdOut, "(/,A5,' ',A6,' ',A6,' ',A7,' ',A11,' ',A11)") "Spin", "KPoint", "State",& & "Action", "Norm", "W. Occup." end if @@ -184,7 +189,7 @@ program waveplot call invert33(invBoxVecs, wp%opt%boxVecs) call invert33(recVecs2p, wp%input%geo%latVecs) recVecs2p = reshape(recVecs2p, [3, 3], order=[2, 1]) - call wp%boundaryCond%foldCoordsToCell(wp%input%geo%coords(:,:), wp%input%geo%latVecs) + call wp%boundaryCond%foldCoordsToCell(wp%input%geo%coords, wp%input%geo%latVecs) end if ! Fill the box with atoms @@ -192,22 +197,20 @@ program waveplot ! Shifting plotted region by integer lattice vectors, to have its center as close to the center ! of the lattice unit cell as possible. cellMiddle(:) = 0.5_dp * sum(wp%input%geo%latVecs, dim=2) - boxMiddle(:) = wp%opt%origin(:) + 0.5_dp * sum(wp%opt%boxVecs, dim=2) + boxMiddle(:) = wp%opt%origin + 0.5_dp * sum(wp%opt%boxVecs, dim=2) ! Workaround for intel 2021 ICE, replacing matmul(boxMiddle - cellMiddle, recVecs2p) shift(:) = boxMiddle - cellMiddle frac(:) = matmul(shift, recVecs2p) - wp%opt%origin(:) = wp%opt%origin(:) - matmul(wp%input%geo%latVecs, real(anint(frac), dp)) - wp%opt%gridOrigin(:) = wp%opt%gridOrigin(:)& - & - matmul(wp%input%geo%latVecs, real(anint(frac), dp)) + wp%opt%origin(:) = wp%opt%origin - matmul(wp%input%geo%latVecs, real(anint(frac), dp)) + wp%opt%gridOrigin(:) = wp%opt%gridOrigin - matmul(wp%input%geo%latVecs, real(anint(frac), dp)) ! We need all cells around, which could contain atoms in the sphere, drawn from the center of ! the unit cell, containing the entire plotted region mDist = 0.0_dp do i1 = 0, 1 do i2 = 0, 1 do i3 = 0, 1 - cubeCorner(:) = wp%opt%origin(:) & - &+ matmul(wp%opt%boxVecs, real([i1, i2, i3], dp)) - dist = sqrt(sum((cubeCorner - cellMiddle)**2)) + cubeCorner(:) = wp%opt%origin + matmul(wp%opt%boxVecs, real([i1, i2, i3], dp)) + dist = norm2(cubeCorner - cellMiddle) mDist = max(dist, mDist) end do end do @@ -238,8 +241,8 @@ program waveplot end if ! Calculate the molecular orbitals and write them to the disk - iSpinOld = 1 tFinished = .false. + hasIoError = .false. lpStates: do while (.not. tFinished) ! Get the next grid and its parameters if (wp%input%tRealHam) then @@ -251,129 +254,159 @@ program waveplot iKPoint = levelIndex(2) iSpin = levelIndex(3) - ! Save spin up density before processing first level for spin down - if (wp%opt%tPlotTotSpin .and. (iSpinOld /= iSpin)) then - iSpinOld = iSpin - spinUp(:,:,:) = totChrg(:,:,:) - end if - ! Build charge if needed for total charge or was explicitely required tPlotLevel = any(wp%opt%plottedSpins == iSpin) & &.and. any(wp%opt%plottedKPoints == iKPoint) .and. any(wp%opt%plottedLevels == iLevel) if (wp%opt%tCalcTotChrg .or. (tPlotLevel .and. (wp%opt%tPlotChrg .or. wp%opt%tPlotChrgDiff)))& & then if (wp%input%tRealHam) then - buffer = gridValReal**2 + buffer(:,:,:) = gridValReal**2 else - buffer = abs(gridValCmpl)**2 + buffer(:,:,:) = abs(gridValCmpl)**2 end if if (wp%opt%tCalcTotChrg) then - totChrg(:,:,:) = totChrg(:,:,:) + wp%input%occupations(iLevel, iKPoint, iSpin)& - & * buffer(:,:,:) + totChrg(:,:,:) = totChrg + wp%input%occupations(iLevel, iKPoint, iSpin) * buffer end if sumChrg = sum(buffer) * wp%loc%gridVol if (wp%opt%tVerbose) then - write(stdout, "(I5,I7,I7,A8,F12.6,F12.6)") iSpin, iKPoint, iLevel, "calc", sumChrg,& + write(*, "(I5,I7,I7,A8,F12.6,F12.6)") iSpin, iKPoint, iLevel, "calc", sumChrg,& & wp%input%occupations(iLevel, iKPoint, iSpin) end if end if + ! Save spin up density before processing first level for spin down + if (wp%opt%tPlotTotSpin .and. (iSpin == 1)) then + spinUp(:,:,:) = spinUp + wp%input%occupations(iLevel, iKPoint, iSpin) * buffer + end if + ! Build and dump desired properties of the current level if (tPlotLevel) then if (wp%opt%tPlotChrg) then - write (comments(2), 9990) wp%input%identity, iSpin, iKPoint, iLevel + write(comments(2), 9990) wp%input%identity, iSpin, iKPoint, iLevel 9990 format('Calc-Id:',I11,', Spin:',I2,', K-Point:',I6,', State:',I6, ', abs2') fileName = "wp-" // i2c(iSpin) // "-" // i2c(iKPoint) // "-" // i2c(iLevel) // "-abs2.cube" call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + hasIoError = hasIoError .or. ioStat /= 0 + if (ioStat == 0) write(*, "(A)") "File '" // trim(fileName) // "' written" end if if (wp%opt%tPlotChrgDiff) then - buffer = buffer - (sumChrg / sumAtomicChrg) * atomicChrg(:,:,:,1) - write (comments(2), 9995) wp%input%identity, iSpin, iKPoint, iLevel + buffer(:,:,:) = buffer - (sumChrg / sumAtomicChrg) * atomicChrg(:,:,:,1) + write(comments(2), 9995) wp%input%identity, iSpin, iKPoint, iLevel 9995 format('Calc-Id:',I11,', Spin:',I2,', K-Point:',I6,', State:',I6, ', abs2diff') fileName = "wp-" // i2c(iSpin) // "-" // i2c(iKPoint) // "-" // i2c(iLevel) //& & "-abs2diff.cube" call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + hasIoError = hasIoError .or. ioStat /= 0 + if (ioStat == 0) write(*, "(A)") "File '" // trim(fileName) // "' written" end if if (wp%opt%tPlotReal) then if (wp%input%tRealHam) then - buffer = gridValReal + buffer(:,:,:) = gridValReal else - buffer = real(gridValCmpl, dp) + buffer(:,:,:) = real(gridValCmpl, dp) end if - write (comments(2), 9991) wp%input%identity, iSpin, iKPoint, iLevel + write(comments(2), 9991) wp%input%identity, iSpin, iKPoint, iLevel 9991 format('Calc-Id:',I11,', Spin:',I2,', K-Point:',I6,', State:',I6, ', real') fileName = "wp-" // i2c(iSpin) // "-" // i2c(iKPoint) // "-" // i2c(iLevel) // "-real.cube" call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + hasIoError = hasIoError .or. ioStat /= 0 + if (ioStat == 0) write(*, "(A)") "File '" // trim(fileName) // "' written" end if if (wp%opt%tPlotImag) then - buffer = aimag(gridValCmpl) - write (comments(2), 9992) wp%input%identity, iSpin, iKPoint, iLevel + buffer(:,:,:) = aimag(gridValCmpl) + write(comments(2), 9992) wp%input%identity, iSpin, iKPoint, iLevel 9992 format('Calc-Id:',I11,', Spin:',I2,', K-Point:',I6,', State:',I6, ', imag') fileName = "wp-" // i2c(iSpin) // "-" // i2c(iKPoint) // "-" // i2c(iLevel) // "-imag.cube" call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + hasIoError = hasIoError .or. ioStat /= 0 + if (ioStat == 0) write(*, "(A)") "File '" // trim(fileName) // "' written" end if end if end do lpStates +#:if WITH_MPI + call mpifx_allreduceip(env%mpi%globalComm, hasIoError, MPI_LOR) +#:endif + + if (env%tGlobalLead .and. hasIoError) then + call error("At least one of the processes encountered an I/O error.") + end if + ! Dump total charge, if required if (wp%opt%tCalcTotChrg) then + #:if WITH_MPI + call mpifx_allreduceip(env%mpi%globalComm, totChrg, MPI_SUM) + #:endif sumTotChrg = sum(totChrg) * wp%loc%gridVol end if - if (wp%opt%tPlotTotChrg) then - write (comments(2), 9993) wp%input%identity + if (env%tGlobalLead .and. wp%opt%tPlotTotChrg) then + write(comments(2), 9993) wp%input%identity 9993 format('Calc-Id:',I11,', abs2') fileName = "wp-abs2.cube" call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & totChrg, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & totChrg, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + if (ioStat /= 0) then + call error("Error while writing file '" // trim(fileName) // "'.") + else + write(stdOut, "(A)") "File '" // trim(fileName) // "' written" + end if if (wp%opt%tVerbose) then - write(stdout, "(/,'Total charge:',F12.6,/)") sumTotChrg + write(stdOut, "(/,'Total charge:',F12.6,/)") sumTotChrg end if end if ! Dump total charge difference - if (wp%opt%tPlotTotDiff) then - buffer = totChrg - (sumTotChrg / sumAtomicChrg) * atomicChrg(:,:,:,1) - write (comments(2), 9994) wp%input%identity + if (env%tGlobalLead .and. wp%opt%tPlotTotDiff) then + buffer(:,:,:) = totChrg - (sumTotChrg / sumAtomicChrg) * atomicChrg(:,:,:,1) + write(comments(2), 9994) wp%input%identity 9994 format('Calc-Id:',I11,', abs2diff') fileName = 'wp-abs2diff.cube' call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + if (ioStat /= 0) then + call error("Error while writing file '" // trim(fileName) // "'.") + else + write(stdOut, "(A)") "File '" // trim(fileName) // "' written" + end if end if - ! Dump spin polarisation +#:if WITH_MPI + ! Collect spin polarisation if (wp%opt%tPlotTotSpin) then - buffer = 2.0_dp * spinUp - totChrg - write (comments(2), 9996) wp%input%identity + call mpifx_allreduceip(env%mpi%globalComm, spinUp, MPI_SUM) + end if +#:endif + + if (env%tGlobalLead .and. wp%opt%tPlotTotSpin) then + buffer(:,:,:) = 2.0_dp * spinUp - totChrg + write(comments(2), 9996) wp%input%identity 9996 format('Calc-Id:',I11,', spinpol') fileName = 'wp-spinpol.cube' call writeCubeFile(wp%input%geo, wp%aNr%atomicNumbers, wp%loc%gridVec, wp%opt%gridOrigin,& - & buffer, fileName, comments, wp%opt%repeatBox) - write(stdout, "(A)") "File '" // trim(fileName) // "' written" + & buffer, fileName, comments=comments, repeatBox=wp%opt%repeatBox, ioStat=ioStat) + if (ioStat /= 0) then + call error("Error while writing file '" // trim(fileName) // "'.") + else + write(stdOut, "(A)") "File '" // trim(fileName) // "' written" + end if end if -#:if WITH_MPI - call mpifx_finalize() -#:endif + call env%destruct() + call destructGlobalEnv() contains !> Writes a 3D function as cube file. subroutine writeCubeFile(geo, atomicNumbers, gridVecs, origin, gridVal, fileName, comments,& - & repeatBox) + & repeatBox, ioStat) !> Geometry information about the structure type(TGeometry), intent(in) :: geo @@ -399,12 +432,16 @@ subroutine writeCubeFile(geo, atomicNumbers, gridVecs, origin, gridVal, fileName !> How often the grid should be repeated along the direction of the grid vectors integer, intent(in), optional :: repeatBox(:) + !> Error status + integer, intent(out), optional :: ioStat + integer, parameter :: bufferSize = 6 real(dp) :: buffer(bufferSize) character(len=*), parameter :: formBuffer = "(6E16.8)" integer :: rep(3) integer :: ii, i1, i2, i3, ir1, ir2, ir3 type(TFileDescr) :: fd + integer :: ioStat_ @:ASSERT(size(atomicNumbers) == size(geo%speciesNames)) @:ASSERT(all(shape(gridVecs) == [3, 3])) @@ -422,26 +459,33 @@ subroutine writeCubeFile(geo, atomicNumbers, gridVecs, origin, gridVal, fileName #:endblock DEBUG_CODE if (present(repeatBox)) then - rep(:) = repeatBox(:) + rep(:) = repeatBox else rep(:) = [1, 1, 1] end if - call openFile(fd, fileName, mode="w") + call openFile(fd, fileName, mode="w", iostat=ioStat_, parallelWriting=.true.) + + ! hand over the error status, if asked for + if (present(ioStat)) ioStat = ioStat_ + + if (ioStat_ /= 0) then + call warning("Error while opening file '" // trim(fileName) // "'.") + return + end if if (present(comments)) then - write (fd%unit, "(A)") trim(comments(1)) - write (fd%unit, "(A)") trim(comments(2)) + write(fd%unit, "(A)") trim(comments(1)) + write(fd%unit, "(A)") trim(comments(2)) else - write (fd%unit, "(A)") "Made by waveplot" - write (fd%unit, *) + write(fd%unit, "(A)") "Made by waveplot" + write(fd%unit, *) end if - write (fd%unit,"(I5,3F12.6)") geo%nAtom, origin(:) - write (fd%unit,"(I5,3F12.6)") rep(1) * size(gridVal, dim=1), gridVecs(:,1) - write (fd%unit,"(I5,3F12.6)") rep(2) * size(gridVal, dim=2), gridVecs(:,2) - write (fd%unit,"(I5,3F12.6)") rep(3) * size(gridVal, dim=3), gridVecs(:,3) + write(fd%unit,"(I5,3F12.6)") geo%nAtom, origin + write(fd%unit,"(I5,3F12.6)") rep(1) * size(gridVal, dim=1), gridVecs(:,1) + write(fd%unit,"(I5,3F12.6)") rep(2) * size(gridVal, dim=2), gridVecs(:,2) + write(fd%unit,"(I5,3F12.6)") rep(3) * size(gridVal, dim=3), gridVecs(:,3) do ii = 1, geo%nAtom - write (fd%unit, "(I5,4F12.6)") atomicNumbers(geo%species(ii)), 0.0_dp, & - &geo%coords(:, ii) + write(fd%unit, "(I5,4F12.6)") atomicNumbers(geo%species(ii)), 0.0_dp, geo%coords(:, ii) end do do ir1 = 1, rep(1) @@ -453,17 +497,18 @@ subroutine writeCubeFile(geo, atomicNumbers, gridVecs, origin, gridVal, fileName ii = mod(i3 - 1, bufferSize) + 1 buffer(ii) = gridVal(i1, i2, i3) if (ii == bufferSize) then - write (fd%unit,formBuffer) real(buffer) + write(fd%unit, formBuffer) real(buffer) end if end do if (ii /= bufferSize) then - write (fd%unit, "(" // i2c(ii) // "E16.8)") real(buffer(:ii)) + write(fd%unit, "(" // i2c(ii) // "E16.8)") real(buffer(:ii)) end if end do end do end do end do end do + call closeFile(fd) end subroutine writeCubeFile