.gitlab-ci.yml

@@ -19,3 +19,9 @@ style_check:
   image: gitlab-registry.mpcdf.mpg.de/tabriz/findent:latest
   script:
     - findent_batch --dir=./src --dry-run
+  artifacts:
+    paths:
+        - findent.patch
+    when: on_failure
+    expire_in: 24 hrs
+

src/electrons/v_ks.F90

@@ -418,10 +418,6 @@ contains
           end if
           ks%oep%level = OEP_LEVEL_NONE
         case default
-          if (kpoints%reduced%npoints > 1) then
-            call messages_not_implemented("OEP exchange with k-points", namespace=namespace)
-          end if
-
           if(oep_type == -1) then ! Else we have a MGGA
             if(ks%has_photons) then
               oep_type = OEP_TYPE_PHOTONS

src/electrons/xc_kli_inc.F90

 !! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
 !! Copyright (C) 2012-2013 M. Gruning, P. Melo, M. Oliveira
+!! Copyright (C) 2021 N. Tancogne-Dejean
 !!
 !! This program is free software; you can redistribute it and/or modify
 !! it under the terms of the GNU General Public License as published by
@@ -18,75 +19,112 @@
 !!
 
 ! ---------------------------------------------------------
-subroutine X(xc_KLI_solve) (mesh, st, is, oep, first)
+subroutine X(xc_KLI_solve) (space, mesh, st, oep, rcell_volume)
+  type(space_t),            intent(in)    :: space
   class(mesh_t),            intent(in)    :: mesh
   type(states_elec_t),      intent(in)    :: st
-  integer,                  intent(in)    :: is
   type(xc_oep_t),           intent(inout) :: oep
-  logical,                  intent(in)    :: first
+  FLOAT,                    intent(in)    :: rcell_volume
 
-  integer :: ist, ip, jst, eigen_n, kssi, kssj, proc
-  FLOAT, allocatable :: rho_sigma(:), v_bar_S(:), sqphi(:, :, :), dd(:)
-  FLOAT, allocatable :: Ma(:,:), xx(:,:), yy(:,:)
+  integer :: ist, ip, jst, eigen_n, kssi, kssj, ispin, ik
+  FLOAT, allocatable :: rho_sigma(:,:), v_bar_S(:), sqphi(:, :, :), dd(:)
+  FLOAT, allocatable :: Ma(:,:), xx(:,:), yy(:,:), kli(:,:)
   R_TYPE, allocatable :: psi(:, :)
+  FLOAT :: cnst(st%d%nspin)
+
+  PUSH_SUB(X(xc_KLI_solve))
 
   ASSERT(oep%type /= OEP_TYPE_PHOTONS)
+  ASSERT(st%d%ispin /= SPINORS)
 
   call profiling_in(C_PROFILING_XC_KLI, TOSTRING(X(XC_KLI)))
 
-  PUSH_SUB(X(xc_KLI_solve))
+  oep%vxc = M_ZERO
+
   ! some intermediate quantities
   ! vxc contains the Slater part!
-  SAFE_ALLOCATE(rho_sigma(1:mesh%np))
-  SAFE_ALLOCATE(sqphi(1:mesh%np, 1:st%d%dim, 1:st%nst))
+  SAFE_ALLOCATE(rho_sigma(1:mesh%np, 1:st%d%spin_channels))
+  rho_sigma = M_ZERO
+
+  SAFE_ALLOCATE(sqphi(1:mesh%np, 1:st%nst, st%d%kpt%start:st%d%kpt%end))
   SAFE_ALLOCATE(psi(1:mesh%np, 1:st%d%dim))
 
+  do ik = st%d%kpt%start, st%d%kpt%end
+    ispin = st%d%get_spin_index(ik)
+
     do ist = st%st_start, st%st_end
-    call states_elec_get_state(st, mesh, ist, is, psi)
-    sqphi(1:mesh%np, 1:st%d%dim, ist) = abs(psi(1:mesh%np, 1:st%d%dim))**2
+      call states_elec_get_state(st, mesh, ist, ik, psi)
+      sqphi(1:mesh%np, ist, ik) = R_REAL(R_CONJ(psi(1:mesh%np, 1))*psi(1:mesh%np, 1))
     end do
 
     do ip = 1, mesh%np
-    rho_sigma(ip) = max(sum(oep%socc * st%occ(st%st_start:st%st_end, is) * sqphi(ip, 1, st%st_start:st%st_end)), CNST(1e-20))
+      rho_sigma(ip, ispin) = rho_sigma(ip, ispin) + st%d%kweights(ik) &
+        * sum(oep%socc * st%occ(st%st_start:st%st_end, ik) * sqphi(ip, st%st_start:st%st_end, ik))
+    end do
   end do
 
-  if (st%parallel_in_states) then
-    SAFE_ALLOCATE(dd(1:mesh%np))
-    call st%mpi_grp%allreduce(rho_sigma(1), dd(1), mesh%np, MPI_FLOAT, MPI_SUM)
-    rho_sigma(1:mesh%np) = dd(1:mesh%np)
+  ! Comparing to KLI paper 1990, oep%vxc corresponds to V_{x \sigma}^S in Eq. 8
+  ! The n_{i \sigma} in Eq. 8 is partitioned in this code into \psi^{*} (included in lxc) and \psi (explicitly below)
+  if (st%parallel_in_states .or. st%d%kpt%parallel) then
+    call comm_allreduce(st%st_kpt_mpi_grp, rho_sigma)
   end if
 
+  ! Add a lower bound to the density to avoid numerical issues
+  do ispin = 1, st%d%spin_channels
+    do ip = 1, mesh%np
+      rho_sigma(ip, ispin) = max(rho_sigma(ip, ispin), CNST(1e-20))
+    end do
+  end do
+
   ! Comparing to KLI paper 1990, oep%vxc corresponds to V_{x \sigma}^S in Eq. 8
   ! The n_{i \sigma} in Eq. 8 is partitioned in this code into \psi^{*} (included in lxc) and \psi (explicitly below)
 
