dmftproj: add option to specify band indices

This adds another option and a mode flag in dmftproj:
 * third value in window line defines now the mode
 * new option to provide an energy window where all
   bands which are within the window (at least at one k-point)
   are taken into account for the projectors.
 * updates to documention to reflects those changes

doc/ChangeLog.md

@@ -8,6 +8,7 @@ DFTTools Version 3.0.0 is a major release that
 * is now aligned with the general [app4triqs](https://github.com/TRIQS/app4triqs) application skeleton
 * brings a major rework of the VASP interface, including thorough documentation, tutorials, a new Hamiltonian mode, the option to select bands instead of an energy window, and many small bugfixes.
 * brings a major update of the block structure functionalities especially for SOC calculations, with detailed documentation and tutorials. Allows more control over the block structure coming from DFT, cutting out certain orbitals or throwing away off-diagonal elements when preparing input for the solver.
+* New option in dmftproj to select the projection window using band indices instead of energie
 
 Restructuring
 -------------
@@ -65,7 +66,7 @@ other changes:
 
 
 Thanks to all commit-contributors (in alphabetical order):
-Markus Aichhorn, Alexander Hampel, Oleg Peil, Hermann Schnait, Malte Schueler, Nils Wentzell, Manuel Zingl
+Markus Aichhorn, Alexander Hampel, Gernot Kraberger, Oleg Peil, Hermann Schnait, Malte Schueler, Nils Wentzell, Manuel Zingl
 
 
 Version 2.2.1

doc/guide/conv_wien2k.rst

@@ -68,9 +68,30 @@ is given by the following 3 to 5 lines:
 
 These lines have to be repeated for each inequivalent atom.
 
-The last line gives the energy window, relative to the Fermi energy,
+The last line gives the lower and upper limit of the energy window,
-that is used for the projective Wannier functions. Note that, in
+relative to the Fermi energy, which is used for the projective Wannier functions.
-accordance with Wien2k, we give energies in Rydberg units!
+Note that, in accordance with Wien2k, we give energies in Rydberg units!
+
+The third number is an optional flag to switch between different modes:
+
+#. 0: The projectors are constructed for the given energy window. The number
+   of bands within the window is usually different at each k-point which
+   will be reflected by the projectors, too. This is the default mode
+   which is also used if no mode flag is provided.
+#. 1: The lowest and highest band indices within the given energy window
+   are calculated. The resulting indices are used at all k-points.
+   Bands which fall within the window only in some parts of the Brillouin zone
+   are fully taken into account. Keep in mind that a different set of k-points
+   or the -band option can change the lower or upper index. This can be avoided
+   by using mode 2.
+#. 2: In this mode the first two values of the line are interpreted as lower
+   and upper band indices to be included in the projective subspace. For example,
+   if the line reads `21 23 2`, bands number 21, 22 and 23 are included at all
+   k-points. For SrVO3 this corresponds to the t2g bands around the Fermi energy.
+   The lowest possible index is 1. Note that the indices need to be provided as integer.
+
+In all modes the used energy range, i.e. band range, is printed to the
+:program:`dmftproj` output.
 
 After setting up the :file:`case.indmftpr` input file, you run:
 

doc/tutorials/images_scripts/SrVO3.indmftpr

@@ -12,4 +12,4 @@ cubic            ! choice of angular harmonics
 complex          ! choice of angular harmonics
 1 1 0 0          ! l included for each sort
 0 0 0 0          ! If split into ireps, gives number of ireps. for a given orbital (otherwise 0)
--0.11 0.14
+-0.11 0.14  0    ! lower energy window limit, upper energy window limit, mode (0,1,2)

fortran/dmftproj/dmftproj.f

@@ -48,7 +48,7 @@ C
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: lnreps
         INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: correps
         INTEGER :: isrt, ie, l, m, isym, jatom
-        INTEGER :: lm, ik, ilo, ib, iatom, imu
+        INTEGER :: lm, ik, ilo, ib, iatom, imu, io
         INTEGER :: idum, i1, i2
         INTEGER :: m1, m2, lm1, lm2
         INTEGER :: is, irep, nbrep
@@ -57,6 +57,7 @@ C
         LOGICAL :: ifcorr
         REAL(KIND=8) :: fdum, rtetr
         REAL(KIND=8),PARAMETER :: Elarge=1d6
+        character(len=120) line
 C ================================
 C Processing of the command line :
 C ================================
@@ -146,7 +147,7 @@ C With SO, the number of ireps must not exceed 2*(2*l+1).
                WRITE(buf,'(a)')'END OF THE PRGM'
                CALL printout(0)
                STOP
-              END IF
+              ENDIF
              ELSE
 C Without SO, the number of ireps must not exceed (2*l+1).
               IF(lnreps(l,isrt).gt.(2*l+1)) THEN
@@ -157,7 +158,7 @@ C Without SO, the number of ireps must not exceed (2*l+1).
                WRITE(buf,'(a)')'END OF THE PRGM'
                CALL printout(0)
                STOP
-              END IF
+              ENDIF
              ENDIF
 C ---------------------------------------------------------------------------------------
 C
@@ -177,24 +178,83 @@ C lsort = index for each orbital (0 : not include / 1 : include / 2 : correlated
 C lnreps = number of irreducible representations for each orbital, table from 0 to lmax, from 1 to nsort (temporary variables)
 C correps = index for each irreducible representations of the correlated orbital, table from 1 to lnreps(l,isrt), from 0 to lmax, from 1 to nsort (temporary variable)
 C ifSOflag = table of correspondance isort -> optionSO (1 or 0). Only used for isort with correlated orbitals
-        READ(iuinp,*) e_bot,e_top
+
-C e_bot, e_top : lower and upper limits of the energy window
+        READ(iuinp,'(A)',iostat=io) line
+C Try reading the energies/bandindices and the proj_mode
+        READ(line,*,iostat=io) e_bot, e_top, proj_mode
+C If it fails we know that we are dealing with an older version of the indmftpr file
+C with only 2 values on the window line. proj_mode = 0.
+        IF(io.ne.0) THEN
+         proj_mode = 0
+         READ(line,*,iostat=io) e_bot, e_top
+         IF(io.ne.0) THEN
+          WRITE(buf,'(a,a)')' The energy window line',
+     &     ' is ill-defined.'
+          CALL printout(0)
+          STOP
+          WRITE(buf,'(a)')'END OF THE PRGM'
+         ENDIF
+        ENDIF
+
+C ---------------------------------------------------------------------------------------
+C proj_mode:
+C 0: use energy window for projection
+C 1: use all band indices present in the given energy window
+C    (same number of bands at all kpoints)
+C 2: use given band indices (same number of bands at all kpoints)
+C ---------------------------------------------------------------------------------------
+
+C ---------------------------------------------------------------------------------------
+C e_bot, e_top : lower/upper energy limits of window (used in mode 0)
+C b_bot, b_top : lower/upper band index of window (used in mode 2)
+C In mode 1 e_bot/e_top are provided in the input file and then
+C translated into b_bot/b_top
+C ---------------------------------------------------------------------------------------
 C
 C ---------------------------------------------------------------------------------------
-C Interruption of the prgm if the energy window is not well-defined.
+C Interruption of the prgm if the energy/band window or proj_mode is not well-defined.
-C -------------------------
+C ---------------------------------------------------------------------------------------
 C
-        IF(e_bot.gt.e_top) THEN
+        IF((proj_mode.lt.0) .or. (proj_mode.gt.2)) THEN
+         WRITE(buf,'(a,a)')' The energy window mode (3rd value)',
+     &     ' must be 0,1 or 2.'
+         CALL printout(0)
+         WRITE(buf,'(a)')'END OF THE PRGM'
+         CALL printout(0)
+         STOP
+        ENDIF
+        
+        IF(proj_mode==0) THEN
+                b_bot = 0
+                b_top = 0
+        ELSEIF(proj_mode==1) THEN
+                b_bot = 1e3
+                b_top = 1
+        ELSEIF(proj_mode==2) THEN
+                b_bot = INT(e_bot)
+                b_top = INT(e_top)
+                e_bot = 0.0
+                e_top = 0.0
+        ENDIF
+
+        IF((proj_mode.lt.2) .and. (e_bot.gt.e_top)) THEN
          WRITE(buf,'(a,a)')' The energy window ',
      &     ' is ill-defined.'
          CALL printout(0)
          WRITE(buf,'(a)')'END OF THE PRGM'
          CALL printout(0)
          STOP
-        END IF
+        ENDIF
-C ---------------------------------------------------------------------------------------
+
-C
+        IF((proj_mode==2) .and. (b_bot.gt.b_top)) THEN
-C =====================================================================
+          WRITE(buf,'(a,a)')' The k-point index window ',
+     &     ' is ill-defined.'
+          CALL printout(0)
+          WRITE(buf,'(a)')'END OF THE PRGM'
+          CALL printout(0)
+          STOP
+        ENDIF
+        
 C Writing in the output file case.outdmftpr the previous informations :
 C =====================================================================
         WRITE(buf,'(a,a)')'Welcome in DMFTPROJ: ',
@@ -371,13 +431,6 @@ C the field crorb%ifSOflag = boolean (if SO are used or not)
           ENDDO      ! End of the l loop
         ENDDO        ! End of the isrt loop
 C
-C Printing the size of the Energy window
-        CALL printout(0)
-        WRITE(buf,'(2(a,f10.5),a)')
-     &     'The Eigenstates are projected in an energy window from ',
-     &        e_bot,' Ry  to  ',e_top,' Ry around the Fermi level.'
-        CALL printout(1)
-C
 C =======================================================================================
 C Reading of the transformation matrices from the complex to the required angular basis :
 C =======================================================================================
@@ -556,7 +609,7 @@ C
                WRITE(buf,'(a)')'END OF THE PRGM'
                CALL printout(0)
                STOP
-              END IF
+              ENDIF
 C ---------------------------------------------------------------------------------------
 C
 C It is assumed in the following that nLO is 0 or 1.
@@ -644,6 +697,49 @@ C Printing in the file case.outdmftpr the fermi level (in Rydberg)
      &    'with respect to this value. (E_Fermi is now 0 Ry)'
         CALL printout(1)
 C
+C ---------------------------------------------------------------------------------------
+C If proj_mode=1 find now the lowest and highes band index
+C ---------------------------------------------------------------------------------------
+C
+        IF(proj_mode==1) THEN
+         DO is=1,ns
+           DO ik=1,nk
+             DO ib=kp(ik,is)%nbmin,kp(ik,is)%nbmax
+               IF(kp(ik,is)%eband(ib) > e_bot.AND.
+     &          kp(ik,is)%eband(ib).LE.e_top) THEN
+                IF(ib.gt.b_top) THEN
+                 b_top = ib
+                ENDIF
+                IF(ib.lt.b_bot) THEN
+                 b_bot = ib
+                ENDIF             
+               ENDIF
+             ENDDO    ! End of the ib loop
+           ENDDO      ! End of the ik loop
+         ENDDO        ! End of the is loop
+         e_top = 0.0
+         e_bot = 0.0
+        ENDIF
+
+C ---------------------------------------------------------------------------------------
+C Printing the size of the Energy window
+C ---------------------------------------------------------------------------------------
+
+        CALL printout(0)
+        IF(proj_mode==0) THEN
+         WRITE(buf,'(2(a,f10.5),a)')
+     &     'The Eigenstates are projected in an energy window from ',
+     &        e_bot,' Ry  to  ',e_top,' Ry around the Fermi level.'
+        ELSEIF(proj_mode==1) THEN
+         WRITE(buf,'(a,2(a,i3),a)')
+     &     'The Eigenstates are projected for the band indices from ',
+     &     'band Nr. ', b_bot,' to  ',b_top,'.'
+        ELSEIF(proj_mode==2) THEN         
+         WRITE(buf,'(a,2(a,i3),a)')
+     &     'The Eigenstates are projected for the band indices from ',
+     &     'band Nr. ',b_bot,' to  ',b_top,'.'
+        ENDIF
+        CALL printout(1)
 C
 C ==============================================================
 C Computation of the density matrices up to the Fermi level Ef :
@@ -657,7 +753,11 @@ C
 C ----------------------------------------
 C Setting up the projections for all bands
 C ----------------------------------------
-        CALL set_projections(-Elarge,Elarge)
+        IF(proj_mode==0) THEN
+         CALL set_projections(-Elarge,Elarge)
+        ELSE
+         CALL set_projections(1d0,1d6)
+        ENDIF
 
 
 C Elarge is an energy variable equal to 1.d6 Rydberg (very large !!!)
@@ -675,21 +775,31 @@ C
 C The calculation of Wannier projectors is performed only if correlated orbitals are included.
         IF(ncrorb.NE.0) THEN
 C 
-C =====================================================================
+C ==========================================================================
-C Computation of the charge below the lower limit e_bot of the window :
+C Computation of the charge below the lower limit e_bot/b_bot of the window :
-C =====================================================================
+C ==========================================================================
 C
          WRITE(buf,'(a)')'======================================='
          CALL printout(0)
-         WRITE(buf,'(a,a,f10.5,a)')'Computation of the total ',
+         IF(proj_mode==0) THEN
+          WRITE(buf,'(a,a,f10.5,a)')'Computation of the total ',
      &     'Charge below the lower limit of the energy window :',
      &     e_bot,' Ry'  
+         ELSE
+          WRITE(buf,'(a,a,i3)')'Computation of the total ',
+     &     'Charge below the lower band index Nr. ', b_bot  
+         ENDIF
          CALL printout(1)
 C
 C ----------------------------------------
 C Setting up the projections for all bands
 C ----------------------------------------
-         CALL set_projections(-Elarge,e_bot)
+         IF(proj_mode==0) THEN
+          CALL set_projections(-Elarge,e_bot)
+         ELSE
+C set_projections expects REAL(8)
+          CALL set_projections(1d0,REAL(b_bot-1,8))
+         ENDIF
 C
 C ---------------------------------------------------------
 C Computation of the density matrices and the total charges
@@ -698,7 +808,7 @@ C
          IF(.NOT.ifBAND) CALL density(.FAlSE.,.FALSE.,qtot,.FALSE.)
 C A simple point integration is used. 
 C The computation is performed for all the included orbitals.
-C qtot is the total charge density below e_bot. 
+C qtot is the total charge density below e_bot/b_bot. 
 C Nothing will be printed in the file case.outdmftpr apart from the total charge qtot.
 C
 C
@@ -708,15 +818,25 @@ C ============================================================
 C
          WRITE(buf,'(a)')'======================================='
          CALL printout(0)
-         WRITE(buf,'(a,a,a,f10.5,a,f10.5,a)')'Computation of the ',
+         IF(proj_mode==0) THEN
+          WRITE(buf,'(a,a,a,f10.5,a,f10.5,a)')'Computation of the ',
      &     'Occupancies and Density Matrices in the desired ',
      &     'energy window [ ',e_bot,'; ',e_top,']' 
+         ELSE
+          WRITE(buf,'(a,a,a,i3,a,i3,a)')'Computation of the ',
+     &     'Occupancies and Density Matrices in the desired ',
+     &     'band range [ ',b_bot,'; ',b_top,']' 
+         ENDIF
          CALL printout(1)
 C
 C ----------------------------------------
 C Setting up the projections for all bands
 C ----------------------------------------
-         CALL set_projections(e_bot,e_top)
+         IF(proj_mode==0) THEN
+          CALL set_projections(e_bot,e_top)
+         ELSE
+          CALL set_projections(REAL(b_bot,8),REAL(b_top,8))
+         ENDIF
 C
 C ------------------------------------------------------------------------------
 C Orthonormalization of the projectors for correlated orbitals P(icrorb,ik,is) :

fortran/dmftproj/modules.f

@@ -68,10 +68,12 @@ C     &    '/workpmc/martins/DMFTprojectors/newDMFTproj'
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: lsort
         INTEGER, DIMENSION(:), ALLOCATABLE :: ifSOflag
         INTEGER, DIMENSION(:), ALLOCATABLE :: timeflag
+        INTEGER :: b_bot, b_top, proj_mode
         LOGICAL :: ifSO, ifSP, ifBAND
         LOGICAL, DIMENSION(:), ALLOCATABLE :: notinclude
         REAL(KIND=8) :: eferm
         REAL(KIND=8) :: e_bot, e_top
+
         REAL(KIND=8), PARAMETER :: PI=3.1415926535898d0
 C New type structure basistrans
         TYPE deftrans

fortran/dmftproj/set_projections.f

@@ -24,7 +24,9 @@ c *****************************************************************************/
 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 C %%                                                                 %%
 C %% This subroutine sets up the projection matrices in the energy   %%
-C %% window [e1,e2].Two types of projection can be defined :         %%
+C %% window [e1,e2]. If proj_mode is 1 or 2 then e1 and e2 are not   %%
+C %% energies, but directly used as band indices.                    %%
+C %% Two types of projection can be defined :                        %%
 C %%   - The projectors <u_orb|ik,ib,is> for the correlated orbital  %%
 C %%     only. (orb = iatom,is,m)                                    %%
 C %%     (They are stored in the table pr_crorb)                     %%
@@ -60,48 +62,71 @@ C
 C
 C
 C ======================================================================
-C Selection of the bands which lie in the chosen energy window [e1;e2] :
+C Selection of the bands which lie in the chosen energy window [e1;e2]
+C or if proj_mode = [1,2] e1 and e2 are directly used as band indices:       
 C ======================================================================
 C
-        kp(:,:)%included=.FALSE.
+        
+        IF(proj_mode.gt.0) THEN
+C The same number of bands are included at each k-point
+         kp(:,:)%included=.TRUE.
+C Use directly e1 and e2 as band indices. Set to nbmin or
+C nbmax if too small or too large
+         DO is=1,ns
+           DO ik=1,nk
+            IF(e1 > kp(ik,is)%nbmin) THEN
+             kp(ik,is)%nb_bot=INT(e1)
+            ELSE
+             kp(ik,is)%nb_bot=kp(ik,is)%nbmin
+            ENDIF
+            IF(e2 < kp(ik,is)%nbmax) THEN
+             kp(ik,is)%nb_top=INT(e2)
+            ELSE
+             kp(ik,is)%nb_top=kp(ik,is)%nbmax
+            ENDIF
+           ENDDO
+         ENDDO
+        ELSE
+         kp(:,:)%included=.FALSE.
 C the field kp%included = boolean which is .TRUE. when there is at least one band 
 C at this k-point whose energy eignevalue is in the energy window.
-        DO is=1,ns
+         DO is=1,ns
-          DO ik=1,nk
+           DO ik=1,nk
-            included=.FALSE.
+             included=.FALSE.
-            DO ib=kp(ik,is)%nbmin,kp(ik,is)%nbmax
+             DO ib=kp(ik,is)%nbmin,kp(ik,is)%nbmax
-              IF(.NOT.included.AND.kp(ik,is)%eband(ib) > e1.AND.
+               IF(.NOT.included.AND.kp(ik,is)%eband(ib) > e1.AND.
      &          kp(ik,is)%eband(ib).LE.e2) THEN
 C If the energy eigenvalue E of the band ib at the k-point ik is such that e1 < E =< e2, 
 C then all the band with ib1>ib must be "included" in the computation and kp%nb_bot is initialized at the value ib.
-               included=.TRUE.
+                included=.TRUE.
-               kp(ik,is)%nb_bot=ib
+                kp(ik,is)%nb_bot=ib
-              ELSEIF(included.AND.kp(ik,is)%eband(ib) > e2) THEN
+               ELSEIF(included.AND.kp(ik,is)%eband(ib) > e2) THEN
 C If the energy eigenvalue E of the current band ib at the k-point ik is such that E > e2 and all the previous 
 C band are "included", then the field kp%included = .TRUE. and kp%nb_top = ib-1 (the index of the previous band)
-               kp(ik,is)%nb_top=ib-1
+                kp(ik,is)%nb_top=ib-1
-               kp(ik,is)%included=.TRUE.
+                kp(ik,is)%included=.TRUE.
-               EXIT
+                EXIT
 C The loop on the band ib is stopped, since all the bands after ib have an energy > that of ib. 
-              ELSEIF(ib==kp(ik,is)%nbmax.AND.kp(ik,is)%eband(ib)
+               ELSEIF(ib==kp(ik,is)%nbmax.AND.kp(ik,is)%eband(ib)
      &          > e1.AND.kp(ik,is)%eband(ib).LE.e2) THEN
 C If the energy eigenvalue E of the last band ib=kp%nbmax at the k-point ik is such that e1 < E =< e2 and all the 
 C previous bands are "included", then the band ib must be "included" and kp%nb_bot is initialized at the value kp%nbmax.
-               kp(ik,is)%nb_top=ib
+                kp(ik,is)%nb_top=ib
-               kp(ik,is)%included=.TRUE.
+                kp(ik,is)%included=.TRUE.
-              ENDIF
+               ENDIF
 C If the eigenvalues of the bands at the k-point ik are < e1 and included=.FALSE. 
 C of if the eigenvalues of the bands at the k-point ik are in the energy window [e1,e2] and included=.TRUE.,
 C nothing is done...
-            ENDDO    ! End of the ib loop
+             ENDDO    ! End of the ib loop
 C If all the eigenvalues of the bands at the k-point ik are not in the window, 
-C then kp%included remains at the value .FALSE. and the field kp%nb_top and kp%nb_bot are set to 0.
+C then kp%included remains at the value .FALSE. and the field kp%nb_top and kp%nb_bot are set to 0
-	    IF (.not.kp(ik,is)%included) THEN
+             IF (.not.kp(ik,is)%included) THEN
-             kp(ik,is)%nb_bot=0
+              kp(ik,is)%nb_bot=0
-             kp(ik,is)%nb_top=0 
+              kp(ik,is)%nb_top=0 
-	    ENDIF
+             ENDIF
-          ENDDO      ! End of the ik loop
+           ENDDO      ! End of the ik loop
-        ENDDO        ! End of the is loop
+         ENDDO        ! End of the is loop
+        ENDIF
 C ---------------------------------------------------------------------------------------
 C Checking of the input files if spin-polarized inputs and SO is taken into account:
 C There should not be any difference between up and dn limits for each k-point.