Revisiting the evolution eigensolver.

This commit contains few bugfixes:
 - The Taylor expansion for imaginary time was not correct.
 - The code was wrongly mixing the density between timesteps of the evolution
 - The code was doing a subspace diagonalization which is not correct
 - The orthonormalization routine was not stable and is replaced by the standard orthogonalization.
 - The code was freezing states once converged, which is not correct unless independent particles

Besides this, the code now check that the exponential is stable. The default time-step is set to a valid
close to the one shown in the litterature.
parent be86811a
Pipeline #148087950 failed with stage
in 0 seconds
......@@ -408,6 +408,23 @@ contains
! now the eigensolver stuff
call eigensolver_init(scf%eigens, namespace, gr, st, geo, mc)
!The evolution operator is a very specific propagation that requires a specific
!setting to work in the current framework
if(scf%eigens%es_type == RS_EVO) then
if(scf%mix_field /= OPTION__MIXFIELD__DENSITY) then
message(1) = "Evolution eigensolver is only compatible with MixingField = density."
call messages_fatal(1)
end if
if(mix_coefficient(scf%smix) /= M_ONE) then
message(1) = "Evolution eigensolver is only compatible with Mixing = 1."
call messages_fatal(1)
end if
if(mix_scheme(scf%smix) /= OPTION__MIXINGSCHEME__LINEAR) then
message(1) = "Evolution eigensolver is only compatible with MixingScheme = linear."
call messages_fatal(1)
end if
end if
!%Variable SCFinLCAO
!%Type logical
!%Default no
......
......@@ -32,7 +32,7 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
integer :: ib, minst, maxst, ist, iter, maxiter, conv, convb, matvec, i
R_TYPE, allocatable :: c(:, :), zeig(:), res(:)
FLOAT, allocatable :: eig(:)
FLOAT, allocatable :: eig(:), norm(:)
type(wfs_elec_t) :: hpsib
#if defined(R_TREAL)
type(wfs_elec_t) :: zpsib
......@@ -49,6 +49,7 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
SAFE_ALLOCATE(eig(1:st%nst))
SAFE_ALLOCATE(zeig(1:st%nst))
SAFE_ALLOCATE(res(1:st%nst))
SAFE_ALLOCATE(norm(1:st%nst))
#if defined(R_TREAL)
SAFE_ALLOCATE(zpsi(1:mesh%np_part, 1:st%d%dim))
SAFE_ALLOCATE(psi(1:mesh%np_part, 1:st%d%dim))
......@@ -58,15 +59,21 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
! whose convergence should not be checked.
convb = st%group%block_start - 1
conv = 0
do ib = st%group%block_start, st%group%block_end
if (conv + st%group%psib(ib, ik)%nst <= converged) then
conv = conv + st%group%psib(ib, ik)%nst
convb = convb + 1
else
exit
end if
end do
!In case of independent particles, the Hamiltonian does not depends upon the density and we can freeze
!the lowest converged states
if(hm%theory_level == INDEPENDENT_PARTICLES) then
do ib = st%group%block_start, st%group%block_end
if (conv + st%group%psib(ib, ik)%nst <= converged) then
conv = conv + st%group%psib(ib, ik)%nst
convb = convb + 1
else
exit
end if
end do
end if
!I am leaving the loop as it is in case we want to change the logic of the code.
! However, we would need to compute the density and v_ks inside the loop if we do so.
do iter = 1, maxiter
#if defined(R_TREAL)
! The application of the exponential for the real case is still done one state at a time, as we need
......@@ -83,12 +90,21 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
call batch_get_state(zpsib, 1, mesh%np, zpsi)
psi(1:mesh%np, 1:st%d%dim) = R_TOTYPE(zpsi(1:mesh%np, 1:st%d%dim))
norm(ist) = X(mf_nrm2)(mesh, st%d%dim, psi)
if(norm(ist) > CNST(1.1)) then
message(1) = "Evolution eigensolver: the time evolution seems to be unstable."
message(2) = "Please reduce the value of EigensolverImaginaryTime."
call messages_fatal(2)
end if
call states_elec_set_state(st, mesh, ist, ik, psi)
call zpsib%end()
matvec = matvec + te%exp_order
end do
#else
do ib = convb + 1, st%group%block_end
minst = states_elec_block_min(st, ib)
maxst = states_elec_block_max(st, ib)
if (hamiltonian_elec_apply_packed(hm)) then
call st%group%psib(ib, ik)%do_pack()
end if
......@@ -98,24 +114,25 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
call hamiltonian_elec_base_unset_phase_corr(hm%hm_base, mesh, st%group%psib(ib, ik))
matvec = matvec + te%exp_order*(states_elec_block_max(st, ib) - states_elec_block_min(st, ib) + 1)
call mesh_batch_nrm2(mesh, st%group%psib(ib, ik), norm(minst:maxst))
if(any(norm(minst:maxst) > CNST(1.1))) then
message(1) = "Evolution eigensolver: the time evolution seems to be unstable."
message(2) = "Please reduce the value of EigensolverImaginaryTime."
call messages_fatal(2)
end if
if (hamiltonian_elec_apply_packed(hm)) then
call st%group%psib(ib, ik)%do_unpack()
end if
end do
#endif
! This is the orthonormalization suggested by Aichinger and Krotschek
! This is not the orthonormalization suggested by Aichinger and Krotschek
! [Comp. Mat. Science 34, 188 (2005)]
call X(states_elec_calc_overlap)(st, mesh, ik, c)
! However, the current scheme seem to be more robbust (Cholesky).
! Otherwise the norms after orthonormalization are not exactly 1, which create a problem
call X(states_elec_orthogonalization_full)(st, namespace, mesh, ik)
call lalg_eigensolve(st%nst, c, eig)
do i = 1, st%nst
c(1:st%nst, i) = c(1:st%nst, i)/sqrt(eig(i))
end do
call states_elec_rotate(st, namespace, mesh, c, ik)
! Get the eigenvalues and the residues.
do ib = convb + 1, st%group%block_end
if (hamiltonian_elec_apply_packed(hm)) then
......@@ -130,7 +147,7 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
call X(hamiltonian_elec_apply_batch)(hm, namespace, mesh, st%group%psib(ib, ik), hpsib)
call X(mesh_batch_dotp_vector)(mesh, st%group%psib(ib, ik), hpsib, zeig(minst:maxst))
st%eigenval(minst:maxst, ik) = R_REAL(zeig(minst:maxst))
call batch_axpy(mesh%np, -st%eigenval(minst:maxst, ik), st%group%psib(ib, ik), hpsib)
call batch_axpy(mesh%np, -st%eigenval(:,ik), st%group%psib(ib, ik), hpsib)
call mesh_batch_nrm2(mesh, hpsib, diff(minst:maxst))
call hpsib%end()
......@@ -148,15 +165,17 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
end do
! And check for convergence. Note that they must be converged *in order*, so that they can be frozen.
do ib = convb + 1, st%group%block_end
minst = states_elec_block_min(st, ib)
maxst = states_elec_block_max(st, ib)
if (all(diff(minst:maxst) < tol) .and. ib == convb + 1) then
convb = convb + 1
conv = conv + maxst - minst + 1
end if
end do
if (convb == st%group%block_end) exit
if(hm%theory_level == INDEPENDENT_PARTICLES) then
do ib = convb + 1, st%group%block_end
minst = states_elec_block_min(st, ib)
maxst = states_elec_block_max(st, ib)
if (all(diff(minst:maxst) < tol) .and. ib == convb + 1) then
convb = convb + 1
conv = conv + maxst - minst + 1
end if
end do
if (convb == st%group%block_end) exit
end if
end do
converged = conv
......@@ -169,6 +188,7 @@ subroutine X(eigensolver_evolution)(namespace, mesh, st, hm, te, tol, niter, con
#endif
SAFE_DEALLOCATE_A(c)
SAFE_DEALLOCATE_A(eig)
SAFE_DEALLOCATE_A(norm)
SAFE_DEALLOCATE_A(zeig)
SAFE_DEALLOCATE_A(res)
......
......@@ -264,14 +264,15 @@ contains
!%Variable EigensolverImaginaryTime
!%Type float
!%Default 10.0
!%Default 0.1
!%Section SCF::Eigensolver
!%Description
!% The imaginary-time step that is used in the imaginary-time evolution
!% method (<tt>Eigensolver = evolution</tt>) to obtain the lowest eigenvalues/eigenvectors.
!% It must satisfy <tt>EigensolverImaginaryTime > 0</tt>.
!% Increasing this value can make the propagation faster, but could lead to unstable propagations.
!%End
call parse_variable(namespace, 'EigensolverImaginaryTime', CNST(10.0), eigens%imag_time)
call parse_variable(namespace, 'EigensolverImaginaryTime', CNST(0.1), eigens%imag_time)
if(eigens%imag_time <= M_ZERO) call messages_input_error(namespace, 'EigensolverImaginaryTime')
call exponential_init(eigens%exponential_operator, namespace)
......@@ -337,6 +338,7 @@ contains
!% except for <tt>rmdiis</tt>, which performs only 3 iterations.
!% Increasing this value for <tt>rmdiis</tt> increases the convergence speed,
!% at the cost of an increased memory footprint.
!% In the case of imaginary time propatation, this variable is not used.
!%End
call parse_variable(namespace, 'EigensolverMaxIter', default_iter, eigens%es_maxiter)
if(eigens%es_maxiter < 1) call messages_input_error(namespace, 'EigensolverMaxIter')
......@@ -367,7 +369,14 @@ contains
! FEAST: subspace diagonalization or not? I guess not.
! But perhaps something could be gained by changing this.
call subspace_init(eigens%sdiag, namespace, st, no_sd = .false.)
!
! In case of the evolution eigensolver, this makes no sense to use subspace diagonalization
! as orthogonalization is done internally at each time-step
if(eigens%es_type == RS_EVO) then
call subspace_init(eigens%sdiag, namespace, st, no_sd = .true.)
else
call subspace_init(eigens%sdiag, namespace, st, no_sd = .false.)
end if
! print memory requirements
select case(eigens%es_type)
......@@ -504,6 +513,7 @@ contains
call deigensolver_plan(namespace, gr, st, hm, eigens%pre, eigens%tolerance, maxiter, eigens%converged(ik), ik, &
eigens%diff(:, ik))
case(RS_EVO)
maxiter = 1
call deigensolver_evolution(namespace, gr%mesh, st, hm, eigens%exponential_operator, eigens%tolerance, maxiter, &
eigens%converged(ik), ik, eigens%diff(:, ik), tau = eigens%imag_time)
case(RS_LOBPCG)
......@@ -550,6 +560,7 @@ contains
call zeigensolver_plan(namespace, gr, st, hm, eigens%pre, eigens%tolerance, maxiter, &
eigens%converged(ik), ik, eigens%diff(:, ik))
case(RS_EVO)
maxiter = 1
call zeigensolver_evolution(namespace, gr%mesh, st, hm, eigens%exponential_operator, eigens%tolerance, maxiter, &
eigens%converged(ik), ik, eigens%diff(:, ik), tau = eigens%imag_time)
case(RS_LOBPCG)
......
......@@ -343,7 +343,7 @@ contains
SAFE_ALLOCATE(hzpsi1(1:mesh%np, 1:hm%d%dim))
zfact = M_z1
zfact_is_real = .true.
zfact_is_real = abs(deltat-real(deltat)) < M_EPSILON
do idim = 1, hm%d%dim
call lalg_copy(mesh%np, zpsi(:, idim), zpsi1(:, idim))
......@@ -711,7 +711,7 @@ contains
zfact = M_z1
zfact2 = M_z1
zfact_is_real = .true.
zfact_is_real = abs(deltat-real(deltat)) < M_EPSILON
if(present(psib2)) call psib%copy_data_to(mesh%np, psib2)
......
......@@ -94,6 +94,7 @@ contains
!% Selects the method to perform subspace diagonalization. The
!% default is <tt>standard</tt>, unless states parallelization is used,
!% when the default is <tt>scalapack</tt>.
!% Note that this variable is not parsed in the case of the evolution eigensolver.
!%Option none 0
!% No subspace diagonalization. WARNING: this will generally give incorrect results.
!%Option standard 1
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment