34 use,
intrinsic :: iso_fortran_env
62 integer,
public,
parameter :: &
69 integer,
public :: exp_method
70 real(real64) :: lanczos_tol
71 real(real64) :: chebyshev_tol
72 integer,
public :: exp_order
74 logical,
public :: full_batch = .false.
85 type(exponential_t),
intent(out) :: te
86 type(namespace_t),
intent(in) :: namespace
87 logical,
optional,
intent(in) :: full_batch
146 select case (te%exp_method)
159 call parse_variable(namespace,
'TDChebyshevTol', 1e-10_real64, te%chebyshev_tol)
190 call parse_variable(namespace,
'TDExpOrder', default__tdexporder, te%exp_order)
194 message(1) =
"TDExpOrder is only relevant for TDExponentialMethod = taylor"
199 te%arnoldi_gs = option__arnoldiorthogonalization__cgs
200 if (te%exp_method == exp_lanczos)
then
215 call parse_variable(namespace,
'ArnoldiOrthogonalization', option__arnoldiorthogonalization__cgs, &
232 teo%exp_method = tei%exp_method
233 teo%lanczos_tol = tei%lanczos_tol
234 teo%exp_order = tei%exp_order
235 teo%arnoldi_gs = tei%arnoldi_gs
245 class(
mesh_t),
intent(in) :: mesh
247 integer,
intent(in) :: ist
248 integer,
intent(in) :: ik
249 complex(real64),
contiguous,
intent(inout) :: zpsi(:, :)
250 real(real64),
intent(in) :: deltat
251 real(real64),
optional,
intent(in) :: vmagnus(mesh%np, hm%d%nspin, 2)
252 logical,
optional,
intent(in) :: imag_time
255 complex(real64),
allocatable :: zpsi_inh(:, :)
261 if (hm%phase%is_allocated())
then
262 call hm%phase%apply_to_single(zpsi, mesh%np, hm%d%dim, ik, .false.)
268 safe_allocate(zpsi_inh(1:mesh%np_part, 1:hm%d%dim))
270 call wfs_elec_init(inh_psib, hm%d%dim, ist, ist, zpsi_inh, ik)
271 call te%apply_batch(namespace, mesh, hm, psib, deltat, &
272 vmagnus=vmagnus, imag_time=imag_time, inh_psib=inh_psib)
274 safe_deallocate_a(zpsi_inh)
276 call te%apply_batch(namespace, mesh, hm, psib, deltat, &
277 vmagnus=vmagnus, imag_time=imag_time)
282 if (hm%phase%is_allocated())
then
283 call hm%phase%apply_to_single(zpsi, mesh%np, hm%d%dim, ik, .
true.)
307 subroutine exponential_apply_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, imag_time, inh_psib)
310 class(
mesh_t),
intent(in) :: mesh
312 class(
batch_t),
intent(inout) :: psib
313 real(real64),
intent(in) :: deltat
314 class(
batch_t),
optional,
intent(inout) :: psib2
315 real(real64),
optional,
intent(in) :: deltat2
316 real(real64),
optional,
intent(in) :: vmagnus(:,:,:)
317 logical,
optional,
intent(in) :: imag_time
318 class(
batch_t),
optional,
intent(inout) :: inh_psib
320 complex(real64) :: deltat_, deltat2_
322 logical :: imag_time_
329 assert(
present(psib2) .eqv.
present(deltat2))
330 if (
present(inh_psib))
then
331 assert(inh_psib%nst == psib%nst)
334 if (
present(vmagnus))
then
337 assert(
size(vmagnus, 1) >= mesh%np)
338 assert(
size(vmagnus, 2) == hm%d%nspin)
339 assert(
size(vmagnus, 3) == 2)
341 write(
message(1),
'(a)')
'Magnus operators only implemented for electrons at the moment'
350 deltat_ = -
m_zi*deltat
351 if (
present(deltat2)) deltat2_ =
m_zi*deltat2
353 deltat_ = cmplx(deltat,
m_zero, real64)
354 if (
present(deltat2)) deltat2_ = cmplx(deltat2,
m_zero, real64)
357 if (.not. hm%is_hermitian() .and. te%exp_method ==
exp_chebyshev)
then
358 write(
message(1),
'(a)')
'The Chebyshev expansion cannot be used for non-Hermitian operators.'
359 write(
message(2),
'(a)')
'Please use the Lanczos exponentiation scheme ("TDExponentialMethod = lanczos")'
360 write(
message(3),
'(a)')
'or the Taylor expansion ("TDExponentialMethod = taylor") method.'
364 select case (te%exp_method)
367 if (
present(deltat2))
then
369 psib2, deltat2_, vmagnus)
374 if (
present(inh_psib))
then
375 if (
present(deltat2))
then
377 psib2, deltat2_, vmagnus, inh_psib)
380 vmagnus=vmagnus, inh_psib=inh_psib)
385 if (
present(psib2))
call psib%copy_data_to(mesh%np, psib2)
387 if (
present(inh_psib))
then
390 if (
present(psib2))
then
392 if (
present(inh_psib))
then
398 if (
present(inh_psib))
then
399 write(
message(1),
'(a)')
'Chebyshev exponential ("TDExponentialMethod = chebyshev")'
400 write(
message(2),
'(a)')
'with inhomogeneous term is not implemented'
406 if (
present(psib2))
then
410 chebyshev_function =>
chebyshev_exp_t(hm%spectral_half_span, hm%spectral_middle_point, deltat)
411 if (
present(psib2))
then
412 chebyshev_function_dt2 =>
chebyshev_exp_t(hm%spectral_half_span, hm%spectral_middle_point, deltat2)
415 if (
present(psib2))
call psib%copy_data_to(mesh%np, psib2)
417 deallocate(chebyshev_function)
418 if (
present(psib2))
then
420 deallocate(chebyshev_function_dt2)
430 subroutine exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, inh_psib, phik_shift)
433 class(
mesh_t),
intent(in) :: mesh
435 class(
batch_t),
intent(inout) :: psib
436 complex(real64),
intent(in) :: deltat
437 class(
batch_t),
optional,
intent(inout) :: psib2
438 complex(real64),
optional,
intent(in) :: deltat2
439 real(real64),
optional,
intent(in) :: vmagnus(:,:,:)
440 class(
batch_t),
optional,
intent(inout) :: inh_psib
441 integer,
optional,
intent(in) :: phik_shift
443 complex(real64) :: zfact, zfact2
444 integer :: iter, denom, phik_shift_
445 logical :: zfact_is_real
446 class(
batch_t),
allocatable :: psi1b, hpsi1b
451 call psib%clone_to(psi1b)
452 call psib%clone_to(hpsi1b)
456 zfact_is_real = abs(deltat-real(deltat)) <
m_epsilon
458 if (
present(psib2))
call psib%copy_data_to(mesh%np, psib2)
460 if (
present(inh_psib))
then
462 call batch_axpy(mesh%np, real(zfact), inh_psib, psib)
464 if (
present(psib2))
then
465 zfact2 = zfact2*deltat2
466 call batch_axpy(mesh%np, real(zfact2), inh_psib, psib2)
473 do iter = 1, te%exp_order
474 denom = iter+phik_shift_
475 if (
present(inh_psib)) denom = denom + 1
476 zfact = zfact*(-
m_zi*deltat)/denom
477 if (
present(deltat2)) zfact2 = zfact2*(-
m_zi*deltat2)/denom
478 zfact_is_real = .not. zfact_is_real
485 call operate_batch(hm, namespace, mesh, psi1b, hpsi1b, vmagnus)
487 if (
present(inh_psib))
then
488 call operate_batch(hm, namespace, mesh, inh_psib, hpsi1b, vmagnus)
490 call operate_batch(hm, namespace, mesh, psib, hpsi1b, vmagnus)
494 if (zfact_is_real)
then
495 call batch_axpy(mesh%np, real(zfact), hpsi1b, psib)
496 if (
present(psib2))
call batch_axpy(mesh%np, real(zfact2), hpsi1b, psib2)
499 if (
present(psib2))
call batch_axpy(mesh%np, zfact2, hpsi1b, psib2)
502 if (iter /= te%exp_order)
call hpsi1b%copy_data_to(mesh%np, psi1b)
508 safe_deallocate_a(psi1b)
509 safe_deallocate_a(hpsi1b)
518 class(
mesh_t),
intent(in) :: mesh
520 class(
batch_t),
intent(inout) :: psib
521 complex(real64),
intent(in) :: deltat
522 real(real64),
optional,
intent(in) :: vmagnus(:,:,:)
523 class(
batch_t),
optional,
intent(in) :: inh_psib
529 if (
present(inh_psib))
then
530 call inh_psib%clone_to(tmpb, copy_data=.
true.)
533 call batch_axpy(mesh%np, real(deltat, real64), tmpb, psib)
557 class(
mesh_t),
intent(in) :: mesh
559 class(
batch_t),
intent(inout) :: psib
560 complex(real64),
intent(in) :: deltat
562 complex(real64) function fun(z)
564 complex(real64),
intent(in) :: z
567 real(real64),
optional,
intent(in) :: vmagnus(:,:,:)
569 integer :: iter, l, ii, ist, max_initialized
570 complex(real64),
allocatable :: hamilt(:,:,:), expo(:,:,:)
571 real(real64),
allocatable :: beta(:), res(:), norm(:)
572 integer,
parameter :: max_order = 200
578 if (te%exp_method /= exp_lanczos)
then
579 message(1) =
"The exponential method needs to be set to Lanzcos (TDExponentialMethod=lanczos)."
583 safe_allocate(beta(1:psib%nst))
584 safe_allocate(res(1:psib%nst))
585 safe_allocate(norm(1:psib%nst))
586 safe_allocate(vb(1:max_order))
587 call psib%clone_to(vb(1)%p)
590 call psib%copy_data_to(mesh%np, vb(1)%p, async=.
true.)
593 if (te%full_batch) beta = norm2(beta)
596 if (all(abs(beta) <= 1.0e-12_real64))
then
597 safe_deallocate_a(beta)
598 safe_deallocate_a(res)
599 safe_deallocate_a(norm)
601 safe_deallocate_a(vb)
609 safe_allocate(hamilt(1:max_order+1, 1:max_order+1, 1:psib%nst))
610 safe_allocate( expo(1:max_order+1, 1:max_order+1, 1:psib%nst))
613 do iter = 1, max_order-1
614 call psib%clone_to(vb(iter + 1)%p)
615 max_initialized = iter + 1
618 call operate_batch(hm, namespace, mesh, vb(iter)%p, vb(iter+1)%p, vmagnus)
621 if (hm%is_hermitian())
then
623 hamilt(1:max(l-1, 1), iter, 1:psib%nst) =
m_zero
627 hamilt(iter, 1:iter-2, 1:psib%nst) =
m_zero
633 normalize = .false., overlap = hamilt(l:iter, iter, 1:psib%nst), norm = hamilt(iter + 1, iter, 1:psib%nst), &
634 gs_scheme = te%arnoldi_gs, full_batch = te%full_batch)
644 hm%is_hermitian(), tridiagonal=hm%is_hermitian())
645 res(ii) = abs(hamilt(iter + 1, iter, ii) * abs(expo(iter, 1, ii)))
649 if (all(abs(hamilt(iter + 1, iter, :)) < 1.0e4_real64*
m_epsilon))
exit
654 if (abs(hamilt(iter + 1, iter, ist)) >= 1.0e4_real64 *
m_epsilon)
then
655 norm(ist) =
m_one / abs(hamilt(iter + 1, iter, ist))
659 call batch_scal(mesh%np, norm, vb(iter+1)%p, a_full = .false.)
661 if (iter > 3 .and. all(res < te%lanczos_tol))
exit
665 if (iter == max_order)
then
666 write(
message(1),
'(a,i5,a,es9.2)')
'Lanczos exponential expansion did not converge after ', iter, &
667 ' iterations. Residual: ', maxval(res)
670 write(
message(1),
'(a,i5)')
'Debug: Lanczos exponential iterations: ', iter
676 call batch_scal(mesh%np, expo(1,1,1:psib%nst), psib, a_full = .false.)
679 call batch_axpy(mesh%np, beta(1:psib%nst)*expo(ii,1,1:psib%nst), vb(ii)%p, psib, a_full = .false.)
684 do ii = 1, max_initialized
688 safe_deallocate_a(vb)
689 safe_deallocate_a(hamilt)
690 safe_deallocate_a(expo)
691 safe_deallocate_a(beta)
692 safe_deallocate_a(res)
693 safe_deallocate_a(norm)
718 class(
mesh_t),
intent(in) :: mesh
720 class(
batch_t),
intent(inout) :: psib
722 real(real64),
optional,
intent(in) :: vmagnus(:,:,:)
724 integer :: j, order_needed
725 complex(real64) :: coefficient
726 complex(real64),
allocatable :: coefficients(:)
727 real(real64) :: error
728 class(
batch_t),
allocatable,
target :: psi0, psi1, psi2
729 class(
batch_t),
pointer :: psi_n, psi_n1, psi_n2
730 integer,
parameter :: max_order = 200
735 call psib%clone_to(psi0)
736 call psib%clone_to(psi1)
737 call psib%clone_to(psi2)
738 call psib%copy_data_to(mesh%np, psi0)
740 order_needed = max_order
742 error = chebyshev_function%get_error(j)
743 if (error >
m_zero .and. error < te%chebyshev_tol)
then
749 call chebyshev_function%get_coefficients(j, coefficients)
752 call batch_scal(mesh%np, coefficients(0), psib)
756 call batch_axpy(mesh%np, -hm%spectral_middle_point, psi0, psi1)
759 call batch_axpy(mesh%np, coefficients(1), psi1, psib)
766 do j = 2, order_needed
768 call operate_batch(hm, namespace, mesh, psi_n1, psi_n, vmagnus)
769 call batch_axpy(mesh%np, -hm%spectral_middle_point, psi_n1, psi_n)
774 call batch_axpy(mesh%np, coefficients(j), psi_n, psib)
777 if (mod(j, 3) == 2)
then
781 else if (mod(j, 3) == 1)
then
792 if (order_needed == max_order)
then
793 write(
message(1),
'(a,i5,a,es9.2)')
'Chebyshev exponential expansion did not converge after ', j, &
794 ' iterations. Coefficient: ', coefficient
797 write(
message(1),
'(a,i5)')
'Debug: Chebyshev exponential iterations: ', j
804 safe_deallocate_a(psi0)
805 safe_deallocate_a(psi1)
806 safe_deallocate_a(psi2)
808 safe_deallocate_a(coefficients)
816 subroutine operate_batch(hm, namespace, mesh, psib, hpsib, vmagnus)
819 class(
mesh_t),
intent(in) :: mesh
820 class(
batch_t),
intent(inout) :: psib
821 class(
batch_t),
intent(inout) :: hpsib
822 real(real64),
optional,
intent(in) :: vmagnus(:, :, :)
826 if (
present(vmagnus))
then
827 call hm%zmagnus_apply(namespace, mesh, psib, hpsib, vmagnus)
829 call hm%zapply(namespace, mesh, psib, hpsib)
842 type(
grid_t),
intent(inout) :: gr
845 real(real64),
intent(in) :: deltat
846 integer,
optional,
intent(inout) :: order
849 real(real64) :: zfact
861 do i = 1, te%exp_order
862 zfact = zfact * deltat / i
870 do ik = st%d%kpt%start, st%d%kpt%end
871 do ib = st%group%block_start, st%group%block_end
873 call batch_axpy(gr%np, -
m_zi, hst1%group%psib(ib, ik), st1%group%psib(ib, ik))
874 call batch_axpy(gr%np, zfact, st1%group%psib(ib, ik), st%group%psib(ib, ik))
892 do ik = st%d%kpt%start, st%d%kpt%end
893 do ib = st%group%block_start, st%group%block_end
894 call batch_axpy(gr%np, deltat, st1%group%psib(ib, ik), st%group%psib(ib, ik))
899 do i = 1, te%exp_order
900 zfact = zfact * deltat / (i+1)
908 do ik = st%d%kpt%start, st%d%kpt%end
909 do ib = st%group%block_start, st%group%block_end
912 call batch_axpy(gr%np, deltat * zfact, st1%group%psib(ib, ik), st%group%psib(ib, ik))
923 if (
present(order)) order = te%exp_order*st%nik*st%nst
931 class(
mesh_t),
intent(in) :: mesh
933 class(
batch_t),
intent(inout) :: psib
934 real(real64),
intent(in) :: deltat
935 integer,
intent(in) :: k
936 real(real64),
optional,
intent(in) :: vmagnus(:,:,:)
939 complex(real64) :: deltat_
944 if (
present(vmagnus))
then
947 assert(
size(vmagnus, 1) >= mesh%np)
948 assert(
size(vmagnus, 2) == hm%d%nspin)
949 assert(
size(vmagnus, 3) == 2)
951 write(
message(1),
'(a)')
'Magnus operators only implemented for electrons at the moment'
956 if (.not. hm%is_hermitian() .and. te%exp_method ==
exp_chebyshev)
then
957 write(
message(1),
'(a)')
'The Chebyshev expansion for the exponential will only converge if the imaginary'
958 write(
message(2),
'(a)')
'eigenvalues are small enough compared to the span of the real eigenvalues,'
959 write(
message(3),
'(a)')
'i.e., for ratios smaller than about 1e-3.'
960 write(
message(4),
'(a)')
'The Lanczos method ("TDExponentialMethod = lanczos") is guaranteed to'
961 write(
message(5),
'(a)')
'always converge in this case.'
965 deltat_ = cmplx(deltat,
m_zero, real64)
967 select case (te%exp_method)
974 else if (k == 2)
then
977 write(
message(1),
'(a)')
'Lanczos expansion not implemented for phi_k, k > 2'
984 else if (k == 2)
then
987 write(
message(1),
'(a)')
'Chebyshev expansion not implemented for phi_k, k > 2'
991 deallocate(chebyshev_function)
batchified version of the BLAS axpy routine:
scale a batch by a constant or vector
subroutine, public accel_finish()
This module implements batches of mesh functions.
This module implements common operations on batches of mesh functions.
subroutine, public batch_set_zero(this, np, async)
fill all mesh functions of the batch with zero
This module contains interfaces for BLAS routines You should not use these routines directly....
subroutine, public exponential_copy(teo, tei)
subroutine, public exponential_apply_all(te, namespace, gr, hm, st, deltat, order)
Note that this routine not only computes the exponential, but also an extra term if there is a inhomo...
subroutine operate_batch(hm, namespace, mesh, psib, hpsib, vmagnus)
subroutine exponential_apply_phi_batch(te, namespace, mesh, hm, psib, deltat, k, vmagnus)
subroutine exponential_cheby_batch(te, namespace, mesh, hm, psib, chebyshev_function, vmagnus)
Calculates the exponential of the Hamiltonian through an expansion in Chebyshev polynomials.
subroutine, public exponential_init(te, namespace, full_batch)
integer, parameter, public exp_taylor
subroutine exponential_lanczos_batch(te, namespace, mesh, hm, psib, deltat, vmagnus, inh_psib)
subroutine exponential_apply_single(te, namespace, mesh, hm, zpsi, ist, ik, deltat, vmagnus, imag_time)
Wrapper to batchified routine for applying exponential to an array.
subroutine exponential_taylor_series_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, inh_psib, phik_shift)
subroutine exponential_apply_batch(te, namespace, mesh, hm, psib, deltat, psib2, deltat2, vmagnus, imag_time, inh_psib)
This routine performs the operation:
subroutine, public exponential_lanczos_function_batch(te, namespace, mesh, hm, psib, deltat, fun, vmagnus)
Compute fun(H) psib, i.e. the application of a function of the Hamiltonian to a batch.
integer, parameter, public exp_chebyshev
real(real64), parameter, public m_two
real(real64), parameter, public m_zero
complex(real64), parameter, public m_zi
real(real64), parameter, public m_epsilon
complex(real64), parameter, public m_z1
real(real64), parameter, public m_one
This module implements the underlying real-space grid.
This module defines an abstract class for Hamiltonians.
subroutine, public zhamiltonian_elec_apply_all(hm, namespace, gr, st, hst)
pure logical function, public hamiltonian_elec_inh_term(hm)
subroutine, public zlalg_matrix_function(n, factor, a, fun_a, fun, hermitian, tridiagonal)
This routine calculates a function of a matrix by using an eigenvalue decomposition.
This module is intended to contain "only mathematical" functions and procedures.
complex(real64) pure function, public phi2(z)
Compute phi2(z) = (phi1(z)-1)/z = (exp(z) - z - 1)/z^2.
complex(real64) pure function, public exponential(z)
Wrapper for exponential.
complex(real64) pure function, public phi1(z)
Compute phi1(z) = (exp(z)-1)/z.
This module defines functions over batches of mesh functions.
subroutine, public mesh_batch_nrm2(mesh, aa, nrm2, reduce)
Calculate the norms (norm2, not the square!) of a batch of mesh functions.
subroutine, public zmesh_batch_orthogonalization(mesh, nst, psib, phib, normalize, overlap, norm, gs_scheme, full_batch)
Orthonormalizes states of phib to the orbitals of nst batches of psi.
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_warning(no_lines, all_nodes, namespace)
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
logical function, public parse_is_defined(namespace, name)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
type(type_t), public type_cmplx
Class defining batches of mesh functions.
Description of the grid, containing information on derivatives, stencil, and symmetries.
The abstract Hamiltonian class defines a skeleton for specific implementations.
Describes mesh distribution to nodes.
The states_elec_t class contains all electronic wave functions.
batches of electronic states