Fill some more arrays correctly

(in the same way as the previous two commits)

cdn/cdnval.F90

@@ -207,7 +207,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
          CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
                     jspin,lapw,noccbd,ev_list,we,eig,den,results,force%f_b8,zMat,dos)
          ! charge of each valence state in this k-point of the SBZ in the layer interstitial region of the film
-         IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
+         IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,ev_list,lapw,slab,oneD,zMat)
          ! valence density in the vacuum region
          IF (input%film) THEN
             CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,&
@@ -234,7 +234,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
 
          ! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film
          IF (l_dosNdir) THEN
-            IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
+            IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ev_list,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
 
             IF (banddos%l_orb.AND.ANY((/banddos%alpha,banddos%beta,banddos%gamma/).NE.0.0)) THEN
                CALL abcrot2(atoms,banddos,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs

cdn/eparas.f90

@@ -99,7 +99,7 @@ CONTAINS
                    DO icore = 1, mcd%ncore(n)
                       DO ipol = 1, 3
                          index = 3*(n-1) + ipol
-                         mcd%mcd(index,icore,i,ikpt,jsp)=mcd%mcd(index,icore,i,ikpt,jsp) + fac*(&
+                         mcd%mcd(index,icore,ev_list(i),ikpt,jsp)=mcd%mcd(index,icore,ev_list(i),ikpt,jsp) + fac*(&
                               suma * CONJG(mcd%m_mcd(icore,lm+1,index,1))*mcd%m_mcd(icore,lm+1,index,1)  +&
                               sumb * CONJG(mcd%m_mcd(icore,lm+1,index,2))*mcd%m_mcd(icore,lm+1,index,2)  +&
                               sumab* CONJG(mcd%m_mcd(icore,lm+1,index,2))*mcd%m_mcd(icore,lm+1,index,1)  +&

cdn/q_int_sl.F90

 MODULE m_qintsl
   USE m_juDFT
 CONTAINS
-  SUBROUTINE q_int_sl(isp,ikpt,stars,atoms,sym,cell,ne,lapw,slab,oneD,zMat)          
+  SUBROUTINE q_int_sl(isp,ikpt,stars,atoms,sym,cell,ne,ev_list,lapw,slab,oneD,zMat)          
     !     *******************************************************
     !     calculate the charge of the En(k) state 
     !     in the interstitial region of each leyer
@@ -24,6 +24,9 @@ CONTAINS
     !
     !     .. Scalar Arguments ..
     INTEGER, INTENT (IN) :: ne,isp,ikpt
+
+    INTEGER, INTENT (IN) :: ev_list(ne)
+
     !     ..
     !     .. Local Scalars ..
     REAL q1,zsl1,zsl2,qi,volsli,volintsli
@@ -105,7 +108,7 @@ CONTAINS
           DO j = 1,stars%ng3
              qi = qi + z_z(j)*stfunint(j,i)
           ENDDO
-          slab%qintsl(i,n,ikpt,isp) = qi 
+          slab%qintsl(i,ev_list(n),ikpt,isp) = qi 
        ENDDO    ! over vacuum%layers         
 
     ENDDO ! over states

cdn/q_mt_sl.f90

@@ -8,7 +8,7 @@ CONTAINS
   !
   !***********************************************************************
   !
-  SUBROUTINE q_mt_sl(jsp,atoms,nobd,ikpt,ne,skip_t,noccbd,eigVecCoeffs,usdus,slab)
+  SUBROUTINE q_mt_sl(jsp,atoms,nobd,ev_list,ikpt,ne,skip_t,noccbd,eigVecCoeffs,usdus,slab)
     USE m_types_setup
     USE m_types_usdus
     USE m_types_cdnval, ONLY: t_eigVecCoeffs, t_slab
@@ -21,6 +21,9 @@ CONTAINS
     !     .. Scalar Arguments ..
     INTEGER, INTENT (IN) :: nobd,jsp      
     INTEGER, INTENT (IN) :: ne,ikpt ,skip_t,noccbd
+
+    INTEGER, INTENT (IN) :: ev_list(nobd)
+
     !     ..
     !     .. Local Scalars ..
     INTEGER i,l,lo ,natom,nn,ntyp,nt1,nt2,m
@@ -136,7 +139,7 @@ CONTAINS
           DO ntyp = 1,atoms%ntype
              qq = qq + qmttot(ntyp,i)*slab%nmtsl(ntyp,nl)
           ENDDO
-          slab%qmtsl(nl,i,ikpt,jsp) = qq
+          slab%qmtsl(nl,ev_list(i),ikpt,jsp) = qq
        ENDDO
     ENDDO
     !        DO ntyp = 1,ntype