Resolve "Enabling free Maxwell propagation in multisystem mode, part 2/4"

Description

Allow Maxwell propagation with external currents.

News snippet

Allow Maxwell propagation with external currents.

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.

Closes #206 (closed)

src/system/xc_oep_inc.F90

@@ -37,8 +37,11 @@ subroutine X(xc_oep_calc)(oep, xcs, apply_sic_pz, gr, hm, st, ex, ec, vxc)
 
   FLOAT :: eig
   integer :: is, ist, ixc, nspin_, isp, idm, jdm
+  integer, save :: itera = 0
   logical, save :: first = .true.
   R_TYPE, allocatable :: psi(:)
+  if ((.not.st%fromscratch).and.(itera==0)) &
+  itera = itera + 1
 
   if(oep%level == XC_OEP_NONE) return
 
@@ -125,9 +128,13 @@ subroutine X(xc_oep_calc)(oep, xcs, apply_sic_pz, gr, hm, st, ex, ec, vxc)
         if(oep%level == XC_OEP_FULL .and. (.not. first)) then
           call X(xc_oep_solve)(gr, hm, st, is, vxc(:,is), oep)
         end if
+        if (is == nspin_) &
+        itera = itera+1
+        if (is < nspin_) &
+        ec = 0
 
         first = .false.
-        vxc(1:gr%mesh%np, is) = vxc(1:gr%mesh%np, is) + oep%vxc(1:gr%mesh%np,1)
+        vxc(1:gr%mesh%np, is) = vxc(1:gr%mesh%np, is) + oep%vxc(1:gr%mesh%np,is)
       end if
     end do spin2
   end if
@@ -147,13 +154,15 @@ subroutine X(xc_oep_solve) (gr, hm, st, is, vxc, oep)
   type(hamiltonian_t), intent(in)    :: hm
   type(states_t),      intent(in)    :: st
   integer,             intent(in)    :: is
-  FLOAT,               intent(inout) :: vxc(:) !< (gr%mesh%np)
+  FLOAT,               intent(inout) :: vxc(:) !< (gr%mesh%np, given for the spin is)
   type(xc_oep_t),      intent(inout) :: oep
 
   integer :: iter, ist, iter_used
   FLOAT :: vxc_bar, ff, residue
   FLOAT, allocatable :: ss(:), vxc_old(:)
   R_TYPE, allocatable :: bb(:,:), psi(:, :)
+  R_TYPE, allocatable :: phi1(:,:,:)
+  logical, allocatable :: orthogonal(:)
 
   call profiling_in(C_PROFILING_XC_OEP_FULL, 'XC_OEP_FULL')
   PUSH_SUB(X(xc_oep_solve))
@@ -165,49 +174,110 @@ subroutine X(xc_oep_solve) (gr, hm, st, is, vxc, oep)
   SAFE_ALLOCATE(     ss(1:gr%mesh%np))
   SAFE_ALLOCATE(vxc_old(1:gr%mesh%np))
   SAFE_ALLOCATE(psi(1:gr%mesh%np, 1:st%d%dim))
+  SAFE_ALLOCATE(orthogonal(1:st%nst))
+
+  if (oep%has_photons) then
+    SAFE_ALLOCATE(phi1(1:gr%mesh%np,1:st%d%dim,1:st%nst))
+  end if
 
   vxc_old(1:gr%mesh%np) = vxc(1:gr%mesh%np)
 
   if(.not. lr_is_allocated(oep%lr)) then
     call lr_allocate(oep%lr, st, gr%mesh)
     ! initialize to something non-zero
-    oep%lr%X(dl_psi)(:,:, :, :) = M_ONE
+    oep%lr%X(dl_psi)(:,:, :, :) = M_ZERO
+  end if
+
+  if (oep%has_photons) then
+    if(.not. lr_is_allocated(oep%pt%lr)) then
+      call lr_allocate(oep%pt%lr, st, gr%mesh)
+      ! initialize to something non-zero
+      oep%pt%lr%X(dl_psi)(:,:, :, :) = M_ZERO
+    end if
+  !write(*,*) oep%has_photons, oep%mixing_density, oep%pt%omega, oep%pt%lambda
+    call X(xc_oep_pt_phi)(gr, hm, st, is, oep, phi1)
   end if
 
   ! fix xc potential (needed for Hpsi)
