Bug fixes to make DOS=t, ndir=-3 mode work

cdn/slab_dim.f90

@@ -47,16 +47,16 @@ CONTAINS
     nz = 1
     na = 0
     DO i=1,atoms%ntype
-       DO j=1,atoms%neq(i)
+       equivAtomsLoop: DO j=1,atoms%neq(i)
           na = na + 1
           zs = atoms%pos(3,na)
           DO iz=1,nz
-             IF(ABS(zs-znz(iz)).LT.epsz) CYCLE
-          ENDDO
+             IF(ABS(zs-znz(iz)).LT.epsz) CYCLE equivAtomsLoop
+          END DO
           nz = nz+1
           znz(nz) = zs
-       ENDDO
-    ENDDO
+       END DO equivAtomsLoop
+    END DO
     nsld = nz
     IF(nsld>atoms%nat)   CALL juDFT_error("nsld.GT.atoms%nat ",calledby="slab_dim")
     !

dos/evaldos.f90

@@ -73,10 +73,9 @@
       DATA spin12/'.1' , '.2'/
       DATA ch_mcd/'.+' , '.-' , '.0'/
 !
+      qdim = lmax*atoms%ntype+3
       l_orbcomp = .false.
-      IF  (banddos%ndir.NE.-3) THEN
-        qdim = lmax*atoms%ntype+3
-      ELSE
+      IF (banddos%ndir.EQ.-3) THEN
         qdim = 2*nsld 
         n_orb = 0
         INQUIRE(file='orbcomp',exist=l_orbcomp)
@@ -170,12 +169,12 @@
             ALLOCATE( qis(dimension%neigd),qvlay(dimension%neigd,vacuum%layerd,2))
             CALL read_eig(eig_id,k,jspin,bk=vk(:,k),wk=wt(k),neig=nevk(k),eig=ev(:,k))
             CALL read_dos(eig_id,k,jspin,qal_tmp,qvac,qis,qvlay,qstars,ksym,jsym,mcd,qintsl,qmtsl,qmtp,orbcomp)
-            qal(1:lmax*atoms%ntype,:,k)=reshape(qal_tmp,(/lmax*atoms%ntype,size(qal_tmp,3)/))
-            qal(lmax*atoms%ntype+2,:,k)=qvac(:,1)         ! vacuum 1
-            qal(lmax*atoms%ntype+3,:,k)=qvac(:,2)         ! vacuum 2
-            qal(lmax*atoms%ntype+1,:,k)=qis              ! interstitial
-            DEALLOCATE( ksym,jsym )
-            IF (l_orbcomp)THEN
+            IF (.NOT.l_orbcomp) THEN
+               qal(1:lmax*atoms%ntype,:,k)=reshape(qal_tmp,(/lmax*atoms%ntype,size(qal_tmp,3)/))
+               qal(lmax*atoms%ntype+2,:,k)=qvac(:,1)         ! vacuum 1
+               qal(lmax*atoms%ntype+3,:,k)=qvac(:,2)         ! vacuum 2
+               qal(lmax*atoms%ntype+1,:,k)=qis              ! interstitial
+            ELSE
                IF (n_orb == 0) THEN
                   qal(1:nsld,:,k)        = qintsl(:,:)
                   qal(nsld+1:2*nsld,:,k) = qmtsl(:,:)
@@ -183,13 +182,14 @@
                   DO i = 1, 23
                      DO l = 1, nevk(k)
                         qal(i,l,k) = orbcomp(l,i,n_orb)*qmtp(l,n_orb)/10000.
-                     ENDDO
+                     END DO
                      DO l = nevk(k)+1, dimension%neigd
                         qal(i,l,k) = 0.0
-                     ENDDO
-                  ENDDO
-               ENDIF
-            ENDIF
+                     END DO
+                  END DO
+               END IF
+            END IF
+            DEALLOCATE( ksym,jsym )
             DEALLOCATE( orbcomp,qintsl,qmtsl,qmtp,qvac,qis,qal_tmp,qvlay)
             ntb = max(ntb,nevk(k))
 !

main/fleur.F90

@@ -519,8 +519,14 @@
                                ! so transform the eig file such that each kpt is ro
                                ! to its symmetry equivalent ones
                                IF ( banddos%dos .AND. banddos%ndir == -3 ) THEN
-                                  ! generate pointer from gpts at one kpt to a diffe
 #ifdef CPP_NEVER
+                                  ! This code section has to be reimplemented to enable
+                                  ! calculations of the orbital decomposed DOS with
+                                  ! eigenvectors up to this point only calculated for the
+                                  ! IBZ. Here we reconstruct the eigenvectors for the full
+                                  ! BZ from this subset.
+
+                                  ! generate pointer from gpts at one kpt to a diffe
                                   ALLOCATE( kpts%pntgptd(3) )
                                   CALL generate_pntgpt(&
                                        &                      dimension,obsolete,input,&
@@ -529,12 +535,13 @@
                                   CALL rotate_eig(&
                                        &                      kpts,dimension,atoms,&
                                        &                      sym,mpi)
-#endif
+
                                   DEALLOCATE( kpts%pntgptd, kpts%pntgpt )
                                   ! this change is sufficient to modify fermie and c
                                   ! the enlarged kpt mesh
                                   kpts%nkpt    = kpts%nkptf
                                   kpts%nkpt   = kpts%nkptf
+#endif
                                END IF
 
                             ENDIF

main/fleur_init.F90

@@ -508,7 +508,26 @@
              END DO
           ELSE
              IF ( banddos%dos .AND. banddos%ndir == -3 ) THEN
+                WRITE(*,*) 'Recalculating k point grid to cover the full BZ.'
                 CALL gen_bz(kpts,sym)
+                kpts%nkpt = kpts%nkptf
+                DEALLOCATE(kpts%bk,kpts%wtkpt)
+                ALLOCATE(kpts%bk(3,kpts%nkptf),kpts%wtkpt(kpts%nkptf))
+                kpts%bk(:,:) = kpts%bkf(:,:)
+                IF (kpts%nkpt3(1)*kpts%nkpt3(2)*kpts%nkpt3(3).NE.kpts%nkptf) THEN
+                   IF(kpts%l_gamma) THEN
+                      kpts%wtkpt = 1.0 / (kpts%nkptf-1)
+                      DO i = 1, kpts%nkptf
+                         IF(ALL(kpts%bk(:,i).EQ.0.0)) THEN
+                            kpts%wtkpt(i) = 0.0
+                         END IF
+                      END DO
+                   ELSE
+                      CALL juDFT_error("nkptf does not match product of nkpt3(i).",calledby="fleur_init")
+                   END IF
+                ELSE
+                   kpts%wtkpt = 1.0 / kpts%nkptf
+                END IF
              END IF
              ALLOCATE(hybrid%map(0,0),hybrid%tvec(0,0,0),hybrid%d_wgn2(0,0,0,0))
              hybrid%l_calhf   = .FALSE.