-  oep%vxc(1:mesh%np, is) = CNST(0.0)
-
+  do ik = st%d%kpt%start, st%d%kpt%end
+    ispin = st%d%get_spin_index(ik)
     do ist = st%st_start, st%st_end
-    call states_elec_get_state(st, mesh, ist, is, psi)
+      call states_elec_get_state(st, mesh, ist, ik, psi)
       do ip = 1, mesh%np
-      oep%vxc(ip, is) = oep%vxc(ip, is) + oep%socc*st%occ(ist, is)*R_REAL(oep%X(lxc)(ip, ist, is)*psi(ip, 1))
+        oep%vxc(ip, ispin) = oep%vxc(ip, ispin) + oep%socc * st%occ(ist, ik) * st%d%kweights(ik) &
+          * R_REAL(oep%X(lxc)(ip, 1, ist, ik)*psi(ip, 1))
       end do
     end do
+  end do
+  SAFE_DEALLOCATE_A(psi)
+
+  if(st%parallel_in_states .or. st%d%kpt%parallel) then
+    call comm_allreduce(st%st_kpt_mpi_grp, oep%vxc)
+  end if
+
 
+  do ispin = 1, st%d%spin_channels
     do ip = 1, mesh%np
-    oep%vxc(ip, is) = oep%vxc(ip, is)/rho_sigma(ip)
+      oep%vxc(ip, ispin) = oep%vxc(ip, ispin)/rho_sigma(ip, ispin)
+    end do
   end do
 
-  SAFE_DEALLOCATE_A(psi)
+  ! We have now computed the Slater part. We are adding the explicit KLI part now
 
-  if (st%parallel_in_states) then
-    call st%mpi_grp%allreduce(oep%vxc(1,is), dd(1), mesh%np, MPI_FLOAT, MPI_SUM)
-    oep%vxc(1:mesh%np,is) = dd(1:mesh%np)
-    SAFE_DEALLOCATE_A(dd)
-  end if
 
+  SAFE_ALLOCATE(kli(1:mesh%np, 1:st%d%nspin))
+  kli = M_ZERO
+
+  do ik = st%d%kpt%start, st%d%kpt%end
+    ispin = st%d%get_spin_index(ik)
+
+    ! In the case of isolated system, the OEP equation is simply solved for N-1 electrons
+    ! whereas for solids we remove a constant, as given in PRB 93, 205205 (2016)
+    if(space%is_periodic()) then
+      eigen_n = st%nst
+    else
+      call xc_oep_AnalyzeEigen(oep, st, ik)
       eigen_n = oep%eigen_n
+    end if
 
     SAFE_ALLOCATE(v_bar_S(1:st%nst))
     v_bar_S = M_ZERO
     do ist = st%st_start, st%st_end
-    if (st%occ(ist, is) > M_EPSILON) then
-      v_bar_S(ist) = dmf_dotp(mesh, sqphi(:, 1, ist) , oep%vxc(:,is), reduce = .false.)
+      if (st%occ(ist, ik) > M_EPSILON) then
+        v_bar_S(ist) = dmf_dotp(mesh, sqphi(:, ist, ik), oep%vxc(:,ispin), reduce = .false.)
       end if
     end do
     if (mesh%parallel_in_domains) call mesh%allreduce(v_bar_S, dim = st%st_end)
@@ -98,11 +136,12 @@ subroutine X(xc_KLI_solve) (mesh, st, is, oep, first)
       end do
       do ist = 1, eigen_n
         kssi = oep%eigen_index(ist)
-      call st%mpi_grp%bcast(sqphi(1, 1, kssi), mesh%np, MPI_FLOAT, st%node(kssi))
+        call st%mpi_grp%bcast(sqphi(1, kssi, ik), mesh%np, MPI_FLOAT, st%node(kssi))
       end do
     end if
+
     ! If there is more than one state, then solve linear equation.
-  linear_equation: if (eigen_n > 0) then
+    if (eigen_n > 0) then
       SAFE_ALLOCATE(dd(1:mesh%np))
       SAFE_ALLOCATE(xx(1:eigen_n, 1))
       SAFE_ALLOCATE(Ma(1:eigen_n, 1:eigen_n))
@@ -111,33 +150,49 @@ subroutine X(xc_KLI_solve) (mesh, st, is, oep, first)
       yy = M_ZERO
       Ma = M_ZERO
       dd = M_ZERO
-    proc = st%mpi_grp%rank
 
       i_loop: do ist = 1, eigen_n
+        if(space%is_periodic()) then
+          kssi = ist
+        else
           kssi = oep%eigen_index(ist)
-      if (proc == st%node(kssi)) then
-        dd(1:mesh%np) = sqphi(1:mesh%np, 1, kssi) / rho_sigma(1:mesh%np)
+        end if
+        if (st%mpi_grp%rank == st%node(kssi)) then
+          dd(1:mesh%np) = sqphi(1:mesh%np, kssi, ik) / rho_sigma(1:mesh%np, ispin)
+
           j_loop: do jst = ist, eigen_n
+            if(space%is_periodic()) then
+              kssj = jst
+            else
               kssj = oep%eigen_index(jst)
-          Ma(ist, jst) = - dmf_dotp(mesh, dd, sqphi(:, 1, kssj))
+            end if
+            Ma(ist, jst) = - dmf_dotp(mesh, dd, sqphi(:, kssj, ik))
           end do j_loop
           Ma(ist, ist) = M_ONE + Ma(ist, ist)
-        yy(ist, 1) = v_bar_S(kssi) - oep%uxc_bar(kssi, is)
+
+          yy(ist, 1) = v_bar_S(kssi) - oep%uxc_bar(1, kssi, ik)
+
         end if
       end do i_loop
 
