Fixed setup of k_grid in the gamma_only case. Subroutines chi_invert_real and...

Fixed setup of k_grid in the gamma_only case. Subroutines chi_invert_real and chi_invert_complex placed in a module, to make use of optional arguments.

Tools/class_bz_grid.f90

@@ -62,6 +62,7 @@ MODULE class_bz_grid
       USE klist,            ONLY : xk, wk, nkstot
       USE start_k,          ONLY : nk1, nk2, nk3
       USE pwcom,            ONLY : nspin
+      USE control_flags,    ONLY : gamma_only
       USE westcom,          ONLY : nq
       USE constants,        ONLY : eps8
       !
@@ -80,7 +81,8 @@ MODULE class_bz_grid
       SELECT CASE( grid_type )
       CASE ( 'K', 'k')
          !
-         this%ngrid(1:3) = (/ nk1, nk2, nk3 /) 
+         ! This is a workaround to prevent ngrid(:) to be set to (/ 0, 0, 0 /) in the gamma_only case (espresso default)
+         IF ( .NOT. gamma_only ) this%ngrid(1:3) = (/ nk1, nk2, nk3 /) 
          this%np = this%ngrid(1) * this%ngrid(2) * this%ngrid(3) 
          this%ns = nspin 
          this%nps = nkstot   !    =   np * ns  

Tools/do_setup.f90

@@ -65,11 +65,18 @@ SUBROUTINE do_setup
   !
   k_grid = bz_grid()
   CALL k_grid%init('K')
+  print *, 'K'
+  print *, k_grid%p_cart
+  print *, k_grid%ngrid
+  print *, 'nks = ', nks
   !
   q_grid = bz_grid()
   CALL q_grid%init('Q')
+  print *, 'Q'
+  print *, q_grid%p_cart
+  print *, q_grid%ngrid
   !
-  IF ( ANY ( q_grid%ngrid(:) - k_grid%ngrid(:) /= 0   ) ) THEN
+  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
   !

Wfreq/chi_invert.f90

@@ -11,261 +11,276 @@
 ! Marco Govoni
 !
 !-----------------------------------------------------------------------