-  vxc(1:gr%mesh%np) = vxc_old(1:gr%mesh%np) + oep%vxc(1:gr%mesh%np,1)
+  vxc(1:gr%mesh%np) = vxc_old(1:gr%mesh%np) + oep%vxc(1:gr%mesh%np,is)
 
   do iter = 1, oep%scftol%max_iter
     ! iteration over all states
     ss = M_ZERO
-    do ist = 1, st%nst
+    oep%pt%ex = M_ZERO
+    do ist = 1, (oep%eigen_n + 1)
 
       call states_get_state(st, gr%mesh, ist, is, psi)
       
       ! evaluate right-hand side
-      vxc_bar = dmf_dotp(gr%mesh, (R_ABS(psi(:, 1)))**2, oep%vxc(1:gr%mesh%np, 1))
-      bb(1:gr%mesh%np, 1) = -(oep%vxc(1:gr%mesh%np, 1) - (vxc_bar - oep%uxc_bar(ist, is)))* &
+      vxc_bar = dmf_dotp(gr%mesh, (R_ABS(psi(:, 1)))**2, oep%vxc(1:gr%mesh%np, is))
+      bb(1:gr%mesh%np, 1) = -(oep%vxc(1:gr%mesh%np, is) - (vxc_bar - oep%uxc_bar(ist, is)))* &
         R_CONJ(psi(:, 1)) + oep%X(lxc)(1:gr%mesh%np, ist, is)
 
-      call X(lr_orth_vector) (gr%mesh, st, bb, ist, is, R_TOTYPE(M_ZERO))
+      if ((oep%eigen_n == 0) .AND. (st%occ(st%st_start, is) == 1)) then
+        if (min(st%d%nspin, 2) == 2) then
+          if ((st%occ(st%st_start, 1)+st%occ(st%st_start, 2)) == 1) &
+          write(*,*) ' OEP is not self-interaction corrected - please use KLI or more than one electron'
+        else 
+          write(*,*) ' OEP is not self-interaction corrected - please use KLI or more than one electron'
+        end if
+        ! single_electron_correction = 2.0
+        ! bb(1:gr%mesh%np, 1) = bb(1:gr%mesh%np, 1)*single_electron_correction
+      end if
+
+      if (oep%has_photons) &
+        call X(xc_oep_pt_rhs)(gr, st, is, oep, phi1, ist, bb)
+
+      if (.not.oep%has_photons) then
+        call X(lr_orth_vector) (gr%mesh, st, bb, ist, is, R_TOTYPE(M_ZERO))
+      else
+       orthogonal = .true.
+       orthogonal(ist) = .false.
+       call X(states_orthogonalize_single)(st, gr%mesh, st%nst, is, bb, normalize = .false., mask = orthogonal)
+      end if
 
       call X(linear_solver_solve_HXeY)(oep%solver, hm, gr, st, ist, is, oep%lr%X(dl_psi)(:,:, ist, is), bb, &
            R_TOTYPE(-st%eigenval(ist, is)), oep%scftol%final_tol, residue, iter_used)
-      
-      call X(lr_orth_vector) (gr%mesh, st, oep%lr%X(dl_psi)(:,:, ist, is), ist, is, R_TOTYPE(M_ZERO))
+
+      if (.not.oep%has_photons) then
+        call X(lr_orth_vector) (gr%mesh, st, oep%lr%X(dl_psi)(:,:, ist, is), ist, is, R_TOTYPE(M_ZERO))
+      else
+        orthogonal = .true.
+        orthogonal(ist) = .false.    
+        call X(states_orthogonalize_single)(st, gr%mesh, st%nst, is, oep%lr%X(dl_psi)(:,:, ist, is), normalize = .false., mask = orthogonal)
+      end if
 
       ! calculate this funny function ss
       ss(1:gr%mesh%np) = ss(1:gr%mesh%np) + M_TWO*R_REAL(oep%lr%X(dl_psi)(1:gr%mesh%np, 1, ist, is)*psi(:, 1))
