From b0f78970415307634d405af59cb8e2393ad8035d Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 09:29:04 +0100 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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 From fbd2f3bb608d076b9818d9bbc4e47d35c66cc603 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 14:12:01 +0100 Subject: [PATCH 06/11] Simplify gauge field propagation call --- src/electrons/electrons.F90 | 7 ++----- src/td/propagation_ops_elec.F90 | 22 +++++++++++++--------- src/td/propagator_etrs.F90 | 30 ++++-------------------------- src/td/propagator_expmid.F90 | 31 ++++++++----------------------- 4 files changed, 27 insertions(+), 63 deletions(-) diff --git a/src/electrons/electrons.F90 b/src/electrons/electrons.F90 index 05fda6a82a..4b06338c42 100644 --- a/src/electrons/electrons.F90 +++ b/src/electrons/electrons.F90 @@ -663,11 +663,8 @@ contains this%mc, time+algo%dt, algo%dt) !Propagate gauge field - gfield => list_get_gauge_field(this%ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(this%td%tr%propagation_ops_elec, gfield, & - algo%dt, time+algo%dt) - end if + call propagation_ops_elec_propagate_gauge_field(this%td%tr%propagation_ops_elec, this%ext_partners, & + algo%dt, time+algo%dt) !Update Hamiltonian and current call get_fields_from_interaction(this, time+algo%dt) diff --git a/src/td/propagation_ops_elec.F90 b/src/td/propagation_ops_elec.F90 index 691c573b52..c1122dc920 100644 --- a/src/td/propagation_ops_elec.F90 +++ b/src/td/propagation_ops_elec.F90 @@ -194,26 +194,30 @@ contains end subroutine propagation_ops_elec_restore_ions ! --------------------------------------------------------- - subroutine propagation_ops_elec_propagate_gauge_field(wo, gfield, dt, time, save_gf) + subroutine propagation_ops_elec_propagate_gauge_field(wo, ext_partners, dt, time, save_gf) class(propagation_ops_elec_t), intent(inout) :: wo - type(gauge_field_t), intent(inout) :: gfield + type(partner_list_t), intent(in) :: ext_partners real(real64), intent(in) :: dt real(real64), intent(in) :: time logical, optional, intent(in) :: save_gf + type(gauge_field_t), pointer :: gfield PUSH_SUB(propagation_ops_elec_propagate_gauge_field) call profiling_in('ELEC_MOVE_GAUGE') - if (gauge_field_is_propagated(gfield)) then - if (optional_default(save_gf, .false.)) then - SAFE_ALLOCATE(wo%vecpot(1:gfield%space%dim)) - SAFE_ALLOCATE(wo%vecpot_vel(1:gfield%space%dim)) - call gauge_field_get_vec_pot(gfield, wo%vecpot) - call gauge_field_get_vec_pot_vel(gfield, wo%vecpot_vel) + gfield => list_get_gauge_field(ext_partners) + if (associated(gfield)) then + if (gauge_field_is_propagated(gfield)) then + if (optional_default(save_gf, .false.)) then + SAFE_ALLOCATE(wo%vecpot(1:gfield%space%dim)) + SAFE_ALLOCATE(wo%vecpot_vel(1:gfield%space%dim)) + call gauge_field_get_vec_pot(gfield, wo%vecpot) + call gauge_field_get_vec_pot_vel(gfield, wo%vecpot_vel) + end if + call gauge_field_do_algorithmic_operation(gfield, OP_VERLET_COMPUTE_ACC, dt, time) end if - call gauge_field_do_algorithmic_operation(gfield, OP_VERLET_COMPUTE_ACC, dt, time) end if call profiling_out('ELEC_MOVE_GAUGE') diff --git a/src/td/propagator_etrs.F90 b/src/td/propagator_etrs.F90 index 5392e4f94c..d2d80c37c2 100644 --- a/src/td/propagator_etrs.F90 +++ b/src/td/propagator_etrs.F90 @@ -86,7 +86,6 @@ contains type(multicomm_t), intent(inout) :: mc !< index and domain communicators type(xc_copied_potentials_t) :: vhxc_t1, vhxc_t2 - type(gauge_field_t), pointer :: gfield PUSH_SUB(td_etrs) @@ -114,11 +113,7 @@ contains ! first move the ions to time t call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, & ext_partners, mc, time, dt) - - gfield => list_get_gauge_field(ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time) - end if + call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, ext_partners, dt, time) call hm%ks_pot%restore_potentials(vhxc_t2) @@ -155,7 +150,6 @@ contains class(wfs_elec_t), allocatable :: psi2(:, :) ! these are hardcoded for the moment integer, parameter :: niter = 10 - type(gauge_field_t), pointer :: gfield type(xc_copied_potentials_t) :: vhxc_t1, vhxc_t2 PUSH_SUB(td_etrs_sc) @@ -185,11 +179,7 @@ contains ! first move the ions to time t call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, & ext_partners, mc, time, dt) - - gfield => list_get_gauge_field(ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time) - end if + call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, ext_partners, dt, time) call hm%ks_pot%restore_potentials(vhxc_t2) @@ -273,8 +263,6 @@ contains type(partner_list_t), intent(in) :: ext_partners type(multicomm_t), intent(inout) :: mc !< index and domain communicators - type(gauge_field_t), pointer :: gfield - PUSH_SUB(td_aetrs) ! propagate half of the time step with H(time - dt) @@ -286,12 +274,7 @@ contains ! move the ions to time t call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, & ions, ext_partners, mc, time, dt) - - !Propagate gauge field - gfield => list_get_gauge_field(ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time) - end if + call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, ext_partners, dt, time) !Update Hamiltonian call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time) @@ -326,7 +309,6 @@ contains type(density_calc_t) :: dens_calc integer(int64) :: pnp, wgsize, dim2, dim3 type(accel_mem_t) :: phase_buff - type(gauge_field_t), pointer :: gfield type(xc_copied_potentials_t) :: vold PUSH_SUB(td_caetrs) @@ -369,11 +351,7 @@ contains ! move the ions to time t call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, & ions, ext_partners, mc, time, dt) - - gfield => list_get_gauge_field(ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, dt, time) - end if + call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, ext_partners, dt, time) call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time) diff --git a/src/td/propagator_expmid.F90 b/src/td/propagator_expmid.F90 index 7b4db2b4a0..14bb37d06a 100644 --- a/src/td/propagator_expmid.F90 +++ b/src/td/propagator_expmid.F90 @@ -75,15 +75,11 @@ contains call hm%ks_pot%interpolate_potentials(tr%vks_old, 3, time, dt, time - dt/M_TWO) - !move the ions to time 'time - dt/2' + !move the ions and gauge field to time 'time - dt/2' call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, & ext_partners, mc, time - M_HALF*dt, M_HALF*dt, save_pos = .true.) - - gfield => list_get_gauge_field(ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, & - M_HALF*dt, time, save_gf = .true.) - end if + call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, ext_partners, & + M_HALF*dt, time, save_gf = .true.) call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt*M_HALF) @@ -91,10 +87,7 @@ contains !restore to time 'time - dt' call propagation_ops_elec_restore_ions(tr%propagation_ops_elec, ions_dyn, ions) - - if(associated(gfield)) then - call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners) - end if + call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners) POP_SUB(propagator_dt.exponential_midpoint) end subroutine exponential_midpoint @@ -119,7 +112,6 @@ contains type(partner_list_t), intent(in) :: ext_partners type(multicomm_t), intent(inout) :: mc !< index and domain communicators - type(gauge_field_t), pointer :: gfield type(states_elec_t) :: st0 PUSH_SUB(propagator_dt.exponential_midpoint_predictor) @@ -141,15 +133,11 @@ contains call states_elec_end(st0) end if - ! move the ions to time time - dt/2 + ! move the ions and gauge field to time time - dt/2 call propagation_ops_elec_move_ions(tr%propagation_ops_elec, gr, hm, st, namespace, space, ions_dyn, ions, & ext_partners, mc, time - M_HALF*dt, M_HALF*dt, save_pos = .true.) - - gfield => list_get_gauge_field(ext_partners) - if(associated(gfield)) then - call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, gfield, & - M_HALF*dt, time, save_gf = .true.) - end if + call propagation_ops_elec_propagate_gauge_field(tr%propagation_ops_elec, ext_partners, & + M_HALF*dt, time, save_gf = .true.) call propagation_ops_elec_update_hamiltonian(namespace, space, st, gr, hm, ext_partners, time - dt*M_HALF) @@ -158,10 +146,7 @@ contains !restore to time 'time - dt' call propagation_ops_elec_restore_ions(tr%propagation_ops_elec, ions_dyn, ions) - - if(associated(gfield)) then - call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners) - end if + call propagation_ops_elec_restore_gauge_field(tr%propagation_ops_elec, namespace, space, hm, gr, ext_partners) POP_SUB(propagator_dt.exponential_midpoint_predictor) end subroutine exponential_midpoint_predictor -- GitLab From bd2afe2385604a165951455af0230cfc23e0ffef Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 15:37:27 +0100 Subject: [PATCH 07/11] Add CFM4 propagator with moving ions and gauge field Apply Hamiltonian separately for both times inside Hamiltonian to capture all features (ions, gauge field, phases, ...). --- src/td/propagator_base.F90 | 3 +- src/td/propagator_elec.F90 | 12 ++- src/td/propagator_magnus.F90 | 170 ++++++++++++++++++++++++++++++++++- 3 files changed, 180 insertions(+), 5 deletions(-) diff --git a/src/td/propagator_base.F90 b/src/td/propagator_base.F90 index 0a537c1324..45393e3102 100644 --- a/src/td/propagator_base.F90 +++ b/src/td/propagator_base.F90 @@ -44,7 +44,8 @@ module propagator_base_oct_m PROP_RUNGE_KUTTA2 = 14, & PROP_EXPLICIT_RUNGE_KUTTA4 = 15, & PROP_CFMAGNUS4 = 16, & - PROP_EXPONENTIAL_MIDPOINT_PREDICTOR = 17 + PROP_EXPONENTIAL_MIDPOINT_PREDICTOR = 17, & + PROP_CFMAGNUS4_FULL = 18 type propagator_base_t ! Components are public by default diff --git a/src/td/propagator_elec.F90 b/src/td/propagator_elec.F90 index 9371eb571d..96fe95f1f9 100644 --- a/src/td/propagator_elec.F90 +++ b/src/td/propagator_elec.F90 @@ -237,6 +237,8 @@ contains !% WARNING: EXPERIMENTAL. Explicit RK4 method. !%Option cfmagnus4 16 !% WARNING EXPERIMENTAL. Commutator-free magnus expansion. + !%Option cfmagnus4_full 18 + !% WARNING EXPERIMENTAL. Commutator-free magnus expansion, including support for ions and gauge field. Might be slower. !%Option exp_mid_predictor 17 !% Exponential Midpoint Rule (EM). !% This is maybe the simplest method, but it is very well grounded theoretically: @@ -324,6 +326,8 @@ contains call messages_experimental("explicit Runge-Kutta 4 propagator", namespace=namespace) case (PROP_CFMAGNUS4) call messages_experimental("Commutator-Free Magnus propagator", namespace=namespace) + case (PROP_CFMAGNUS4_FULL) + call messages_experimental("Commutator-Free Magnus propagator", namespace=namespace) case default call messages_input_error(namespace, 'TDPropagator') end select @@ -339,9 +343,10 @@ contains tr%method /= PROP_RUNGE_KUTTA4 .and. & tr%method /= PROP_EXPLICIT_RUNGE_KUTTA4 .and. & tr%method /= PROP_RUNGE_KUTTA2 .and. & + tr%method /= PROP_CFMAGNUS4_FULL .and. & tr%method /= PROP_CRANK_NICOLSON_SPARSKIT) then message(1) = "To move the ions or put in a gauge field, use one of the following propagators:" - message(2) = "etrs, aetrs, crank-nicolson, rk2, rk4, or exp_mid." + message(2) = "etrs, aetrs, crank-nicolson, rk2, rk4, exp_mid, or cfmagnus4_full." call messages_fatal(2, namespace=namespace) end if end if @@ -349,6 +354,7 @@ contains if (relax_cell) then if (tr%method /= PROP_ETRS .and. & tr%method /= PROP_AETRS .and. & + tr%method /= PROP_CFMAGNUS4_FULL .and. & tr%method /= PROP_EXPONENTIAL_MIDPOINT .and. & tr%method /= PROP_EXPONENTIAL_MIDPOINT_PREDICTOR .and. & tr%method /= PROP_CRANK_NICOLSON .and. & @@ -362,6 +368,8 @@ contains select case (tr%method) case (PROP_CFMAGNUS4) call ks_pot%init_interpolation(tr%vks_old, order = 4) + case (PROP_CFMAGNUS4_FULL) + call ks_pot%init_interpolation(tr%vks_old, order = 6) case default call ks_pot%init_interpolation(tr%vks_old) end select @@ -553,6 +561,8 @@ contains end if case (PROP_CFMAGNUS4) call td_cfmagnus4(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, nt) + case (PROP_CFMAGNUS4_FULL) + call td_cfmagnus4_full(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc) end select generate = .false. diff --git a/src/td/propagator_magnus.F90 b/src/td/propagator_magnus.F90 index 7d2d124713..99dd8a259d 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -41,6 +41,7 @@ module propagator_magnus_oct_m use lasers_oct_m use messages_oct_m use mesh_oct_m + use multicomm_oct_m use namespace_oct_m use parser_oct_m use potential_interpolation_oct_m @@ -58,9 +59,10 @@ module propagator_magnus_oct_m private - public :: & - td_magnus, & - td_cfmagnus4 + public :: & + td_magnus, & + td_cfmagnus4, & + td_cfmagnus4_full !< @brief Class to encapsulate the operation for the Magnus 4 propagator type, extends(hamiltonian_operator_t) :: magnus4_operator_t @@ -74,6 +76,28 @@ module propagator_magnus_oct_m procedure magnus4_operator_constructor end interface magnus4_operator_t + !< @brief Class to encapsulate the operation for the commutator-free Magnus 4 propagator + type, extends(operator_t) :: cfm4_operator_t + private + type(v_ks_t), pointer :: ks + type(electron_space_t), pointer :: space + type(hamiltonian_elec_t), pointer :: hm + type(grid_t), pointer :: gr + type(states_elec_t), pointer :: st + type(propagator_base_t), pointer :: tr + type(ion_dynamics_t), pointer :: ions_dyn + type(ions_t), pointer :: ions + type(partner_list_t), pointer :: ext_partners + type(multicomm_t), pointer :: mc + real(real64) :: time, dt, alpha1, alpha2 + contains + procedure :: apply => cfm4_operator_apply + end type cfm4_operator_t + + interface cfm4_operator_t + procedure cfm4_operator_constructor + end interface cfm4_operator_t + contains @@ -304,6 +328,146 @@ contains POP_SUB(propagator_dt.td_cfmagnus4) end subroutine td_cfmagnus4 + ! --------------------------------------------------------- + !> Commutator-free Magnus propagator of order 4. + subroutine td_cfmagnus4_full(ks, namespace, space, hm, gr, st, tr, time, dt, ions_dyn, ions, ext_partners, mc) + type(v_ks_t), target, intent(inout) :: ks + type(namespace_t), intent(in) :: namespace + type(electron_space_t), intent(in) :: space + type(hamiltonian_elec_t), target, intent(inout) :: hm + type(grid_t), target, intent(in) :: gr + type(states_elec_t), target, intent(inout) :: st + type(propagator_base_t), target, intent(inout) :: tr + real(real64), intent(in) :: time + real(real64), intent(in) :: dt + type(ion_dynamics_t), intent(inout) :: ions_dyn + type(ions_t), intent(inout) :: ions + type(partner_list_t), intent(in) :: ext_partners + type(multicomm_t), intent(inout) :: mc !< index and domain communicators + + real(real64) :: alpha1, alpha2 + class(operator_t), pointer :: op + + PUSH_SUB(propagator_dt.td_cfmagnus4_full) + + alpha1 = (M_THREE - M_TWO * sqrt(M_THREE))/12.0_real64 + alpha2 = (M_THREE + M_TWO * sqrt(M_THREE))/12.0_real64 + + ! first exponential + op => cfm4_operator_t(ks, namespace, space, hm, gr, st, tr, ions_dyn, ions, ext_partners, mc, time, dt, alpha2, alpha1) + call propagation_ops_elec_exp_apply(tr%te, namespace, st, gr, hm, dt, op=op) + SAFE_DEALLOCATE_P(op) + + ! second exponential + op => cfm4_operator_t(ks, namespace, space, hm, gr, st, tr, ions_dyn, ions, ext_partners, mc, time, dt, alpha1, alpha2) + call propagation_ops_elec_fuse_density_exp_apply(tr%te, namespace, st, gr, hm, dt, op=op) + SAFE_DEALLOCATE_P(op) + + POP_SUB(propagator_dt.td_cfmagnus4_full) + end subroutine td_cfmagnus4_full + + function cfm4_operator_constructor(ks, namespace, space, hm, gr, st, tr, ions_dyn, ions, ext_partners, mc, & + time, dt, alpha1, alpha2) result(this) + type(v_ks_t), target, intent(inout) :: ks + type(namespace_t), target, intent(in) :: namespace + type(electron_space_t), target, intent(in) :: space + type(hamiltonian_elec_t), target, intent(inout) :: hm + type(grid_t), target, intent(in) :: gr + type(states_elec_t), target, intent(inout) :: st + type(propagator_base_t), target, intent(inout) :: tr + type(ion_dynamics_t), target, intent(inout) :: ions_dyn + type(ions_t), target, intent(inout) :: ions + type(partner_list_t), target, intent(in) :: ext_partners + type(multicomm_t), target, intent(inout) :: mc + real(real64), target, intent(in) :: time, dt, alpha1, alpha2 + type(cfm4_operator_t), pointer :: this + + PUSH_SUB(cfm4_operator_constructor) + + allocate(this) + this%ks => ks + this%namespace => namespace + this%space => space + this%hm => hm + this%gr => gr + this%st => st + this%tr => tr + this%ions_dyn => ions_dyn + this%ions => ions + this%ext_partners => ext_partners + this%mc => mc + this%time = time + this%dt = dt + this%alpha1 = alpha1 + this%alpha2 = alpha2 + + POP_SUB(cfm4_operator_constructor) + end function cfm4_operator_constructor + + subroutine cfm4_operator_apply(this, psib, hpsib) + class(cfm4_operator_t), intent(in) :: this + class(batch_t), intent(inout) :: psib + class(batch_t), intent(inout) :: hpsib + + type(wfs_elec_t) :: auxpsib + real(real64) :: alpha1, alpha2, c1, c2, t1, t2 + + PUSH_SUB(cfm4_operator_apply) + + c1 = M_HALF - sqrt(M_THREE)/6.0_real64 + c2 = M_HALF + sqrt(M_THREE)/6.0_real64 + t1 = this%time - this%dt + c1*this%dt + t2 = this%time - this%dt + c2*this%dt + + select type (psib) + class is (wfs_elec_t) + select type (hpsib) + class is (wfs_elec_t) + call hpsib%copy_to(auxpsib, copy_data=.false.) + + call this%hm%ks_pot%interpolate_potentials(this%tr%vks_old, 6, this%time, this%dt, t1) + !move the ions and gauge field to t1 + call propagation_ops_elec_move_ions(this%tr%propagation_ops_elec, this%gr, this%hm, this%st, this%namespace, & + this%space, this%ions_dyn, this%ions, this%ext_partners, this%mc, t1, c1*this%dt, save_pos = .true.) + call propagation_ops_elec_propagate_gauge_field(this%tr%propagation_ops_elec, this%ext_partners, & + c1*this%dt, this%time, save_gf = .true.) + call propagation_ops_elec_update_hamiltonian(this%namespace, this%space, this%st, this%gr, this%hm, & + this%ext_partners, t1) + call this%hm%zapply(this%namespace, this%gr, psib, hpsib) + call batch_scal(this%gr%np, this%alpha1, hpsib) + + call propagation_ops_elec_restore_ions(this%tr%propagation_ops_elec, this%ions_dyn, this%ions) + call propagation_ops_elec_restore_gauge_field(this%tr%propagation_ops_elec, this%namespace, this%space, & + this%hm, this%gr, this%ext_partners) + + call this%hm%ks_pot%interpolate_potentials(this%tr%vks_old, 6, this%time, this%dt, t2) + !move the ions and gauge field to t2 + call propagation_ops_elec_move_ions(this%tr%propagation_ops_elec, this%gr, this%hm, this%st, this%namespace, & + this%space, this%ions_dyn, this%ions, this%ext_partners, this%mc, t2, c2*this%dt, save_pos = .true.) + call propagation_ops_elec_propagate_gauge_field(this%tr%propagation_ops_elec, this%ext_partners, & + c2*this%dt, this%time, save_gf = .true.) + call propagation_ops_elec_update_hamiltonian(this%namespace, this%space, this%st, this%gr, this%hm, & + this%ext_partners, t2) + call this%hm%zapply(this%namespace, this%gr, psib, auxpsib) + call batch_axpy(this%gr%np, this%alpha2, auxpsib, hpsib) + + call propagation_ops_elec_restore_ions(this%tr%propagation_ops_elec, this%ions_dyn, this%ions) + call propagation_ops_elec_restore_gauge_field(this%tr%propagation_ops_elec, this%namespace, this%space, & + this%hm, this%gr, this%ext_partners) + + class default + message(1) = "Internal error: imcompatible batch_t in cfm4_operator_apply for argument hpsib." + call messages_fatal(1, namespace=this%namespace) + end select + class default + message(1) = "Internal error: imcompatible batch_t in cfm4_operator_apply for argument psib." + call messages_fatal(1, namespace=this%namespace) + end select + + POP_SUB(cfm4_operator_apply) + end subroutine cfm4_operator_apply + + end module propagator_magnus_oct_m !! Local Variables: -- GitLab From e863d22d37481bf668594733724d6b08cb2043f9 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 15:56:03 +0100 Subject: [PATCH 08/11] Add test for cfmagnus4_full propagator --- .../01-propagators.13-cfmagnus4_full.inp | 41 +++++++++++++++++++ testsuite/real_time/01-propagators.test | 10 +++++ 2 files changed, 51 insertions(+) create mode 100644 testsuite/real_time/01-propagators.13-cfmagnus4_full.inp diff --git a/testsuite/real_time/01-propagators.13-cfmagnus4_full.inp b/testsuite/real_time/01-propagators.13-cfmagnus4_full.inp new file mode 100644 index 0000000000..436e3bcdcb --- /dev/null +++ b/testsuite/real_time/01-propagators.13-cfmagnus4_full.inp @@ -0,0 +1,41 @@ +CalculationMode = td +ExperimentalFeatures = yes +FromScratch = yes + +bond_length = 2.8 + +%Coordinates +"C" | -bond_length/2 | 0.0 | 0.0 +"C" | bond_length/2 | 0.0 | 0.0 +% + +ExtraStates = 1 +%Occupations + 2 | 2 | 4/3 | 4/3 | 4/3 +% + +# We use the default values for Carbon, except the nuclear mass, that is +# divided by 1000. +%Species +"C" | species_pseudo | lmax | 1 | lloc| 0 | mass | 0.0120107000 +% + +Spacing = 0.6 + +Radius = 6 + +MoveIons = no + +TDPropagator = cfmagnus4_full + +TDMaxSteps = 20 +TDTimeStep = 0.1 + +TDDeltaStrength = 0.01 +TDPolarizationDirection = 1 + +TDEnergyUpdateIter = 1 + + + +PoissonSolver = isf diff --git a/testsuite/real_time/01-propagators.test b/testsuite/real_time/01-propagators.test index 233358ac01..c37761d64f 100644 --- a/testsuite/real_time/01-propagators.test +++ b/testsuite/real_time/01-propagators.test @@ -153,3 +153,13 @@ Precision: 3.00e-14 match ; Forces [step 0] ; LINEFIELD(td.general/coordinates, -21, 15) ; 1.4215273371425141e-03 Precision: 2.08e-14 match ; Forces [step 20] ; LINEFIELD(td.general/coordinates, -1, 15) ; 1.6743533403078992e-02 + +Processors : 4 +Input : 01-propagators.13-cfmagnus4_full.inp +Precision: 1.20e-13 +match ; Energy [step 0] ; LINEFIELD(td.general/energy, -21, 3) ; -1.0610604457225079e+01 +match ; Energy [step 20] ; LINEFIELD(td.general/energy, -1, 3) ; -1.0610210266499825e+01 +Precision: 2.80e-14 +match ; Multipoles [step 0] ; LINEFIELD(td.general/multipoles, -21, 4) ; 0.0 +Precision: 3.25e-12 +match ; Multipoles [step 20] ; LINEFIELD(td.general/multipoles, -1, 4) ; -1.0902162672451297e-01 -- GitLab From 4ad39933122804d6b060d80b714dcf23c088d7d8 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 16:02:47 +0100 Subject: [PATCH 09/11] Fix indentation --- src/td/propagator_magnus.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/td/propagator_magnus.F90 b/src/td/propagator_magnus.F90 index 99dd8a259d..8b67d8abc1 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -367,7 +367,7 @@ contains end subroutine td_cfmagnus4_full function cfm4_operator_constructor(ks, namespace, space, hm, gr, st, tr, ions_dyn, ions, ext_partners, mc, & - time, dt, alpha1, alpha2) result(this) + time, dt, alpha1, alpha2) result(this) type(v_ks_t), target, intent(inout) :: ks type(namespace_t), target, intent(in) :: namespace type(electron_space_t), target, intent(in) :: space -- GitLab From 66ff304ffc48139be00f46b6d61b3296ac8c2737 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 16:11:53 +0100 Subject: [PATCH 10/11] Remove unused variable --- src/td/propagator_magnus.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/td/propagator_magnus.F90 b/src/td/propagator_magnus.F90 index 8b67d8abc1..1a65540b48 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -410,7 +410,7 @@ contains class(batch_t), intent(inout) :: hpsib type(wfs_elec_t) :: auxpsib - real(real64) :: alpha1, alpha2, c1, c2, t1, t2 + real(real64) :: c1, c2, t1, t2 PUSH_SUB(cfm4_operator_apply) -- GitLab From 8c9185a84a0db2676abe85da482bf0ccd6aac6d1 Mon Sep 17 00:00:00 2001 From: Sebastian Ohlmann Date: Thu, 18 Dec 2025 16:34:07 +0100 Subject: [PATCH 11/11] Fix memory leak --- src/td/propagator_magnus.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/td/propagator_magnus.F90 b/src/td/propagator_magnus.F90 index 1a65540b48..2848623037 100644 --- a/src/td/propagator_magnus.F90 +++ b/src/td/propagator_magnus.F90 @@ -455,6 +455,7 @@ contains call propagation_ops_elec_restore_gauge_field(this%tr%propagation_ops_elec, this%namespace, this%space, & this%hm, this%gr, this%ext_partners) + call auxpsib%end() class default message(1) = "Internal error: imcompatible batch_t in cfm4_operator_apply for argument hpsib." call messages_fatal(1, namespace=this%namespace) -- GitLab