+#if defined(HAVE_MPI)
       if (st%parallel_in_states) then
         do ist = 1, eigen_n
+          if(space%is_periodic()) then
+            kssi = ist
+          else
             kssi = oep%eigen_index(ist)
+          end if
           call  st%mpi_grp%bcast(yy(ist, 1), 1, MPI_FLOAT, st%node(kssi))
           do jst = 1, eigen_n
             call st%mpi_grp%bcast(Ma(ist, jst), 1, MPI_FLOAT, st%node(kssi))
           end do
         end do
       end if
+#endif
 
       do ist = 1, eigen_n
-      do jst = ist, eigen_n
+        do jst = ist+1, eigen_n
           Ma(jst, ist) = Ma(ist, jst)
         end do
       end do
@@ -145,9 +200,14 @@ subroutine X(xc_KLI_solve) (mesh, st, is, oep, first)
       call lalg_linsyssolve(eigen_n, 1, Ma, yy, xx)
 
       do ist = 1, eigen_n
+        if(space%is_periodic()) then
+          kssi = ist
+        else
           kssi = oep%eigen_index(ist)
-      oep%vxc(1:mesh%np,is) = oep%vxc(1:mesh%np,is) + &
-        oep%socc * st%occ(kssi, is) * xx(ist, 1) * sqphi(1:mesh%np, 1, kssi) / rho_sigma(1:mesh%np)
+        end if
+
+        kli(1:mesh%np,ispin) = kli(1:mesh%np,ispin) + st%d%kweights(ik) * &
+          oep%socc * st%occ(kssi, ik) * xx(ist, 1) * sqphi(1:mesh%np, kssi, ik) / rho_sigma(1:mesh%np, ispin)
       end do
 
       SAFE_DEALLOCATE_A(dd)
@@ -155,10 +215,46 @@ subroutine X(xc_KLI_solve) (mesh, st, is, oep, first)
       SAFE_DEALLOCATE_A(Ma)
       SAFE_DEALLOCATE_A(yy)
 
-  end if linear_equation
+    end if
     ! The previous stuff is only needed if eigen_n>0.
-
     SAFE_DEALLOCATE_A(v_bar_S)
+
+  end do
+
+  if(st%d%kpt%parallel) then
+    call comm_allreduce(st%d%kpt%mpi_grp, kli)
+  end if
+
+  !Adding the KLI part
+  call lalg_axpy(mesh%np, st%d%nspin, M_ONE, kli, oep%vxc)
+  SAFE_DEALLOCATE_A(kli)
+
+  ! Substracting the constant, see PRB 93, 205205 (2016)
+  if(space%is_periodic()) then
+    SAFE_ALLOCATE(dd(1:mesh%np))
+    cnst = M_ZERO
+    do ik = st%d%kpt%start, st%d%kpt%end
+      ispin = st%d%get_spin_index(ik)
+      do ist = st%st_start, st%st_end
+        if (st%occ(ist, ik) > M_EPSILON) then
+          dd(:) = sqphi(:, ist, ik) / rho_sigma(:, ispin)
+          cnst(ispin) = cnst(ispin) + oep%socc * st%occ(ist, ik) * st%d%kweights(ik) &
+            * dmf_dotp(mesh, sqphi(:, ist, ik), oep%vxc(:,ispin)) * dmf_integrate(mesh, dd)
+        end if
+      end do
+    end do
+
+    if (st%parallel_in_states .or. st%d%kpt%parallel) then
+      call comm_allreduce(st%st_kpt_mpi_grp, cnst)
+    end if
+
+    do ispin = 1, st%d%spin_channels
+      oep%vxc(:, ispin) = oep%vxc(:, ispin) - cnst(ispin)/rcell_volume
+    end do
+    SAFE_DEALLOCATE_A(dd)
+  end if
+
+
   SAFE_DEALLOCATE_A(rho_sigma)
   SAFE_DEALLOCATE_A(sqphi)
   POP_SUB(X(xc_KLI_solve))
@@ -240,8 +336,8 @@ subroutine X(xc_KLI_solve_photon) (namespace, mesh, hm, st, is, oep, first)
     call states_elec_get_state(st, mesh, ist, is, psi)
 #ifdef R_TREAL
     if (ist>(oep%eigen_n + 1)) exit ! included to guarantee that the photonic KLI finishes correctly but the parallel in states feature of the normal KLI works still
-    bb(:,1) = oep%X(lxc)(1:mesh%np, ist, is)
-    if (ist /= oep%eigen_n + 1) bb(:,1) = bb(:,1) - oep%uxc_bar(ist, is)*R_CONJ(psi(:, 1))
+    bb(:,1) = oep%X(lxc)(1:mesh%np, 1, ist, is)
+    if (ist /= oep%eigen_n + 1) bb(:,1) = bb(:,1) - oep%uxc_bar(1, ist, is)*R_CONJ(psi(:, 1))
     if (.not. first) then
       call xc_oep_pt_rhs(mesh, st, is, oep, phi1, ist, bb)
     end if

src/electrons/xc_kli_pauli_inc.F90

@@ -54,9 +54,9 @@ subroutine xc_kli_pauli_solve(mesh, st, oep)
       call states_elec_get_state(st, mesh, ist, ik, psi)
       !Here we accumulate the result for the potential
       do ip = 1, mesh%np
