Extended use of k_grid object to gamma_only routines.

IO_kernel/wfreq_db.f90

@@ -41,6 +41,7 @@ MODULE wfreq_db
       USE io_push,              ONLY : io_push_bar
       USE json_module,          ONLY : json_file 
       USE constants,            ONLY : rytoev
+      USE types_bz_grid,        ONLY : k_grid
       !
       IMPLICIT NONE
       !
@@ -87,7 +88,7 @@ MODULE wfreq_db
          DEALLOCATE(ilist)
          IF( l_generate_plot ) CALL json%add('output.P.freqlist',sigma_freq(1:n_spectralf)*rytoev)
          !
-         DO iks = 1, nks 
+         DO iks = 1, k_grid%nps
             !
             WRITE(my_label_k,'(i6.6)') iks_l2g(iks)
             !

Tools/do_setup.f90

@@ -55,18 +55,17 @@ SUBROUTINE do_setup
   ! INIT PW
   !
   CALL init_pw_arrays(nbnd)
-  CALL set_iks_l2g()
   !
   CALL set_dirs()
   !
   ! INIT K, Q GRIDS
   !
-  !k_grid = bz_grid()
   CALL k_grid%init('K')
   !
-  !q_grid = bz_grid()
   CALL q_grid%init('Q')
   !
+  CALL set_iks_l2g()
+  !
   IF ( ANY ( (q_grid%ngrid(:) - k_grid%ngrid(:)) /= 0   ) ) THEN
      CALL errore( 'do_setup','q-point grid must be the same as k-point grid ',1)
   ENDIF
@@ -79,7 +78,7 @@ SUBROUTINE do_setup
   IF ( lsda ) THEN
      IF ( INT( nelup ) == 0 .AND. INT( neldw ) == 0 ) THEN
      !IF ( .NOT. two_fermi_energies ) THEN
-        DO iks = 1, nks
+        DO iks = 1, k_grid%nps
            spin = isk(iks)
            !
            SELECT CASE(spin)

Tools/get_alpha_pv.f90

@@ -20,6 +20,7 @@ FUNCTION get_alpha_pv()
   USE mp_global,              ONLY : inter_pool_comm
   USE klist,                  ONLY : nks
   USE westcom,                ONLY : nbnd_occ
+  USE types_bz_grid,          ONLY : k_grid
   !
   IMPLICIT NONE
   !
@@ -34,7 +35,7 @@ FUNCTION get_alpha_pv()
   ! Calculate ALPHA_PV
   !
   emin = et (1, 1)
-  DO iks = 1, nks
+  DO iks = 1, k_grid%nps
      DO ibnd = 1, nbnd
         emin = MIN (emin, et (ibnd, iks) )
      ENDDO
@@ -43,7 +44,7 @@ FUNCTION get_alpha_pv()
   CALL mp_min( emin, inter_pool_comm) 
   !
   emax = et (1, 1)
-  DO iks = 1, nks 
+  DO iks = 1, k_grid%nps 
      DO ibnd = 1, nbnd_occ(iks)
         emax = MAX (emax, et (ibnd, iks) )
      ENDDO

Tools/set_iks_l2g.f90

@@ -19,6 +19,7 @@ SUBROUTINE set_iks_l2g
   USE westcom,                ONLY : iks_l2g
   USE mp,                     ONLY : mp_sum
   USE mp_global,              ONLY : inter_pool_comm,npool,my_pool_id
+  USE types_bz_grid,          ONLY : k_grid
   !
   IMPLICIT NONE
   !
@@ -28,7 +29,7 @@ SUBROUTINE set_iks_l2g
   ALLOCATE( iks_l2g(nkstot) )
   !
   my_nks=0
-  my_nks(my_pool_id) = nks
+  my_nks(my_pool_id) = k_grid%nps
   CALL mp_sum( my_nks, inter_pool_comm )
   !
   my_offset = 0 
@@ -36,7 +37,7 @@ SUBROUTINE set_iks_l2g
      my_offset = my_offset + my_nks(ip)
   ENDDO
   !
-  DO iks = 1, nks
+  DO iks = 1, k_grid%nps
      iks_l2g(iks) = my_offset + iks 
   ENDDO
   !

Tools/set_nbndocc.f90

@@ -22,6 +22,7 @@ SUBROUTINE set_nbndocc( )
   USE noncollin_module,       ONLY : noncolin,npol
   USE io_global,              ONLY : stdout
   USE westcom,                ONLY : nbnd_occ
+  USE types_bz_grid,          ONLY : k_grid
   !
   IMPLICIT NONE
   !
@@ -30,11 +31,11 @@ SUBROUTINE set_nbndocc( )
   ! Calculate NBNDVAL
   !
   IF(ALLOCATED(nbnd_occ)) DEALLOCATE(nbnd_occ)
-  ALLOCATE( nbnd_occ(nks) )
+  ALLOCATE( nbnd_occ(k_grid%nps) )
   !
   IF(lsda) THEN  
      !
-     DO iks = 1, nks
+     DO iks = 1, k_grid%nps
         spin = isk(iks)
         !
         SELECT CASE(spin)

Westpp/do_rho.f90

@@ -90,12 +90,12 @@ SUBROUTINE do_rho ( )
         IF( gamma_only ) THEN 
            CALL single_invfft_gamma(dffts,npw,npwx,evc(1,global_ib),psic,'Wave')
            DO ir = 1, dffts%nnr
-              auxr(ir) = auxr(ir) + REAL( psic(ir), KIND=DP) *  REAL( psic(ir), KIND=DP) * wk(iks)  
+              auxr(ir) = auxr(ir) + REAL( psic(ir), KIND=DP) *  REAL( psic(ir), KIND=DP) * k_grid%weight(iks)  
            ENDDO 
         ELSE
            CALL single_invfft_k(dffts,npw,npwx,evc(1,global_ib),psic,'Wave',igk_k(1,current_k))
            DO ir = 1, dffts%nnr
-              auxr(ir) = auxr(ir) + REAL( CONJG( psic(ir) ) * psic(ir), KIND=DP) * wk(iks)
+              auxr(ir) = auxr(ir) + REAL( CONJG( psic(ir) ) * psic(ir), KIND=DP) * k_grid%weight(iks)
            ENDDO 
         ENDIF
         !

Westpp/do_sxx.f90

@@ -66,8 +66,8 @@ SUBROUTINE do_sxx ( )
   !
   CALL io_push_title("(S)creened eXact eXchange")
   !
-  ALLOCATE( sigma_exx( westpp_range(1):westpp_range(2), nks) )
-  ALLOCATE( sigma_sxx( westpp_range(1):westpp_range(2), nks) )
+  ALLOCATE( sigma_exx( westpp_range(1):westpp_range(2), k_grid%nps) )
+  ALLOCATE( sigma_sxx( westpp_range(1):westpp_range(2), k_grid%nps) )
   !
   sigma_exx = 0._DP
   sigma_sxx = 0._DP

Wfreq/calc_corr.f90

@@ -30,13 +30,14 @@ SUBROUTINE calc_corr_gamma( sigma_corr, energy, l_verbose)
   USE io_push,              ONLY : io_push_bar,io_push_value,io_push_title
   USE distribution_center,  ONLY : pert,ifr,rfr,aband
   USE types_coulomb,        ONLY : pot3D
+  USE types_bz_grid,        ONLY : k_grid
   !
   IMPLICIT NONE
   !
   ! I/O
   !
