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