Fix OEP full - spin polarized case

Description

The OEP full was not taking the right occupations for the spin-polarized case.

News snippet

The OEP full was not taking the right occupations for the spin-polarized case.

Checklist

• I have checked that my code follows the Octopus coding standards
• I have added tests for all the new features added in this request.

src/electrons/exponential.F90

@@ -421,7 +421,7 @@ contains
     ! ---------------------------------------------------------
 
     ! ---------------------------------------------------------
-    !TODO: Add a reference
+    ! See the batch routine for the reference
     subroutine lanczos()
       integer ::  iter, l, idim
       CMPLX, allocatable :: hamilt(:,:), v(:,:,:), expo(:,:), psi(:, :)
@@ -756,7 +756,12 @@ contains
   end subroutine exponential_taylor_series_batch
 
   ! ---------------------------------------------------------
-  !TODO: Add a reference
+  ! Some details of the implementation can be understood from
+  ! Saad, Y. (1992). Analysis of some Krylov subspace approximations to the matrix exponential operator.
+  ! SIAM Journal on Numerical Analysis, 29(1), 209-228.
+  ! 
+  ! A pdf can be accessed here https://www-users.cse.umn.edu/~saad/PDF/RIACS-90-ExpTh.pdf
+  ! Equation numbers below refer to this paper
   subroutine exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat, vmagnus, inh_psib)
     type(exponential_t),                intent(inout) :: te
     type(namespace_t),                  intent(in)    :: namespace
@@ -812,20 +817,24 @@ contains
     ! This is the Lanczos loop...
     do iter = 1, te%exp_order
 
-      !to apply the Hamiltonian
+      ! to apply the Hamiltonian
       call operate_batch(hm, namespace, mesh, vb(iter), vb(iter+1), vmagnus)
 
+      ! We use either the Lanczos method (Hermitian case) or the Arnoldi method
       if (hm%is_hermitian()) then
         l = max(1, iter - 1)
       else
         l = 1
       end if
 
-      !orthogonalize against previous vectors
+      ! Orthogonalize against previous vectors
       call zmesh_batch_orthogonalization(mesh, iter - l + 1, vb(l:iter), vb(iter+1), &
         normalize = .false., overlap = hamilt(l:iter, iter, 1:psib%nst), norm = hamilt(iter + 1, iter, 1:psib%nst), &
         gs_scheme = te%arnoldi_gs)
 
+      ! We now need to compute exp(Hm), where Hm is the projection of the linear transformation
+      ! of the Hamiltonian onto the Krylov subspace Km
+      ! See Eq. 4
       do ii = 1, psib%nst
         if (present(inh_psib)) then
           call zlalg_phi(iter, -M_zI*deltat, hamilt(:,:,ii), expo(:,:,ii), hm%is_hermitian())
@@ -836,8 +845,9 @@ contains
         res(ii) = abs(hamilt(iter + 1, iter, ii) * abs(expo(iter, 1, ii)))
       end do !ii
 
+      ! We now estimate the error we made. This is given by the formula denoted Er2 in Sec. 5.2
       if (all(abs(hamilt(iter + 1, iter, :)) < CNST(1.0e4)*M_EPSILON)) exit ! "Happy breakdown"
-      !We normalize only if the norm is non-zero
+      ! We normalize only if the norm is non-zero
       ! see http://www.netlib.org/utk/people/JackDongarra/etemplates/node216.html#alg:arn0
       norm = M_ONE
       do ist = 1, psib%nst
@@ -862,13 +872,14 @@ contains
         call batch_axpy(mesh%np, TOFLOAT(deltat)*beta(1:psib%nst)*expo(ii,1,1:psib%nst), vb(ii), psib, a_full = .false.)
       end do
     else
+      ! See Eq. 4 for the expression here
       ! zpsi = nrm * V * expo(1:iter, 1) = nrm * V * expo * V^(T) * zpsi
       call batch_scal(mesh%np, expo(1,1,1:psib%nst), psib, a_full = .false.)
-      !TODO: We should have a routine batch_gemv fro improve performances
+      ! TODO: We should have a routine batch_gemv fro improve performances
       do ii = 2, iter
         call batch_axpy(mesh%np, beta(1:psib%nst)*expo(ii,1,1:psib%nst), vb(ii), psib, a_full = .false.)
-        !In order to apply the two exponentials, we mush store the eigenvales and eigenvectors given by zlalg_exp
-        !And to recontruct here the exp(i*dt*H) for deltat2
+        ! In order to apply the two exponentials, we must store the eigenvalues and eigenvectors given by zlalg_exp
+        ! And to recontruct here the exp(i*dt*H) for deltat2
       end do
     end if