From b0f78970415307634d405af59cb8e2393ad8035d Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 09:29:04 +0100 Subject: [PATCH 1/5] Encapsulate operation inside exponential in a class As a first step, this allows to remove the vmagnus argument to many functions. As a next step, this allows to implement more complex propagation schemes that involve combinations of Hamiltonians at different times. --- src/electrons/exponential.F90 | 242 ++++++++++++++++++++++++++-------- 1 file changed, 184 insertions(+), 58 deletions(-) diff --git a/src/electrons/exponential.F90 b/src/electrons/exponential.F90 index 53ac0b7d46..52900f416f 100644 --- a/src/electrons/exponential.F90 +++ b/src/electrons/exponential.F90 @@ -78,6 +78,50 @@ 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 + private + 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 + private + 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 + + !< @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 ! --------------------------------------------------------- @@ -304,7 +348,7 @@ 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) + subroutine exponential_apply_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, imag_time, inh_psib, op) class(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh @@ -316,10 +360,12 @@ contains 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 + 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") @@ -343,6 +389,16 @@ contains end select end if + if (present(op)) then + op_ => op + else + if (present(vmagnus)) then + op_ => magnus4_operator_t(namespace, mesh, hm, vmagnus) + else + op_ => hamiltonian_operator_t(namespace, mesh, hm) + end if + end if + deltat2_ = cmplx(optional_default(deltat2, M_ZERO), M_ZERO, real64) imag_time_ = optional_default(imag_time, .false.) @@ -365,32 +421,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,10 +468,10 @@ 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 @@ -427,16 +482,16 @@ contains 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 +537,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 +567,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 +584,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 +606,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 +619,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 +669,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 +767,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 +807,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 +820,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 +867,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 +960,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, vmagnus, op) class(exponential_t), intent(inout) :: te type(namespace_t), intent(in) :: namespace class(mesh_t), intent(in) :: mesh @@ -934,9 +969,11 @@ contains 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) @@ -962,17 +999,27 @@ contains call messages_warning(5, namespace=namespace) end if + if (present(op)) then + op_ => op + else + if (present(vmagnus)) then + op_ => magnus4_operator_t(namespace, mesh, hm, vmagnus) + else + op_ => hamiltonian_operator_t(namespace, mesh, hm) + end if + 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 +1034,92 @@ 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 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 + + 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 + + PUSH_SUB(magnus4_operator_apply) + + call this%hm%zmagnus_apply(this%namespace, this%mesh, psib, hpsib, this%vmagnus) + + POP_SUB(magnus4_operator_apply) + end subroutine magnus4_operator_apply + end module exponential_oct_m !! Local Variables: -- GitLab From ac16b68fe6dba7ad6ce3993a8714e4bdc057ec09 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 10:26:04 +0100 Subject: [PATCH 2/5] Move definition of Magnus operator to propagator module With this, the vmagnus argument can be completely removed as it is stored inside the operator class. --- src/electrons/exponential.F90 | 130 ++++++-------------------------- src/td/propagation_ops_elec.F90 | 47 ++++++------ src/td/propagator_magnus.F90 | 57 +++++++++++++- 3 files changed, 105 insertions(+), 129 deletions(-) diff --git a/src/electrons/exponential.F90 b/src/electrons/exponential.F90 index 52900f416f..5eb19ebf60 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, & @@ -80,7 +82,6 @@ module exponential_oct_m !< @brief Class to encapsulate the operation inside the exponential type, abstract :: operator_t - private type(namespace_t), pointer :: namespace class(mesh_t), pointer :: mesh contains @@ -99,7 +100,6 @@ module exponential_oct_m !< @brief Class to encapsulate the Hamiltonian inside the exponential type, extends(operator_t) :: hamiltonian_operator_t - private class(hamiltonian_abst_t), pointer :: hm contains procedure :: apply => hamiltonian_operator_apply @@ -110,18 +110,6 @@ module exponential_oct_m procedure hamiltonian_operator_constructor end interface hamiltonian_operator_t - !< @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 ! --------------------------------------------------------- @@ -283,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 @@ -292,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 @@ -313,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() @@ -348,19 +335,18 @@ 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, 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 - 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 - class(operator_t), target,optional, intent(in) :: op !< operator to be applied in the exponential + 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 @@ -377,26 +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 - end if - if (present(op)) then op_ => op else - if (present(vmagnus)) then - op_ => magnus4_operator_t(namespace, mesh, hm, vmagnus) - else - op_ => hamiltonian_operator_t(namespace, mesh, hm) - end if + op_ => hamiltonian_operator_t(namespace, mesh, hm) end if deltat2_ = cmplx(optional_default(deltat2, M_ZERO), M_ZERO, real64) @@ -477,6 +447,10 @@ contains 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 @@ -960,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, op) + 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 @@ -968,7 +942,6 @@ 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 @@ -978,17 +951,6 @@ contains 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' @@ -1002,11 +964,7 @@ contains if (present(op)) then op_ => op else - if (present(vmagnus)) then - op_ => magnus4_operator_t(namespace, mesh, hm, vmagnus) - else - op_ => hamiltonian_operator_t(namespace, mesh, hm) - end if + op_ => hamiltonian_operator_t(namespace, mesh, hm) end if deltat_ = cmplx(deltat, M_ZERO, real64) @@ -1082,44 +1040,6 @@ contains POP_SUB(hamiltonian_operator_apply) end subroutine hamiltonian_operator_apply - 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 - - PUSH_SUB(magnus4_operator_apply) - - call this%hm%zmagnus_apply(this%namespace, this%mesh, psib, hpsib, this%vmagnus) - - POP_SUB(magnus4_operator_apply) - end subroutine magnus4_operator_apply - end module exponential_oct_m !! Local Variables: diff --git a/src/td/propagation_ops_elec.F90 b/src/td/propagation_ops_elec.F90 index 1babb5cc78..691c573b52 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_magnus.F90 b/src/td/propagator_magnus.F90 index 150530bd2a..f29a84f524 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -28,6 +28,7 @@ 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 interaction_partner_oct_m use ion_dynamics_oct_m @@ -36,6 +37,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 @@ -56,6 +58,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 ! --------------------------------------------------------- @@ -74,6 +89,7 @@ contains real(real64) :: atime(2) real(real64), allocatable :: vaux(:, :, :), pot(:) type(lasers_t), pointer :: lasers + class(operator_t), pointer :: op PUSH_SUB(propagator_dt.td_magnus) @@ -118,7 +134,9 @@ contains tr%vmagnus(:, :, 2) = M_HALF*(vaux(:, :, 1) + vaux(:, :, 2)) tr%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, tr%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) POP_SUB(propagator_dt.td_magnus) @@ -190,6 +208,43 @@ contains POP_SUB(propagator_dt.td_cfmagnus4) end subroutine td_cfmagnus4 + 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 + + PUSH_SUB(magnus4_operator_apply) + + call this%hm%zmagnus_apply(this%namespace, this%mesh, psib, hpsib, this%vmagnus) + + POP_SUB(magnus4_operator_apply) + end subroutine magnus4_operator_apply end module propagator_magnus_oct_m !! Local Variables: -- GitLab From 0490c9f001d00d17fb8321ef6a248cefb8612981 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 12:04:05 +0100 Subject: [PATCH 3/5] Move all magnus-related code to the propagator file It didn't make much sense that the abstract Hamiltonian base class exposes an interface for this functionality that is only used for electrons. With the new operator classes, all the code related to the propagation is in the corresponding propagator module. --- src/hamiltonian/hamiltonian_abst.F90 | 22 ---- src/hamiltonian/hamiltonian_elec.F90 | 78 +------------ src/hamiltonian/hamiltonian_elec_inc.F90 | 80 ------------- src/maxwell/hamiltonian_mxll.F90 | 32 ------ src/td/propagator_base.F90 | 1 - src/td/propagator_elec.F90 | 8 -- src/td/propagator_magnus.F90 | 139 ++++++++++++++++------- 7 files changed, 99 insertions(+), 261 deletions(-) diff --git a/src/hamiltonian/hamiltonian_abst.F90 b/src/hamiltonian/hamiltonian_abst.F90 index 7ae72767d5..33d9c5b87b 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 b9f812defd..b24de19390 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 4cad89e373..1d595d8a02 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 f4a38251a4..592ae66654 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/propagator_base.F90 b/src/td/propagator_base.F90 index bf0ad77a9f..0a537c1324 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 a53093d43c..9371eb571d 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 f29a84f524..78606e8aee 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 @@ -30,6 +32,7 @@ module propagator_magnus_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 @@ -48,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 @@ -87,7 +91,7 @@ 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 @@ -96,6 +100,7 @@ contains 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 @@ -131,17 +136,106 @@ 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)) - op => magnus4_operator_t(namespace, gr, hm, 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. @@ -208,43 +302,6 @@ contains POP_SUB(propagator_dt.td_cfmagnus4) end subroutine td_cfmagnus4 - 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 - - PUSH_SUB(magnus4_operator_apply) - - call this%hm%zmagnus_apply(this%namespace, this%mesh, psib, hpsib, this%vmagnus) - - POP_SUB(magnus4_operator_apply) - end subroutine magnus4_operator_apply end module propagator_magnus_oct_m !! Local Variables: -- GitLab From f87cf7b95f7b1f10923c1b37237987e97116923e Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 12:10:26 +0100 Subject: [PATCH 4/5] Fix line length --- src/td/propagator_magnus.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/td/propagator_magnus.F90 b/src/td/propagator_magnus.F90 index 78606e8aee..7d2d124713 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -201,7 +201,8 @@ contains 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) + 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) @@ -219,7 +220,8 @@ contains ! 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 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() -- GitLab From adc4d531108d74797587650d2731500e3845c267 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 14:08:00 +0100 Subject: [PATCH 5/5] Fix memory leak for phi functions --- src/electrons/exponential.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/electrons/exponential.F90 b/src/electrons/exponential.F90 index 5eb19ebf60..c872caadbe 100644 --- a/src/electrons/exponential.F90 +++ b/src/electrons/exponential.F90 @@ -994,8 +994,12 @@ contains end if 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 -- GitLab