-  COMPLEX(DP),INTENT(OUT) :: sigma_corr( qp_bandrange(1):qp_bandrange(2), nks )  ! The correlation self-energy, imaginary part is lifetime.  
-  REAL(DP),INTENT(IN) :: energy( qp_bandrange(1):qp_bandrange(2), nks )          ! The energy variable
+  COMPLEX(DP),INTENT(OUT) :: sigma_corr( qp_bandrange(1):qp_bandrange(2), k_grid%nps )  ! The correlation self-energy, imaginary part is lifetime.  
+  REAL(DP),INTENT(IN) :: energy( qp_bandrange(1):qp_bandrange(2), k_grid%nps )          ! The energy variable
   LOGICAL,INTENT(IN) :: l_verbose
   !
   ! Workspace
@@ -73,12 +74,12 @@ SUBROUTINE calc_corr_gamma( sigma_corr, energy, l_verbose)
   IF(l_verbose) WRITE(stdout,'(5x,"Integrating along the IM axis...")') 
   IF(l_verbose) CALL io_push_bar
   !
-  IF(l_verbose) barra_load = nks * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
+  IF(l_verbose) barra_load = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
   IF(l_verbose) CALL start_bar_type( barra, 'sigmac_i', barra_load )
   !
   ! LOOP 
   !
-  DO iks = 1, nks   ! KPOINT-SPIN
+  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
      !
      nbndval = nbnd_occ(iks)
      !
@@ -145,12 +146,12 @@ SUBROUTINE calc_corr_gamma( sigma_corr, energy, l_verbose)
   IF(l_verbose) WRITE(stdout,'(5x,"Residues along the RE axis...")') 
   IF(l_verbose) CALL io_push_bar
   !
-  IF(l_verbose) barra_load = nks * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
+  IF(l_verbose) barra_load = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
   IF(l_verbose) CALL start_bar_type( barra, 'sigmac_r', barra_load )
   !
   ! LOOP 
   !
-  DO iks = 1, nks
+  DO iks = 1, k_grid%nps
      !
      nbndval = nbnd_occ(iks)
      !
@@ -298,9 +299,6 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
            ikk = k_grid%ip(ikks)
            !
            CALL k_grid%find( k_grid%p_cart(:,ik) - k_grid%p_cart(:,ikk), 1, 'cart', iq, g0 )
-           !CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' ) 
-           !ikqs = k_grid%find( kmq, 'cart' )
-           !ikqs = kmq_grid%index_kq(iks,iq)
            l_gammaq = q_grid%l_pIsGamma(iq)
            nbndval = nbnd_occ(ikks)
            !
@@ -387,9 +385,6 @@ SUBROUTINE calc_corr_k( sigma_corr, energy, l_verbose)
            ikk = k_grid%ip(ikks)
            !
            CALL k_grid%find( k_grid%p_cart(:,ik) - k_grid%p_cart(:,ikk), 1, 'cart', iq, g0 )
-           !CALL k_grid%add( k_grid%p_cart(:,ik), -q_grid%p_cart(:,iq), kmq, g0, 'cart' ) 
-           !ikqs = k_grid%find( kmq, 'cart' )
-           !!ikqs = kmq_grid%index_kq(iks,iq)
            l_gammaq = q_grid%l_pIsGamma(iq)
            nbndval = nbnd_occ(ikks)
            !

Wfreq/calc_exx2.f90

@@ -53,7 +53,7 @@ SUBROUTINE calc_exx2( sigma_exx, nb1, nb2 )
   ! I/O
   !
   INTEGER, INTENT(IN) :: nb1, nb2
-  REAL(DP),INTENT(OUT) :: sigma_exx( nb1:nb2, nks) 
+  REAL(DP),INTENT(OUT) :: sigma_exx( nb1:nb2, k_grid%nps) 
   !
   ! Workspace
   !

Wfreq/calc_vxc.f90

@@ -42,13 +42,14 @@ SUBROUTINE calc_vxc( sigma_vxcl, sigma_vxcnl )
   USE funct,                ONLY : dft_is_hybrid,get_exx_fraction
   USE class_idistribute,    ONLY : idistribute
   USE exx,                  ONLY : vexx,exxalfa
