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