+      if (oep%has_photons) &
+        call X(xc_oep_pt_inhomog)(gr, st, is, oep, phi1, ist, ss)
     end do
 
-    oep%vxc(1:gr%mesh%np,1) = oep%vxc(1:gr%mesh%np,1) + oep%mixing*ss(1:gr%mesh%np)
 
+    if ((oep%mixing_scheme == OEP_MIXING_SCHEME_CONST)) then
+      oep%vxc(1:gr%mesh%np,is) = oep%vxc(1:gr%mesh%np,is) + oep%mixing*ss(1:gr%mesh%np)
+    else if (oep%mixing_scheme == OEP_MIXING_SCHEME_DENS) then
+      oep%vxc(1:gr%mesh%np,is) = oep%vxc(1:gr%mesh%np,is) + oep%mixing*ss(1:gr%mesh%np)/st%rho(1:gr%mesh%np,is)
+    else if (oep%mixing_scheme == OEP_MIXING_SCHEME_BB) then
+      if (dmf_nrm2(gr%mesh, oep%vxc_old(1:gr%mesh%np,is)) > M_EPSILON ) then ! do not do it for the first run
+        oep%mixing = -dmf_dotp(gr%mesh, oep%vxc(1:gr%mesh%np,is) - oep%vxc_old(1:gr%mesh%np,is), ss - oep%ss_old) &
+          / dmf_dotp(gr%mesh, ss - oep%ss_old, ss - oep%ss_old)
+      end if
+      write(message(1), '(a,es14.6,a,es14.8)') "Info: oep%mixing:", oep%mixing, " norm2ss: ", dmf_nrm2(gr%mesh, ss)
+      call messages_info(1)
+
+      oep%vxc_old(1:gr%mesh%np,is) = oep%vxc(1:gr%mesh%np,is)
+      oep%ss_old(1:gr%mesh%np) = ss(1:gr%mesh%np)
+      oep%vxc(1:gr%mesh%np,is) = oep%vxc(1:gr%mesh%np,is) + oep%mixing*ss(1:gr%mesh%np)
+    end if
 
     do ist = 1, st%nst
       if(oep%eigen_type(ist) == 2) then
         call states_get_state(st, gr%mesh, ist, is, psi)
-        vxc_bar = dmf_dotp(gr%mesh, (R_ABS(psi(:, 1)))**2, oep%vxc(1:gr%mesh%np,1))
-        oep%vxc(1:gr%mesh%np,1) = oep%vxc(1:gr%mesh%np,1) - (vxc_bar - oep%uxc_bar(ist,is))
+        vxc_bar = dmf_dotp(gr%mesh, (R_ABS(psi(:, 1)))**2, oep%vxc(1:gr%mesh%np,is))
+        if (oep%has_photons) &
+          call X(xc_oep_pt_uxcbar)(gr, st, is, oep, phi1(:,:,ist), ist, vxc_bar)
+        oep%vxc(1:gr%mesh%np,is) = oep%vxc(1:gr%mesh%np,is) - (vxc_bar - oep%uxc_bar(ist,is))
       end if
     end do
 
@@ -215,6 +285,8 @@ subroutine X(xc_oep_solve) (gr, hm, st, is, vxc, oep)
     if(ff < oep%scftol%conv_abs_dens) exit
   end do
 
+  oep%norm2ss = ff
+
   if(ff > oep%scftol%conv_abs_dens) then
     write(message(1), '(a)') "OEP did not converge."
     call messages_warning(1)
@@ -232,6 +304,10 @@ subroutine X(xc_oep_solve) (gr, hm, st, is, vxc, oep)
   SAFE_DEALLOCATE_A(ss)
   SAFE_DEALLOCATE_A(vxc_old)
   SAFE_DEALLOCATE_A(psi)
+  SAFE_DEALLOCATE_A(orthogonal)
+  if (oep%has_photons) then
+    SAFE_DEALLOCATE_A(phi1)
+  end if
 
   POP_SUB(X(xc_oep_solve))
   call profiling_out(C_PROFILING_XC_OEP_FULL)