+  USE types_bz_grid,        ONLY : k_grid
   !
   IMPLICIT NONE
   !
   ! I/O
   !
-  REAL(DP),INTENT(OUT) :: sigma_vxcl( qp_bandrange(1):qp_bandrange(2), nks )
-  REAL(DP),INTENT(OUT) :: sigma_vxcnl( qp_bandrange(1):qp_bandrange(2), nks )
+  REAL(DP),INTENT(OUT) :: sigma_vxcl( qp_bandrange(1):qp_bandrange(2), k_grid%nps )
+  REAL(DP),INTENT(OUT) :: sigma_vxcnl( qp_bandrange(1):qp_bandrange(2), k_grid%nps )
   !
   ! Workspace
   !
@@ -84,12 +85,12 @@ SUBROUTINE calc_vxc( sigma_vxcl, sigma_vxcnl )
   !
   nnr = REAL( dfftp%nr1*dfftp%nr2*dfftp%nr3, KIND=DP )
   !
-  barra_load = nks 
+  barra_load = k_grid%nps
   CALL start_bar_type( barra, 'sigmavxc', barra_load )
   !
   ! LOOP 
   !
-  DO iks = 1, nks   ! KPOINT-SPIN
+  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
      !
      ! ... Set k-point, spin, kinetic energy, needed by Hpsi
      !
@@ -104,7 +105,7 @@ SUBROUTINE calc_vxc( sigma_vxcl, sigma_vxcnl )
      !
      ! ... read in wavefunctions from the previous iteration
      !
-     IF(nks>1) THEN
+     IF(k_grid%nps>1) THEN
         !iuwfc = 20
         !lrwfc = nbnd * npwx * npol 
         !!CALL get_buffer( evc, nwordwfc, iunwfc, iks )

Wfreq/solve_gfreq.f90

@@ -63,6 +63,7 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
   USE wfreq_restart,        ONLY : solvegfreq_restart_write,solvegfreq_restart_read,bks_type
   USE wfreq_io,             ONLY : writeout_overlap,writeout_solvegfreq,preallocate_solvegfreq
   USE types_coulomb,        ONLY : pot3D
+  USE types_bz_grid,        ONLY : k_grid
   !
   IMPLICIT NONE
   !
@@ -106,12 +107,12 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
      bks%lastdone_band = 0 
      bks%old_ks        = 0 
      bks%old_band      = 0 
-     bks%max_ks        = nks 
+     bks%max_ks        = k_grid%nps
      bks%min_ks        = 1 
   ENDIF
   !
   barra_load = 0
-  DO iks = 1, nks
+  DO iks = 1, k_grid%nps
      IF(iks<bks%lastdone_ks) CYCLE
      DO ib = qp_bandrange(1), qp_bandrange(2)
         IF(iks==bks%lastdone_ks .AND. ib <= bks%lastdone_band ) CYCLE
@@ -128,7 +129,7 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
   !
   ! LOOP 
   !
-  DO iks = 1, nks   ! KPOINT-SPIN
+  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
      IF(iks<bks%lastdone_ks) CYCLE
      !
      ! ... Set k-point, spin, kinetic energy, needed by Hpsi
@@ -144,7 +145,7 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
      !
      ! ... read in wavefunctions from the previous iteration
      !
-     IF(nks>1) THEN
+     IF(k_grid%nps>1) THEN
         !iuwfc = 20
         !lrwfc = nbnd * npwx * npol 
         !!CALL get_buffer( evc, nwordwfc, iunwfc, iks )
@@ -187,14 +188,7 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
         !
         ! PSIC
         !