-SUBROUTINE chi_invert_real(matilda,head,lambda,nma)
-  !-----------------------------------------------------------------------
-  !
-  ! For each frequency I calculate X, ky, head and lambda
-  !
-  ! X = (1-B)^{-1}
-  ! temph = X * wh
-  ! templ = wl * X
-  ! tempt = wl * temph
-  !
-  ! ky = 1 - f - Tr ( tempt ) / 3
-  !
-  ! Head = (1-ky) / ky
-  !
-  ! Lamba = X * B + temph * templ / (3*ky)
-  !
-  USE kinds,                 ONLY : DP 
-  USE linear_algebra_kernel, ONLY : matinvrs_dge 
-  USE westcom,               ONLY : west_prefix,n_pdep_eigen_to_use,l_macropol
-  USE io_files,              ONLY : tmp_dir
-  !
-  ! I/O
-  !
-  REAL(DP),INTENT(IN) :: matilda(nma,nma)
-  REAL(DP),INTENT(OUT) :: head,lambda(n_pdep_eigen_to_use,n_pdep_eigen_to_use) 
-  INTEGER,INTENT(IN) :: nma
-  !
-  ! Workspace
-  ! 
-  REAL(DP) :: f
-  REAL(DP), ALLOCATABLE :: body(:,:)
-  REAL(DP),ALLOCATABLE :: wh(:,:)
-  REAL(DP),ALLOCATABLE :: wl(:,:)
-  INTEGER :: i1,i2
-  REAL(DP),ALLOCATABLE :: x(:,:)
-  REAL(DP),ALLOCATABLE :: temph(:,:)
-  REAL(DP),ALLOCATABLE :: templ(:,:)
-  REAL(DP) :: tempt(3,3)
-  REAL(DP) :: ky
-  !
-  ALLOCATE( body(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
-  ALLOCATE( x(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
-  !
-  DO i2 = 1, n_pdep_eigen_to_use
-     DO i1 = 1, n_pdep_eigen_to_use
-        body(i1,i2) = matilda(i1,i2) 
-     ENDDO
-  ENDDO
-  !
-  IF(l_macropol) THEN
-     !
-     f = 0._DP
-     DO i1 = 1, 3
-        f = f + matilda(n_pdep_eigen_to_use+i1,n_pdep_eigen_to_use+i1) / 3._DP
-     ENDDO
-     !
-     ALLOCATE( wh( n_pdep_eigen_to_use, 3)) 
-     DO i2 = 1, 3
-        DO i1 = 1, n_pdep_eigen_to_use
-           wh(i1,i2) = matilda(i1,n_pdep_eigen_to_use+i2)
-        ENDDO
-     ENDDO
-     !
-     ALLOCATE( wl( 3, n_pdep_eigen_to_use)) 
-     DO i2 = 1, n_pdep_eigen_to_use
-        DO i1 = 1, 3
-           wl(i1,i2) = matilda(n_pdep_eigen_to_use+i1,i2)
-        ENDDO
-     ENDDO
-     !
-  ENDIF
-  !
-  ! X = (1-B)^{-1}
-  x(:,:) = -body(:,:)
-  DO i1=1,n_pdep_eigen_to_use
-     x(i1,i1) = x(i1,i1) + 1._DP
-  ENDDO
-  !
-  CALL matinvrs_dge(n_pdep_eigen_to_use,x)
-  !
-  IF( l_macropol) THEN
-     !
-     ! temph = X * wh
-     ALLOCATE( temph(n_pdep_eigen_to_use,3) )
-     CALL DGEMM( 'N', 'N', n_pdep_eigen_to_use, 3, n_pdep_eigen_to_use, 1._DP, x, n_pdep_eigen_to_use, &
-     & wh, n_pdep_eigen_to_use, 0._DP, temph, n_pdep_eigen_to_use )
-     DEALLOCATE(wh)
-     !
-     ! templ = wl * X
-     ALLOCATE( templ(3,n_pdep_eigen_to_use) )
-     CALL DGEMM( 'N', 'N', 3, n_pdep_eigen_to_use, n_pdep_eigen_to_use, 1._DP, wl, 3, x, &
-     & n_pdep_eigen_to_use, 0._DP, templ, 3 )
-     !
-     ! tempt = wl * temph
-     CALL DGEMM( 'N', 'N', 3, 3, n_pdep_eigen_to_use, 1._DP, wl, 3, temph, n_pdep_eigen_to_use, 0._DP, tempt, 3 )
-     DEALLOCATE(wl)
-     !
-     ! ky = 1 - f - Tr ( tempt ) / 3
-     ky = 1._DP - f
-     DO i1=1,3
-        ky = ky - tempt(i1,i1) / 3._DP
-     ENDDO
-     !
-     head = (1._DP - ky) / ky
-     !
-  ELSE
-     head = 0._DP
-  ENDIF
-  !
-  CALL DGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, n_pdep_eigen_to_use, 1._DP, x, n_pdep_eigen_to_use, &
-  & body, n_pdep_eigen_to_use, 0._DP, lambda, n_pdep_eigen_to_use )
-  !
-  IF( l_macropol ) THEN
-     CALL DGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, 3, 1._DP/(3._DP*ky), temph, &
-     & n_pdep_eigen_to_use, templ, 3, 1._DP, lambda, n_pdep_eigen_to_use )
-     DEALLOCATE( temph )
-     DEALLOCATE( templ )
-  ENDIF
-  !
-  DEALLOCATE( body, x) 
-  !
-END SUBROUTINE
-!
-!
-!-----------------------------------------------------------------------
-SUBROUTINE chi_invert_complex(matilda,head,lambda,nma,l_gammaq)
-  !-----------------------------------------------------------------------
-  !
-  ! For each frequency and q-point I calculate X, ky, head and lambda
-  !
-  ! X = (1-B)^{-1}
-  ! temph = X * wh
-  ! templ = wl * X
-  ! tempt = wl * temph
-  !
-  ! ky = 1 - f - Tr ( tempt ) / 3
-  !
-  ! Head = (1-ky) / ky
-  !
-  ! Lamba = X * B + temph * templ / (3*ky)
-  !
-  USE kinds,                 ONLY : DP 
-  USE linear_algebra_kernel, ONLY : matinvrs_zge 
-  USE westcom,               ONLY : west_prefix,n_pdep_eigen_to_use,l_macropol
-  USE io_files,              ONLY : tmp_dir
-  !
-  ! I/O
-  !
-  COMPLEX(DP),INTENT(IN) :: matilda(nma,nma)
-  COMPLEX(DP),INTENT(OUT) :: head,lambda(n_pdep_eigen_to_use,n_pdep_eigen_to_use) 
-  INTEGER,INTENT(IN) :: nma
-  LOGICAL,OPTIONAL :: l_gammaq
-  !
-  ! Workspace
-  ! 
-  COMPLEX(DP) :: f
-  COMPLEX(DP),ALLOCATABLE :: body(:,:)
-  COMPLEX(DP),ALLOCATABLE :: wh(:,:)
-  COMPLEX(DP),ALLOCATABLE :: wl(:,:)
-  INTEGER :: i1,i2
-  COMPLEX(DP),ALLOCATABLE :: x(:,:)
-  COMPLEX(DP),ALLOCATABLE :: temph(:,:)
-  COMPLEX(DP),ALLOCATABLE :: templ(:,:)
-  COMPLEX(DP) :: tempt(3,3)
-  COMPLEX(DP) :: ky,Zone,Zzero
-  LOGICAL :: l_dohead
-  !
-  IF( PRESENT(l_gammaq) ) THEN 
-     l_dohead = l_macropol .AND. l_gammaq 
-  ELSE 
-     l_dohead = l_macropol
-  ENDIF 
-  !
-  ALLOCATE( body(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
-  ALLOCATE( x(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
-  !
-  Zone = CMPLX(1._DP,0._DP,KIND=DP)
-  Zzero = CMPLX(0._DP,0._DP,KIND=DP)
-  !
-  DO i2 = 1, n_pdep_eigen_to_use
-     DO i1 = 1, n_pdep_eigen_to_use
-        body(i1,i2) = matilda(i1,i2) 
-     ENDDO
-  ENDDO
-  !
-  IF(l_dohead) THEN
-     !
-     f = Zzero
-     DO i1 = 1, 3
-        f = f + matilda(n_pdep_eigen_to_use+i1,n_pdep_eigen_to_use+i1) / 3._DP
-     ENDDO
-     !
-     ALLOCATE( wh( n_pdep_eigen_to_use, 3)) 
-     DO i2 = 1, 3
-        DO i1 = 1, n_pdep_eigen_to_use
-           wh(i1,i2) = matilda(i1,n_pdep_eigen_to_use+i2)
-        ENDDO
-     ENDDO
-     !
-     ALLOCATE( wl( 3, n_pdep_eigen_to_use)) 
-     DO i2 = 1, n_pdep_eigen_to_use
-        DO i1 = 1, 3
-           wl(i1,i2) = matilda(n_pdep_eigen_to_use+i1,i2)
-        ENDDO
-     ENDDO
-     !
-  ENDIF
-  !
-  ! X = (1-B)^{-1}
-  x(:,:) = -body(:,:)
-  DO i1=1,n_pdep_eigen_to_use
-     x(i1,i1) = x(i1,i1) + Zone
-  ENDDO
-  !
-  CALL matinvrs_zge(n_pdep_eigen_to_use,x)
-  !
-  IF(l_dohead) THEN
-     !
-     ! temph = X * wh
-     ALLOCATE( temph(n_pdep_eigen_to_use,3) )
-     CALL ZGEMM( 'N', 'N', n_pdep_eigen_to_use, 3, n_pdep_eigen_to_use, Zone, x, n_pdep_eigen_to_use, &
-     & wh, n_pdep_eigen_to_use, Zzero, temph, n_pdep_eigen_to_use )
-     DEALLOCATE(wh)
-     !
-     ! templ = wl * X
-     ALLOCATE( templ(3,n_pdep_eigen_to_use) )
-     CALL ZGEMM( 'N', 'N', 3, n_pdep_eigen_to_use, n_pdep_eigen_to_use, Zone, wl, 3, x, &
-     & n_pdep_eigen_to_use, Zzero, templ, 3 )
-     !
-     ! tempt = wl * temph
-     CALL ZGEMM( 'N', 'N', 3, 3, n_pdep_eigen_to_use, Zone, wl, 3, temph, n_pdep_eigen_to_use, Zzero, tempt, 3 )
-     DEALLOCATE(wl)
-     !
-     ! ky = 1 - f - Tr ( tempt ) / 3
-     ky = Zone - f
-     DO i1=1,3
-        ky = ky - tempt(i1,i1) / 3._DP
-     ENDDO
-     !
-     head = (Zone - ky) / ky
-     !
-  ELSE
-     head = Zzero
-  ENDIF
-  !
-  CALL ZGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, n_pdep_eigen_to_use, Zone, x, n_pdep_eigen_to_use, &
-  & body, n_pdep_eigen_to_use, Zzero, lambda, n_pdep_eigen_to_use )
-  !
-  IF( l_dohead ) THEN
-     CALL ZGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, 3, Zone/(3._DP*ky), temph, &
-     & n_pdep_eigen_to_use, templ, 3, Zone, lambda, n_pdep_eigen_to_use )
-     DEALLOCATE( temph )
-     DEALLOCATE( templ )
-  ENDIF
-  !
-  DEALLOCATE( body, x )
-  !
-END SUBROUTINE
+MODULE chi_invert
+   !-----------------------------------------------------------------------
+   !
+   IMPLICIT NONE
+   !
+   PRIVATE
+   !
+   PUBLIC :: chi_invert_real, chi_invert_complex
+   !
+   CONTAINS
+   ! 
+   ! 
+   !-----------------------------------------------------------------------
+   SUBROUTINE chi_invert_real(matilda,head,lambda,nma)
+      !-----------------------------------------------------------------------
+      !
+      ! For each frequency I calculate X, ky, head and lambda
+      !
+      ! X = (1-B)^{-1}
+      ! temph = X * wh
+      ! templ = wl * X
+      ! tempt = wl * temph
+      !
+      ! ky = 1 - f - Tr ( tempt ) / 3
+      !
+      ! Head = (1-ky) / ky
+      !
+      ! Lamba = X * B + temph * templ / (3*ky)
+      !
+      USE kinds,                 ONLY : DP 
+      USE linear_algebra_kernel, ONLY : matinvrs_dge 
+      USE westcom,               ONLY : west_prefix,n_pdep_eigen_to_use,l_macropol
+      USE io_files,              ONLY : tmp_dir
+      !
+      ! I/O
+      !
+      REAL(DP),INTENT(IN) :: matilda(nma,nma)
+      REAL(DP),INTENT(OUT) :: head,lambda(n_pdep_eigen_to_use,n_pdep_eigen_to_use) 
+      INTEGER,INTENT(IN) :: nma
+      !
+      ! Workspace
+      ! 
+      REAL(DP) :: f
+      REAL(DP), ALLOCATABLE :: body(:,:)
+      REAL(DP),ALLOCATABLE :: wh(:,:)
+      REAL(DP),ALLOCATABLE :: wl(:,:)
+      INTEGER :: i1,i2
+      REAL(DP),ALLOCATABLE :: x(:,:)
+      REAL(DP),ALLOCATABLE :: temph(:,:)
+      REAL(DP),ALLOCATABLE :: templ(:,:)
+      REAL(DP) :: tempt(3,3)
+      REAL(DP) :: ky
+      !
+      ALLOCATE( body(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
+      ALLOCATE( x(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
+      !
+      DO i2 = 1, n_pdep_eigen_to_use
+         DO i1 = 1, n_pdep_eigen_to_use
+            body(i1,i2) = matilda(i1,i2) 
+         ENDDO
+      ENDDO
+      !
+      IF(l_macropol) THEN
+         !
+         f = 0._DP
+         DO i1 = 1, 3
+            f = f + matilda(n_pdep_eigen_to_use+i1,n_pdep_eigen_to_use+i1) / 3._DP
+         ENDDO
+         !
+         ALLOCATE( wh( n_pdep_eigen_to_use, 3)) 
+         DO i2 = 1, 3
+            DO i1 = 1, n_pdep_eigen_to_use
+               wh(i1,i2) = matilda(i1,n_pdep_eigen_to_use+i2)
+            ENDDO
+         ENDDO
+         !
+         ALLOCATE( wl( 3, n_pdep_eigen_to_use)) 
+         DO i2 = 1, n_pdep_eigen_to_use
+            DO i1 = 1, 3
+               wl(i1,i2) = matilda(n_pdep_eigen_to_use+i1,i2)
+            ENDDO
+         ENDDO
+         !
+      ENDIF
+      !
+      ! X = (1-B)^{-1}
+      x(:,:) = -body(:,:)
+      DO i1=1,n_pdep_eigen_to_use
+         x(i1,i1) = x(i1,i1) + 1._DP
+      ENDDO
+      !
+      CALL matinvrs_dge(n_pdep_eigen_to_use,x)
+      !
+      IF( l_macropol) THEN
+         !
+         ! temph = X * wh
+         ALLOCATE( temph(n_pdep_eigen_to_use,3) )
+         CALL DGEMM( 'N', 'N', n_pdep_eigen_to_use, 3, n_pdep_eigen_to_use, 1._DP, x, n_pdep_eigen_to_use, &
+         & wh, n_pdep_eigen_to_use, 0._DP, temph, n_pdep_eigen_to_use )
+         DEALLOCATE(wh)
+         !
+         ! templ = wl * X
+         ALLOCATE( templ(3,n_pdep_eigen_to_use) )
+         CALL DGEMM( 'N', 'N', 3, n_pdep_eigen_to_use, n_pdep_eigen_to_use, 1._DP, wl, 3, x, &
+         & n_pdep_eigen_to_use, 0._DP, templ, 3 )
+         !
+         ! tempt = wl * temph
+         CALL DGEMM( 'N', 'N', 3, 3, n_pdep_eigen_to_use, 1._DP, wl, 3, temph, n_pdep_eigen_to_use, 0._DP, tempt, 3 )
+         DEALLOCATE(wl)
+         !
+         ! ky = 1 - f - Tr ( tempt ) / 3
+         ky = 1._DP - f
+         DO i1=1,3
+            ky = ky - tempt(i1,i1) / 3._DP
+         ENDDO
+         !
+         head = (1._DP - ky) / ky
+         !
+      ELSE
+         head = 0._DP
+      ENDIF
+      !
+      CALL DGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, n_pdep_eigen_to_use, 1._DP, x, n_pdep_eigen_to_use, &
+      & body, n_pdep_eigen_to_use, 0._DP, lambda, n_pdep_eigen_to_use )
+      !
+      IF( l_macropol ) THEN
+         CALL DGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, 3, 1._DP/(3._DP*ky), temph, &
+         & n_pdep_eigen_to_use, templ, 3, 1._DP, lambda, n_pdep_eigen_to_use )
+         DEALLOCATE( temph )
+         DEALLOCATE( templ )
+      ENDIF
+      !
+      DEALLOCATE( body, x) 
+      !
+    END SUBROUTINE
+    !
+    !
+    !-----------------------------------------------------------------------
+    SUBROUTINE chi_invert_complex(matilda,head,lambda,nma,l_gammaq)
+      !-----------------------------------------------------------------------
+      !
+      ! For each frequency and q-point I calculate X, ky, head and lambda
+      !
+      ! X = (1-B)^{-1}
+      ! temph = X * wh
+      ! templ = wl * X
+      ! tempt = wl * temph
+      !
+      ! ky = 1 - f - Tr ( tempt ) / 3
+      !
+      ! Head = (1-ky) / ky
+      !
+      ! Lamba = X * B + temph * templ / (3*ky)
+      !
+      USE kinds,                 ONLY : DP 
+      USE linear_algebra_kernel, ONLY : matinvrs_zge 
+      USE westcom,               ONLY : west_prefix,n_pdep_eigen_to_use,l_macropol
+      USE io_files,              ONLY : tmp_dir
+      !
+      ! I/O
+      !
+      COMPLEX(DP),INTENT(IN) :: matilda(nma,nma)
+      COMPLEX(DP),INTENT(OUT) :: head,lambda(n_pdep_eigen_to_use,n_pdep_eigen_to_use) 
+      INTEGER,INTENT(IN) :: nma
+      LOGICAL,INTENT(IN) :: l_gammaq
+      !
+      ! Workspace
+      ! 
+      COMPLEX(DP) :: f
+      COMPLEX(DP),ALLOCATABLE :: body(:,:)
+      COMPLEX(DP),ALLOCATABLE :: wh(:,:)
+      COMPLEX(DP),ALLOCATABLE :: wl(:,:)
+      INTEGER :: i1,i2
+      COMPLEX(DP),ALLOCATABLE :: x(:,:)
+      COMPLEX(DP),ALLOCATABLE :: temph(:,:)
+      COMPLEX(DP),ALLOCATABLE :: templ(:,:)
+      COMPLEX(DP) :: tempt(3,3)
+      COMPLEX(DP) :: ky,Zone,Zzero
+      LOGICAL :: l_dohead
+      !
+      IF( PRESENT(l_gammaq) ) THEN 
+         l_dohead = l_macropol .AND. l_gammaq 
+      ELSE 
+         l_dohead = l_macropol
+      ENDIF 
+      !
+      ALLOCATE( body(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
+      ALLOCATE( x(n_pdep_eigen_to_use,n_pdep_eigen_to_use) )
+      !
+      Zone = CMPLX(1._DP,0._DP,KIND=DP)
+      Zzero = CMPLX(0._DP,0._DP,KIND=DP)
+      !
+      DO i2 = 1, n_pdep_eigen_to_use
+         DO i1 = 1, n_pdep_eigen_to_use
+            body(i1,i2) = matilda(i1,i2) 
+         ENDDO
+      ENDDO
+      !
+      IF(l_dohead) THEN
+         !
+         f = Zzero
+         DO i1 = 1, 3
+            f = f + matilda(n_pdep_eigen_to_use+i1,n_pdep_eigen_to_use+i1) / 3._DP
+         ENDDO
+         !
+         ALLOCATE( wh( n_pdep_eigen_to_use, 3)) 
+         DO i2 = 1, 3
+            DO i1 = 1, n_pdep_eigen_to_use
+               wh(i1,i2) = matilda(i1,n_pdep_eigen_to_use+i2)
+            ENDDO
+         ENDDO
+         !
+         ALLOCATE( wl( 3, n_pdep_eigen_to_use)) 
+         DO i2 = 1, n_pdep_eigen_to_use
+            DO i1 = 1, 3
+               wl(i1,i2) = matilda(n_pdep_eigen_to_use+i1,i2)
+            ENDDO
+         ENDDO
+         !
+      ENDIF
+      !
+      ! X = (1-B)^{-1}
+      x(:,:) = -body(:,:)
+      DO i1=1,n_pdep_eigen_to_use
+         x(i1,i1) = x(i1,i1) + Zone
+      ENDDO
+      !
+      CALL matinvrs_zge(n_pdep_eigen_to_use,x)
+      !
+      IF(l_dohead) THEN
+         !
+         ! temph = X * wh
+         ALLOCATE( temph(n_pdep_eigen_to_use,3) )
+         CALL ZGEMM( 'N', 'N', n_pdep_eigen_to_use, 3, n_pdep_eigen_to_use, Zone, x, n_pdep_eigen_to_use, &
+         & wh, n_pdep_eigen_to_use, Zzero, temph, n_pdep_eigen_to_use )
+         DEALLOCATE(wh)
+         !
+         ! templ = wl * X
+         ALLOCATE( templ(3,n_pdep_eigen_to_use) )
+         CALL ZGEMM( 'N', 'N', 3, n_pdep_eigen_to_use, n_pdep_eigen_to_use, Zone, wl, 3, x, &
+         & n_pdep_eigen_to_use, Zzero, templ, 3 )
+         !
+         ! tempt = wl * temph
+         CALL ZGEMM( 'N', 'N', 3, 3, n_pdep_eigen_to_use, Zone, wl, 3, temph, n_pdep_eigen_to_use, Zzero, tempt, 3 )
+         DEALLOCATE(wl)
+         !
+         ! ky = 1 - f - Tr ( tempt ) / 3
+         ky = Zone - f
+         DO i1=1,3
+            ky = ky - tempt(i1,i1) / 3._DP
+         ENDDO
+         !
+         head = (Zone - ky) / ky
+         !
+      ELSE
+         head = Zzero
+      ENDIF
+      !
+      CALL ZGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, n_pdep_eigen_to_use, Zone, x, n_pdep_eigen_to_use, &
+      & body, n_pdep_eigen_to_use, Zzero, lambda, n_pdep_eigen_to_use )
+      !
+      IF( l_dohead ) THEN
+         CALL ZGEMM( 'N', 'N', n_pdep_eigen_to_use, n_pdep_eigen_to_use, 3, Zone/(3._DP*ky), temph, &
+         & n_pdep_eigen_to_use, templ, 3, Zone, lambda, n_pdep_eigen_to_use )
+         DEALLOCATE( temph )
+         DEALLOCATE( templ )
+      ENDIF
+      !
+      DEALLOCATE( body, x )
+      !
+   END SUBROUTINE
+   !
+END MODULE

Wfreq/solve_wfreq.f90

@@ -66,6 +66,7 @@ SUBROUTINE solve_wfreq_gamma(l_read_restart,l_generate_plot)
   USE class_bz_grid,        ONLY : bz_grid
   USE types_bz_grid,        ONLY : k_grid, q_grid
   USE coulomb,              ONLY : store_sqvc
+  USE chi_invert,           ONLY : chi_invert_real, chi_invert_complex
   !
   IMPLICIT NONE
   !
@@ -644,6 +645,7 @@ SUBROUTINE solve_wfreq_k(l_read_restart,l_generate_plot)
   USE class_bz_grid,        ONLY : bz_grid
   USE types_bz_grid,        ONLY : k_grid, q_grid, compute_phase
   USE coulomb,              ONLY : store_sqvc
+  USE chi_invert,           ONLY : chi_invert_complex