Octopus
eigensolver.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22 use batch_oct_m
24 use debug_oct_m
32 use global_oct_m
33 use grid_oct_m
37 use, intrinsic :: iso_fortran_env
41 use loct_oct_m
42 use mesh_oct_m
46 use mpi_oct_m
50 use parser_oct_m
54 use smear_oct_m
55 use space_oct_m
62 use types_oct_m
63 use unit_oct_m
66 use xc_oct_m
67
68 implicit none
69
70 private
71 public :: &
75
76 type eigensolver_t
77 private
78 integer, public :: es_type
79
80 real(real64), public :: tolerance
81 integer, public :: es_maxiter
82
83 real(real64) :: tau
84 real(real64) :: tau0
85 real(real64) :: time
86 integer :: it_propagator
87 type(potential_interpolation_t) :: vks_old
88 real(real64), allocatable :: normalization_energies(:, :)
89 ! !! evolution eigensolver
90 real(real64), allocatable :: normalization_energies_prev(:, :)
92 logical :: variable_timestep
93
95 real(real64), allocatable, public :: diff(:, :)
96 integer, public :: matvec
97 integer, allocatable, public :: converged(:)
98
100 type(preconditioner_t), public :: pre
101
103 type(subspace_t) :: sdiag
104
105 integer :: rmmdiis_minimization_iter
106
107 logical, public :: folded_spectrum
108
109 ! cg options
110 logical, public :: orthogonalize_to_all
111 integer, public :: conjugate_direction
112 real(real64), public :: energy_change_threshold
113
114 ! Chebyshev filtering options
115 type(eigen_chebyshev_t), public :: cheby_params
117 type(exponential_t) :: exponential_operator
118 contains
119 procedure :: run => eigensolver_run
120 end type eigensolver_t
121
122
123 integer, public, parameter :: &
124 RS_PLAN = 11, &
125 rs_cg = 5, &
126 rs_evo = 9, &
127 rs_rmmdiis = 10, &
128 rs_chebyshev = 12
129
130contains
131
132 ! ---------------------------------------------------------
133 subroutine eigensolver_init(eigens, namespace, gr, st, hm, mc, space, deactivate_oracle)
134 type(eigensolver_t), intent(out) :: eigens
135 type(namespace_t), intent(in) :: namespace
136 type(grid_t), intent(in) :: gr
137 type(states_elec_t), intent(in) :: st
138 type(hamiltonian_elec_t), intent(in) :: hm
139 type(multicomm_t), intent(in) :: mc
140 class(space_t), intent(in) :: space
141 logical, optional, intent(in) :: deactivate_oracle
142
143 integer :: default_iter, default_es
144 real(real64) :: default_tol
145 real(real64) :: mem
146
147 push_sub(eigensolver_init)
148
149 !%Variable Eigensolver
150 !%Type integer
151 !%Section SCF::Eigensolver
152 !%Description
153 !% Which eigensolver to use to obtain the lowest eigenvalues and
154 !% eigenfunctions of the Kohn-Sham Hamiltonian. The default is
155 !% conjugate gradients (<tt>cg</tt>), except that when parallelization in states is
156 !% enabled, the default is <tt>chebyshev_filter</tt>.
157 !%Option cg 5
158 !% Conjugate-gradients algorithm.
159 !%Option plan 11
160 !% Preconditioned Lanczos scheme. Ref: Y. Saad, A. Stathopoulos, J. Chelikowsky, K. Wu and S. Ogut,
161 !% "Solution of Large Eigenvalue Problems in Electronic Structure Calculations", <i>BIT</i> <b>36</b>, 1 (1996).
162 !%Option evolution 9
163 !% (Experimental) Propagation in imaginary time.
164 !% This eigensolver uses a time propagation in imaginary time (see, e.g., Aichinger, M. & Krotscheck, E.,
165 !% Computational Materials Science 34, 188–212 (2005), 10.1016/j.commatsci.2004.11.002).
166 !% The timestep is set using <tt>EigensolverImaginaryTime</tt>. The convergence typically depends
167 !% on the total imaginary time in a propagation which means that it should be faster with larger
168 !% timesteps. However, simulations typically diverge above a certain timestep that depends on the system.
169 !% In general, it is advisable to use an exponential method that always converges, especially for large
170 !% timesteps, i.e., set <tt>TDExponentialMethod</tt> to <tt>lanczos</tt> or <tt>chebyshev</tt>.
171 !% The propagation method can be set with <tt>ImaginaryTimePropagator</tt> and variable timestepping
172 !% can be enabled with <tt>ImaginaryTimeVariableTimestep</tt>.
173 !% Moreover, this method usually converges quite slowly, so you might need to increase <tt>MaximumIter</tt>.
174 !% It is incompatible with smearing for the occupations.
175 !%Option rmmdiis 10
176 !% Residual minimization scheme, direct inversion in the
177 !% iterative subspace eigensolver, based on the implementation of
178 !% Kresse and Furthm&uuml;ller [<i>Phys. Rev. B</i> <b>54</b>, 11169
179 !% (1996)]. This eigensolver requires almost no orthogonalization
180 !% so it can be considerably faster than the other options for
181 !% large systems. To improve its performance a large number of <tt>ExtraStates</tt>
182 !% are required (around 10-20% of the number of occupied states).
183 !% Note: with <tt>unocc</tt>, you will need to stop the calculation
184 !% by hand, since the highest states will probably never converge.
185 !% Usage with more than one block of states per node is experimental, unfortunately.
186 !%Option chebyshev_filter 12
187 !% Chebyshev polynomials are used to construct a spectral filter that is iteratively applied to a trial subspace,
188 !% amplifying a subset of the lowest eigenvalues of the Hamiltonian, which effectively isolates the invariant
189 !% subspace associated with the these states. The filtered subspace is projected to define a reduced eigenvalue problem,
190 !% which is then solved using dense diagonalisation. Finally, the eigenstates of the original Hamiltonian are recovered
191 !% via subspace rotation. For further details, see [Zhou et. al.](http://dx.doi.org/10.1016/j.jcp.2014.06.056)
192 !%End
193
194 if(st%parallel_in_states) then
195 default_es = rs_chebyshev
196 else
197 default_es = rs_cg
198 end if
199
200 call parse_variable(namespace, 'Eigensolver', default_es, eigens%es_type)
201
202 if(st%parallel_in_states .and. .not. eigensolver_parallel_in_states(eigens)) then
203 message(1) = "The selected eigensolver is not parallel in states."
204 message(2) = "Please use the rmmdiis or Chebyshev filtering eigensolvers."
205 call messages_fatal(2, namespace=namespace)
206 end if
208 call messages_obsolete_variable(namespace, 'EigensolverVerbose')
209 call messages_obsolete_variable(namespace, 'EigensolverSubspaceDiag', 'SubspaceDiagonalization')
211 default_iter = 25
213 select case(eigens%es_type)
214 case(rs_cg)
215 !%Variable CGOrthogonalizeAll
216 !%Type logical
217 !%Default no
218 !%Section SCF::Eigensolver
219 !%Description
220 !% Used by the cg solver only.
221 !% During the cg iterations, the current band can be orthogonalized
222 !% against all other bands or only against the lower bands. Orthogonalizing
223 !% against all other bands can improve convergence properties, whereas
224 !% orthogonalizing against lower bands needs less operations.
225 !% Moreover, orthogonalizing against all bands can make converging
226 !% the highest band or unoccupied bands more difficult.
227 !%End
228 call parse_variable(namespace, 'CGOrthogonalizeAll', .false., eigens%orthogonalize_to_all)
229
230 !%Variable CGDirection
231 !%Type integer
232 !%Section SCF::Eigensolver
233 !%Description
234 !% Used by the cg solver only.
235 !% The conjugate direction is updated using a certain coefficient to the previous
236 !% direction. This coeffiction can be computed in different ways. The default is
237 !% to use Fletcher-Reeves (FR), an alternative is Polak-Ribiere (PR).
238 !%Option fletcher 1
239 !% The coefficient for Fletcher-Reeves consists of the current norm of the
240 !% steepest descent vector divided by that of the previous iteration.
241 !%Option polak 2
242 !% For the Polak-Ribiere scheme, a product of the current with the previous
243 !% steepest descent vector is subtracted in the nominator.
244 !%End
245 call parse_variable(namespace, 'CGDirection', option__cgdirection__fletcher, eigens%conjugate_direction)
246
247 !%Variable CGEnergyChangeThreshold
248 !%Type float
249 !%Section SCF::Eigensolver
250 !%Default 0.1
251 !%Description
252 !% Used by the cg solver only.
253 !% For each band, the CG iterations are stopped when the change in energy is smaller than the
254 !% change in the first iteration multiplied by this factor. This limits the number of CG
255 !% iterations for each band, while still showing good convergence for the SCF cycle. The criterion
256 !% is discussed in Sec. V.B.6 of Payne et al. (1992), Rev. Mod. Phys. 64, 4.
257 !% The default value is 0.1, which is usually a good choice for LDA and GGA potentials. If you
258 !% are solving the OEP equation, you might want to set this value to 1e-3 or smaller. In general,
259 !% smaller values might help if you experience convergence problems.
260 !% For very small convergence tolerances, choose 0 to disable this criterion.
261 !%End
262 call parse_variable(namespace, 'CGEnergyChangeThreshold', 0.1_real64, eigens%energy_change_threshold)
263
264 case(rs_plan)
265 case(rs_evo)
266 call messages_experimental("imaginary-time evolution eigensolver")
267
268 !%Variable EigensolverImaginaryTime
269 !%Type float
270 !%Default 0.1
271 !%Section SCF::Eigensolver
272 !%Description
273 !% The imaginary-time step that is used in the imaginary-time evolution
274 !% method (<tt>Eigensolver = evolution</tt>) to obtain the lowest eigenvalues/eigenvectors.
275 !% It must satisfy <tt>EigensolverImaginaryTime > 0</tt>.
276 !% Increasing this value can make the propagation faster, but could lead to unstable propagations.
277 !%End
278 call parse_variable(namespace, 'EigensolverImaginaryTime', 0.1_real64, eigens%tau)
279 if(eigens%tau <= m_zero) call messages_input_error(namespace, 'EigensolverImaginaryTime')
280 ! save initial timestep
281 eigens%tau0 = eigens%tau
282
283 call exponential_init(eigens%exponential_operator, namespace)
284
285 if(st%smear%method /= smear_semiconductor .and. st%smear%method /= smear_fixed_occ) then
286 message(1) = "Smearing of occupations is incompatible with imaginary time evolution."
287 call messages_fatal(1, namespace=namespace)
288 end if
289
290 !%Variable ImaginaryTimePropagator
291 !%Type integer
292 !%Section SCF::Eigensolver
293 !%Default it_forward_euler
294 !%Description
295 !% The propagation method for imaginary-time evolution.
296 !%Option it_forward_euler 1
297 !% Approximate the evolution operator by one forward Euler step with the exponential:
298 !% <math>
299 !% U_{\rm FE}(t+\delta t, t) = \exp \left( -\tau H_{t}\right)\,.
300 !% </math>
301 !% This scheme is first order.
302 !%Option it_expmid 2
303 !% Similar to forward Euler, but evaluate the Hamiltonian at the half-timestep by
304 !% extrapolation:
305 !% <math>
306 !% U_{\rm EM}(t+\delta t, t) = \exp \left( -\tau H_{t+\delta t/2}\right)\,.
307 !% </math>
308 !% This scheme is second order.
309 !%End
310 call parse_variable(namespace, 'ImaginaryTimePropagator', it_forward_euler, eigens%it_propagator)
311 ! initialize potential interpolation
312 if (eigens%it_propagator == it_expmid) then
313 call hm%ks_pot%init_interpolation(eigens%vks_old)
314 call hm%ks_pot%run_zero_iter(eigens%vks_old)
315 end if
316
317 !%Variable ImaginaryTimeVariableTimestep
318 !%Type logical
319 !%Default no
320 !%Section SCF::Eigensolver
321 !%Description
322 !% Enable variable timesteps for the imaginary time propagation.
323 !%End
324 call parse_variable(namespace, 'ImaginaryTimeVariableTimestep', .false., eigens%variable_timestep)
325
326 eigens%time = m_zero
327
328 safe_allocate(eigens%normalization_energies(1:st%nst, 1:st%nik))
329 safe_allocate(eigens%normalization_energies_prev(1:st%nst, 1:st%nik))
330
331 case(rs_rmmdiis)
332 default_iter = 5
333
334 !%Variable EigensolverMinimizationIter
335 !%Type integer
336 !%Default 0
337 !%Section SCF::Eigensolver
338 !%Description
339 !% During the first iterations, the RMMDIIS eigensolver requires
340 !% some steepest-descent minimizations to improve
341 !% convergence. This variable determines the number of those
342 !% minimizations.
343 !%End
344
345 call parse_variable(namespace, 'EigensolverMinimizationIter', 0, eigens%rmmdiis_minimization_iter)
346
347 if(gr%use_curvilinear) call messages_experimental("RMMDIIS eigensolver for curvilinear coordinates")
348
349 case(rs_chebyshev)
350 !%Variable ChebyshevFilterLanczosOrder
351 !%Type integer
352 !%Default -1
353 !%Section SCF::Eigensolver
354 !%Description
355 !% Used by the Chebyshev filter only.
356 !% The number of Lanczos iterations used to construct the tridiagonal matrix,
357 !% from which the upper bound of H is estimated.
358 !% A value in the range 4 <= ChebyshevFilterLanczosOrder <= 10 is reasonable.
359 !% Values greater than 30 will raise an assertion.
360 !%
361 !% The default value, -1, indicate that the order is selected automatically.
362 !% In this case, the trial vector is also reused to make the procedure faster.
363 !%End
364 call parse_variable(namespace, 'ChebyshevFilterLanczosOrder', default_chebyshev_params%n_lanczos, &
365 eigens%cheby_params%n_lanczos)
366
367 ! TODO Rewrite this, having changed the behaviour
368 !%Variable ChebyshevFilterDegree
369 !%Type integer
370 !%Default 25
371 !%Section SCF::Eigensolver
372 !%Description
373 !% Used by the Chebyshev filter only.
374 !% The degree of the Chebyshev polynomial used to dampen the interval of eigenstates one is not interested in.
375 !% If used in conjunction with "OptimizeChebyshevFilterDegree", which is the default, "ChebyshevFilterDegree" defines the
376 !% the maximum Chebyshev polynomial degree Octopus will consider when determining an approximate, optimal degree.
377 !%
378 !% If not used in conjunction with "OptimizeChebyshevFilterDegree", this defines the polynomial degree used for all SCF steps.
379 !% The larger the degree, the larger the separation between the subspaces, which will reduce the number of SCF iterations
380 !% required to reach convergence. However, ChebyshevFilterDegree also directly correlates with the number of matrix-vector
381 !% products performed on the Hamiltonian per SCF step, and so larger values result in slower SCF cycles.
382 !% A value in the range 8 <= ChebyshevFilterDegree <= 20 is reasonable.
383 !%End
384 call parse_variable(namespace, 'ChebyshevFilterDegree', default_chebyshev_params%degree, eigens%cheby_params%degree)
385
386 !%Variable OptimizeChebyshevFilterDegree
387 !%Type logical
388 !%Default yes
389 !%Section SCF::Eigensolver
390 !%Description
391 !% Used by the Chebyshev filter only.
392 !% Octopus generates a best-estimate for the degree of the Chebyshev polynomial used to filter the subspace.
393 !% A separate estimate is generated for each state block, per SCF iteration. Note that if running parallelism
394 !% over states, the block/batch containing the largest eigenstates will converge more slowly and will
395 !% therefore use a larger degree relative to all other batches. One should therefore avoid setting "ChebyshevFilterDegree"
396 !% to an excessive value (for example > 50). For more details regarding how the degree is estimated, one can
397 !% refer to Section 5.4 in "Computer Physics Communications 187 (2015) 98–105"
398 !% (http://dx.doi.org/10.1016/j.cpc.2014.10.015).
399 !%
400 !% This is deactivated by default for spinors as this is not found to be advantageous in test systems.
401 !%End
402 if (st%d%ispin == spinors .or. optional_default(deactivate_oracle, .false.)) then
403 default_chebyshev_params%optimize_degree = .false.
404 end if
405 call parse_variable(namespace, 'OptimizeChebyshevFilterDegree', default_chebyshev_params%optimize_degree, &
406 eigens%cheby_params%optimize_degree)
407
408 !%Variable ChebyshevFilterBoundMixing
409 !%Type float
410 !%Default 0.5
411 !%Section SCF::Eigensolver
412 !%Description
413 !% In the first application of the filter, for the first SCF step, the initial lower bound estimate
414 !% is defined as a linear combination of the smallest and largest eigenvalues of the Hamiltonian.
415 !% The bound mixing determines the proportion of linear mixing, beta:
416 !% $b_{lower} = beta min{\lambda} + (\beta - 1) max{\lambda}$
417 !% of the smallest and largest eigenvalues.
418 !%End
419 call parse_variable(namespace, 'ChebyshevFilterBoundMixing', default_chebyshev_params%bound_mixing, &
420 eigens%cheby_params%bound_mixing)
421
422 !%Variable ChebyshevFilterNIter
423 !%Type integer
424 !%Default 5
425 !%Section SCF::Eigensolver
426 !%Description
427 !% The max number of iterations in the first SCF step of the Chebyshev solver. In practice this
428 !% does not need to exceed 10, as the eigenstates will vary alot between the first and second
429 !% SCF steps.
430 !%End
431 call parse_variable(namespace, 'ChebyshevFilterNIter', default_chebyshev_params%n_iter, &
432 eigens%cheby_params%n_iter)
433
434 case default
435 call messages_input_error(namespace, 'Eigensolver')
436 end select
437
438 call messages_print_with_emphasis(msg='Eigensolver', namespace=namespace)
439
440 call messages_print_var_option("Eigensolver", eigens%es_type, namespace=namespace)
441
442 call messages_obsolete_variable(namespace, 'EigensolverInitTolerance', 'EigensolverTolerance')
443 call messages_obsolete_variable(namespace, 'EigensolverFinalTolerance', 'EigensolverTolerance')
444 call messages_obsolete_variable(namespace, 'EigensolverFinalToleranceIteration')
445
446 ! this is an internal option that makes the solver use the
447 ! folded operator (H-shift)^2 to converge first eigenvalues around
448 ! the values of shift
449 ! c.f. L. W. Wang and A. Zunger
450 ! JCP 100, 2394 (1994); doi: http://dx.doi.org/10.1063/1.466486
451 eigens%folded_spectrum = .false.
452
453 !%Variable EigensolverTolerance
454 !%Type float
455 !%Section SCF::Eigensolver
456 !%Description
457 !% This is the tolerance for the eigenvectors. The default is 1e-7.
458 !% However, if <tt>ConvRelDens</tt> is defined, the default is set to a tenth of its value.
459 !%End
460 default_tol = 1e-7_real64
461 if (parse_is_defined(namespace, 'ConvRelDens')) then
462 call parse_variable(namespace, 'ConvRelDens', 1e-6_real64, default_tol)
463 default_tol = default_tol / 10.0_real64
464 end if
465 call parse_variable(namespace, 'EigensolverTolerance', default_tol, eigens%tolerance)
466
467 !%Variable EigensolverMaxIter
468 !%Type integer
469 !%Section SCF::Eigensolver
470 !%Description
471 !% Determines the maximum number of iterations that the
472 !% eigensolver will perform if the desired tolerance is not
473 !% achieved. The default is 25 iterations for all eigensolvers
474 !% except for <tt>rmdiis</tt>, which performs only 5 iterations.
475 !% Increasing this value for <tt>rmdiis</tt> increases the convergence speed,
476 !% at the cost of an increased memory footprint.
477 !%End
478 call parse_variable(namespace, 'EigensolverMaxIter', default_iter, eigens%es_maxiter)
479 if(eigens%es_maxiter < 1) call messages_input_error(namespace, 'EigensolverMaxIter')
480
481 if(eigens%es_maxiter > default_iter) then
482 call messages_write('You have specified a large number of eigensolver iterations (')
483 call messages_write(eigens%es_maxiter)
484 call messages_write(').', new_line = .true.)
485 call messages_write('This is not a good idea as it might slow down convergence, even for', new_line = .true.)
486 call messages_write('independent particles, as subspace diagonalization will not be used', new_line = .true.)
487 call messages_write('often enough.')
488 call messages_warning(namespace=namespace)
489 end if
490
491 if (any(eigens%es_type == (/rs_plan, rs_cg, rs_rmmdiis/))) then
492 call preconditioner_init(eigens%pre, namespace, gr, mc, space)
493 end if
494
495 safe_allocate(eigens%diff(1:st%nst, 1:st%nik))
496 eigens%diff(1:st%nst, 1:st%nik) = 0
497
498 safe_allocate(eigens%converged(1:st%nik))
499 eigens%converged(1:st%nik) = 0
500 eigens%matvec = 0
501
502 call eigens%sdiag%init(namespace, st)
503
504 ! print memory requirements
505 select case(eigens%es_type)
506 case(rs_rmmdiis)
507 call messages_write('Info: The rmmdiis eigensolver requires ')
508 mem = (m_two*eigens%es_maxiter - m_one)*st%block_size*real(gr%np_part, real64)
509 if(states_are_real(st)) then
510 mem = mem*8.0_real64
511 else
512 mem = mem*16.0_real64
513 end if
514 call messages_write(mem, units = unit_megabytes, fmt = '(f9.1)')
515 call messages_write(' of additional')
516 call messages_new_line()
517 call messages_write(' memory. This amount can be reduced by decreasing the value')
518 call messages_new_line()
519 call messages_write(' of the variable StatesBlockSize (currently set to ')
520 call messages_write(st%block_size)
521 call messages_write(').')
522 call messages_info()
523 end select
524
525 call messages_print_with_emphasis(namespace=namespace)
526
527 pop_sub(eigensolver_init)
528 end subroutine eigensolver_init
529
530
531 ! ---------------------------------------------------------
532 subroutine eigensolver_end(eigens)
533 type(eigensolver_t), intent(inout) :: eigens
534
535 push_sub(eigensolver_end)
536
537 select case(eigens%es_type)
538 case(rs_plan, rs_cg, rs_rmmdiis)
539 call preconditioner_end(eigens%pre)
540 end select
541
542 safe_deallocate_a(eigens%converged)
543 safe_deallocate_a(eigens%diff)
544
545 safe_deallocate_a(eigens%normalization_energies)
546 safe_deallocate_a(eigens%normalization_energies_prev)
547
548 pop_sub(eigensolver_end)
549 end subroutine eigensolver_end
550
551
552 ! ---------------------------------------------------------
553 subroutine eigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter, conv, nstconv)
554 class(eigensolver_t), intent(inout) :: eigens
555 type(namespace_t), intent(in) :: namespace
556 type(grid_t), intent(in) :: gr
557 type(states_elec_t), intent(inout) :: st
558 type(hamiltonian_elec_t), intent(inout) :: hm
559 class(space_t), intent(in) :: space
560 type(partner_list_t), intent(in) :: ext_partners
561 integer, intent(in) :: iter
562 logical, optional, intent(out) :: conv
563 integer, optional, intent(in) :: nstconv
564 ! !< the convergence criteria
565
566 integer :: nstconv_
567#ifdef HAVE_MPI
568 logical :: conv_reduced
569 integer :: ist, outcount, lmatvec
570 real(real64), allocatable :: ldiff(:), leigenval(:)
571 real(real64), allocatable :: ldiff_out(:), leigenval_out(:)
572 integer, allocatable :: lconv(:)
573#endif
574
575 call profiling_in("EIGEN_SOLVER")
576 push_sub(eigensolver_run)
577
578 if(present(conv)) conv = .false.
579 if(present(nstconv)) then
580 nstconv_ = nstconv
581 else
582 nstconv_ = st%nst
583 end if
584
585 eigens%matvec = 0
586
587 if(mpi_world%is_root() .and. eigensolver_has_progress_bar(eigens) .and. .not. debug%info) then
588 call loct_progress_bar(-1, st%lnst*st%d%kpt%nlocal)
589 end if
590
591 if (states_are_real(st)) then
592 call deigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter)
593 else
594 call zeigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter)
595 end if
596
597 if(mpi_world%is_root() .and. eigensolver_has_progress_bar(eigens) .and. .not. debug%info) then
598 write(stdout, '(1x)')
599 end if
600
601 if(present(conv)) conv = all(eigens%converged(st%d%kpt%start:st%d%kpt%end) >= nstconv_)
602
603#ifdef HAVE_MPI
604 if (st%d%kpt%parallel) then
605 if (present(conv)) then
606 call st%d%kpt%mpi_grp%allreduce(conv, conv_reduced, 1, mpi_logical, mpi_land)
607 conv = conv_reduced
608 end if
609
610 lmatvec = eigens%matvec
611 call st%d%kpt%mpi_grp%allreduce(lmatvec, eigens%matvec, 1, mpi_integer, mpi_sum)
612
613 safe_allocate(lconv(1:st%d%kpt%nlocal))
614 lconv(1:st%d%kpt%nlocal) = eigens%converged(st%d%kpt%start:st%d%kpt%end)
615 call lmpi_gen_allgatherv(st%d%kpt%nlocal, lconv, outcount, eigens%converged, st%d%kpt%mpi_grp)
616 assert(outcount == st%nik)
617 safe_deallocate_a(lconv)
618
619 ! every node needs to know all eigenvalues (and diff)
620 safe_allocate(ldiff(1:st%d%kpt%nlocal))
621 safe_allocate(leigenval(1:st%d%kpt%nlocal))
622 safe_allocate(ldiff_out(1:st%nik))
623 safe_allocate(leigenval_out(1:st%nik))
624 do ist = st%st_start, st%st_end
625 ldiff(1:st%d%kpt%nlocal) = eigens%diff(ist, st%d%kpt%start:st%d%kpt%end)
626 leigenval(1:st%d%kpt%nlocal) = st%eigenval(ist, st%d%kpt%start:st%d%kpt%end)
627 call lmpi_gen_allgatherv(st%d%kpt%nlocal, ldiff, outcount, ldiff_out, st%d%kpt%mpi_grp)
628 eigens%diff(ist, :) = ldiff_out
629 assert(outcount == st%nik)
630 call lmpi_gen_allgatherv(st%d%kpt%nlocal, leigenval, outcount, leigenval_out, st%d%kpt%mpi_grp)
631 st%eigenval(ist, :) = leigenval_out
632 assert(outcount == st%nik)
633 end do
634 safe_deallocate_a(ldiff)
635 safe_deallocate_a(ldiff_out)
636 safe_deallocate_a(leigenval)
637 safe_deallocate_a(leigenval_out)
638 end if
639#endif
640
641 pop_sub(eigensolver_run)
642 call profiling_out("EIGEN_SOLVER")
643 end subroutine eigensolver_run
644
645
646 ! ---------------------------------------------------------
647 logical function eigensolver_parallel_in_states(this) result(par_stat)
648 type(eigensolver_t), intent(in) :: this
649
651
652 par_stat = .false.
653
654 select case(this%es_type)
656 par_stat = .true.
657 end select
658
661
662
663 ! ---------------------------------------------------------
664 logical function eigensolver_has_progress_bar(this) result(has)
665 type(eigensolver_t), intent(in) :: this
666
668
669 has = .false.
670
671 select case(this%es_type)
672 case(rs_rmmdiis, rs_cg)
673 has = .true.
674 end select
675
678
679
680#include "undef.F90"
681#include "real.F90"
682#include "eigensolver_inc.F90"
683#include "eigen_plan_inc.F90"
684
685#include "undef.F90"
686#include "complex.F90"
687#include "eigensolver_inc.F90"
688#include "eigen_plan_inc.F90"
689
690end module eigensolver_oct_m
691
692
693!! Local Variables:
694!! mode: f90
695!! coding: utf-8
696!! End:
This module implements batches of mesh functions.
Definition: batch.F90:135
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:118
type(debug_t), save, public debug
Definition: debug.F90:158
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
type(eigen_chebyshev_t), public default_chebyshev_params
Default Chebyshev input parameters Arguments 1 and 2 taken from 10.1016/j.jcp.2006....
integer, parameter, public it_expmid
integer, parameter, public it_forward_euler
subroutine zeigensolver_run(eigens, namespace, mesh, st, hm, space, ext_partners, iter)
integer, parameter, public rs_evo
subroutine, public eigensolver_init(eigens, namespace, gr, st, hm, mc, space, deactivate_oracle)
integer, parameter, public rs_cg
subroutine deigensolver_run(eigens, namespace, mesh, st, hm, space, ext_partners, iter)
subroutine eigensolver_run(eigens, namespace, gr, st, hm, space, ext_partners, iter, conv, nstconv)
logical function eigensolver_parallel_in_states(this)
integer, parameter, public rs_chebyshev
logical function eigensolver_has_progress_bar(this)
subroutine, public eigensolver_end(eigens)
integer, parameter, public rs_rmmdiis
integer, parameter, public spinors
subroutine, public exponential_init(te, namespace, full_batch)
real(real64), parameter, public m_two
Definition: global.F90:192
real(real64), parameter, public m_zero
Definition: global.F90:190
real(real64), parameter, public m_one
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:119
This module defines classes and functions for interaction partners.
A module to handle KS potential, without the external potential.
This module defines functions over batches of mesh functions.
Definition: mesh_batch.F90:118
This module defines various routines, operating on mesh functions.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:120
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:904
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:531
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1029
subroutine, public messages_new_line()
Definition: messages.F90:1118
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:416
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:697
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1069
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
This module contains some common usage patterns of MPI routines.
Definition: mpi_lib.F90:117
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:272
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:147
logical function, public parse_is_defined(namespace, name)
Definition: parser.F90:505
subroutine, public preconditioner_end(this)
subroutine, public preconditioner_init(this, namespace, gr, mc, space)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:625
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
integer, parameter, public smear_semiconductor
Definition: smear.F90:173
integer, parameter, public smear_fixed_occ
Definition: smear.F90:173
pure logical function, public states_are_real(st)
This module handles spin dimensions of the states and the k-point distribution.
This module provides routines for communicating states when using states parallelization.
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:134
This module defines the unit system, used for input and output.
type(unit_t), public unit_megabytes
For large amounts of data (natural code units are bytes)
Definition: xc.F90:116
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:171
The states_elec_t class contains all electronic wave functions.
int true(void)