-!        IF(gamma_only) THEN
-           CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ib),psic,'Wave') 
-!        ELSEIF(noncolin) THEN
-!           CALL SINGLEBAND_invfft_k(npw,evc(1     ,ib),npwx,psic_nc(1,1),dffts%nnr,.TRUE.)
-!           CALL SINGLEBAND_invfft_k(npw,evc(1+npwx,ib),npwx,psic_nc(1,2),dffts%nnr,.TRUE.)
-!        ELSE
-!           CALL SINGLEBAND_invfft_k(npw,evc(1,ib),npwx,psic,dffts%nnr,.TRUE.)
-!        ENDIF
+        CALL single_invfft_gamma(dffts,npw,npwx,evc(1,ib),psic,'Wave') 
         !
         ! ZEROS
         !
@@ -221,30 +215,11 @@ SUBROUTINE solve_gfreq_gamma(l_read_restart)
            ENDDO
            !
            ! Bring it to R-space
-!           IF(gamma_only) THEN
-              CALL single_invfft_gamma(dffts,npwq,npwqx,pertg(1),pertr,TRIM(fftdriver))
-              DO ir=1,dffts%nnr 
-                 pertr(ir)=psic(ir)*pertr(ir)
-              ENDDO
-              CALL single_fwfft_gamma(dffts,npw,npwx,pertr,dvpsi(1,ip),'Wave')
-!           ELSEIF(noncolin) THEN
-!              CALL SINGLEBAND_invfft_k(npwq,pertg(1),npwx,pertr,dffts%nnr,.FALSE.)
-!              DO ir=1,dffts%nnr 
-!                 pertr(ir)=psic_nc(ir,1)*DCONJG(pertr(ir))
-!              ENDDO
-!              CALL SINGLEBAND_fwfft_k(npw,pertr,dffts%nnr,dvpsi(1,ip),npwx,.TRUE.)
-!              CALL SINGLEBAND_invfft_k(npwq,pertg(1),npwx,pertr,dffts%nnr,.FALSE.)
-!              DO ir=1,dffts%nnr 
-!                 pertr(ir)=psic_nc(ir,2)*DCONJG(pertr(ir))
-!              ENDDO
-!              CALL SINGLEBAND_fwfft_k(npw,pertr,dffts%nnr,dvpsi(1+npwx,ip),npwx,.TRUE.)
-!           ELSE
-!              CALL SINGLEBAND_invfft_k(npwq,pertg(1),npwx,pertr,dffts%nnr,.FALSE.)
-!              DO ir=1,dffts%nnr 
-!                 pertr(ir)=psic(ir)*DCONJG(pertr(ir))
-!              ENDDO
-!              CALL SINGLEBAND_fwfft_k(npw,pertr,dffts%nnr,dvpsi(1,ip),npwx,.TRUE.)
-!           ENDIF 
+           CALL single_invfft_gamma(dffts,npwq,npwqx,pertg(1),pertr,TRIM(fftdriver))
+           DO ir=1,dffts%nnr 
+              pertr(ir)=psic(ir)*pertr(ir)
+           ENDDO
+           CALL single_fwfft_gamma(dffts,npw,npwx,pertr,dvpsi(1,ip),'Wave')
            !
            !
         ENDDO ! pert
@@ -451,8 +426,8 @@ SUBROUTINE solve_gfreq_k(l_read_restart)
   !
   ! LOOP 
   !
