Octopus
scf.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2014 M. Marques, A. Castro, A. Rubio, G. Bertsch, M. Oliveira
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
21module scf_oct_m
23 use berry_oct_m
26 use debug_oct_m
34 use forces_oct_m
35 use global_oct_m
36 use grid_oct_m
39 use io_oct_m
40 use ions_oct_m
41 use, intrinsic :: iso_fortran_env
44 use lcao_oct_m
45 use lda_u_oct_m
49 use loct_oct_m
51 use math_oct_m
52 use mesh_oct_m
55 use mix_oct_m
57 use mpi_oct_m
60 use output_oct_m
63 use parser_oct_m
67 use smear_oct_m
68 use space_oct_m
74 use stress_oct_m
76 use types_oct_m
77 use unit_oct_m
79 use utils_oct_m
80 use v_ks_oct_m
82 use vdw_ts_oct_m
86 use xc_oct_m
87 use xc_f03_lib_m
88 use xc_functional_oct_m
90 use xc_oep_oct_m
92
93 implicit none
94
95 private
96 public :: &
97 scf_t, &
98 scf_init, &
100 scf_load, &
101 scf_start, &
102 scf_run, &
103 scf_iter, &
105 scf_finish, &
106 scf_end, &
109
110 integer, public, parameter :: &
111 VERB_NO = 0, &
112 verb_compact = 1, &
113 verb_full = 3
116 type scf_t
117 private
118 integer, public :: max_iter
119
120 real(real64), public :: lmm_r
121
122 ! several convergence criteria
123 logical :: conv_eigen_error
124 logical :: check_conv
125
126 integer :: mix_field
127 logical :: lcao_restricted
128 logical :: calc_force
129 logical, public :: calc_stress
130 logical :: calc_dipole
131 logical :: calc_partial_charges
132 type(mix_t) :: smix
133 type(mixfield_t), pointer :: mixfield
134 type(eigensolver_t) :: eigens
135 integer :: mixdim1
136 logical :: forced_finish
137 type(lda_u_mixer_t) :: lda_u_mix
138 type(vtau_mixer_t) :: vtau_mix
139 type(berry_t) :: berry
140 integer :: matvec
141
142 type(restart_t), public :: restart_load, restart_dump
143
144 type(criterion_list_t), public :: criterion_list
145 real(real64) :: energy_in, energy_diff, abs_dens_diff, evsum_in, evsum_out, evsum_diff
146
147 ! Variables needed to store information accross scf_start, scf_run, and scf_finish
148 logical :: converged_current, converged_last
149 integer :: verbosity_
150 type(lcao_t) :: lcao
151 real(real64), allocatable :: rhoout(:,:), rhoin(:,:)
152 real(real64), allocatable :: vhxc_old(:,:)
153 class(wfs_elec_t), allocatable :: psioutb(:, :)
154 logical :: output_forces, calc_current, output_during_scf, finish
155 end type scf_t
156
157contains
158
159 ! ---------------------------------------------------------
160 subroutine scf_init(scf, namespace, gr, ions, st, mc, hm, space)
161 type(scf_t), intent(inout) :: scf
162 type(grid_t), intent(in) :: gr
163 type(namespace_t), intent(in) :: namespace
164 type(ions_t), intent(in) :: ions
165 type(states_elec_t), intent(in) :: st
166 type(multicomm_t), intent(in) :: mc
167 type(hamiltonian_elec_t), intent(inout) :: hm
168 class(space_t), intent(in) :: space
169
170 real(real64) :: rmin
171 integer :: mixdefault
172 type(type_t) :: mix_type
173 class(convergence_criterion_t), pointer :: crit
174 type(criterion_iterator_t) :: iter
175
176 push_sub(scf_init)
177
178 !%Variable MaximumIter
179 !%Type integer
180 !%Default 200
181 !%Section SCF::Convergence
182 !%Description
183 !% Maximum number of SCF iterations. The code will stop even if convergence
184 !% has not been achieved. -1 means unlimited.
185 !% 0 means just do LCAO (or read from restart), compute the eigenvalues and energy,
186 !% and stop, without updating the wavefunctions or density.
187 !%
188 !% If convergence criteria are set, the SCF loop will only stop once the criteria
189 !% are fulfilled for two consecutive iterations.
190 !%
191 !% Note that this variable is also used in the section Calculation Modes::Unoccupied States,
192 !% where it denotes the maximum number of calls of the eigensolver. In this context, the
193 !% default value is 50.
194 !%End
195 call parse_variable(namespace, 'MaximumIter', 200, scf%max_iter)
196
197 if (allocated(hm%vberry)) then
198 call berry_init(scf%berry, namespace)
199 end if
200
201 !Create the list of convergence criteria
202 call criteria_factory_init(scf%criterion_list, namespace, scf%check_conv)
203 !Setting the pointers
204 call iter%start(scf%criterion_list)
205 do while (iter%has_next())
206 crit => iter%get_next()
207 select type (crit)
208 type is (energy_criterion_t)
209 call crit%set_pointers(scf%energy_diff, scf%energy_in)
210 type is (density_criterion_t)
211 call crit%set_pointers(scf%abs_dens_diff, st%qtot)
212 type is (eigenval_criterion_t)
213 call crit%set_pointers(scf%evsum_diff, scf%evsum_out)
214 class default
215 assert(.false.)
216 end select
217 end do
218
220 if(.not. scf%check_conv .and. scf%max_iter < 0) then
221 call messages_write("All convergence criteria are disabled. Octopus is cowardly refusing")
223 call messages_write("to enter an infinite loop.")
226 call messages_write("Please set one of the following variables to a positive value:")
229 call messages_write(" | MaximumIter | ConvEnergy | ConvAbsDens | ConvRelDens |")
231 call messages_write(" | ConvAbsEv | ConvRelEv |")
233 call messages_fatal(namespace=namespace)
234 end if
236 !%Variable ConvEigenError
237 !%Type logical
238 !%Default false
239 !%Section SCF::Convergence
240 !%Description
241 !% If true, the calculation will not be considered converged unless all states have
242 !% individual errors less than <tt>EigensolverTolerance</tt>.
243 !%
244 !% If this criterion is used, the SCF loop will only stop once it is
245 !% fulfilled for two consecutive iterations.
246 !%End
247 call parse_variable(namespace, 'ConvEigenError', .false., scf%conv_eigen_error)
248
249 if(scf%max_iter < 0) scf%max_iter = huge(scf%max_iter)
250
251 call messages_obsolete_variable(namespace, 'What2Mix', 'MixField')
252
253 ! now the eigensolver stuff
254 call eigensolver_init(scf%eigens, namespace, gr, st, mc, space)
255
256 if(scf%eigens%es_type /= rs_evo) then
257 !%Variable MixField
258 !%Type integer
259 !%Section SCF::Mixing
260 !%Description
261 !% Selects what should be mixed during the SCF cycle. Note that
262 !% currently the exact-exchange part of hybrid functionals is not
263 !% mixed at all, which would require wavefunction-mixing, not yet
264 !% implemented. This may lead to instabilities in the SCF cycle,
265 !% so starting from a converged LDA/GGA calculation is recommended
266 !% for hybrid functionals. The default depends on the <tt>TheoryLevel</tt>
267 !% and the exchange-correlation potential used.
268 !% This is not used in case of imaginary-time evolution.
269 !%Option none 0
270 !% No mixing is done. This is the default for independent
271 !% particles.
272 !%Option potential 1
273 !% The Kohn-Sham potential is mixed. This is the default for other cases.
274 !%Option density 2
275 !% Mix the density.
276 !%Option states 3
277 !% (Experimental) Mix the states. In this case, the mixing is always linear.
278 !%End
279
280 mixdefault = option__mixfield__potential
281 if(hm%theory_level == independent_particles) mixdefault = option__mixfield__none
282
283 call parse_variable(namespace, 'MixField', mixdefault, scf%mix_field)
284 if(.not.varinfo_valid_option('MixField', scf%mix_field)) call messages_input_error(namespace, 'MixField')
285 call messages_print_var_option('MixField', scf%mix_field, "what to mix during SCF cycles", namespace=namespace)
286
287 if (scf%mix_field == option__mixfield__potential .and. hm%theory_level == independent_particles) then
288 call messages_write('Input: Cannot mix the potential for non-interacting particles.')
289 call messages_fatal(namespace=namespace)
290 end if
291
292 if (scf%mix_field == option__mixfield__potential .and. hm%pcm%run_pcm) then
293 call messages_write('Input: You have selected to mix the potential.', new_line = .true.)
294 call messages_write(' This might produce convergence problems for solvated systems.', new_line = .true.)
295 call messages_write(' Mix the Density instead.')
296 call messages_warning(namespace=namespace)
297 end if
298
299 if(scf%mix_field == option__mixfield__density &
300 .and. bitand(hm%xc%family, xc_family_oep + xc_family_mgga + xc_family_hyb_mgga + xc_family_nc_mgga) /= 0) then
301
302 call messages_write('Input: You have selected to mix the density with OEP or MGGA XC functionals.', new_line = .true.)
303 call messages_write(' This might produce convergence problems. Mix the potential instead.')
304 call messages_warning(namespace=namespace)
305 end if
306
307 if(scf%mix_field == option__mixfield__states) then
308 call messages_experimental('MixField = states', namespace=namespace)
309 end if
310
311 ! Handle mixing now...
312 select case(scf%mix_field)
313 case (option__mixfield__potential, option__mixfield__density)
314 scf%mixdim1 = gr%np
315 case(option__mixfield__states)
316 ! we do not really need the mixer, except for the value of the mixing coefficient
317 scf%mixdim1 = 1
318 end select
319
320 mix_type = type_float
321
322 if (scf%mix_field /= option__mixfield__none) then
323 call mix_init(scf%smix, namespace, space, gr%der, scf%mixdim1, st%d%nspin, func_type_ = mix_type)
324 end if
325
326 ! If we use DFT+U, we also have do mix it
327 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none ) then
328 call lda_u_mixer_init(hm%lda_u, scf%lda_u_mix, st)
329 call lda_u_mixer_init_auxmixer(hm%lda_u, namespace, scf%lda_u_mix, scf%smix, st)
330 end if
331
332 ! If we use tau-dependent MGGA, we need to mix vtau
333 if(scf%mix_field == option__mixfield__potential) then
334 call vtau_mixer_init_auxmixer(namespace, scf%vtau_mix, scf%smix, hm, gr%np, st%d%nspin)
335 end if
336
337 call mix_get_field(scf%smix, scf%mixfield)
338 else
339 scf%mix_field = option__mixfield__none
340 end if
341
342 !%Variable SCFinLCAO
343 !%Type logical
344 !%Default no
345 !%Section SCF
346 !%Description
347 !% Performs the SCF cycle with the calculation restricted to the LCAO subspace.
348 !% This may be useful for systems with convergence problems (first do a
349 !% calculation within the LCAO subspace, then restart from that point for
350 !% an unrestricted calculation).
351 !%End
352 call parse_variable(namespace, 'SCFinLCAO', .false., scf%lcao_restricted)
353 if(scf%lcao_restricted) then
354 call messages_experimental('SCFinLCAO', namespace=namespace)
355 message(1) = 'Info: SCF restricted to LCAO subspace.'
356 call messages_info(1, namespace=namespace)
357
358 if(scf%conv_eigen_error) then
359 message(1) = "ConvEigenError cannot be used with SCFinLCAO, since error is unknown."
360 call messages_fatal(1, namespace=namespace)
361 end if
362 end if
363
364
365 !%Variable SCFCalculateForces
366 !%Type logical
367 !%Section SCF
368 !%Description
369 !% This variable controls whether the forces on the ions are
370 !% calculated at the end of a self-consistent iteration. The
371 !% default is yes, unless the system only has user-defined
372 !% species.
373 !%End
374 call parse_variable(namespace, 'SCFCalculateForces', .not. ions%only_user_def, scf%calc_force)
375
376 if(scf%calc_force .and. gr%der%boundaries%spiralBC) then
377 message(1) = 'Forces cannot be calculated when using spiral boundary conditions.'
378 write(message(2),'(a)') 'Please use SCFCalculateForces = no.'
379 call messages_fatal(2, namespace=namespace)
380 end if
381 if(scf%calc_force) then
382 if (allocated(hm%ep%b_field) .or. allocated(hm%ep%a_static)) then
383 write(message(1),'(a)') 'The forces are currently not properly calculated if static'
384 write(message(2),'(a)') 'magnetic fields or static vector potentials are present.'
385 write(message(3),'(a)') 'Please use SCFCalculateForces = no.'
386 call messages_fatal(3, namespace=namespace)
387 end if
388 end if
389
390 !%Variable SCFCalculateStress
391 !%Type logical
392 !%Section SCF
393 !%Description
394 !% This variable controls whether the stress on the lattice is
395 !% calculated at the end of a self-consistent iteration. The
396 !% default is no.
397 !%End
398 call parse_variable(namespace, 'SCFCalculateStress', .false. , scf%calc_stress)
399
400 !%Variable SCFCalculateDipole
401 !%Type logical
402 !%Section SCF
403 !%Description
404 !% This variable controls whether the dipole is calculated at the
405 !% end of a self-consistent iteration. For finite systems the
406 !% default is yes. For periodic systems the default is no, unless
407 !% an electric field is being applied in a periodic direction.
408 !% The single-point Berry`s phase approximation is used for
409 !% periodic directions. Ref:
410 !% E Yaschenko, L Fu, L Resca, and R Resta, <i>Phys. Rev. B</i> <b>58</b>, 1222-1229 (1998).
411 !%End
412 call parse_variable(namespace, 'SCFCalculateDipole', .not. space%is_periodic(), scf%calc_dipole)
413 if (allocated(hm%vberry)) scf%calc_dipole = .true.
414
415 !%Variable SCFCalculatePartialCharges
416 !%Type logical
417 !%Default no
418 !%Section SCF
419 !%Description
420 !% (Experimental) This variable controls whether partial charges
421 !% are calculated at the end of a self-consistent iteration.
422 !%End
423 call parse_variable(namespace, 'SCFCalculatePartialCharges', .false., scf%calc_partial_charges)
424 if (scf%calc_partial_charges) call messages_experimental('SCFCalculatePartialCharges', namespace=namespace)
425
426 rmin = ions%min_distance()
427
428 !%Variable LocalMagneticMomentsSphereRadius
429 !%Type float
430 !%Section Output
431 !%Description
432 !% The local magnetic moments are calculated by integrating the
433 !% magnetization density in spheres centered around each atom.
434 !% This variable controls the radius of the spheres.
435 !% The default is half the minimum distance between two atoms
436 !% in the input coordinates, or 100 a.u. if there is only one atom (for isolated systems).
437 !%End
438 call parse_variable(namespace, 'LocalMagneticMomentsSphereRadius', min(m_half*rmin, lmm_r_single_atom), scf%lmm_r, &
439 unit=units_inp%length)
440 ! this variable is also used in td/td_write.F90
441
442 scf%forced_finish = .false.
443
444 pop_sub(scf_init)
445 end subroutine scf_init
446
447
448 ! ---------------------------------------------------------
449 subroutine scf_end(scf)
450 type(scf_t), intent(inout) :: scf
451
452 class(convergence_criterion_t), pointer :: crit
453 type(criterion_iterator_t) :: iter
454
455 push_sub(scf_end)
456
457 call eigensolver_end(scf%eigens)
458
459 if(scf%mix_field /= option__mixfield__none) call mix_end(scf%smix)
460
461 nullify(scf%mixfield)
462
463 if (scf%mix_field /= option__mixfield__states) then
464 call lda_u_mixer_end(scf%lda_u_mix, scf%smix)
465 call vtau_mixer_end(scf%vtau_mix, scf%smix)
466 end if
467
468 call iter%start(scf%criterion_list)
469 do while (iter%has_next())
470 crit => iter%get_next()
471 safe_deallocate_p(crit)
472 end do
473
474 pop_sub(scf_end)
475 end subroutine scf_end
476
477
478 ! ---------------------------------------------------------
479 subroutine scf_mix_clear(scf)
480 type(scf_t), intent(inout) :: scf
481
482 push_sub(scf_mix_clear)
483
484 call mix_clear(scf%smix)
485
486 if (scf%mix_field /= option__mixfield__states) then
487 call lda_u_mixer_clear(scf%lda_u_mix, scf%smix)
488 call vtau_mixer_clear(scf%vtau_mix, scf%smix)
489 end if
490
491 pop_sub(scf_mix_clear)
492 end subroutine scf_mix_clear
493
494 ! ---------------------------------------------------------
496 subroutine scf_load(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, restart_load)
497 type(scf_t), intent(inout) :: scf
498 type(namespace_t), intent(in) :: namespace
499 type(electron_space_t), intent(in) :: space
500 type(multicomm_t), intent(in) :: mc
501 type(grid_t), intent(inout) :: gr
502 type(ions_t), intent(in) :: ions
503 type(partner_list_t), intent(in) :: ext_partners
504 type(states_elec_t), intent(inout) :: st
505 type(v_ks_t), intent(inout) :: ks
506 type(hamiltonian_elec_t), intent(inout) :: hm
507 type(restart_t), intent(in) :: restart_load
508
509 integer :: ierr, is
510
511 push_sub(scf_load)
512
513 if (restart_has_flag(restart_load, restart_flag_rho)) then
514 ! Load density and used it to recalculated the KS potential.
515 call states_elec_load_rho(restart_load, space, st, gr, ierr)
516 if (ierr /= 0) then
517 message(1) = 'Unable to read density. Density will be calculated from states.'
518 call messages_warning(1, namespace=namespace)
519 else
520 if (bitand(ks%xc_family, xc_family_oep) == 0) then
521 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
522 else
523 if (.not. restart_has_flag(restart_load, restart_flag_vhxc) .and. ks%oep%level /= oep_level_full) then
524 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
525 end if
526 end if
527 end if
528 end if
529
530 if (restart_has_flag(restart_load, restart_flag_vhxc)) then
531 call hm%ks_pot%load(restart_load, space, gr, ierr)
532 if (ierr /= 0) then
533 message(1) = 'Unable to read Vhxc. Vhxc will be calculated from states.'
534 call messages_warning(1, namespace=namespace)
535 else
536 call hm%update(gr, namespace, space, ext_partners)
537 if (bitand(ks%xc_family, xc_family_oep) /= 0) then
538 if (ks%oep%level == oep_level_full) then
539 do is = 1, st%d%nspin
540 ks%oep%vxc(1:gr%np, is) = hm%ks_pot%vhxc(1:gr%np, is) - hm%ks_pot%vhartree(1:gr%np)
541 end do
542 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
543 end if
544 end if
545 end if
546 end if
547
548 if (restart_has_flag(restart_load, restart_flag_mix)) then
549 if (scf%mix_field == option__mixfield__density .or. scf%mix_field == option__mixfield__potential) then
550 call mix_load(namespace, restart_load, scf%smix, space, gr, ierr)
551 end if
552 if (ierr /= 0) then
553 message(1) = "Unable to read mixing information. Mixing will start from scratch."
554 call messages_warning(1, namespace=namespace)
555 end if
556 end if
557
558 if(hm%lda_u_level /= dft_u_none) then
559 call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr)
560 if (ierr /= 0) then
561 message(1) = "Unable to read DFT+U information. DFT+U data will be calculated from states."
562 call messages_warning(1, namespace=namespace)
563 end if
564 end if
565
566 !TODO: Create a dedicated routine and call it from the initialize
567
568! if (present(outp) .and. mpi_grp_is_root(mpi_world)) then
569! call io_rm(STATIC_DIR //"info")
570! end if
571! end if
573 pop_sub(scf_load)
574 end subroutine scf_load
575
576 ! ---------------------------------------------------------
578 subroutine scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
579 type(scf_t), intent(inout) :: scf
580 type(namespace_t), intent(in) :: namespace
581 type(electron_space_t), intent(in) :: space
582 type(grid_t), intent(inout) :: gr
583 type(ions_t), intent(inout) :: ions
584 type(states_elec_t), intent(inout) :: st
585 type(v_ks_t), intent(inout) :: ks
586 type(hamiltonian_elec_t), intent(inout) :: hm
587 type(output_t), optional, intent(in) :: outp
588 integer, optional, intent(in) :: verbosity
590 integer :: ib, iqn
591
592 push_sub(scf_start)
593
594 if(scf%forced_finish) then
595 message(1) = "Previous clean stop, not doing SCF and quitting."
596 call messages_fatal(1, only_root_writes = .true., namespace=namespace)
597 end if
598
599 scf%verbosity_ = optional_default(verbosity, verb_full)
600
601 scf%output_during_scf = .false.
602 scf%output_forces = .false.
603 scf%calc_current = .false.
604
605 if (present(outp)) then
606 ! if the user has activated output=stress but not SCFCalculateStress,
607 ! we assume that is implied
608 if (outp%what(option__output__stress)) then
609 scf%calc_stress = .true.
610 end if
611
612 scf%output_during_scf = outp%duringscf
613 scf%calc_current = output_needs_current(outp, states_are_real(st))
614
615 if (outp%duringscf .and. outp%what(option__output__forces)) then
616 scf%output_forces = .true.
617 end if
618 end if
619
620 if(scf%lcao_restricted) then
621 call lcao_init(scf%lcao, namespace, space, gr, ions, st)
622 if(.not. lcao_is_available(scf%lcao)) then
623 message(1) = 'LCAO is not available. Cannot do SCF in LCAO.'
624 call messages_fatal(1, namespace=namespace)
625 end if
626 end if
627
628 safe_allocate(scf%rhoout(1:gr%np, 1:st%d%nspin))
629 safe_allocate(scf%rhoin (1:gr%np, 1:st%d%nspin))
630
631 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
632 scf%rhoout = m_zero
633
634 if (scf%calc_force .or. scf%output_forces) then
635 !We store the Hxc potential for the contribution to the forces
636 safe_allocate(scf%vhxc_old(1:gr%np, 1:st%d%nspin))
637 call lalg_copy(gr%np, st%d%nspin, hm%ks_pot%vhxc, scf%vhxc_old)
638 end if
639
640
641 select case(scf%mix_field)
642 case(option__mixfield__potential)
643 call mixfield_set_vin(scf%mixfield, hm%ks_pot%vhxc)
644 call vtau_mixer_set_vin(scf%vtau_mix, hm)
645 case(option__mixfield__density)
646 call mixfield_set_vin(scf%mixfield, scf%rhoin)
647
648 case(option__mixfield__states)
649
650 ! There is a ICE with foss2022a-serial. I am changing to allocate - NTD
651 allocate(wfs_elec_t::scf%psioutb (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
652
653 do iqn = st%d%kpt%start, st%d%kpt%end
654 do ib = st%group%block_start, st%group%block_end
655 call st%group%psib(ib, iqn)%copy_to(scf%psioutb(ib, iqn))
656 end do
657 end do
658
659 end select
660
661 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
662 ! If we use DFT+U, we also have do mix it
663 if (scf%mix_field /= option__mixfield__states) call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
664
665 call create_convergence_file(static_dir, "convergence")
666
667 if ( scf%verbosity_ /= verb_no ) then
668 if(scf%max_iter > 0) then
669 write(message(1),'(a)') 'Info: Starting SCF iteration.'
670 else
671 write(message(1),'(a)') 'Info: No SCF iterations will be done.'
672 ! we cannot tell whether it is converged.
673 scf%finish = .false.
674 end if
675 call messages_info(1, namespace=namespace)
676 end if
677
678 scf%converged_current = .false.
679 scf%matvec = 0
680
681 pop_sub(scf_start)
682
683 contains
684
685 ! -----------------------------------------------------
686
687 subroutine create_convergence_file(dir, fname)
688 character(len=*), intent(in) :: dir
689 character(len=*), intent(in) :: fname
690
691 integer :: iunit
692 character(len=12) :: label
693 if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes
694 call io_mkdir(dir, namespace)
695 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
696 write(iunit, '(a)', advance = 'no') '#iter energy '
697 label = 'energy_diff'
698 write(iunit, '(1x,a)', advance = 'no') label
699 label = 'abs_dens'
700 write(iunit, '(1x,a)', advance = 'no') label
701 label = 'rel_dens'
702 write(iunit, '(1x,a)', advance = 'no') label
703 label = 'abs_ev'
704 write(iunit, '(1x,a)', advance = 'no') label
705 label = 'rel_ev'
706 write(iunit, '(1x,a)', advance = 'no') label
707 if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
708 .and. ks%theory_level /= generalized_kohn_sham_dft) then
709 if (ks%oep%level == oep_level_full) then
710 label = 'OEP norm2ss'
711 write(iunit, '(1x,a)', advance = 'no') label
712 end if
713 end if
714 write(iunit,'(a)') ''
715 call io_close(iunit)
716 end if
717
718 end subroutine create_convergence_file
719
720 end subroutine scf_start
721
722 ! ---------------------------------------------------------
724 subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
725 verbosity, iters_done, restart_dump)
726 type(scf_t), intent(inout) :: scf
727 type(namespace_t), intent(in) :: namespace
728 type(electron_space_t), intent(in) :: space
729 type(multicomm_t), intent(in) :: mc
730 type(grid_t), intent(inout) :: gr
731 type(ions_t), intent(inout) :: ions
732 type(partner_list_t), intent(in) :: ext_partners
733 type(states_elec_t), intent(inout) :: st
734 type(v_ks_t), intent(inout) :: ks
735 type(hamiltonian_elec_t), intent(inout) :: hm
736 type(output_t), optional, intent(in) :: outp
737 integer, optional, intent(in) :: verbosity
738 integer, optional, intent(out) :: iters_done
739 type(restart_t), optional, intent(in) :: restart_dump
740
741 integer :: iter
742 logical :: completed
743
744 push_sub(scf_run)
745
746 call scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
747
748 ! SCF cycle
749 do iter = 1, scf%max_iter
750
751 call scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
752 verbosity, iters_done, restart_dump)
753
754 completed = scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
755
756 if(scf%forced_finish .or. completed) then
757 exit
758 end if
759 end do
760
761 call scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
762 verbosity, iters_done, restart_dump)
763
764 pop_sub(scf_run)
765 end subroutine scf_run
766
767 ! ---------------------------------------------------------
768 subroutine scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
769 verbosity, iters_done, restart_dump)
770 type(scf_t), intent(inout) :: scf
771 type(namespace_t), intent(in) :: namespace
772 type(electron_space_t), intent(in) :: space
773 type(multicomm_t), intent(in) :: mc
774 type(grid_t), intent(inout) :: gr
775 type(ions_t), intent(inout) :: ions
776 type(partner_list_t), intent(in) :: ext_partners
777 type(states_elec_t), intent(inout) :: st
778 type(v_ks_t), intent(inout) :: ks
779 type(hamiltonian_elec_t), intent(inout) :: hm
780 integer, intent(in) :: iter
781 type(output_t), optional, intent(in) :: outp
782 integer, optional, intent(in) :: verbosity
783 integer, optional, intent(out) :: iters_done
784 type(restart_t), optional, intent(in) :: restart_dump
785
786 integer :: iqn, ib, ierr
787 class(convergence_criterion_t), pointer :: crit
788 type(criterion_iterator_t) :: iterator
789 logical :: is_crit_conv
790 real(real64) :: etime, itime
791
792 push_sub(scf_iter)
793
794 call profiling_in("SCF_CYCLE")
795
796 itime = loct_clock()
797
798 ! this initialization seems redundant but avoids improper optimization at -O3 by PGI 7 on chum,
799 ! which would cause a failure of testsuite/linear_response/04-vib_modes.03-vib_modes_fd.inp
800 scf%eigens%converged = 0
801
802 ! keep the information about the spectrum up to date, needed e.g. for Chebyshev expansion for imaginary time
803 call hm%update_span(gr%spacing(1:space%dim), minval(st%eigenval(:, :)), namespace)
804
805 !We update the quantities at the begining of the scf cycle
806 if (iter == 1) then
807 scf%evsum_in = states_elec_eigenvalues_sum(st)
808 end if
809 call iterator%start(scf%criterion_list)
810 do while (iterator%has_next())
811 crit => iterator%get_next()
812 call scf_update_initial_quantity(scf, hm, crit)
813 end do
814
815 if (scf%calc_force .or. scf%output_forces) then
816 !Used for computing the imperfect convegence contribution to the forces
817 scf%vhxc_old(1:gr%np, 1:st%d%nspin) = hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin)
818 end if
819
820 if(scf%lcao_restricted) then
821 call lcao_init_orbitals(scf%lcao, namespace, st, gr, ions)
822 call lcao_wf(scf%lcao, st, gr, ions, hm, namespace)
823 else
824
825 !We check if the system is coupled with a partner that requires self-consistency
826 ! if(hamiltonian_has_scf_partner(hm)) then
827 if (allocated(hm%vberry)) then
828 !In this case, v_Hxc is frozen and we do an internal SCF loop over the
829 ! partners that require SCF
830 ks%frozen_hxc = .true.
831 ! call perform_scf_partners()
832 call berry_perform_internal_scf(scf%berry, namespace, space, scf%eigens, gr, st, hm, iter, ks, ions, ext_partners)
833 !and we unfreeze the potential once finished
834 ks%frozen_hxc = .false.
835 else
836 scf%eigens%converged = 0
837 call scf%eigens%run(namespace, gr, st, hm, iter)
838 end if
839 end if
840
841 scf%matvec = scf%matvec + scf%eigens%matvec
842
843 ! occupations
844 call states_elec_fermi(st, namespace, gr)
845 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
846
847 ! compute output density, potential (if needed) and eigenvalues sum
848 call density_calc(st, gr, st%rho)
849
850 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoout)
851
852 select case (scf%mix_field)
853 case (option__mixfield__potential)
854 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
855 call mixfield_set_vout(scf%mixfield, hm%ks_pot%vhxc)
856 call vtau_mixer_set_vout(scf%vtau_mix, hm)
857 case (option__mixfield__density)
858 call mixfield_set_vout(scf%mixfield, scf%rhoout)
859 case(option__mixfield__states)
860
861 do iqn = st%d%kpt%start, st%d%kpt%end
862 do ib = st%group%block_start, st%group%block_end
863 call st%group%psib(ib, iqn)%copy_data_to(gr%np, scf%psioutb(ib, iqn))
864 end do
865 end do
866 end select
867
868 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
869 call lda_u_mixer_set_vout(hm%lda_u, scf%lda_u_mix)
870 endif
871
872 ! recalculate total energy
873 call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit = -1)
874
875 if (present(outp)) then
876 ! compute forces only if requested
877 if (outp%duringscf .and. outp%what_now(option__output__forces, iter)) then
878 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
879 end if
880 end if
881
882 !We update the quantities at the end of the scf cycle
883 call iterator%start(scf%criterion_list)
884 do while (iterator%has_next())
885 crit => iterator%get_next()
886 call scf_update_diff_quantity(scf, hm, st, gr, scf%rhoout, scf%rhoin, crit)
887 end do
888
889 ! are we finished?
890 scf%converged_last = scf%converged_current
891
892 scf%converged_current = scf%check_conv .and. &
893 (.not. scf%conv_eigen_error .or. all(scf%eigens%converged == st%nst))
894 !Loop over the different criteria
895 call iterator%start(scf%criterion_list)
896 do while (iterator%has_next())
897 crit => iterator%get_next()
898 call crit%is_converged(is_crit_conv)
899 scf%converged_current = scf%converged_current .and. is_crit_conv
900 end do
901
902 ! only finish if the convergence criteria are fulfilled in two
903 ! consecutive iterations
904 scf%finish = scf%converged_last .and. scf%converged_current
905
906 etime = loct_clock() - itime
907 call scf_write_iter(namespace)
908
909 ! mixing
910 select case (scf%mix_field)
911 case (option__mixfield__density)
912 ! mix input and output densities and compute new potential
913 call mixing(namespace, scf%smix)
914 call mixfield_get_vnew(scf%mixfield, st%rho)
915 ! for spinors, having components 3 or 4 be negative is not unphysical
916 if (minval(st%rho(1:gr%np, 1:st%d%spin_channels)) < -1e-6_real64) then
917 write(message(1),*) 'Negative density after mixing. Minimum value = ', &
918 minval(st%rho(1:gr%np, 1:st%d%spin_channels))
919 call messages_warning(1, namespace=namespace)
920 end if
921 call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
922 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
923 case (option__mixfield__potential)
924 ! mix input and output potentials
925 call mixing(namespace, scf%smix)
926 call mixfield_get_vnew(scf%mixfield, hm%ks_pot%vhxc)
927 call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
928 call vtau_mixer_get_vnew(scf%vtau_mix, hm)
929 call hamiltonian_elec_update_pot(hm, gr)
930
931 case(option__mixfield__states)
932
933 do iqn = st%d%kpt%start, st%d%kpt%end
934 do ib = st%group%block_start, st%group%block_end
935 call batch_scal(gr%np, m_one - mix_coefficient(scf%smix), st%group%psib(ib, iqn))
936 call batch_axpy(gr%np, mix_coefficient(scf%smix), scf%psioutb(ib, iqn), st%group%psib(ib, iqn))
937 end do
938 end do
939
940 call density_calc(st, gr, st%rho)
941 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
942
943 case (option__mixfield__none)
944 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
945 end select
946
947
948 ! Are we asked to stop? (Whenever Fortran is ready for signals, this should go away)
949 scf%forced_finish = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
950
951 if (scf%finish .and. st%modelmbparticles%nparticle > 0) then
952 call modelmb_sym_all_states(space, gr, st)
953 end if
954
955 if (present(outp) .and. present(restart_dump)) then
956 ! save restart information
957
958 if ( (scf%finish .or. (modulo(iter, outp%restart_write_interval) == 0) &
959 .or. iter == scf%max_iter .or. scf%forced_finish) ) then
960
961 call states_elec_dump(scf%restart_dump, space, st, gr, hm%kpoints, ierr, iter=iter)
962 if (ierr /= 0) then
963 message(1) = 'Unable to write states wavefunctions.'
964 call messages_warning(1, namespace=namespace)
965 end if
966
967 call states_elec_dump_rho(scf%restart_dump, space, st, gr, ierr, iter=iter)
968 if (ierr /= 0) then
969 message(1) = 'Unable to write density.'
970 call messages_warning(1, namespace=namespace)
971 end if
972
973 if(hm%lda_u_level /= dft_u_none) then
974 call lda_u_dump(scf%restart_dump, namespace, hm%lda_u, st, gr, ierr)
975 if (ierr /= 0) then
976 message(1) = 'Unable to write DFT+U information.'
977 call messages_warning(1, namespace=namespace)
978 end if
979 end if
980
981 select case (scf%mix_field)
982 case (option__mixfield__density)
983 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
984 if (ierr /= 0) then
985 message(1) = 'Unable to write mixing information.'
986 call messages_warning(1, namespace=namespace)
987 end if
988 case (option__mixfield__potential)
989 call hm%ks_pot%dump(scf%restart_dump, space, gr, ierr)
990 if (ierr /= 0) then
991 message(1) = 'Unable to write Vhxc.'
992 call messages_warning(1, namespace=namespace)
993 end if
994
995 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
996 if (ierr /= 0) then
997 message(1) = 'Unable to write mixing information.'
998 call messages_warning(1, namespace=namespace)
999 end if
1000 end select
1001 end if
1002 end if
1003
1004 call write_convergence_file(static_dir, "convergence")
1005
1006 call profiling_out("SCF_CYCLE")
1007
1008 pop_sub(scf_iter)
1009 contains
1010
1011 ! ---------------------------------------------------------
1012 subroutine scf_write_iter(namespace)
1013 type(namespace_t), intent(in) :: namespace
1014
1015 character(len=50) :: str
1016 real(real64) :: dipole(1:space%dim)
1017
1018 push_sub(scf_run.scf_write_iter)
1019
1020 if ( scf%verbosity_ == verb_full ) then
1021
1022 write(str, '(a,i5)') 'SCF CYCLE ITER #' ,iter
1023 call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
1024 write(message(1),'(a,es15.8,2(a,es9.2))') ' etot = ', units_from_atomic(units_out%energy, hm%energy%total), &
1025 ' abs_ev = ', units_from_atomic(units_out%energy, scf%evsum_diff), &
1026 ' rel_ev = ', scf%evsum_diff/(abs(scf%evsum_out)+1e-20)
1027 write(message(2),'(a,es15.2,2(a,es9.2))') &
1028 ' ediff = ', scf%energy_diff, ' abs_dens = ', scf%abs_dens_diff, &
1029 ' rel_dens = ', scf%abs_dens_diff/st%qtot
1030 call messages_info(2, namespace=namespace)
1031
1032 if(.not.scf%lcao_restricted) then
1033 write(message(1),'(a,i6)') 'Matrix vector products: ', scf%eigens%matvec
1034 write(message(2),'(a,i6)') 'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%nik))
1035 call messages_info(2, namespace=namespace)
1036 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, scf%eigens%diff, compact = .true., namespace=namespace)
1037 else
1038 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, compact = .true., namespace=namespace)
1039 end if
1040
1041 if (allocated(hm%vberry)) then
1042 call calc_dipole(dipole, space, gr, st, ions)
1043 call write_dipole(st, hm, space, dipole, namespace=namespace)
1044 end if
1045
1046 if(st%d%ispin > unpolarized) then
1047 call write_magnetic_moments(gr, st, ions, gr%der%boundaries, scf%lmm_r, namespace=namespace)
1048 end if
1049
1050 if(hm%lda_u_level == dft_u_acbn0) then
1051 call lda_u_write_u(hm%lda_u, namespace=namespace)
1052 call lda_u_write_v(hm%lda_u, namespace=namespace)
1053 end if
1054
1055 write(message(1),'(a)') ''
1056 write(message(2),'(a,i5,a,f14.2)') 'Elapsed time for SCF step ', iter,':', etime
1057 call messages_info(2, namespace=namespace)
1058
1059 call scf_print_mem_use(namespace)
1060
1061 call messages_print_with_emphasis(namespace=namespace)
1062
1063 end if
1064
1065 if ( scf%verbosity_ == verb_compact ) then
1066 write(message(1),'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') &
1067 'iter ', iter, &
1068 ' : etot ', units_from_atomic(units_out%energy, hm%energy%total), &
1069 ' : abs_dens', scf%abs_dens_diff, &
1070 ' : etime ', etime, 's'
1071 call messages_info(1, namespace=namespace)
1072 end if
1073
1074 pop_sub(scf_run.scf_write_iter)
1075 end subroutine scf_write_iter
1076
1077
1078 ! -----------------------------------------------------
1079 subroutine write_convergence_file(dir, fname)
1080 character(len=*), intent(in) :: dir
1081 character(len=*), intent(in) :: fname
1082
1083 integer :: iunit
1084
1085 if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes
1086 call io_mkdir(dir, namespace)
1087 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write', position='append')
1088 write(iunit, '(i5,es18.8)', advance = 'no') iter, units_from_atomic(units_out%energy, hm%energy%total)
1089 call iterator%start(scf%criterion_list)
1090 do while (iterator%has_next())
1091 crit => iterator%get_next()
1092 select type (crit)
1093 type is (energy_criterion_t)
1094 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, crit%val_abs)
1095 type is (density_criterion_t)
1096 write(iunit, '(2es13.5)', advance = 'no') crit%val_abs, crit%val_rel
1097 type is (eigenval_criterion_t)
1098 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, crit%val_abs)
1099 write(iunit, '(es13.5)', advance = 'no') crit%val_rel
1100 class default
1101 assert(.false.)
1102 end select
1103 end do
1104 if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
1105 .and. ks%theory_level /= generalized_kohn_sham_dft) then
1106 if (ks%oep%level == oep_level_full) then
1107 write(iunit, '(es13.5)', advance = 'no') ks%oep%norm2ss
1108 end if
1109 end if
1110 write(iunit,'(a)') ''
1111 call io_close(iunit)
1112 end if
1113 end subroutine write_convergence_file
1114
1115 end subroutine scf_iter
1116
1117 logical function scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, &
1118 verbosity, iters_done) result(completed)
1119 type(scf_t), intent(inout) :: scf
1120 type(namespace_t), intent(in) :: namespace
1121 type(electron_space_t), intent(in) :: space
1122 type(grid_t), intent(inout) :: gr
1123 type(ions_t), intent(inout) :: ions
1124 type(states_elec_t), intent(inout) :: st
1125 type(v_ks_t), intent(inout) :: ks
1126 type(hamiltonian_elec_t), intent(inout) :: hm
1127 integer, intent(in) :: iter
1128 type(output_t), optional, intent(in) :: outp
1129 integer, optional, intent(in) :: verbosity
1130 integer, optional, intent(out) :: iters_done
1131
1132 character(len=MAX_PATH_LEN) :: dirname
1133 integer(int64) :: what_i
1134
1135 push_sub(scf_iter_finish)
1136
1137 completed = .false.
1138
1139 if(scf%finish) then
1140 if(present(iters_done)) iters_done = iter
1141 if(scf%verbosity_ >= verb_compact) then
1142 write(message(1), '(a, i4, a)') 'Info: SCF converged in ', iter, ' iterations'
1143 write(message(2), '(a)') ''
1144 call messages_info(2, namespace=namespace)
1145 end if
1146 completed = .true.
1147 pop_sub(scf_iter_finish)
1148 return
1149 end if
1150 if (present(outp)) then
1151 if (any(outp%what) .and. outp%duringscf) then
1152 do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
1153 if (outp%what_now(what_i, iter)) then
1154 write(dirname,'(a,a,i4.4)') trim(outp%iter_dir),"scf.", iter
1155 call output_all(outp, namespace, space, dirname, gr, ions, iter, st, hm, ks)
1156 call output_modelmb(outp, namespace, space, dirname, gr, ions, iter, st)
1157 exit
1158 end if
1159 end do
1160 end if
1161 end if
1162
1163 ! save information for the next iteration
1164 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
1165
1166 ! restart mixing
1167 if (scf%mix_field /= option__mixfield__none) then
1168 if (scf%smix%ns_restart > 0) then
1169 if (mod(iter, scf%smix%ns_restart) == 0) then
1170 message(1) = "Info: restarting mixing."
1171 call messages_info(1, namespace=namespace)
1173 end if
1174 end if
1175 end if
1176
1177 select case(scf%mix_field)
1178 case(option__mixfield__potential)
1179 call mixfield_set_vin(scf%mixfield, hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin))
1180 call vtau_mixer_set_vin(scf%vtau_mix, hm)
1181 case (option__mixfield__density)
1182 call mixfield_set_vin(scf%mixfield, scf%rhoin)
1183 end select
1184
1185 !If we use LDA+U, we also have do mix it
1186 if (scf%mix_field /= option__mixfield__states) then
1187 call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
1188 end if
1189
1190 ! check if debug mode should be enabled or disabled on the fly
1191 call io_debug_on_the_fly(namespace)
1192
1193 pop_sub(scf_iter_finish)
1194 end function scf_iter_finish
1195
1196 ! ---------------------------------------------------------
1197 subroutine scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
1198 verbosity, iters_done, restart_dump)
1199 type(scf_t), intent(inout) :: scf
1200 type(namespace_t), intent(in) :: namespace
1201 type(electron_space_t), intent(in) :: space
1202 type(multicomm_t), intent(in) :: mc
1203 type(grid_t), intent(inout) :: gr
1204 type(ions_t), intent(inout) :: ions
1205 type(partner_list_t), intent(in) :: ext_partners
1206 type(states_elec_t), intent(inout) :: st
1207 type(v_ks_t), intent(inout) :: ks
1208 type(hamiltonian_elec_t), intent(inout) :: hm
1209 integer, intent(in) :: iter
1210 type(output_t), optional, intent(in) :: outp
1211 integer, optional, intent(in) :: verbosity
1212 integer, optional, intent(out) :: iters_done
1213 type(restart_t), optional, intent(in) :: restart_dump
1214
1215 integer :: iqn, ib
1216 class(convergence_criterion_t), pointer :: crit
1217 type(criterion_iterator_t) :: iterator
1218
1219
1220 push_sub(scf_finish)
1221
1222 if(scf%lcao_restricted) call lcao_end(scf%lcao)
1223
1224 if ((scf%max_iter > 0 .and. scf%mix_field == option__mixfield__potential) .or. scf%calc_current) then
1225 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
1226 calc_current=scf%calc_current)
1227 end if
1228
1229 select case(scf%mix_field)
1230 case(option__mixfield__states)
1231
1232 do iqn = st%d%kpt%start, st%d%kpt%end
1233 do ib = st%group%block_start, st%group%block_end
1234 call scf%psioutb(ib, iqn)%end()
1235 end do
1236 end do
1237
1238 ! There is a ICE with foss2022a-serial. I am changing to deallocate - NTD
1239 deallocate(scf%psioutb)
1240 end select
1241
1242 safe_deallocate_a(scf%rhoout)
1243 safe_deallocate_a(scf%rhoin)
1244
1245 if (scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then
1246 write(message(1),'(a)') 'Some of the states are not fully converged!'
1247 call messages_warning(1, namespace=namespace)
1248 end if
1249
1250 if (.not.scf%finish) then
1251 write(message(1), '(a,i4,a)') 'SCF *not* converged after ', iter - 1, ' iterations.'
1252 call messages_warning(1, namespace=namespace)
1253 end if
1254
1255 write(message(1), '(a,i10)') 'Info: Number of matrix-vector products: ', scf%matvec
1256 call messages_info(1)
1257
1258 if (scf%calc_force) then
1259 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
1260 end if
1261
1262 if (scf%calc_stress) call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
1263
1264 if(scf%max_iter == 0) then
1265 call energy_calc_eigenvalues(namespace, hm, gr%der, st)
1266 call states_elec_fermi(st, namespace, gr)
1267 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, namespace=namespace)
1268 end if
1269
1270 if(present(outp)) then
1271 ! output final information
1272 call scf_write_static(static_dir, "info")
1273 call output_all(outp, namespace, space, static_dir, gr, ions, -1, st, hm, ks)
1274 call output_modelmb(outp, namespace, space, static_dir, gr, ions, -1, st)
1275 end if
1276
1277 if (space%is_periodic() .and. st%nik > st%d%nspin) then
1278 if (bitand(hm%kpoints%method, kpoints_path) /= 0) then
1279 call states_elec_write_bandstructure(static_dir, namespace, st%nst, st, &
1280 ions, gr, hm%kpoints, hm%phase, vec_pot = hm%hm_base%uniform_vector_potential, &
1281 vec_pot_var = hm%hm_base%vector_potential)
1282 end if
1283 end if
1284
1285 if (ks%vdw%vdw_correction == option__vdwcorrection__vdw_ts) then
1286 call vdw_ts_write_c6ab(ks%vdw%vdw_ts, ions, static_dir, 'c6ab_eff', namespace)
1287 end if
1288
1289 safe_deallocate_a(scf%vhxc_old)
1291 pop_sub(scf_finish)
1292
1293 contains
1294
1295 ! ---------------------------------------------------------
1296 subroutine scf_write_static(dir, fname)
1297 character(len=*), intent(in) :: dir, fname
1298
1299 integer :: iunit, iatom
1300 real(real64), allocatable :: hirshfeld_charges(:)
1301 real(real64) :: dipole(1:space%dim)
1302 real(real64) :: ex_virial
1303
1304 push_sub(scf_run.scf_write_static)
1305
1306 if(mpi_grp_is_root(mpi_world)) then ! this the absolute master writes
1307 call io_mkdir(dir, namespace)
1308 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
1309
1310 call grid_write_info(gr, iunit=iunit)
1311
1312 call symmetries_write_info(gr%symm, space, iunit=iunit)
1313
1314 if (space%is_periodic()) then
1315 call hm%kpoints%write_info(iunit=iunit)
1316 write(iunit,'(1x)')
1317 end if
1318
1319 call v_ks_write_info(ks, iunit=iunit)
1320
1321 ! scf information
1322 if(scf%finish) then
1323 write(iunit, '(a, i4, a)')'SCF converged in ', iter, ' iterations'
1324 else
1325 write(iunit, '(a)') 'SCF *not* converged!'
1326 end if
1327 write(iunit, '(1x)')
1328
1329 if(any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then
1330 write(iunit,'(a)') 'Some of the states are not fully converged!'
1331 end if
1332
1333 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, iunit=iunit)
1334 write(iunit, '(1x)')
1335
1336 if (space%is_periodic()) then
1337 call states_elec_write_gaps(iunit, st, space)
1338 write(iunit, '(1x)')
1339 end if
1340
1341 write(iunit, '(3a)') 'Energy [', trim(units_abbrev(units_out%energy)), ']:'
1342 else
1343 iunit = -1
1344 end if
1345
1346 call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full = .true.)
1347
1348 if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
1349 if(st%d%ispin > unpolarized) then
1350 call write_magnetic_moments(gr, st, ions, gr%der%boundaries, scf%lmm_r, iunit=iunit)
1351 if (mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
1352 end if
1353
1354 if(st%d%ispin == spinors .and. space%dim == 3 .and. &
1355 (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) ) then
1356 call write_total_xc_torque(iunit, gr, hm%ks_pot%vxc, st)
1357 if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
1358 end if
1359
1360 if(hm%lda_u_level == dft_u_acbn0) then
1361 call lda_u_write_u(hm%lda_u, iunit=iunit)
1362 call lda_u_write_v(hm%lda_u, iunit=iunit)
1363 if(mpi_grp_is_root(mpi_world)) write(iunit, '(1x)')
1364 end if
1365
1366 if(scf%calc_dipole) then
1367 call calc_dipole(dipole, space, gr, st, ions)
1368 call write_dipole(st, hm, space, dipole, iunit=iunit)
1369 end if
1370
1371 ! This only works when we do not have a correlation part
1372 if(ks%theory_level == kohn_sham_dft .and. &
1373 hm%xc%functional(func_c,1)%family == xc_family_none .and. st%d%ispin /= spinors) then
1374 call energy_calc_virial_ex(gr%der, hm%ks_pot%vxc, st, ex_virial)
1375
1376 if (mpi_grp_is_root(mpi_world)) then
1377 write(iunit, '(3a)') 'Virial relation for exchange [', trim(units_abbrev(units_out%energy)), ']:'
1378 write(iunit,'(a,es14.6)') "Energy from the orbitals ", units_from_atomic(units_out%energy, hm%energy%exchange)
1379 write(iunit,'(a,es14.6)') "Energy from the potential (virial) ", units_from_atomic(units_out%energy, ex_virial)
1380 write(iunit, '(1x)')
1381 end if
1382 end if
1383
1384 if(mpi_grp_is_root(mpi_world)) then
1385 if(scf%max_iter > 0) then
1386 write(iunit, '(a)') 'Convergence:'
1387 call iterator%start(scf%criterion_list)
1388 do while (iterator%has_next())
1389 crit => iterator%get_next()
1390 call crit%write_info(iunit)
1391 end do
1392 write(iunit,'(1x)')
1393 end if
1394 ! otherwise, these values are uninitialized, and unknown.
1395
1396 if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
1397 .and. ks%theory_level /= generalized_kohn_sham_dft) then
1398 if ((ks%oep_photon%level == oep_level_full) .or. (ks%oep_photon%level == oep_level_kli)) then
1399 write(iunit, '(a)') 'Photon observables:'
1400 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon number = ', ks%oep_photon%pt%number(1)
1401 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon ex. = ', ks%oep_photon%pt%ex
1402 write(iunit,'(1x)')
1403 end if
1404 end if
1405
1406 if (scf%calc_force) call forces_write_info(iunit, ions, dir, namespace)
1407
1408 if (scf%calc_stress) then
1409 call output_stress(iunit, space%periodic_dim, st%stress_tensors, all_terms=.false.)
1410 call output_pressure(iunit, space%periodic_dim, st%stress_tensors%total)
1411 end if
1412
1413 end if
1414
1415 if(scf%calc_partial_charges) then
1416 safe_allocate(hirshfeld_charges(1:ions%natoms))
1417
1418 call partial_charges_calculate(gr, st, ions, hirshfeld_charges)
1419
1420 if(mpi_grp_is_root(mpi_world)) then
1421
1422 write(iunit,'(a)') 'Partial ionic charges'
1423 write(iunit,'(a)') ' Ion Hirshfeld'
1424
1425 do iatom = 1, ions%natoms
1426 write(iunit,'(i4,a10,f16.3)') iatom, trim(ions%atom(iatom)%species%get_label()), hirshfeld_charges(iatom)
1427
1428 end do
1429
1430 end if
1431
1432 safe_deallocate_a(hirshfeld_charges)
1433
1434 end if
1435
1436 if(mpi_grp_is_root(mpi_world)) then
1437 call io_close(iunit)
1438 end if
1439
1440 pop_sub(scf_run.scf_write_static)
1441 end subroutine scf_write_static
1442
1443 end subroutine scf_finish
1444
1445 ! ---------------------------------------------------------
1446 subroutine scf_state_info(namespace, st)
1447 type(namespace_t), intent(in) :: namespace
1448 class(states_abst_t), intent(in) :: st
1449
1450 push_sub(scf_state_info)
1451
1452 if (states_are_real(st)) then
1453 call messages_write('Info: SCF using real wavefunctions.')
1454 else
1455 call messages_write('Info: SCF using complex wavefunctions.')
1456 end if
1457 call messages_info(namespace=namespace)
1458
1459 pop_sub(scf_state_info)
1460
1461 end subroutine scf_state_info
1462
1463 ! ---------------------------------------------------------
1464 subroutine scf_print_mem_use(namespace)
1465 type(namespace_t), intent(in) :: namespace
1466 real(real64) :: mem
1467 real(real64) :: mem_tmp
1468
1469 push_sub(scf_print_mem_use)
1470
1471 if(conf%report_memory) then
1472 mem = loct_get_memory_usage()/(1024.0_real64**2)
1473 call mpi_world%allreduce(mem, mem_tmp, 1, mpi_double_precision, mpi_sum)
1474 mem = mem_tmp
1475 write(message(1),'(a,f14.2)') 'Memory usage [Mbytes] :', mem
1476 call messages_info(1, namespace=namespace)
1477 end if
1478
1479 pop_sub(scf_print_mem_use)
1480 end subroutine scf_print_mem_use
1481
1482 ! --------------------------------------------------------
1484 subroutine scf_update_initial_quantity(scf, hm, criterion)
1485 type(scf_t), intent(inout) :: scf
1486 type(hamiltonian_elec_t), intent(in) :: hm
1487 class(convergence_criterion_t), intent(in) :: criterion
1488
1490
1491 select type (criterion)
1492 type is (energy_criterion_t)
1493 scf%energy_in = hm%energy%total
1494 type is (density_criterion_t)
1495 !Do nothing here
1496 type is (eigenval_criterion_t)
1497 !Setting of the value is done in the scf_update_diff_quantity routine
1498 class default
1499 assert(.false.)
1500 end select
1501
1503 end subroutine scf_update_initial_quantity
1504
1505 ! --------------------------------------------------------
1507 subroutine scf_update_diff_quantity(scf, hm, st, gr, rhoout, rhoin, criterion)
1508 type(scf_t), intent(inout) :: scf
1509 type(hamiltonian_elec_t), intent(in) :: hm
1510 type(states_elec_t), intent(in) :: st
1511 type(grid_t), intent(in) :: gr
1512 real(real64), intent(in) :: rhoout(:,:), rhoin(:,:)
1513 class(convergence_criterion_t), intent(in) :: criterion
1514
1515 integer :: is
1516 real(real64), allocatable :: tmp(:)
1517
1518 push_sub(scf_update_diff_quantity)
1519
1520 select type (criterion)
1521 type is (energy_criterion_t)
1522 scf%energy_diff = abs(hm%energy%total - scf%energy_in)
1523
1524 type is (density_criterion_t)
1525 scf%abs_dens_diff = m_zero
1526 safe_allocate(tmp(1:gr%np))
1527 do is = 1, st%d%nspin
1528 tmp = abs(rhoin(1:gr%np, is) - rhoout(1:gr%np, is))
1529 scf%abs_dens_diff = scf%abs_dens_diff + dmf_integrate(gr, tmp)
1530 end do
1531 safe_deallocate_a(tmp)
1532
1533 type is (eigenval_criterion_t)
1534 scf%evsum_out = states_elec_eigenvalues_sum(st)
1535 scf%evsum_diff = abs(scf%evsum_out - scf%evsum_in)
1536 scf%evsum_in = scf%evsum_out
1537
1538 class default
1539 assert(.false.)
1540 end select
1541
1543 end subroutine scf_update_diff_quantity
1544
1545 ! ---------------------------------------------------------
1546 subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
1547 type(states_elec_t), intent(in) :: st
1548 type(hamiltonian_elec_t), intent(in) :: hm
1549 type(electron_space_t), intent(in) :: space
1550 real(real64), intent(in) :: dipole(:)
1551 integer, optional, intent(in) :: iunit
1552 type(namespace_t), optional, intent(in) :: namespace
1553
1554 push_sub(write_dipole)
1555
1556 if(mpi_grp_is_root(mpi_world)) then
1557 call output_dipole(dipole, space%dim, iunit=iunit, namespace=namespace)
1558
1559 if (space%is_periodic()) then
1560 message(1) = "Defined only up to quantum of polarization (e * lattice vector)."
1561 message(2) = "Single-point Berry's phase method only accurate for large supercells."
1562 call messages_info(2, iunit=iunit, namespace=namespace)
1563
1564 if (hm%kpoints%full%npoints > 1) then
1565 message(1) = &
1566 "WARNING: Single-point Berry's phase method for dipole should not be used when there is more than one k-point."
1567 message(2) = "Instead, finite differences on k-points (not yet implemented) are needed."
1568 call messages_info(2, iunit=iunit, namespace=namespace)
1569 end if
1570
1571 if(.not. smear_is_semiconducting(st%smear)) then
1572 message(1) = "Single-point Berry's phase dipole calculation not correct without integer occupations."
1573 call messages_info(1, iunit=iunit, namespace=namespace)
1574 end if
1575 end if
1576
1577 call messages_info(iunit=iunit, namespace=namespace)
1578 end if
1579
1580 pop_sub(write_dipole)
1581 end subroutine write_dipole
1582
1583
1584end module scf_oct_m
1585
1586
1587!! Local Variables:
1588!! mode: f90
1589!! coding: utf-8
1590!! End:
batchified version of the BLAS axpy routine:
Definition: batch_ops.F90:154
scale a batch by a constant or vector
Definition: batch_ops.F90:162
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:186
This module implements common operations on batches of mesh functions.
Definition: batch_ops.F90:116
subroutine, public berry_perform_internal_scf(this, namespace, space, eigensolver, gr, st, hm, iter, ks, ions, ext_partners)
Definition: berry.F90:184
subroutine, public berry_init(this, namespace)
Definition: berry.F90:159
subroutine, public calc_dipole(dipole, space, mesh, st, ions)
Definition: berry.F90:250
subroutine, public criteria_factory_init(list, namespace, check_conv)
This module implements a calculator for the density and defines related functions.
Definition: density.F90:120
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:609
subroutine, public eigensolver_init(eigens, namespace, gr, st, mc, space)
integer, parameter, public rs_evo
subroutine, public eigensolver_end(eigens)
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
subroutine, public energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
This subroutine calculates the total energy of the system. Basically, it adds up the KS eigenvalues,...
subroutine, public energy_calc_virial_ex(der, vxc, st, ex)
subroutine, public energy_calc_eigenvalues(namespace, hm, der, st)
subroutine, public forces_write_info(iunit, ions, dir, namespace)
Definition: forces.F90:593
subroutine, public forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old, t, dt)
Definition: forces.F90:339
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
Definition: global.F90:216
type(conf_t), public conf
Global instance of Octopus configuration.
Definition: global.F90:178
character(len= *), parameter, public static_dir
Definition: global.F90:252
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
subroutine, public grid_write_info(gr, iunit, namespace)
Definition: grid.F90:522
subroutine, public hamiltonian_elec_update_pot(this, mesh, accumulate)
Update the KS potential of the electronic Hamiltonian.
This module defines classes and functions for interaction partners.
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
subroutine, public io_debug_on_the_fly(namespace)
check if debug mode should be enabled or disabled on the fly
Definition: io.F90:486
subroutine, public io_mkdir(fname, namespace, parents)
Definition: io.F90:311
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
integer, parameter, public kpoints_path
Definition: kpoints.F90:209
A module to handle KS potential, without the external potential.
integer, parameter, public hartree_fock
integer, parameter, public independent_particles
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public kohn_sham_dft
subroutine, public lcao_init_orbitals(this, namespace, st, gr, ions, start)
Definition: lcao.F90:1533
subroutine, public lcao_wf(this, st, gr, ions, hm, namespace, start)
Definition: lcao.F90:990
subroutine, public lcao_end(this)
Definition: lcao.F90:960
subroutine, public lcao_init(this, namespace, space, gr, ions, st)
Definition: lcao.F90:243
logical function, public lcao_is_available(this)
Definition: lcao.F90:1028
subroutine, public lda_u_dump(restart, namespace, this, st, mesh, ierr)
Definition: lda_u_io.F90:646
subroutine, public lda_u_write_u(this, iunit, namespace)
Definition: lda_u_io.F90:530
subroutine, public lda_u_load(restart, this, st, dftu_energy, ierr, occ_only, u_only)
Definition: lda_u_io.F90:727
subroutine, public lda_u_write_v(this, iunit, namespace)
Definition: lda_u_io.F90:579
subroutine, public lda_u_mixer_set_vin(this, mixer)
subroutine, public lda_u_mixer_init(this, mixer, st)
subroutine, public lda_u_mixer_clear(mixer, smix)
subroutine, public lda_u_mixer_init_auxmixer(this, namespace, mixer, smix, st)
subroutine, public lda_u_mixer_get_vnew(this, mixer, st)
subroutine, public lda_u_mixer_set_vout(this, mixer)
subroutine, public lda_u_mixer_end(mixer, smix)
integer, parameter, public dft_u_none
Definition: lda_u.F90:201
integer, parameter, public dft_u_acbn0
Definition: lda_u.F90:201
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, hm_base, phase, energy)
Definition: lda_u.F90:797
subroutine, public write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
Definition: magnetic.F90:211
subroutine, public write_total_xc_torque(iunit, mesh, vxc, st)
Definition: magnetic.F90:433
This module is intended to contain "only mathematical" functions and procedures.
Definition: math.F90:115
This module defines various routines, operating on mesh functions.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
character(len=512), private msg
Definition: messages.F90:165
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
subroutine, public messages_new_line()
Definition: messages.F90:1134
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
real(real64) pure function, public mix_coefficient(this)
Definition: mix.F90:766
subroutine, public mixing(namespace, smix)
Main entry-point to SCF mixer.
Definition: mix.F90:792
subroutine, public mix_get_field(this, mixfield)
Definition: mix.F90:784
subroutine, public mix_dump(namespace, restart, smix, space, mesh, ierr)
Definition: mix.F90:543
subroutine, public mix_load(namespace, restart, smix, space, mesh, ierr)
Definition: mix.F90:641
subroutine, public mix_init(smix, namespace, space, der, d1, d2, def_, func_type_, prefix_)
Initialise mix_t instance.
Definition: mix.F90:265
subroutine, public mix_end(smix)
Definition: mix.F90:522
subroutine, public mix_clear(smix)
Definition: mix.F90:507
subroutine, public modelmb_sym_all_states(space, mesh, st)
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
Definition: mpi.F90:430
type(mpi_grp_t), public mpi_world
Definition: mpi.F90:266
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
this module contains the low-level part of the output system
Definition: output_low.F90:115
subroutine, public output_modelmb(outp, namespace, space, dir, gr, ions, iter, st)
this module contains the output system
Definition: output.F90:115
logical function, public output_needs_current(outp, states_are_real)
Definition: output.F90:969
subroutine, public output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
Definition: output.F90:493
subroutine, public partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:623
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:552
logical function, public clean_stop(comm)
returns true if a file named stop exists
Definition: restart.F90:276
integer, parameter, public restart_flag_mix
Definition: restart.F90:249
integer, parameter, public restart_flag_rho
Definition: restart.F90:249
integer, parameter, public restart_flag_vhxc
Definition: restart.F90:249
logical pure function, public restart_has_flag(restart, flag)
Returns true if...
Definition: restart.F90:979
subroutine, public scf_load(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, restart_load)
Loading of restarting data of the SCF cycle.
Definition: scf.F90:590
logical function, public scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
Definition: scf.F90:1212
subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
Definition: scf.F90:1640
subroutine scf_update_initial_quantity(scf, hm, criterion)
Update the quantity at the begining of a SCF cycle.
Definition: scf.F90:1578
subroutine scf_update_diff_quantity(scf, hm, st, gr, rhoout, rhoin, criterion)
Update the quantity at the begining of a SCF cycle.
Definition: scf.F90:1601
subroutine, public scf_state_info(namespace, st)
Definition: scf.F90:1540
subroutine, public scf_print_mem_use(namespace)
Definition: scf.F90:1558
subroutine, public scf_mix_clear(scf)
Definition: scf.F90:573
integer, parameter, public verb_full
Definition: scf.F90:203
subroutine, public scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, verbosity, iters_done, restart_dump)
Definition: scf.F90:1292
integer, parameter, public verb_compact
Definition: scf.F90:203
subroutine, public scf_init(scf, namespace, gr, ions, st, mc, hm, space)
Definition: scf.F90:254
subroutine, public scf_end(scf)
Definition: scf.F90:543
subroutine, public scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, verbosity, iters_done, restart_dump)
Legacy version of the SCF code.
Definition: scf.F90:819
subroutine, public scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, verbosity, iters_done, restart_dump)
Definition: scf.F90:863
subroutine, public scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
Preparation of the SCF cycle.
Definition: scf.F90:672
logical pure function, public smear_is_semiconducting(this)
Definition: smear.F90:833
pure logical function, public states_are_real(st)
This module defines routines to write information about states.
subroutine, public states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
write the eigenvalues for some states to a file.
subroutine, public states_elec_write_gaps(iunit, st, space)
calculate gaps and write to a file.
subroutine, public states_elec_write_bandstructure(dir, namespace, nst, st, ions, mesh, kpoints, phase, vec_pot, vec_pot_var)
calculate and write the bandstructure
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
real(real64) function, public states_elec_eigenvalues_sum(st, alt_eig)
function to calculate the eigenvalues sum using occupations as weights
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_dump(restart, space, st, mesh, kpoints, ierr, iter, lr, st_start_writing, verbose)
subroutine, public states_elec_load_rho(restart, space, st, mesh, ierr)
subroutine, public states_elec_dump_rho(restart, space, st, mesh, ierr, iter)
This module implements the calculation of the stress tensor.
Definition: stress.F90:118
subroutine, public output_pressure(iunit, space_dim, total_stress_tensor)
Definition: stress.F90:1113
subroutine, public stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
This computes the total stress on the lattice.
Definition: stress.F90:186
subroutine, public output_stress(iunit, space_dim, stress_tensors, all_terms)
Definition: stress.F90:1051
subroutine, public symmetries_write_info(this, space, iunit, namespace)
Definition: symmetries.F90:601
type(type_t), public type_float
Definition: types.F90:133
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
Definition: unit.F90:132
character(len=20) pure function, public units_abbrev(this)
Definition: unit.F90:223
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
type(unit_system_t), public units_inp
the units systems for reading and writing
This module is intended to contain simple general-purpose utility functions and procedures.
Definition: utils.F90:118
subroutine, public output_dipole(dipole, ndim, iunit, namespace)
Definition: utils.F90:278
subroutine, public v_ks_write_info(ks, iunit, namespace)
Definition: v_ks.F90:642
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Definition: v_ks.F90:730
Tkatchenko-Scheffler pairwise method for van der Waals (vdW, dispersion) interactions.
Definition: vdw_ts.F90:119
subroutine, public vdw_ts_write_c6ab(this, ions, dir, fname, namespace)
Definition: vdw_ts.F90:548
subroutine, public vtau_mixer_end(mixer, smix)
Definition: vtau_mixer.F90:187
subroutine, public vtau_mixer_init_auxmixer(namespace, mixer, smix, hm, np, nspin)
Definition: vtau_mixer.F90:150
subroutine, public vtau_mixer_set_vout(mixer, hm)
Definition: vtau_mixer.F90:200
subroutine, public vtau_mixer_get_vnew(mixer, hm)
Definition: vtau_mixer.F90:226
subroutine, public vtau_mixer_clear(mixer, smix)
Definition: vtau_mixer.F90:174
subroutine, public vtau_mixer_set_vin(mixer, hm)
Definition: vtau_mixer.F90:213
This module provices a simple timer class which can be used to trigger the writing of a restart file ...
Definition: walltimer.F90:121
logical function, public walltimer_alarm(comm, print)
indicate whether time is up
Definition: walltimer.F90:303
integer, parameter, public func_c
integer, parameter, public oep_level_full
Definition: xc_oep.F90:172
integer, parameter, public oep_level_kli
Definition: xc_oep.F90:172
subroutine scf_write_static(dir, fname)
Definition: rdmft.F90:581
subroutine create_convergence_file(dir, fname)
Definition: scf.F90:781
subroutine scf_write_iter(namespace)
Definition: scf.F90:1106
subroutine write_convergence_file(dir, fname)
Definition: scf.F90:1173
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Definition: grid.F90:168
Stores all communicators and groups.
Definition: multicomm.F90:206
output handler class
Definition: output_low.F90:164
some variables used for the SCF cycle
Definition: scf.F90:209
abstract class for states
The states_elec_t class contains all electronic wave functions.
batches of electronic states
Definition: wfs_elec.F90:138
int true(void)