diff --git a/src/electrons/exponential.F90 b/src/electrons/exponential.F90 index 53ac0b7d46ad5c155835b1393598e9ac90b8a6a7..c872caadbe00c7b3d8522cb08007bd66254a704d 100644 --- a/src/electrons/exponential.F90 +++ b/src/electrons/exponential.F90 @@ -57,7 +57,9 @@ module exponential_oct_m exponential_init, & exponential_copy, & exponential_apply_all, & - exponential_lanczos_function_batch + exponential_lanczos_function_batch,& + operator_t, & + hamiltonian_operator_t integer, public, parameter :: & EXP_LANCZOS = 2, & @@ -78,6 +80,36 @@ module exponential_oct_m procedure :: apply_phi_batch => exponential_apply_phi_batch end type exponential_t + !< @brief Class to encapsulate the operation inside the exponential + type, abstract :: operator_t + type(namespace_t), pointer :: namespace + class(mesh_t), pointer :: mesh + contains + procedure(operator_apply), deferred :: apply + end type operator_t + + abstract interface + subroutine operator_apply(this, psib, hpsib) + import operator_t + import batch_t + class(operator_t), intent(in) :: this + class(batch_t), intent(inout) :: psib + class(batch_t), intent(inout) :: hpsib + end subroutine operator_apply + end interface + + !< @brief Class to encapsulate the Hamiltonian inside the exponential + type, extends(operator_t) :: hamiltonian_operator_t + class(hamiltonian_abst_t), pointer :: hm + contains + procedure :: apply => hamiltonian_operator_apply + procedure :: init => hamiltonian_operator_init + end type hamiltonian_operator_t + + interface hamiltonian_operator_t + procedure hamiltonian_operator_constructor + end interface hamiltonian_operator_t + contains ! --------------------------------------------------------- @@ -239,7 +271,7 @@ contains ! --------------------------------------------------------- !> Wrapper to batchified routine for applying exponential to an array - subroutine exponential_apply_single(te, namespace, mesh, hm, zpsi, ist, ik, deltat, vmagnus, imag_time) + subroutine exponential_apply_single(te, namespace, mesh, hm, zpsi, ist, ik, deltat, imag_time) class(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh @@ -248,7 +280,6 @@ contains integer, intent(in) :: ik complex(real64), contiguous, intent(inout) :: zpsi(:, :) real(real64), intent(in) :: deltat - real(real64), optional, intent(in) :: vmagnus(mesh%np, hm%d%nspin, 2) logical, optional, intent(in) :: imag_time type(wfs_elec_t) :: psib, inh_psib @@ -269,12 +300,12 @@ contains call states_elec_get_state(hm%inh_st, mesh, ist, ik, zpsi_inh(:, :)) call wfs_elec_init(inh_psib, hm%d%dim, ist, ist, zpsi_inh, ik) call te%apply_batch(namespace, mesh, hm, psib, deltat, & - vmagnus=vmagnus, imag_time=imag_time, inh_psib=inh_psib) + imag_time=imag_time, inh_psib=inh_psib) call inh_psib%end() SAFE_DEALLOCATE_A(zpsi_inh) else call te%apply_batch(namespace, mesh, hm, psib, deltat, & - vmagnus=vmagnus, imag_time=imag_time) + imag_time=imag_time) end if call psib%end() @@ -304,22 +335,23 @@ contains !! \phi_1(x) = (e^x - 1)/x !! \f] ! --------------------------------------------------------- - subroutine exponential_apply_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, imag_time, inh_psib) - class(exponential_t), intent(inout) :: te - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(hamiltonian_abst_t), intent(inout) :: hm - class(batch_t), intent(inout) :: psib - real(real64), intent(in) :: deltat - class(batch_t), optional, intent(inout) :: psib2 - real(real64), optional, intent(in) :: deltat2 - real(real64), optional, intent(in) :: vmagnus(:,:,:) !(mesh%np, hm%d%nspin, 2) - logical, optional, intent(in) :: imag_time - class(batch_t), optional, intent(inout) :: inh_psib !< inhomogeneous term + subroutine exponential_apply_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, imag_time, inh_psib, op) + class(exponential_t), intent(inout) :: te + type(namespace_t), intent(in) :: namespace + class(mesh_t), intent(in) :: mesh + class(hamiltonian_abst_t), intent(inout) :: hm + class(batch_t), intent(inout) :: psib + real(real64), intent(in) :: deltat + class(batch_t), optional, intent(inout) :: psib2 + real(real64), optional, intent(in) :: deltat2 + logical, optional, intent(in) :: imag_time + class(batch_t), optional, intent(inout) :: inh_psib !< inhomogeneous term + class(operator_t), target, optional, intent(in) :: op !< operator to be applied in the exponential complex(real64) :: deltat_, deltat2_ class(chebyshev_function_t), pointer :: chebyshev_function, chebyshev_function_dt2 logical :: imag_time_ + class(operator_t), pointer :: op_ PUSH_SUB(exponential_apply_batch) call profiling_in("EXPONENTIAL_BATCH") @@ -331,16 +363,10 @@ contains ASSERT(inh_psib%nst == psib%nst) end if - if (present(vmagnus)) then - select type(hm) - class is (hamiltonian_elec_t) - ASSERT(size(vmagnus, 1) >= mesh%np) - ASSERT(size(vmagnus, 2) == hm%d%nspin) - ASSERT(size(vmagnus, 3) == 2) - class default - write(message(1), '(a)') 'Magnus operators only implemented for electrons at the moment' - call messages_fatal(1, namespace=namespace) - end select + if (present(op)) then + op_ => op + else + op_ => hamiltonian_operator_t(namespace, mesh, hm) end if deltat2_ = cmplx(optional_default(deltat2, M_ZERO), M_ZERO, real64) @@ -365,32 +391,31 @@ contains case (EXP_TAYLOR) ! Note that delttat2_ is only initialized if deltat2 is present. if (present(deltat2)) then - call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, & - psib2, deltat2_, vmagnus) + call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, op_, & + psib2, deltat2_) else - call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, & - vmagnus=vmagnus) + call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, op_) end if if (present(inh_psib)) then if (present(deltat2)) then - call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, & - psib2, deltat2_, vmagnus, inh_psib) + call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, op_, & + psib2, deltat2_, inh_psib) else - call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, & - vmagnus=vmagnus, inh_psib=inh_psib) + call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, op_, & + inh_psib=inh_psib) end if end if case (EXP_LANCZOS) if (present(psib2)) call psib%copy_data_to(mesh%np, psib2) - call exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat_, vmagnus) + call exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat_, op_) if (present(inh_psib)) then - call exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat_, vmagnus, inh_psib) + call exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat_, op_, inh_psib) end if if (present(psib2)) then - call exponential_lanczos_batch(te, namespace, mesh, hm, psib2, deltat2_, vmagnus) + call exponential_lanczos_batch(te, namespace, mesh, hm, psib2, deltat2_, op_) if (present(inh_psib)) then - call exponential_lanczos_batch(te, namespace, mesh, hm, psib2, deltat2_, vmagnus, inh_psib) + call exponential_lanczos_batch(te, namespace, mesh, hm, psib2, deltat2_, op_, inh_psib) end if end if @@ -413,30 +438,34 @@ contains end if end if if (present(psib2)) call psib%copy_data_to(mesh%np, psib2) - call exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, vmagnus) + call exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, op_) deallocate(chebyshev_function) if (present(psib2)) then - call exponential_cheby_batch(te, namespace, mesh, hm, psib2, chebyshev_function_dt2, vmagnus) + call exponential_cheby_batch(te, namespace, mesh, hm, psib2, chebyshev_function_dt2, op_) deallocate(chebyshev_function_dt2) end if end select + if (.not.present(op)) then + SAFE_DEALLOCATE_P(op_) + end if + call profiling_out("EXPONENTIAL_BATCH") POP_SUB(exponential_apply_batch) end subroutine exponential_apply_batch ! --------------------------------------------------------- - subroutine exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, inh_psib, phik_shift) + subroutine exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat, op, psib2, deltat2, inh_psib, phik_shift) type(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh class(hamiltonian_abst_t), intent(inout) :: hm class(batch_t), intent(inout) :: psib complex(real64), intent(in) :: deltat + class(operator_t), intent(in) :: op !< operator to be applied in the exponential class(batch_t), optional, intent(inout) :: psib2 complex(real64), optional, intent(in) :: deltat2 - real(real64), optional, intent(in) :: vmagnus(:,:,:) !(mesh%np, hm%d%nspin, 2) class(batch_t), optional, intent(inout) :: inh_psib !< inhomogeneous term integer, optional, intent(in) :: phik_shift !< shift in the Taylor expansion coefficients for phi_k @@ -482,12 +511,12 @@ contains ! the code stops in ZAXPY below without saying why. if (iter /= 1) then - call operate_batch(hm, namespace, mesh, psi1b, hpsi1b, vmagnus) + call op%apply(psi1b, hpsi1b) else if (present(inh_psib)) then - call operate_batch(hm, namespace, mesh, inh_psib, hpsi1b, vmagnus) + call op%apply(inh_psib, hpsi1b) else - call operate_batch(hm, namespace, mesh, psib, hpsi1b, vmagnus) + call op%apply(psib, hpsi1b) end if end if @@ -512,15 +541,15 @@ contains POP_SUB(exponential_taylor_series_batch) end subroutine exponential_taylor_series_batch - subroutine exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat, vmagnus, inh_psib) + subroutine exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat, op, inh_psib) type(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh class(hamiltonian_abst_t), intent(inout) :: hm class(batch_t), intent(inout) :: psib complex(real64), intent(in) :: deltat - real(real64), optional, intent(in) :: vmagnus(:,:,:) !(mesh%np, hm%d%nspin, 2) - class(batch_t), optional, intent(in) :: inh_psib !< inhomogeneous term + class(operator_t), intent(in) :: op !< operator to be applied in the exponential + class(batch_t), optional, intent(in) :: inh_psib !< inhomogeneous term class(batch_t), allocatable :: tmpb @@ -529,11 +558,11 @@ contains if (present(inh_psib)) then call inh_psib%clone_to(tmpb, copy_data=.true.) ! psib = psib + deltat * phi1(-i*deltat*H) inh_psib - call exponential_lanczos_function_batch(te, namespace, mesh, hm, tmpb, deltat, phi1, vmagnus) + call exponential_lanczos_function_batch(te, namespace, mesh, hm, tmpb, deltat, phi1, op) call batch_axpy(mesh%np, real(deltat, real64), tmpb, psib) call tmpb%end() else - call exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat, exponential, vmagnus) + call exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat, exponential, op) end if POP_SUB(exponential_lanczos_batch) @@ -551,7 +580,7 @@ contains !! A pdf can be accessed [here](https://www-users.cse.umn.edu/~saad/PDF/RIACS-90-ExpTh.pdf) !! or [via the DOI](https://doi.org/10.1137/0729014). !! Equation numbers below refer to this paper. - subroutine exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat, fun, vmagnus) + subroutine exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat, fun, op) type(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh @@ -564,7 +593,7 @@ contains complex(real64), intent(in) :: z end end interface - real(real64), optional, intent(in) :: vmagnus(:,:,:) !(mesh%np, hm%d%nspin, 2) + class(operator_t), intent(in) :: op !< operator to be applied in the exponential integer :: iter, l, ii, ist, max_initialized complex(real64), allocatable :: hamilt(:,:,:), expo(:,:,:) @@ -614,8 +643,8 @@ contains call psib%clone_to(vb(iter + 1)%p) max_initialized = iter + 1 - ! to apply the Hamiltonian - call operate_batch(hm, namespace, mesh, vb(iter)%p, vb(iter+1)%p, vmagnus) + ! to apply the operator (default is Hamiltonian) + call op%apply(vb(iter)%p, vb(iter+1)%p) ! We use either the Lanczos method (Hermitian case) or the Arnoldi method if (hm%is_hermitian()) then @@ -712,14 +741,14 @@ contains !! Section III.1 and !! [Kosloff, Annu. Rev. Phys. Chem. (1994), 45, 145–178](https://doi.org/10.1146/annurev.pc.45.100194.001045), !! equations 7.1ff and especially figure 2 for the convergence properties of the coefficients. - subroutine exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, vmagnus) + subroutine exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, op) type(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh class(hamiltonian_abst_t), intent(inout) :: hm class(batch_t), intent(inout) :: psib class(chebyshev_function_t), intent(in) :: chebyshev_function - real(real64), optional, intent(in) :: vmagnus(:,:,:) !(mesh%np, hm%d%nspin, 2) + class(operator_t), intent(in) :: op !< operator to be applied in the exponential integer :: j, order_needed complex(real64) :: coefficient @@ -752,7 +781,7 @@ contains call batch_scal(mesh%np, coefficients(0), psib) ! first-order term ! shifted Hamiltonian - call operate_batch(hm, namespace, mesh, psi0, psi1, vmagnus) + call op%apply(psi0, psi1) call batch_axpy(mesh%np, -hm%spectral_middle_point, psi0, psi1) call batch_scal(mesh%np, M_ONE/hm%spectral_half_span, psi1) ! accumulate result @@ -765,7 +794,7 @@ contains do j = 2, order_needed ! compute shifted Hamiltonian and Chebyshev recurrence formula - call operate_batch(hm, namespace, mesh, psi_n1, psi_n, vmagnus) + call op%apply(psi_n1, psi_n) call batch_axpy(mesh%np, -hm%spectral_middle_point, psi_n1, psi_n) call batch_xpay(mesh%np, psi_n2, -M_TWO/hm%spectral_half_span, psi_n) call batch_scal(mesh%np, -M_ONE, psi_n) @@ -812,26 +841,6 @@ contains POP_SUB(exponential_cheby_batch) end subroutine exponential_cheby_batch - - subroutine operate_batch(hm, namespace, mesh, psib, hpsib, vmagnus) - class(hamiltonian_abst_t), intent(inout) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(batch_t), intent(inout) :: psib - class(batch_t), intent(inout) :: hpsib - real(real64), optional, intent(in) :: vmagnus(:, :, :) - - PUSH_SUB(operate_batch) - - if (present(vmagnus)) then - call hm%zmagnus_apply(namespace, mesh, psib, hpsib, vmagnus) - else - call hm%zapply(namespace, mesh, psib, hpsib) - end if - - POP_SUB(operate_batch) - end subroutine operate_batch - ! --------------------------------------------------------- !> Note that this routine not only computes the exponential, but !! also an extra term if there is a inhomogeneous term in the @@ -925,7 +934,7 @@ contains POP_SUB(exponential_apply_all) end subroutine exponential_apply_all - subroutine exponential_apply_phi_batch(te, namespace, mesh, hm, psib, deltat, k, vmagnus) + subroutine exponential_apply_phi_batch(te, namespace, mesh, hm, psib, deltat, k, op) class(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh @@ -933,25 +942,15 @@ contains class(batch_t), intent(inout) :: psib real(real64), intent(in) :: deltat integer, intent(in) :: k - real(real64), optional, intent(in) :: vmagnus(:,:,:) !(mesh%np, hm%d%nspin, 2) + class(operator_t), target,optional, intent(in) :: op class(chebyshev_function_t), pointer :: chebyshev_function complex(real64) :: deltat_ + class(operator_t), pointer :: op_ PUSH_SUB_WITH_PROFILE(exponential_apply_phi_batch) ASSERT(psib%type() == TYPE_CMPLX) - if (present(vmagnus)) then - select type(hm) - class is (hamiltonian_elec_t) - ASSERT(size(vmagnus, 1) >= mesh%np) - ASSERT(size(vmagnus, 2) == hm%d%nspin) - ASSERT(size(vmagnus, 3) == 2) - class default - write(message(1), '(a)') 'Magnus operators only implemented for electrons at the moment' - call messages_fatal(1, namespace=namespace) - end select - end if if (.not. hm%is_hermitian() .and. te%exp_method == EXP_CHEBYSHEV) then write(message(1), '(a)') 'The Chebyshev expansion for the exponential will only converge if the imaginary' @@ -962,17 +961,23 @@ contains call messages_warning(5, namespace=namespace) end if + if (present(op)) then + op_ => op + else + op_ => hamiltonian_operator_t(namespace, mesh, hm) + end if + deltat_ = cmplx(deltat, M_ZERO, real64) select case (te%exp_method) case (EXP_TAYLOR) - call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, vmagnus=vmagnus, phik_shift=k) + call exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat_, op_, phik_shift=k) case (EXP_LANCZOS) if (k == 1) then - call exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat_, phi1, vmagnus) + call exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat_, phi1, op_) else if (k == 2) then - call exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat_, phi2, vmagnus) + call exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat_, phi2, op_) else write(message(1), '(a)') 'Lanczos expansion not implemented for phi_k, k > 2' call messages_fatal(1, namespace=namespace) @@ -987,13 +992,58 @@ contains write(message(1), '(a)') 'Chebyshev expansion not implemented for phi_k, k > 2' call messages_fatal(1, namespace=namespace) end if - call exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, vmagnus) + call exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, op_) deallocate(chebyshev_function) - end select + + if (.not.present(op)) then + SAFE_DEALLOCATE_P(op_) + end if + POP_SUB_WITH_PROFILE(exponential_apply_phi_batch) end subroutine exponential_apply_phi_batch + function hamiltonian_operator_constructor(namespace, mesh, hm) result(this) + type(namespace_t), target, intent(in) :: namespace + class(mesh_t), target, intent(in) :: mesh + class(hamiltonian_abst_t), target, intent(in) :: hm + type(hamiltonian_operator_t), pointer :: this + + PUSH_SUB(hamiltonian_operator_constructor) + + allocate(this) + call this%init(namespace, mesh, hm) + + POP_SUB(hamiltonian_operator_constructor) + end function hamiltonian_operator_constructor + + subroutine hamiltonian_operator_init(this, namespace, mesh, hm) + class(hamiltonian_operator_t), intent(inout) :: this + type(namespace_t), target, intent(in) :: namespace + class(mesh_t), target, intent(in) :: mesh + class(hamiltonian_abst_t), target, intent(in) :: hm + + PUSH_SUB(hamiltonian_operator_init) + + this%mesh => mesh + this%namespace => namespace + this%hm => hm + + POP_SUB(hamiltonian_operator_init) + end subroutine hamiltonian_operator_init + + subroutine hamiltonian_operator_apply(this, psib, hpsib) + class(hamiltonian_operator_t), intent(in) :: this + class(batch_t), intent(inout) :: psib + class(batch_t), intent(inout) :: hpsib + + PUSH_SUB(hamiltonian_operator_apply) + + call this%hm%zapply(this%namespace, this%mesh, psib, hpsib) + + POP_SUB(hamiltonian_operator_apply) + end subroutine hamiltonian_operator_apply + end module exponential_oct_m !! Local Variables: diff --git a/src/hamiltonian/hamiltonian_abst.F90 b/src/hamiltonian/hamiltonian_abst.F90 index 7ae72767d51559f80eb1829d22a0c4aa2a9849fc..33d9c5b87b694e18d96180cb565c9766d7fc8295 100644 --- a/src/hamiltonian/hamiltonian_abst.F90 +++ b/src/hamiltonian/hamiltonian_abst.F90 @@ -45,8 +45,6 @@ module hamiltonian_abst_oct_m procedure(hamiltonian_update_span), deferred :: update_span !< @copydoc hamiltonian_update_span procedure(dhamiltonian_apply), deferred :: dapply !< @copydoc dhamiltonian_apply procedure(zhamiltonian_apply), deferred :: zapply !< @copydoc zhamiltonian_apply - procedure(dhamiltonian_magnus_apply), deferred :: dmagnus_apply !< @copydoc dhamiltonian_magnus_apply - procedure(zhamiltonian_magnus_apply), deferred :: zmagnus_apply !< @copydoc zhamiltonian_magnus_apply end type hamiltonian_abst_t abstract interface @@ -84,26 +82,6 @@ module hamiltonian_abst_oct_m integer, optional, intent(in) :: terms logical, optional, intent(in) :: set_bc end subroutine zhamiltonian_apply - - subroutine dhamiltonian_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus) - import - class(hamiltonian_abst_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(batch_t), intent(inout) :: psib - class(batch_t), intent(inout) :: hpsib - real(real64), intent(in) :: vmagnus(:, :, :) - end subroutine dhamiltonian_magnus_apply - - subroutine zhamiltonian_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus) - import - class(hamiltonian_abst_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(batch_t), intent(inout) :: psib - class(batch_t), intent(inout) :: hpsib - real(real64), intent(in) :: vmagnus(:, :, :) - end subroutine zhamiltonian_magnus_apply end interface end module hamiltonian_abst_oct_m diff --git a/src/hamiltonian/hamiltonian_elec.F90 b/src/hamiltonian/hamiltonian_elec.F90 index b9f812defdb1a2df154c3856ff4804c05b0864e9..b24de1939076cd0257b853ac5387a46b0cd71d62 100644 --- a/src/hamiltonian/hamiltonian_elec.F90 +++ b/src/hamiltonian/hamiltonian_elec.F90 @@ -110,7 +110,6 @@ module hamiltonian_elec_oct_m dhamiltonian_elec_apply_batch, & zhamiltonian_elec_apply_batch, & hamiltonian_elec_diagonal, & - magnus, & dvmask, & zvmask, & hamiltonian_elec_inh_term, & @@ -125,6 +124,7 @@ module hamiltonian_elec_oct_m hamiltonian_elec_get_time, & hamiltonian_elec_apply_packed, & zhamiltonian_elec_apply_atom, & + zhamiltonian_elec_external, & hamiltonian_elec_has_kick, & hamiltonian_elec_copy_and_set_phase @@ -215,8 +215,6 @@ module hamiltonian_elec_oct_m procedure :: update_span => hamiltonian_elec_span procedure :: dapply => dhamiltonian_elec_apply procedure :: zapply => zhamiltonian_elec_apply - procedure :: dmagnus_apply => dhamiltonian_elec_magnus_apply - procedure :: zmagnus_apply => zhamiltonian_elec_magnus_apply procedure :: is_hermitian => hamiltonian_elec_hermitian procedure :: needs_mgga_term => hamiltonian_elec_needs_mgga_term procedure :: set_mass => hamiltonian_elec_set_mass @@ -1573,80 +1571,6 @@ contains POP_SUB(zhamiltonian_elec_apply_all) end subroutine zhamiltonian_elec_apply_all - - ! --------------------------------------------------------- - - subroutine magnus(hm, namespace, mesh, psi, hpsi, ik, vmagnus, set_phase) - type(hamiltonian_elec_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - complex(real64), contiguous, intent(inout) :: psi(:,:) - complex(real64), contiguous, intent(out) :: hpsi(:,:) - integer, intent(in) :: ik - real(real64), intent(in) :: vmagnus(:, :, :) - logical, optional, intent(in) :: set_phase !< If set to .false. the phase will not be added to the states. - - complex(real64), allocatable :: auxpsi(:, :), aux2psi(:, :) - integer :: idim, ispin - - PUSH_SUB(magnus) - - ! We will assume, for the moment, no spinors. - if (hm%d%dim /= 1) then - call messages_not_implemented("Magnus with spinors", namespace=namespace) - end if - - SAFE_ALLOCATE( auxpsi(1:mesh%np_part, 1:hm%d%dim)) - SAFE_ALLOCATE(aux2psi(1:mesh%np, 1:hm%d%dim)) - - ispin = hm%d%get_spin_index(ik) - - ! Compute (T + Vnl)|psi> and store it - call zhamiltonian_elec_apply_single(hm, namespace, mesh, psi, auxpsi, 1, ik, & - terms = TERM_KINETIC + TERM_NON_LOCAL_POTENTIAL, set_phase = set_phase) - - ! H|psi> = (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders - do idim = 1, hm%d%dim - call lalg_copy(mesh%np, auxpsi(:, idim), hpsi(:, idim)) - hpsi(1:mesh%np, idim) = hpsi(1:mesh%np, idim) + hm%ep%Vpsl(1:mesh%np)*psi(1:mesh%np,idim) - call vborders(mesh, hm, psi(:, idim), hpsi(:, idim)) - end do - hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + vmagnus(1:mesh%np, ispin, 2)*psi(1:mesh%np, 1) - - ! Add first term of the commutator: - i Vmagnus(t1) (T + Vnl) |psi> - hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) - M_zI*vmagnus(1:mesh%np, ispin, 1)*auxpsi(1:mesh%np, 1) - - ! Add second term of commutator: i (T + Vnl) Vmagnus(t1) |psi> - auxpsi(1:mesh%np, 1) = vmagnus(1:mesh%np, ispin, 1)*psi(1:mesh%np, 1) - call zhamiltonian_elec_apply_single(hm, namespace, mesh, auxpsi, aux2psi, 1, ik, & - terms = TERM_KINETIC + TERM_NON_LOCAL_POTENTIAL, set_phase = set_phase) - hpsi(1:mesh%np, 1) = hpsi(1:mesh%np, 1) + M_zI*aux2psi(1:mesh%np, 1) - - SAFE_DEALLOCATE_A(auxpsi) - SAFE_DEALLOCATE_A(aux2psi) - POP_SUB(magnus) - end subroutine magnus - - ! --------------------------------------------------------- - subroutine vborders (mesh, hm, psi, hpsi) - class(mesh_t), intent(in) :: mesh - type(hamiltonian_elec_t), intent(in) :: hm - complex(real64), intent(in) :: psi(:) - complex(real64), intent(inout) :: hpsi(:) - - integer :: ip - - PUSH_SUB(vborders) - - if (hm%abs_boundaries%abtype == IMAGINARY_ABSORBING) then - do ip = 1, mesh%np - hpsi(ip) = hpsi(ip) + M_zI*hm%abs_boundaries%mf(ip)*psi(ip) - end do - end if - - POP_SUB(vborders) - end subroutine vborders - ! --------------------------------------------------------- logical function hamiltonian_elec_has_kick(hm) type(hamiltonian_elec_t), intent(in) :: hm diff --git a/src/hamiltonian/hamiltonian_elec_inc.F90 b/src/hamiltonian/hamiltonian_elec_inc.F90 index 4cad89e373c59f94e78c329e1b0bbfcfee610cb8..1d595d8a0270f89fc39cd1fb7185936900236d71 100644 --- a/src/hamiltonian/hamiltonian_elec_inc.F90 +++ b/src/hamiltonian/hamiltonian_elec_inc.F90 @@ -45,34 +45,6 @@ subroutine X(hamiltonian_elec_apply) (hm, namespace, mesh, psib, hpsib, terms, s POP_SUB(X(hamiltonian_elec_apply)) end subroutine X(hamiltonian_elec_apply) -! --------------------------------------------------------- -subroutine X(hamiltonian_elec_magnus_apply) (hm, namespace, mesh, psib, hpsib, vmagnus) - class(hamiltonian_elec_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(batch_t), intent(inout) :: psib - class(batch_t), intent(inout) :: hpsib - real(real64), intent(in) :: vmagnus(:, :, :) - - PUSH_SUB(X(hamiltonian_elec_magnus_apply)) - - select type (psib) - class is (wfs_elec_t) - select type (hpsib) - class is (wfs_elec_t) - call X(hamiltonian_elec_magnus_apply_batch) (hm, namespace, mesh, psib, hpsib, vmagnus) - class default - message(1) = "Internal error: imcompatible batch_t in hamiltonian_elec_magnus_apply for argument hpsib." - call messages_fatal(1, namespace=namespace) - end select - class default - message(1) = "Internal error: imcompatible batch_t in hamiltonian_elec_magnus_apply for argument psib." - call messages_fatal(1, namespace=namespace) - end select - - POP_SUB(X(hamiltonian_elec_magnus_apply)) -end subroutine X(hamiltonian_elec_magnus_apply) - ! --------------------------------------------------------- subroutine X(hamiltonian_elec_apply_batch) (hm, namespace, mesh, psib, hpsib, terms, set_bc) type(hamiltonian_elec_t), intent(in) :: hm @@ -371,58 +343,6 @@ subroutine X(hamiltonian_elec_apply_single) (hm, namespace, mesh, psi, hpsi, ist POP_SUB(X(hamiltonian_elec_apply_single)) end subroutine X(hamiltonian_elec_apply_single) - -subroutine X(hamiltonian_elec_magnus_apply_batch) (hm, namespace, mesh, psib, hpsib, vmagnus) - type(hamiltonian_elec_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - type(wfs_elec_t), intent(inout) :: psib - type(wfs_elec_t), intent(inout) :: hpsib - real(real64), intent(in) :: vmagnus(:, :, :) - - integer :: ispin - type(wfs_elec_t) :: auxpsib, aux2psib - - PUSH_SUB(X(hamiltonian_elec_magnus_apply_batch)) - - ! We will assume, for the moment, no spinors. - if (hm%d%dim /= 1) call messages_not_implemented("Magnus with spinors", namespace=namespace) - - ASSERT(psib%nst == hpsib%nst) - - ispin = hm%d%get_spin_index(psib%ik) - - call hpsib%copy_to(auxpsib, copy_data=.false.) - call hpsib%copy_to(aux2psib, copy_data=.false.) - - ! Compute (T + Vnl)|psi> and store it - call X(hamiltonian_elec_apply_batch)(hm, namespace, mesh, psib, auxpsib, terms = TERM_KINETIC + TERM_NON_LOCAL_POTENTIAL) - - ! H|psi> = (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders|psi> - call auxpsib%copy_data_to(mesh%np, hpsib) - call X(hamiltonian_elec_external)(hm, mesh, psib, hpsib) - if (hm%abs_boundaries%abtype == IMAGINARY_ABSORBING) then - call batch_mul(mesh%np, hm%abs_boundaries%mf(1:mesh%np), psib, aux2psib) - call batch_axpy(mesh%np, M_zI, aux2psib, hpsib) - end if - call batch_mul(mesh%np, vmagnus(1:mesh%np, ispin, 2), psib, aux2psib) - call batch_axpy(mesh%np, M_ONE, aux2psib, hpsib) - - ! Add first term of the commutator: - i Vmagnus(t1) (T + Vnl) |psi> - call batch_mul(mesh%np, vmagnus(1:mesh%np, ispin, 1), auxpsib, aux2psib) - call batch_axpy(mesh%np, -M_zI, aux2psib, hpsib) - - ! Add second term of commutator: i (T + Vnl) Vmagnus(t1) |psi> - call batch_mul(mesh%np, vmagnus(1:mesh%np, ispin, 1), psib, auxpsib) - call X(hamiltonian_elec_apply_batch)(hm, namespace, mesh, auxpsib, aux2psib, terms = TERM_KINETIC + TERM_NON_LOCAL_POTENTIAL) - call batch_axpy(mesh%np, M_zI, aux2psib, hpsib) - - call auxpsib%end() - call aux2psib%end() - - POP_SUB(X(hamiltonian_elec_magnus_apply_batch)) -end subroutine X(hamiltonian_elec_magnus_apply_batch) - ! --------------------------------------------------------- !> Computes the orbital-dependent part of the mGGA potential !! diff --git a/src/maxwell/hamiltonian_mxll.F90 b/src/maxwell/hamiltonian_mxll.F90 index f4a38251a4f3d9985c06b26a1e50e0311dcd45d7..592ae6665427cd7143fa427b1e7354d36b79116a 100644 --- a/src/maxwell/hamiltonian_mxll.F90 +++ b/src/maxwell/hamiltonian_mxll.F90 @@ -56,8 +56,6 @@ module hamiltonian_mxll_oct_m hamiltonian_mxll_end, & dhamiltonian_mxll_apply, & zhamiltonian_mxll_apply, & - dhamiltonian_mxll_magnus_apply, & - zhamiltonian_mxll_magnus_apply, & hamiltonian_mxll_apply_batch, & hamiltonian_mxll_span, & hamiltonian_mxll_adjoint, & @@ -142,8 +140,6 @@ module hamiltonian_mxll_oct_m procedure :: update_span => hamiltonian_mxll_span procedure :: dapply => dhamiltonian_mxll_apply procedure :: zapply => zhamiltonian_mxll_apply - procedure :: dmagnus_apply => dhamiltonian_mxll_magnus_apply - procedure :: zmagnus_apply => zhamiltonian_mxll_magnus_apply procedure :: is_hermitian => hamiltonian_mxll_hermitian end type hamiltonian_mxll_t @@ -1418,34 +1414,6 @@ contains POP_SUB(maxwell_medium_boxes_calculation) end subroutine maxwell_medium_boxes_calculation - ! --------------------------------------------------------- - !> Maxwell hamiltonian Magnus (not implemented) - subroutine dhamiltonian_mxll_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus) - class(hamiltonian_mxll_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(batch_t), intent(inout) :: psib - class(batch_t), intent(inout) :: hpsib - real(real64), intent(in) :: vmagnus(:, :, :) - - call messages_not_implemented ('dhamiltonian_mxll_magnus_apply', namespace=namespace) - - end subroutine dhamiltonian_mxll_magnus_apply - - ! --------------------------------------------------------- - !> Maxwell hamiltonian Magnus (not implemented) - subroutine zhamiltonian_mxll_magnus_apply(hm, namespace, mesh, psib, hpsib, vmagnus) - class(hamiltonian_mxll_t), intent(in) :: hm - type(namespace_t), intent(in) :: namespace - class(mesh_t), intent(in) :: mesh - class(batch_t), intent(inout) :: psib - class(batch_t), intent(inout) :: hpsib - real(real64), intent(in) :: vmagnus(:, :, :) - - call messages_not_implemented ('zhamiltonian_mxll_magnus_apply', namespace=namespace) - - end subroutine zhamiltonian_mxll_magnus_apply - end module hamiltonian_mxll_oct_m !! Local Variables: diff --git a/src/td/propagation_ops_elec.F90 b/src/td/propagation_ops_elec.F90 index 1babb5cc78660f97d0eb79d6592ba3828944fc1a..691c573b5264b4f58da07127b92df3e9f90ddc2e 100644 --- a/src/td/propagation_ops_elec.F90 +++ b/src/td/propagation_ops_elec.F90 @@ -251,13 +251,14 @@ contains end subroutine propagation_ops_elec_restore_gauge_field ! --------------------------------------------------------- - subroutine propagation_ops_elec_exp_apply(te, namespace, st, mesh, hm, dt) - type(exponential_t), intent(inout) :: te - type(namespace_t), intent(in) :: namespace - type(states_elec_t), intent(inout) :: st - class(mesh_t), intent(in) :: mesh - type(hamiltonian_elec_t), intent(inout) :: hm - real(real64), intent(in) :: dt + subroutine propagation_ops_elec_exp_apply(te, namespace, st, mesh, hm, dt, op) + type(exponential_t), intent(inout) :: te + type(namespace_t), intent(in) :: namespace + type(states_elec_t), intent(inout) :: st + class(mesh_t), intent(in) :: mesh + type(hamiltonian_elec_t), intent(inout) :: hm + real(real64), intent(in) :: dt + class(operator_t), optional, intent(in) :: op !< operator to be applied in the exponential integer :: ik, ib @@ -274,9 +275,9 @@ contains call hm%phase%set_phase_corr(mesh, st%group%psib(ib, ik), async=.true.) if (hamiltonian_elec_inh_term(hm)) then call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt, & - inh_psib = hm%inh_st%group%psib(ib, ik)) + inh_psib = hm%inh_st%group%psib(ib, ik), op=op) else - call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt) + call te%apply_batch(namespace, mesh, hm, st%group%psib(ib, ik), dt, op=op) end if call hm%phase%unset_phase_corr(mesh, st%group%psib(ib, ik)) @@ -295,15 +296,15 @@ contains end subroutine propagation_ops_elec_exp_apply ! --------------------------------------------------------- - subroutine propagation_ops_elec_fuse_density_exp_apply(te, namespace, st, gr, hm, dt, dt2, vmagnus) - type(exponential_t), intent(inout) :: te - type(namespace_t), intent(in) :: namespace - type(states_elec_t), intent(inout) :: st - type(grid_t), intent(in) :: gr - type(hamiltonian_elec_t), intent(inout) :: hm - real(real64), intent(in) :: dt - real(real64), optional, intent(in) :: dt2 - real(real64), optional, intent(in) :: vmagnus(:,:,:) + subroutine propagation_ops_elec_fuse_density_exp_apply(te, namespace, st, gr, hm, dt, dt2, op) + type(exponential_t), intent(inout) :: te + type(namespace_t), intent(in) :: namespace + type(states_elec_t), intent(inout) :: st + type(grid_t), intent(in) :: gr + type(hamiltonian_elec_t), intent(inout) :: hm + real(real64), intent(in) :: dt + real(real64), optional, intent(in) :: dt2 + class(operator_t), optional, intent(in) :: op !< operator to be applied in the exponential integer :: ik, ib type(wfs_elec_t) :: zpsib_dt @@ -329,10 +330,10 @@ contains !propagate the state to dt/2 and dt, simultaneously, with H(time - dt) if (hamiltonian_elec_inh_term(hm)) then call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, psib2 = zpsib_dt, & - deltat2 = dt2, inh_psib = hm%inh_st%group%psib(ib, ik)) + deltat2 = dt2, inh_psib = hm%inh_st%group%psib(ib, ik), op = op) else call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, psib2 = zpsib_dt, & - deltat2 = dt2) + deltat2 = dt2, op = op) end if call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik), async=.true.) @@ -343,10 +344,10 @@ contains else !propagate the state to dt with H(time - dt) if (hamiltonian_elec_inh_term(hm)) then - call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, vmagnus=vmagnus, & - inh_psib = hm%inh_st%group%psib(ib, ik)) + call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, & + inh_psib = hm%inh_st%group%psib(ib, ik), op = op) else - call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, vmagnus=vmagnus) + call te%apply_batch(namespace, gr, hm, st%group%psib(ib, ik), dt, op = op) end if call hm%phase%unset_phase_corr(gr, st%group%psib(ib, ik), async=.true.) diff --git a/src/td/propagator_base.F90 b/src/td/propagator_base.F90 index bf0ad77a9ffba544da0b2e72bfefcf6038345b3a..0a537c1324a51476b1beecef8963672aae210421 100644 --- a/src/td/propagator_base.F90 +++ b/src/td/propagator_base.F90 @@ -53,7 +53,6 @@ module propagator_base_oct_m !> Storage of the KS potential of previous iterations. type(potential_interpolation_t) :: vks_old !> Auxiliary function to store the Magnus potentials. - real(real64), allocatable :: vmagnus(:, :, :) integer :: scf_propagation_steps type(sparskit_solver_t) :: tdsk integer :: tdsk_size diff --git a/src/td/propagator_elec.F90 b/src/td/propagator_elec.F90 index a53093d43c5edd7a8fcd50809235a055987be9ff..9371eb571d2026485f7eb86eeb269a3bed22fdf6 100644 --- a/src/td/propagator_elec.F90 +++ b/src/td/propagator_elec.F90 @@ -89,9 +89,6 @@ contains tro%method = tri%method select case (tro%method) - case (PROP_MAGNUS) - SAFE_ALLOCATE_SOURCE_A(tro%vmagnus, tri%vmagnus) - case (PROP_CRANK_NICOLSON_SPARSKIT) tro%tdsk_size = tri%tdsk_size call sparskit_solver_copy(tro%tdsk, tri%tdsk) @@ -321,7 +318,6 @@ contains message(1) = "Magnus propagator with MGGA" call messages_fatal(1, namespace=namespace) end if - SAFE_ALLOCATE(tr%vmagnus(1:gr%np, 1:st%d%nspin, 1:2)) case (PROP_QOCT_TDDFT_PROPAGATOR) call messages_experimental("QOCT+TDDFT propagator", namespace=namespace) case (PROP_EXPLICIT_RUNGE_KUTTA4) @@ -466,10 +462,6 @@ contains call potential_interpolation_end(tr%vks_old) select case (tr%method) - case (PROP_MAGNUS) - ASSERT(allocated(tr%vmagnus)) - SAFE_DEALLOCATE_A(tr%vmagnus) - case (PROP_RUNGE_KUTTA4, PROP_RUNGE_KUTTA2) call propagator_rk_end() call sparskit_solver_end(tr%tdsk) diff --git a/src/td/propagator_magnus.F90 b/src/td/propagator_magnus.F90 index 150530bd2a4f00adfc21f54583f26cbde5c2cb76..7d2d124713a573791d6199e2f305b38e0107c560 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -19,7 +19,9 @@ #include "global.h" module propagator_magnus_oct_m + use absorbing_boundaries_oct_m use batch_oct_m + use batch_ops_oct_m use debug_oct_m use density_oct_m use electron_space_oct_m @@ -28,7 +30,9 @@ module propagator_magnus_oct_m use gauge_field_oct_m use global_oct_m use grid_oct_m + use hamiltonian_abst_oct_m use hamiltonian_elec_oct_m + use hamiltonian_elec_base_oct_m use interaction_partner_oct_m use ion_dynamics_oct_m use ions_oct_m @@ -36,6 +40,7 @@ module propagator_magnus_oct_m use ks_potential_oct_m use lasers_oct_m use messages_oct_m + use mesh_oct_m use namespace_oct_m use parser_oct_m use potential_interpolation_oct_m @@ -46,6 +51,7 @@ module propagator_magnus_oct_m use space_oct_m use states_elec_oct_m use v_ks_oct_m + use wfs_elec_oct_m use xc_oct_m implicit none @@ -56,6 +62,19 @@ module propagator_magnus_oct_m td_magnus, & td_cfmagnus4 + !< @brief Class to encapsulate the operation for the Magnus 4 propagator + type, extends(hamiltonian_operator_t) :: magnus4_operator_t + private + real(real64), allocatable :: vmagnus(:, :, :) + contains + procedure :: apply => magnus4_operator_apply + end type magnus4_operator_t + + interface magnus4_operator_t + procedure magnus4_operator_constructor + end interface magnus4_operator_t + + contains ! --------------------------------------------------------- @@ -72,14 +91,16 @@ contains integer :: j, is, i real(real64) :: atime(2) - real(real64), allocatable :: vaux(:, :, :), pot(:) + real(real64), allocatable :: vaux(:, :, :), vmagnus(:, :, :), pot(:) type(lasers_t), pointer :: lasers + class(operator_t), pointer :: op PUSH_SUB(propagator_dt.td_magnus) ASSERT(.not. family_is_mgga_with_exc(hm%xc)) SAFE_ALLOCATE(vaux(1:gr%np, 1:st%d%nspin, 1:2)) + SAFE_ALLOCATE(vmagnus(1:gr%np, 1:st%d%nspin, 1:2)) atime(1) = (M_HALF-sqrt(M_THREE)/6.0_real64)*dt atime(2) = (M_HALF+sqrt(M_THREE)/6.0_real64)*dt @@ -115,15 +136,108 @@ contains end if end do - tr%vmagnus(:, :, 2) = M_HALF*(vaux(:, :, 1) + vaux(:, :, 2)) - tr%vmagnus(:, :, 1) = (sqrt(M_THREE)/12.0_real64)*dt*(vaux(:, :, 2) - vaux(:, :, 1)) + vmagnus(:, :, 2) = M_HALF*(vaux(:, :, 1) + vaux(:, :, 2)) + vmagnus(:, :, 1) = (sqrt(M_THREE)/12.0_real64)*dt*(vaux(:, :, 2) - vaux(:, :, 1)) - call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt, vmagnus = tr%vmagnus) + op => magnus4_operator_t(namespace, gr, hm, vmagnus) + call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt, op=op) + SAFE_DEALLOCATE_P(op) SAFE_DEALLOCATE_A(vaux) + SAFE_DEALLOCATE_A(vmagnus) POP_SUB(propagator_dt.td_magnus) end subroutine td_magnus + function magnus4_operator_constructor(namespace, mesh, hm, vmagnus) result(this) + type(namespace_t), target, intent(in) :: namespace + class(mesh_t), target, intent(in) :: mesh + class(hamiltonian_abst_t), target, intent(in) :: hm + real(real64), intent(in) :: vmagnus(:, :, :) + type(magnus4_operator_t), pointer :: this + + PUSH_SUB(magnus4_operator_constructor) + + allocate(this) + call this%init(namespace, mesh, hm) + + select type(hm) + class is (hamiltonian_elec_t) + ASSERT(size(vmagnus, 1) >= mesh%np) + ASSERT(size(vmagnus, 2) == hm%d%nspin) + ASSERT(size(vmagnus, 3) == 2) + class default + write(message(1), '(a)') 'Magnus operators only implemented for electrons at the moment' + call messages_fatal(1, namespace=namespace) + end select + SAFE_ALLOCATE_SOURCE(this%vmagnus, vmagnus) + + POP_SUB(magnus4_operator_constructor) + end function magnus4_operator_constructor + + subroutine magnus4_operator_apply(this, psib, hpsib) + class(magnus4_operator_t), intent(in) :: this + class(batch_t), intent(inout) :: psib + class(batch_t), intent(inout) :: hpsib + + integer :: ispin + type(wfs_elec_t) :: auxpsib, aux2psib + + PUSH_SUB(magnus4_operator_apply) + + select type(hm => this%hm) + class is (hamiltonian_elec_t) + select type (psib) + class is (wfs_elec_t) + select type (hpsib) + class is (wfs_elec_t) + ! We will assume, for the moment, no spinors. + if (hm%d%dim /= 1) call messages_not_implemented("Magnus with spinors", namespace=this%namespace) + + ASSERT(psib%nst == hpsib%nst) + + ispin = hm%d%get_spin_index(psib%ik) + + call hpsib%copy_to(auxpsib, copy_data=.false.) + call hpsib%copy_to(aux2psib, copy_data=.false.) + + ! Compute (T + Vnl)|psi> and store it + call zhamiltonian_elec_apply_batch(hm, this%namespace, this%mesh, psib, auxpsib, & + terms = TERM_KINETIC + TERM_NON_LOCAL_POTENTIAL) + + ! H|psi> = (T + Vnl)|psi> + Vpsl|psi> + Vmagnus(t2)|psi> + Vborders|psi> + call auxpsib%copy_data_to(this%mesh%np, hpsib) + call zhamiltonian_elec_external(hm, this%mesh, psib, hpsib) + if (hm%abs_boundaries%abtype == IMAGINARY_ABSORBING) then + call batch_mul(this%mesh%np, hm%abs_boundaries%mf(1:this%mesh%np), psib, aux2psib) + call batch_axpy(this%mesh%np, M_zI, aux2psib, hpsib) + end if + call batch_mul(this%mesh%np, this%vmagnus(1:this%mesh%np, ispin, 2), psib, aux2psib) + call batch_axpy(this%mesh%np, M_ONE, aux2psib, hpsib) + + ! Add first term of the commutator: - i Vmagnus(t1) (T + Vnl) |psi> + call batch_mul(this%mesh%np, this%vmagnus(1:this%mesh%np, ispin, 1), auxpsib, aux2psib) + call batch_axpy(this%mesh%np, -M_zI, aux2psib, hpsib) + + ! Add second term of commutator: i (T + Vnl) Vmagnus(t1) |psi> + call batch_mul(this%mesh%np, this%vmagnus(1:this%mesh%np, ispin, 1), psib, auxpsib) + call zhamiltonian_elec_apply_batch(hm, this%namespace, this%mesh, auxpsib, aux2psib, & + terms = TERM_KINETIC + TERM_NON_LOCAL_POTENTIAL) + call batch_axpy(this%mesh%np, M_zI, aux2psib, hpsib) + + call auxpsib%end() + call aux2psib%end() + class default + message(1) = "Internal error: imcompatible batch_t in magnus4_operator_apply for argument hpsib." + call messages_fatal(1, namespace=this%namespace) + end select + class default + message(1) = "Internal error: imcompatible batch_t in magnus4_operator_apply for argument psib." + call messages_fatal(1, namespace=this%namespace) + end select + end select + + POP_SUB(magnus4_operator_apply) + end subroutine magnus4_operator_apply ! --------------------------------------------------------- !> Commutator-free Magnus propagator of order 4.