-        bij(ip, 1) = bij(ip, 1) + weight * TOFLOAT(conjg(oep%zlxc(ip, ist, 1))*conjg(psi(ip, 1)))
-        bij(ip, 2) = bij(ip, 2) + weight * TOFLOAT(conjg(oep%zlxc(ip, ist, 2))*conjg(psi(ip, 2)))
-        bij(ip, 3) = bij(ip, 3) + weight * conjg(oep%zlxc(ip, ist, 1)) * conjg(psi(ip, 2))
+        bij(ip, 1) = bij(ip, 1) + weight * TOFLOAT(conjg(oep%zlxc(ip, 1, ist, 1))*conjg(psi(ip, 1)))
+        bij(ip, 2) = bij(ip, 2) + weight * TOFLOAT(conjg(oep%zlxc(ip, 2, ist, 1))*conjg(psi(ip, 2)))
+        bij(ip, 3) = bij(ip, 3) + weight * conjg(oep%zlxc(ip, 1, ist, 1)) * conjg(psi(ip, 2))
         !The last component is simply the complex conjuguate of bij(3), so we do not compute it.
 
         ! We store \phi_{i,\alpha}(r)\phi^*_{i,\beta}(r). Needed for the KLI part
@@ -270,7 +270,7 @@ subroutine xc_kli_pauli_solve(mesh, st, oep)
         Ma(ist, jst) = Ma(ist, jst) - M_FOUR * dmf_dotp(mesh, 2, dd(:,3:4, jst), sqphi(:, 3:4, kssi))
       end do j_loop
       Ma(ist, ist) = M_ONE + Ma(ist, ist)
-      yy(ist, 1) = M_TWO * v_bar_S(kssi) - M_TWO * (oep%uxc_bar(kssi, 1) + oep%uxc_bar(kssi, 2))
+      yy(ist, 1) = M_TWO * v_bar_S(kssi) - M_TWO * (oep%uxc_bar(1, kssi, 1) + oep%uxc_bar(2, kssi, 1))
 
     end do i_loop
 

src/electrons/xc_oep.F90

@@ -96,9 +96,9 @@ module xc_oep_oct_m
     integer                      :: eigen_n
     integer, allocatable         :: eigen_type(:), eigen_index(:)
     FLOAT                        :: socc, sfact
-    FLOAT,   allocatable, public :: vxc(:,:), uxc_bar(:,:)
-    FLOAT,   allocatable         :: dlxc(:, :, :)
-    CMPLX,   allocatable         :: zlxc(:, :, :)
+    FLOAT,   allocatable, public :: vxc(:,:), uxc_bar(:,:,:)
+    FLOAT,   allocatable         :: dlxc(:, :, :, :)
+    CMPLX,   allocatable         :: zlxc(:, :, :, :)
     integer                      :: mixing_scheme
     type(photon_mode_t),  public :: pt
     type(lr_t)                   :: photon_lr     !< to solve the equation H psi = b
@@ -172,6 +172,13 @@ contains
       if (states_are_complex(st)) then
         call messages_not_implemented('Photon OEP with complex wavefunctions')
       end if
+
+      if (st%d%nik > st%d%ispin .and. oep%level == OEP_LEVEL_FULL) then
+        call messages_not_implemented("Full OEP for periodic systems", namespace=namespace)
+      end if
+      if (st%d%nik > st%d%ispin .and. st%d%ispin==SPINORS) then
+        call messages_not_implemented("OEP for periodic systems with spinors", namespace=namespace)
+      end if
     end if
 
     if (oep%level == OEP_LEVEL_FULL) then
@@ -222,11 +229,6 @@ contains
       oep%norm2ss = M_ZERO
     end if
 
-    ! this routine is only prepared for finite systems. (Why not?)
-    if (st%d%nik > st%d%ispin) then
-      call messages_not_implemented("OEP for periodic systems", namespace=namespace)
-    end if
-
     ! obtain the spin factors
     call xc_oep_SpinFactor(oep, st%d%nspin)
 
@@ -264,8 +266,8 @@ contains
       call messages_experimental("OEP with spinors")
     end if
 
-    if (st%d%kpt%parallel) then
-      call messages_not_implemented("OEP parallel in spin/k-points", namespace=namespace)
+    if (st%d%kpt%parallel .and. oep%level == OEP_LEVEL_FULL) then
+      call messages_experimental("OEP parallel in spin/k-points", namespace=namespace)
     end if
 
     POP_SUB(xc_oep_init)

src/electrons/xc_oep_inc.F90

