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