-  ! outer k-point loop (matrix element): ikks, npwk, evck, psick
-  ! inner k-point loop (sum over k'): iks, npw, evc (passed to h_psi)
+  ! ... Outer k-point loop (wfc matrix element): ikks, npwk, evck, psick
+  ! ... Inner k-point loop (wfc summed over k'): iks, npw, evc (passed to h_psi: current_k = iks)
   !
   DO ikks = 1, k_grid%nps   ! KPOINT-SPIN (MATRIX ELEMENT)
      IF(ikks<bksks%lastdone_ks) CYCLE

Wfreq/solve_hf.f90

@@ -51,8 +51,8 @@ SUBROUTINE solve_hf_gamma( )
   USE west_io,              ONLY : serial_table_output
   USE distribution_center,  ONLY : pert
   USE funct,                ONLY : get_exx_fraction,dft_is_hybrid
-  USE klist,                ONLY : wk
   USE wfreq_io,             ONLY : writeout_solvehf 
+  USE types_bz_grid,        ONLY : k_grid
   !
   IMPLICIT NONE
   !
@@ -80,12 +80,12 @@ SUBROUTINE solve_hf_gamma( )
   !
   sigma_hf(:,:) = sigma_exx(:,:) - sigma_vxcl(:,:) - sigma_vxcnl(:,:)
   !
-  CALL writeout_solvehf( sigma_hf(qp_bandrange(1),1), qp_bandrange(2)-qp_bandrange(1)+1, nks  )
+  CALL writeout_solvehf( sigma_hf(qp_bandrange(1),1), qp_bandrange(2)-qp_bandrange(1)+1, k_grid%nps )
  !
  ! Output it per k-point
  !
  ALLOCATE(out_tab(qp_bandrange(2)-qp_bandrange(1)+1,6))
- DO iks=1,nks
+ DO iks=1,k_grid%nps
     DO ib = qp_bandrange(1), qp_bandrange(2)
        out_tab( ib - qp_bandrange(1) + 1, 1) = REAL( ib, KIND=DP) 
        out_tab( ib - qp_bandrange(1) + 1, 2) = et(ib,iks) * rytoev
@@ -107,14 +107,14 @@ SUBROUTINE solve_hf_gamma( )
   IF( l_enable_gwetot) THEN
      !
      nbndval = MIN( MAXVAL( nbnd_occ(:) ), nbnd ) 
-     ALLOCATE(sigma_exx_all_occupied(nbndval,nks))
+     ALLOCATE(sigma_exx_all_occupied(nbndval,k_grid%nps))
      !
      CALL calc_exx2( sigma_exx_all_occupied, 1, nbndval ) 
      !
      exx_etot = 0._DP
-     DO iks = 1, nks 
+     DO iks = 1, k_grid%nps 
         DO ib = 1, nbnd_occ(iks) 
-           exx_etot = exx_etot + sigma_exx_all_occupied( ib, iks) * wk(iks) / 2._DP
+           exx_etot = exx_etot + sigma_exx_all_occupied( ib, iks) * k_grid%weight(iks) / 2._DP
         ENDDO
      ENDDO
      !

Wfreq/solve_qp.f90

@@ -56,6 +56,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
   USE bar,                  ONLY : bar_type,start_bar_type,update_bar_type,stop_bar_type
   USE wfreq_io,             ONLY : readin_overlap,readin_solvegfreq,readin_solvehf
   USE wfreq_db,             ONLY : wfreq_db_write
+  USE types_bz_grid,        ONLY : k_grid
   !
   IMPLICIT NONE
   !
@@ -126,11 +127,11 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
   !
   ! TEMP
   ! 
-  ALLOCATE( d_body1_ifr( aband%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), nks ) )
-  ALLOCATE( z_body_rfr( aband%nloc, rfr%nloc, qp_bandrange(1):qp_bandrange(2), nks ) )
+  ALLOCATE( d_body1_ifr( aband%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
+  ALLOCATE( z_body_rfr( aband%nloc, rfr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
   IF( l_enable_lanczos ) THEN
-     ALLOCATE( d_body2_ifr( n_lanczos, pert%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), nks ) )
-     ALLOCATE( d_diago( n_lanczos, pert%nloc, qp_bandrange(1):qp_bandrange(2), nks ) )
+     ALLOCATE( d_body2_ifr( n_lanczos, pert%nloc, ifr%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
+     ALLOCATE( d_diago( n_lanczos, pert%nloc, qp_bandrange(1):qp_bandrange(2), k_grid%nps ) )
   ENDIF
   !
   d_body1_ifr = 0._DP
@@ -144,12 +145,12 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
   !
   CALL io_push_title("Collecting results from W and G")
   !
-  barra_load = nks * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
+  barra_load = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
   CALL start_bar_type( barra, 'coll_gw', barra_load )
   !
   ! LOOP 
   !
-  DO iks = 1, nks   ! KPOINT-SPIN
+  DO iks = 1, k_grid%nps   ! KPOINT-SPIN
      !
      nbndval = nbnd_occ(iks)
      !
@@ -270,8 +271,8 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
   !
   ! Get Sigma_X
   !
-  ALLOCATE( sigma_hf (qp_bandrange(1):qp_bandrange(2),nks) )
-  CALL readin_solvehf( sigma_hf(qp_bandrange(1),1), qp_bandrange(2)-qp_bandrange(1)+1, nks )
+  ALLOCATE( sigma_hf (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
+  CALL readin_solvehf( sigma_hf(qp_bandrange(1),1), qp_bandrange(2)-qp_bandrange(1)+1, k_grid%nps )
   !
   ! For CORR
   !
@@ -279,24 +280,24 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
      !
      CALL io_push_title("(Q)uasiparticle energies")
      !
-     ALLOCATE( en(qp_bandrange(1):qp_bandrange(2),nks,2) ) 
-     ALLOCATE( sc(qp_bandrange(1):qp_bandrange(2),nks,2) ) 
-     ALLOCATE( l_conv(qp_bandrange(1):qp_bandrange(2),nks) ) 
-     ALLOCATE( sigma_cor_in (qp_bandrange(1):qp_bandrange(2),nks) )
-     ALLOCATE( sigma_cor_out(qp_bandrange(1):qp_bandrange(2),nks) )
-     ALLOCATE( z_in (qp_bandrange(1):qp_bandrange(2),nks) )
-     ALLOCATE( qp_energy (qp_bandrange(1):qp_bandrange(2),nks) )
+     ALLOCATE( en(qp_bandrange(1):qp_bandrange(2),k_grid%nps,2) ) 
+     ALLOCATE( sc(qp_bandrange(1):qp_bandrange(2),k_grid%nps,2) ) 
+     ALLOCATE( l_conv(qp_bandrange(1):qp_bandrange(2),k_grid%nps) ) 
+     ALLOCATE( sigma_cor_in (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
+     ALLOCATE( sigma_cor_out(qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
+     ALLOCATE( z_in (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
+     ALLOCATE( qp_energy (qp_bandrange(1):qp_bandrange(2),k_grid%nps) )
      !
      ! 1st step of secant solver : E_KS - 0.5 * eshift
      !
-     DO iks = 1, nks
+     DO iks = 1, k_grid%nps
         DO ib = qp_bandrange(1), qp_bandrange(2)
            en(ib,iks,1) = et(ib,iks) - eshift*0.5_DP
         ENDDO
      ENDDO
      CALL calc_corr_gamma( sc(:,:,1), en(:,:,1), .TRUE.)  
      !
-     !DO iks = 1, nks
+     !DO iks = 1, k_grid%nps
      !   DO ib = qp_bandrange(1), qp_bandrange(2)
      !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,1) * rytoev, sc(ib,iks,1) * rytoev
      !   ENDDO
@@ -304,13 +305,13 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
      !
      ! 1st step of secant solver : E_KS + 0.5 * eshift
      !
-     DO iks = 1, nks
+     DO iks = 1, k_grid%nps
         DO ib = qp_bandrange(1), qp_bandrange(2)
            en(ib,iks,2) = et(ib,iks) + eshift*0.5_DP
         ENDDO
      ENDDO
      CALL calc_corr_gamma( sc(:,:,2), en(:,:,2), .TRUE.) 
-     !DO iks = 1, nks
+     !DO iks = 1, k_grid%nps
      !   DO ib = qp_bandrange(1), qp_bandrange(2)
      !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,2) * rytoev, sc(ib,iks,2) * rytoev
      !   ENDDO
@@ -318,7 +319,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
      !
      ! Stage sigma_corr_in
      !
-     DO iks = 1, nks
+     DO iks = 1, k_grid%nps
         DO ib = qp_bandrange(1), qp_bandrange(2)
            sigma_cor_in(ib,iks) = ( sc(ib,iks,2) + sc(ib,iks,1) ) * 0.5_DP 
         ENDDO
@@ -326,7 +327,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
      !
      ! Stage z_in
      !
-     DO iks = 1, nks
+     DO iks = 1, k_grid%nps
         DO ib = qp_bandrange(1), qp_bandrange(2)
            z_in(ib,iks) = 1._DP / ( 1._DP-REAL( sc(ib,iks,2)-sc(ib,iks,1), KIND=DP ) / eshift  ) 
         ENDDO
@@ -335,7 +336,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
      ! en 1 = EKS, sc 1 = sigma_corr_in
      ! en 2 = EKS + Z * ( sigmax - vxc + sigmacorrin)
      !
-     DO iks = 1, nks
+     DO iks = 1, k_grid%nps
         DO ib = qp_bandrange(1), qp_bandrange(2)
            en(ib,iks,1) = et(ib,iks) 
            sc(ib,iks,1) = sigma_cor_in(ib,iks)
@@ -350,14 +351,14 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
      ! nth step of the secant solver
      !
      l_conv = .FALSE.
-     notconv = nks * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
+     notconv = k_grid%nps * ( qp_bandrange(2) - qp_bandrange(1) + 1 )
      DO ifixed = 1, n_secant_maxiter
         ! 
         WRITE(cifixed,"(i6.6)") ifixed
         !CALL io_push_title("Iter : "//cifixed)
         !
         CALL calc_corr_gamma( sc(:,:,2), en(:,:,2), .TRUE.)  
-        !DO iks = 1, nks
+        !DO iks = 1, k_grid%nps
         !   DO ib = qp_bandrange(1), qp_bandrange(2)
         !      WRITE(stdout,"(5X,1f14.6,' : ',2f14.6)") en(ib,iks,2) * rytoev, sc(ib,iks,2) * rytoev
         !   ENDDO
@@ -365,7 +366,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
         !
         ! Resulting new energy:
         !
-        DO iks = 1, nks
+        DO iks = 1, k_grid%nps
            DO ib = qp_bandrange(1), qp_bandrange(2)
                IF( .NOT. l_conv(ib,iks) ) THEN
                   qp_energy(ib,iks) = en(ib,iks,2) + &
@@ -377,7 +378,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
         !
         ! Estimate l_conv
         !
-        DO iks = 1, nks
+        DO iks = 1, k_grid%nps
            DO ib = qp_bandrange(1), qp_bandrange(2)
               l_conv(ib,iks) = ( ABS( qp_energy(ib,iks) - en(ib,iks,2) ) < trev_secant ) 
            ENDDO
@@ -388,7 +389,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
        !WRITE(stdout,"(5X,'Iter: ',i6.6)") ifixed
        !WRITE(stdout,"(5X,a,1X,a,1X,a,1X,a)") 'K     ', 'B     ', 'QP energ. [eV]', 'conv'
        !CALL io_push_bar()
-       !DO iks = 1, nks
+       !DO iks = 1, k_grid%nps
        !   DO ib = qp_bandrange(1), qp_bandrange(2)
        !      WRITE(stdout,"(5X,i6.6,1X,i6.6,1X,1f14.6,4X,l)") iks, ib, qp_energy(ib,iks) * rytoev, l_conv(ib,iks) 
        !   ENDDO
@@ -398,7 +399,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)
         ! Count the number of notconverged QP energies
         !
         notconv = 0 
-        DO iks = 1, nks
+        DO iks = 1, k_grid%nps
            DO ib = qp_bandrange(1), qp_bandrange(2)
               IF( .NOT.l_conv(ib,iks) ) notconv = notconv + 1
            ENDDO
@@ -445,7 +446,7 @@ SUBROUTINE solve_qp_gamma(l_secant,l_generate_plot)