@@ -53,9 +53,9 @@ subroutine X(xc_oep_calc)(oep, namespace, xcs, gr, hm, st, space, rcell_volume,
   nspin_ = min(st%d%nspin, 2)
   SAFE_ALLOCATE(oep%eigen_type (1:st%nst))
   SAFE_ALLOCATE(oep%eigen_index(1:st%nst))
-  SAFE_ALLOCATE(oep%X(lxc)(1:gr%np, st%st_start:st%st_end, 1:nspin_))
+  SAFE_ALLOCATE(oep%X(lxc)(1:gr%np, 1:st%d%dim, st%st_start:st%st_end, st%d%kpt%start:st%d%kpt%end))
   oep%X(lxc) = M_ZERO
-  SAFE_ALLOCATE(oep%uxc_bar(1:st%nst, 1:nspin_))
+  SAFE_ALLOCATE(oep%uxc_bar(1:st%d%dim, 1:st%nst, st%d%kpt%start:st%d%kpt%end))
 
   !We first apply the exchange operator to all the states
   call xst%nullify()
@@ -82,7 +82,8 @@ subroutine X(xc_oep_calc)(oep, namespace, xcs, gr, hm, st, space, rcell_volume,
 
     do ik = st%d%kpt%start, st%d%kpt%end
       do ib = st%group%block_start, st%group%block_end
-        ! Here we apply the term by calling hamiltonian_elec_apply_batch, which takes care of setting the        ! phases properly, as well as upating boundary points/ghost points when needed
+        ! Here we apply the term by calling hamiltonian_elec_apply_batch, which takes care of setting the
+        ! phases properly, as well as upating boundary points/ghost points when needed
         call X(hamiltonian_elec_apply_batch)(hm, namespace, gr, st%group%psib(ib, ik), &
           xst%group%psib(ib, ik), terms = TERM_MGGA)
       end do
@@ -92,6 +93,10 @@ subroutine X(xc_oep_calc)(oep, namespace, xcs, gr, hm, st, space, rcell_volume,
     ! this part handles the (pure) orbital functionals
     ! SIC a la PZ is handled here
   case(OEP_TYPE_SIC)
+    ! this routine is only prepared for finite systems. (Why not?)
+    if (st%d%nik > st%d%ispin) then
+      call messages_not_implemented("OEP-SIC for periodic systems", namespace=namespace)
+    end if
     spin: do is = 1, nspin_
       ! distinguish between 'is' being the spin_channel index (collinear)
       ! and being the spinor (noncollinear)
@@ -115,37 +120,29 @@ subroutine X(xc_oep_calc)(oep, namespace, xcs, gr, hm, st, space, rcell_volume,
   SAFE_ALLOCATE(psi(1:gr%np))
   SAFE_ALLOCATE(xpsi(1:gr%np))
 
-  oep%uxc_bar(:, :) = M_ZERO
-  do is = 1, nspin_
-    ! distinguish between 'is' being the spin_channel index (collinear)
-    ! and being the spinor (noncollinear)
-    if (st%d%ispin == SPINORS) then
-      isp = 1
-      idm = is
-    else
-      isp = is
-      idm = 1
-    end if
-
+  oep%uxc_bar(:, :, :) = M_ZERO
+  do ik = st%d%kpt%start, st%d%kpt%end
     do ist = st%st_start, st%st_end
-      call states_elec_get_state(st, gr, idm, ist, isp, psi)
+      do idm = 1, st%d%dim
+        call states_elec_get_state(st, gr, idm, ist, ik, psi)
         if (exx) then
           ! Here we copy the state from xst to X(lxc).
           ! This will be removed in the future, but it allows to keep both EXX and PZ-SIC in the code
-        call states_elec_get_state(xst, gr, idm, ist, isp, xpsi)
+          call states_elec_get_state(xst, gr, idm, ist, ik, xpsi)
           ! There is a complex conjugate here, as the lxc is defined as <\psi|X and
           ! exchange_operator_compute_potentials returns X|\psi>
 #ifndef R_TREAL
           xpsi = R_CONJ(xpsi)
 #endif
-        call lalg_axpy(gr%np, M_ONE, xpsi, oep%X(lxc)(1:gr%np, ist, is))
+          call lalg_axpy(gr%np, M_ONE, xpsi, oep%X(lxc)(1:gr%np, idm, ist, ik))
         end if
-
-      oep%uxc_bar(ist, is) = R_REAL(X(mf_dotp)(gr, psi, oep%X(lxc)(1:gr%np, ist, is), reduce = .false., dotu = .true.))
+        oep%uxc_bar(idm, ist, ik) = R_REAL(X(mf_dotp)(gr, psi, oep%X(lxc)(1:gr%np, idm, ist, ik), reduce = .false., dotu = .true.))
       end do
-
-    if (gr%parallel_in_domains) call gr%allreduce(oep%uxc_bar(1:st%st_end, is), dim = st%st_end)
     end do
+  end do
+  if (gr%parallel_in_domains) then
+    call gr%allreduce(oep%uxc_bar(1:st%d%dim, 1:st%st_end, st%d%kpt%start:st%d%kpt%end))
+  end if
 
   SAFE_DEALLOCATE_A(psi)
   SAFE_DEALLOCATE_A(xpsi)
@@ -154,17 +151,11 @@ subroutine X(xc_oep_calc)(oep, namespace, xcs, gr, hm, st, space, rcell_volume,
 
   if (st%parallel_in_states) then
     call st%mpi_grp%barrier()
-    if (st%d%ispin == SPIN_POLARIZED) then
-      do isp = 1, 2
+    do ik = st%d%kpt%start, st%d%kpt%end
       do ist = 1, st%nst
-          call st%mpi_grp%bcast(oep%uxc_bar(ist, isp), 1, MPI_FLOAT, st%node(ist))
+        call st%mpi_grp%bcast(oep%uxc_bar(1, ist, ik), st%d%dim, MPI_FLOAT, st%node(ist))
       end do
     end do
-    else
-      do ist = 1, st%nst
-        call st%mpi_grp%bcast(oep%uxc_bar(ist, 1), 1, MPI_FLOAT, st%node(ist))
-      end do
-    end if
   end if
 
   if (st%d%ispin == SPINORS) then
@@ -175,39 +166,37 @@ subroutine X(xc_oep_calc)(oep, namespace, xcs, gr, hm, st, space, rcell_volume,
     end if
     ! full OEP not implemented!
   else
-    spin2: do is = 1, nspin_
+    if (oep%level == OEP_LEVEL_FULL .and. (.not. first)) then
+      do is = 1, nspin_
         ! get the HOMO state
         call xc_oep_AnalyzeEigen(oep, st, is)
-      !
-      if (present(vxc)) then
-        ! solve the KLI equation
-        if (oep%level /= OEP_LEVEL_FULL .or. first) then
-          oep%vxc = M_ZERO
-          if (oep%type == OEP_TYPE_PHOTONS) then
-            call X(xc_KLI_solve_photon) (namespace, gr, hm, st, is, oep, first)
-          else
-            call X(xc_KLI_solve) (gr, st, is, oep, first)
-          end if
-          if (present(vxc)) then
-            vxc(1:gr%np, is) = vxc(1:gr%np, is) + oep%vxc(1:gr%np, is)
-          end if
-        end if
-        ! if asked, solve the full OEP equation
-        if (oep%level == OEP_LEVEL_FULL .and. (.not. first)) then
 
+        if (present(vxc)) then
           if (oep%type == OEP_TYPE_PHOTONS) then
             call X(xc_oep_solve_photon)(namespace, gr, hm, st, is, vxc(:,is), oep)
           else
             call X(xc_oep_solve)(namespace, gr, hm, st, is, vxc(:,is), oep)
           end if
 
-          if (present(vxc)) then
           call lalg_axpy(gr%np, M_ONE, oep%vxc(1:gr%np, is), vxc(1:gr%np, is))
         end if
+      end do
+
+    else !KLI
+      if (oep%type == OEP_TYPE_PHOTONS) then
+        oep%vxc = M_ZERO
+        do is = 1, nspin_
+          call xc_oep_AnalyzeEigen(oep, st, is)
+          call X(xc_KLI_solve_photon) (namespace, gr, hm, st, is, oep, first)
+          call lalg_axpy(gr%np, M_ONE, oep%vxc(1:gr%np, is), vxc(1:gr%np, is))
+        end do
+      else
+        call X(xc_KLI_solve) (space, gr, st, oep, rcell_volume)
+        call lalg_axpy(gr%np, st%d%nspin, M_ONE, oep%vxc, vxc)
       end if
-        if (is == nspin_) first = .false.
+
     end if
-    end do spin2
+    first = .false.
   end if
   SAFE_DEALLOCATE_A(oep%eigen_type)
   SAFE_DEALLOCATE_A(oep%eigen_index)
@@ -275,8 +264,8 @@ subroutine X(xc_oep_solve) (namespace, mesh, hm, st, is, vxc, oep)
       vxc_bar = dmf_integrate(mesh, psi2(:, 1)*oep%vxc(1:mesh%np, is))
 
       ! This the right-hand side of Eq. 21
-      bb(1:mesh%np, 1) = -(oep%vxc(1:mesh%np, is) - (vxc_bar - oep%uxc_bar(ist, is)))* &
-        R_CONJ(psi(:, 1)) + oep%X(lxc)(1:mesh%np, ist, is)
+      bb(1:mesh%np, 1) = -(oep%vxc(1:mesh%np, is) - (vxc_bar - oep%uxc_bar(1, ist, is)))* &
+        R_CONJ(psi(:, 1)) + oep%X(lxc)(1:mesh%np, 1, ist, is)
 
       call X(lr_orth_vector) (mesh, st, bb, ist, is, R_TOTYPE(M_ZERO))
 
@@ -310,7 +299,7 @@ subroutine X(xc_oep_solve) (namespace, mesh, hm, st, is, vxc, oep)
         call states_elec_get_state(st, mesh, ist, is, psi)
         psi2(:, 1) = TOFLOAT(R_CONJ(psi(:, 1))*psi(:,1))
         vxc_bar = dmf_integrate(mesh, psi2(:, 1)*oep%vxc(1:mesh%np, is))
-        oep%vxc(1:mesh%np,is) = oep%vxc(1:mesh%np,is) - (vxc_bar - oep%uxc_bar(ist,is))
+        oep%vxc(1:mesh%np,is) = oep%vxc(1:mesh%np,is) - (vxc_bar - oep%uxc_bar(1, ist,is))
       end if
     end do
 
@@ -424,8 +413,8 @@ subroutine X(xc_oep_solve_photon) (namespace, mesh, hm, st, is, vxc, oep)
       vxc_bar = dmf_integrate(mesh, psi2(:, 1)*oep%vxc(1:mesh%np, is))
 
       ! This the right-hand side of Eq. 21
-      bb(1:mesh%np, 1) = -(oep%vxc(1:mesh%np, is) - (vxc_bar - oep%uxc_bar(ist, is)))* &
-        R_CONJ(psi(:, 1)) + oep%X(lxc)(1:mesh%np, ist, is)
+      bb(1:mesh%np, 1) = -(oep%vxc(1:mesh%np, is) - (vxc_bar - oep%uxc_bar(1, ist, is)))* &
+        R_CONJ(psi(:, 1)) + oep%X(lxc)(1:mesh%np, 1, ist, is)
 
 #ifdef R_TREAL
       call xc_oep_pt_rhs(mesh, st, is, oep, phi1, ist, bb)
@@ -475,7 +464,7 @@ subroutine X(xc_oep_solve_photon) (namespace, mesh, hm, st, is, vxc, oep)
 #ifdef R_TREAL
         call xc_oep_pt_uxcbar(mesh, st, is, oep, phi1, ist, vxc_bar)
 #endif
-        oep%vxc(1:mesh%np,is) = oep%vxc(1:mesh%np,is) - (vxc_bar - oep%uxc_bar(ist,is))
+        oep%vxc(1:mesh%np,is) = oep%vxc(1:mesh%np,is) - (vxc_bar - oep%uxc_bar(1, ist,is))
       end if
     end do
 

src/electrons/xc_oep_sic_inc.F90

@@ -69,7 +69,7 @@ subroutine X(oep_sic) (xcs, gr, psolver, namespace, space, rcell_volume, st, kpo
       ex_ = ex_ - oep%sfact*ex2
       ec_ = ec_ - oep%sfact*ec2
 
-      oep%X(lxc)(1:gr%np, ist, is) = oep%X(lxc)(1:gr%np, ist, is) - vxc(1:gr%np, 1)*R_CONJ(psi(1:gr%np, 1))
+      oep%X(lxc)(1:gr%np, 1, ist, is) = oep%X(lxc)(1:gr%np, 1, ist, is) - vxc(1:gr%np, 1)*R_CONJ(psi(1:gr%np, 1))
 
       ! calculate the Hartree contribution using Poisson equation
       vxc(1:gr%np, 1) = M_ZERO
@@ -78,7 +78,7 @@ subroutine X(oep_sic) (xcs, gr, psolver, namespace, space, rcell_volume, st, kpo
       ! The exchange energy.
       ex_ = ex_ - M_HALF*oep%sfact*dmf_dotp(gr, vxc(1:gr%np, 1), rho(1:gr%np, 1))
 
-      oep%X(lxc)(1:gr%np, ist, is) = oep%X(lxc)(1:gr%np, ist, is) - &
+      oep%X(lxc)(1:gr%np, 1, ist, is) = oep%X(lxc)(1:gr%np, 1, ist, is) - &
         vxc(1:gr%np, 1)*R_CONJ(psi(1:gr%np, 1))
     end if
   end do

src/math/fft.F90

@@ -394,14 +394,7 @@ contains
       end if
 
       do ii = 1, fft_dim
-        ! the AMD OpenCL FFT only supports sizes 2, 3 and 5, but size
-        ! 5 gives an fpe error on the Radeon 7970 (APPML 1.8), so we
-        ! only use factors 2 and 3
-#ifdef HAVE_CLFFT
-        nn_temp(ii) = fft_size(nn(ii), (/2, 3/), optimize_parity(ii))
-#else
         nn_temp(ii) = fft_size(nn(ii), (/2, 3, 5, 7/), optimize_parity(ii))
-#endif
         if (fft_optimize .and. optimize(ii)) nn(ii) = nn_temp(ii)
       end do
 

testsuite/finite_systems_3d/12-forces.test

@@ -58,5 +58,5 @@ Precision: 2.75e-07
 match ;  Force Local   ; LINEFIELD(static/forces, 2, 12) ; 5.501625519999999
 Precision: 1.49e-07
 match ;  Force NL      ; LINEFIELD(static/forces, 2, 15) ; 2.97727227
-Precision: 2.36e-14
-match ;  Force SCF     ; LINEFIELD(static/forces, 2, 24) ; -2.0183700199999998e-08
+Precision: 5.00e-14
+match ;  Force SCF     ; LINEFIELD(static/forces, 2, 24) ; -2.0183716e-08

testsuite/functionals/10-vdw_d3_dna.test

@@ -191,7 +191,7 @@ Precision: 2.06e-08
 match ;  Force 50 (z)    ; GREPFIELD(static/info, '50         H', 5) ; 0.006280363215
 Precision: 1.39e-08
 match ;  Force 60 (x)    ; GREPFIELD(static/info, '60         H', 3) ; -0.0187789085
-Precision: 1.37e-09
+Precision: 1.7e-09
 match ;  Force 60 (y)    ; GREPFIELD(static/info, '60         H', 4) ; 0.000414181492
 Precision: 4.84e-09
 match ;  Force 60 (z)    ; GREPFIELD(static/info, '60         H', 5) ; 0.0110064259

testsuite/periodic_systems/28-mgga_kli.01-Si_scan.inp

+CalculationMode = gs
+FromScratch = yes
+PeriodicDimensions = 3
+BoxShape = parallelepiped
+
+a = 10.26120
+%LatticeParameters
+ a | a | a
+%
+
+%LatticeVectors
+ 0.0 | 0.5 | 0.5
+ 0.5 | 0.0 | 0.5
+ 0.5 | 0.5 | 0.0
+%
+
+%ReducedCoordinates
+  "Si" |   0.0       | 0.0       | 0.0
+  "Si" |   1/4       | 1/4       | 1/4
+%
+
+Spacing = 0.45
+%KPointsGrid
+9 | 9 | 9
+%
+
+KPointsUseSymmetries = yes
+ExperimentalFeatures = yes
+
+ExtraStates = 2
+
+TheoryLevel = kohn_sham
+XCFunctional = mgga_x_revscan + mgga_c_revscan
+
+ConvRelEv = 1e-8
+
+FilterPotentials = filter_none
+

testsuite/periodic_systems/28-mgga_kli.test

+# -*- coding: utf-8 mode: shell-script -*-
+
+Test       : MGGA within Kohn-Sham DFT at the KLI level with k-points
+Program    : octopus
+TestGroups : long-run, periodic_systems
+Enabled    : no-GPU-MPI
+
+Input      : 28-mgga_kli.01-Si_scan.inp
+
+match ; SCF convergence ; GREPCOUNT(static/info, 'SCF converged') ; 1
+
+match ; Total k-points   ; GREPFIELD(static/info, 'Total number of k-points', 6) ; 729
+match ; Space group        ; GREPFIELD(out, 'Space group', 4) ; 227
+match ; No. of symmetries  ; GREPFIELD(out, 'symmetries that can be used', 5)  ;  24
+
+Precision: 3.98e-07
+match ;  Total energy        ; GREPFIELD(static/info, 'Total       =', 3) ; -7.95879216
+Precision: 3.91e-07
+match ;  Ion-ion energy      ; GREPFIELD(static/info, 'Ion-ion     =', 3) ; -7.81793472
+Precision: 5.55e-07
+match ;  Eigenvalues sum     ; GREPFIELD(static/info, 'Eigenvalues =', 3) ; -0.663380765
+Precision: 7.64e-07
+match ;  Hartree energy      ; GREPFIELD(static/info, 'Hartree     =', 3) ; 0.565675935
+Precision: 4.35e-07
+match ;  Exchange energy     ; GREPFIELD(static/info, 'Exchange    =', 3) ; -2.151788
+Precision: 1.39e-07
+match ;  Correlation energy  ; GREPFIELD(static/info, 'Correlation =', 3) ; -0.27781913
+Precision: 9.63e-07
+match ;  Kinetic energy      ; GREPFIELD(static/info, 'Kinetic     =', 3) ; 3.0580798399999995
+Precision: 1.39e-06
+match ;  External energy     ; GREPFIELD(static/info, 'External    =', 3) ; -1.33519501
+
+Precision: 1.00e-04
+match ;  k-point 1 (x)  ; GREPFIELD(static/info, '#k =       1', 7) ; 0.0
+match ;  k-point 1 (y)  ; GREPFIELD(static/info, '#k =       1', 8) ; 0.0
+match ;  k-point 1 (z)  ; GREPFIELD(static/info, '#k =       1', 9) ; 0.0
+
+Precision: 1.67e-05
+match ;  Eigenvalue  1  ; GREPFIELD(static/info, '#k =       1', 3, 1) ; -0.334804
+Precision: 4.95e-05
+match ;  Eigenvalue  2  ; GREPFIELD(static/info, '#k =       1', 3, 2) ; 0.098998
+Precision: 5.03e-06
+match ;  Eigenvalue  4  ; GREPFIELD(static/info, '#k =       1', 3, 4) ; 0.100506
+Precision: 1.00e-05
+match ;  Eigenvalue  5  ; GREPFIELD(static/info, '#k =       1', 3, 5) ; 0.199966

testsuite/periodic_systems/Makefile.am

@@ -133,6 +133,8 @@ dist_share_DATA =                                \
         27-Ar.01-gs.inp                          \
         27-Ar.02-kdotp.inp                       \
         27-Ar.03-em_resp_mo.inp                  \
-        27-Ar.test
+        27-Ar.test                               \
+        28-mgga_kli.01-Si_scan.inp               \
+        28-mgga_kli.test
 
 CLEANFILES = *~ *.bak

testsuite/real_time/09-angular_momentum.test

@@ -39,7 +39,7 @@ Precision: 4.13e-15
 match ;  Ly [step  25]      ; LINEFIELD(td.general/angular, -76, 4) ; -0.0718464109480849
 Precision: 4.19e-15
 match ;  Ly [step  50]      ; LINEFIELD(td.general/angular, -51, 4) ; -0.00711776925568968
-Precision: 4.40e-15
+Precision: 5.10e-15
 match ;  Ly [step  75]      ; LINEFIELD(td.general/angular, -26, 4) ; -0.0417208122257626
 Precision: 5.72e-15
 match ;  Ly [step 100]      ; LINEFIELD(td.general/angular, -1, 4) ; -0.057772136201210605
@@ -127,9 +127,9 @@ Precision: 1.00e-04
 match ;  Lx [step   1]      ; LINEFIELD(td.general/angular, -101, 3) ; 0.0
 Precision: 5.11e-15
 match ;  Lx [step  25]      ; LINEFIELD(td.general/angular, -76, 3) ; 0.04378813412511501
-Precision: 4.39e-15
+Precision: 5.40e-15
 match ;  Lx [step  50]      ; LINEFIELD(td.general/angular, -51, 3) ; 0.007223501819495281
-Precision: 5.34e-15
+Precision: 8.00e-15
 match ;  Lx [step  75]      ; LINEFIELD(td.general/angular, -26, 3) ; -0.019482698210286253
 Precision: 6.00e-15
 match ;  Lx [step 100]      ; LINEFIELD(td.general/angular, -1, 3) ; -0.0304720312519288
@@ -140,7 +140,7 @@ Precision: 4.13e-15
 match ;  Ly [step  25]      ; LINEFIELD(td.general/angular, -76, 4) ; -0.0718464109480849
 Precision: 4.19e-15
 match ;  Ly [step  50]      ; LINEFIELD(td.general/angular, -51, 4) ; -0.00711776925568968
-Precision: 4.40e-15
+Precision: 5.10e-15
 match ;  Ly [step  75]      ; LINEFIELD(td.general/angular, -26, 4) ; -0.0417208122257626
 Precision: 5.72e-15
 match ;  Ly [step 100]      ; LINEFIELD(td.general/angular, -1, 4) ; -0.057772136201210605
@@ -153,7 +153,7 @@ Precision: 3.86e-13
 match ;  Lz [step  50]      ; LINEFIELD(td.general/angular, -51, 5) ; 0.00771299300973
 Precision: 5.61e-15
 match ;  Lz [step  75]      ; LINEFIELD(td.general/angular, -26, 5) ; 0.0720814863526959
-Precision: 4.73e-15
+Precision: 5.8-15
 match ;  Lz [step 100]      ; LINEFIELD(td.general/angular, -1, 5) ; 0.0511931067883849
 
 Util : oct-propagation_spectrum

testsuite/real_time/19-td_move_ions.test

@@ -51,5 +51,5 @@ Precision: 4.83e-17
 match ;  X Velocity Atom 1 [step 40]  ; LINEFIELD(td.general/coordinates, -1, 9) ; -0.00966789973133459
 Precision: 7.90e-12
 match ;  X Force Atom 1 [step 30]  ; LINEFIELD(td.general/coordinates, -11, 15) ; -15.803865466118
-Precision: 7.85e-14
-match ;  X Force Atom 1 [step 40]  ; LINEFIELD(td.general/coordinates, -1, 15) ; -15.707794025243349
+Precision: 1.60e-13
+match ;  X Force Atom 1 [step 40]  ; LINEFIELD(td.general/coordinates, -1, 15) ; -15.707794025243398

testsuite/symmetries/09-symmetrization_gga.test

@@ -44,7 +44,7 @@ Precision: 1.00e-01
 match ;  Partial charge  1  ; GREPFIELD(static/info, 'Partial ionic charges', 3, 2) ; 1.0
 Precision: 1.00e-01
 match ;  Partial charge  2  ; GREPFIELD(static/info, 'Partial ionic charges', 3, 3) ; 1.0
-Precision: 6.00e-15
+Precision: 7.00e-15
 match ;  Density value 1    ; LINEFIELD(static/density.y=0\,z=0, 2, 2) ; 0.00974594722143737
 match ;  Density value 2    ; LINEFIELD(static/density.y=0\,z=0, 3, 2) ; 0.00866615894594306
 match ;  Bader value 1      ; LINEFIELD(static/bader.y=0\,z=0, 6, 2) ; 0.00993816715475956