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, mc, 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(multicomm_t), intent(in) :: mc
481 type(grid_t), intent(inout) :: gr
482 type(ions_t), intent(in) :: ions
483 type(partner_list_t), intent(in) :: ext_partners
484 type(states_elec_t), intent(inout) :: st
485 type(v_ks_t), intent(inout) :: ks
486 type(hamiltonian_elec_t), intent(inout) :: hm
487 type(restart_t), intent(in) :: restart_load
488
489 integer :: ierr, is
490
491 push_sub(scf_load)
492
493 if (restart_load%has_flag(restart_flag_rho)) then
494 ! Load density and used it to recalculated the KS potential.
495 call states_elec_load_rho(restart_load, space, st, gr, ierr)
496 if (ierr /= 0) then
497 message(1) = 'Unable to read density. Density will be calculated from states.'
498 call messages_warning(1, namespace=namespace)
499 else
500 if (bitand(ks%xc_family, xc_family_oep) == 0) then
501 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
502 else
503 if (.not. restart_load%has_flag(restart_flag_vhxc) .and. ks%oep%level /= oep_level_full) then
504 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
505 end if
506 end if
507 end if
508 end if
509
510 if (restart_load%has_flag(restart_flag_vhxc)) then
511 call hm%ks_pot%load(restart_load, space, gr, ierr)
512 if (ierr /= 0) then
513 message(1) = 'Unable to read Vhxc. Vhxc will be calculated from states.'
514 call messages_warning(1, namespace=namespace)
515 else
516 call hm%update(gr, namespace, space, ext_partners)
517 if (bitand(ks%xc_family, xc_family_oep) /= 0) then
518 if (ks%oep%level == oep_level_full) then
519 do is = 1, st%d%nspin
520 ks%oep%vxc(1:gr%np, is) = hm%ks_pot%vhxc(1:gr%np, is) - hm%ks_pot%vhartree(1:gr%np)
521 end do
522 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
523 end if
524 end if
525 end if
526 end if
527
528 if (restart_load%has_flag(restart_flag_mix)) then
529 if (scf%mix_field == option__mixfield__density .or. scf%mix_field == option__mixfield__potential) then
530 call mix_load(namespace, restart_load, scf%smix, space, gr, ierr)
531 end if
532 if (ierr /= 0) then
533 message(1) = "Unable to read mixing information. Mixing will start from scratch."
534 call messages_warning(1, namespace=namespace)
535 end if
536 end if
537
538 if(hm%lda_u_level /= dft_u_none) then
539 call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr)
540 if (ierr /= 0) then
541 message(1) = "Unable to read DFT+U information. DFT+U data will be calculated from states."
542 call messages_warning(1, namespace=namespace)
543 end if
544
545 ! As v_ks_calc has already been called, we need to update hm%energy%int_dft_u
546 call v_ks_update_dftu_energy(ks, namespace, hm, st, hm%energy%int_dft_u)
547 end if
548
549 !TODO: Create a dedicated routine and call it from the initialize
550
551! if (present(outp) .and. mpi_world%is_root()) then
552! call io_rm(STATIC_DIR //"info")
553! end if
554! end if
555
556 pop_sub(scf_load)
557 end subroutine scf_load
558
559 ! ---------------------------------------------------------
561 subroutine scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
562 type(scf_t), intent(inout) :: scf
563 type(namespace_t), intent(in) :: namespace
564 type(electron_space_t), intent(in) :: space
565 type(grid_t), intent(inout) :: gr
566 type(ions_t), intent(inout) :: ions
567 type(states_elec_t), intent(inout) :: st
568 type(v_ks_t), intent(inout) :: ks
569 type(hamiltonian_elec_t), intent(inout) :: hm
570 type(output_t), optional, intent(in) :: outp
571 integer, optional, intent(in) :: verbosity
572
573 integer :: ib, iqn
574
575 push_sub(scf_start)
576
577 if(scf%forced_finish) then
578 message(1) = "Previous clean stop, not doing SCF and quitting."
579 call messages_fatal(1, only_root_writes = .true., namespace=namespace)
580 end if
581
582 if (.not. hm%is_hermitian()) then
583 message(1) = "Trying to run a SCF calculation for a non-hermitian Hamiltonian. This is not supported."
584 call messages_fatal(1, namespace=namespace)
585 end if
586
587 scf%verbosity_ = optional_default(verbosity, verb_full)
588
589 scf%output_during_scf = .false.
590 scf%output_forces = .false.
591 scf%calc_current = .false.
592
593 if (present(outp)) then
594 ! if the user has activated output=stress but not SCFCalculateStress,
595 ! we assume that is implied
596 if (outp%what(option__output__stress)) then
597 scf%calc_stress = .true.
598 end if
599
600 scf%output_during_scf = outp%duringscf
601 scf%calc_current = output_needs_current(outp, states_are_real(st))
602
603 if (outp%duringscf .and. outp%what(option__output__forces)) then
604 scf%output_forces = .true.
605 end if
606 end if
607
608 safe_allocate(scf%rhoout(1:gr%np, 1:st%d%nspin))
609 safe_allocate(scf%rhoin (1:gr%np, 1:st%d%nspin))
610
611 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
612 scf%rhoout = m_zero
613
614 if (scf%calc_force .or. scf%output_forces) then
615 !We store the Hxc potential for the contribution to the forces
616 safe_allocate(scf%vhxc_old(1:gr%np, 1:st%d%nspin))
617 call lalg_copy(gr%np, st%d%nspin, hm%ks_pot%vhxc, scf%vhxc_old)
618 end if
619
620
621 select case(scf%mix_field)
622 case(option__mixfield__potential)
623 call mixfield_set_vin(scf%mixfield, hm%ks_pot%vhxc)
624 call vtau_mixer_set_vin(scf%vtau_mix, hm)
625 case(option__mixfield__density)
626 call mixfield_set_vin(scf%mixfield, scf%rhoin)
627
628 case(option__mixfield__states)
629
630 ! There is a ICE with foss2022a-serial. I am changing to allocate - NTD
631 allocate(wfs_elec_t::scf%psioutb (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
632
633 do iqn = st%d%kpt%start, st%d%kpt%end
634 do ib = st%group%block_start, st%group%block_end
635 call st%group%psib(ib, iqn)%copy_to(scf%psioutb(ib, iqn))
636 end do
637 end do
638
639 end select
640
641 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
642 ! If we use DFT+U, we also have do mix it
643 if (scf%mix_field /= option__mixfield__states) call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
644
645 call create_convergence_file(static_dir, "convergence")
646
647 if ( scf%verbosity_ /= verb_no ) then
648 if(scf%max_iter > 0) then
649 write(message(1),'(a)') 'Info: Starting SCF iteration.'
650 else
651 write(message(1),'(a)') 'Info: No SCF iterations will be done.'
652 ! we cannot tell whether it is converged.
653 scf%finish = .false.
654 end if
655 call messages_info(1, namespace=namespace)
656 end if
657
658 scf%converged_current = .false.
659 scf%matvec = 0
660
661 pop_sub(scf_start)
662
663 contains
664
665 ! -----------------------------------------------------
666
667 subroutine create_convergence_file(dir, fname)
668 character(len=*), intent(in) :: dir
669 character(len=*), intent(in) :: fname
670
671 integer :: iunit
672 character(len=12) :: label
673 if(mpi_world%is_root()) then
674 call io_mkdir(dir, namespace)
675 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
676 write(iunit, '(a)', advance = 'no') '#iter energy '
677 label = 'energy_diff'
678 write(iunit, '(1x,a)', advance = 'no') label
679 label = 'abs_dens'
680 write(iunit, '(1x,a)', advance = 'no') label
681 label = 'rel_dens'
682 write(iunit, '(1x,a)', advance = 'no') label
683 label = 'abs_ev'
684 write(iunit, '(1x,a)', advance = 'no') label
685 label = 'rel_ev'
686 write(iunit, '(1x,a)', advance = 'no') label
687 if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
688 .and. ks%theory_level /= generalized_kohn_sham_dft) then
689 if (ks%oep%level == oep_level_full) then
690 label = 'OEP norm2ss'
691 write(iunit, '(1x,a)', advance = 'no') label
692 end if
693 end if
694 write(iunit,'(a)') ''
695 call io_close(iunit)
696 end if
697
698 end subroutine create_convergence_file
699
700 end subroutine scf_start
701
702 ! ---------------------------------------------------------
704 subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
705 verbosity, iters_done, restart_dump)
706 type(scf_t), intent(inout) :: scf
707 type(namespace_t), intent(in) :: namespace
708 type(electron_space_t), intent(in) :: space
709 type(multicomm_t), intent(in) :: mc
710 type(grid_t), intent(inout) :: gr
711 type(ions_t), intent(inout) :: ions
712 type(partner_list_t), intent(in) :: ext_partners
713 type(states_elec_t), intent(inout) :: st
714 type(v_ks_t), intent(inout) :: ks
715 type(hamiltonian_elec_t), intent(inout) :: hm
716 type(output_t), optional, intent(in) :: outp
717 integer, optional, intent(in) :: verbosity
718 integer, optional, intent(out) :: iters_done
719 type(restart_t), optional, intent(in) :: restart_dump
720
721 integer :: iter
722 logical :: completed
723
724 push_sub(scf_run)
725
726 call scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
727
728 ! SCF cycle
729 do iter = 1, scf%max_iter
730
731 call scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
732 verbosity, restart_dump)
733
734 completed = scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
735
736 if(scf%forced_finish .or. completed) then
737 exit
738 end if
739 end do
740
741 call scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
742 verbosity, iters_done, restart_dump)
743
744 pop_sub(scf_run)
745 end subroutine scf_run
746
747 ! ---------------------------------------------------------
748 subroutine scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
749 verbosity, restart_dump)
750 type(scf_t), intent(inout) :: scf
751 type(namespace_t), intent(in) :: namespace
752 type(electron_space_t), intent(in) :: space
753 type(multicomm_t), intent(in) :: mc
754 type(grid_t), intent(inout) :: gr
755 type(ions_t), intent(inout) :: ions
756 type(partner_list_t), intent(in) :: ext_partners
757 type(states_elec_t), intent(inout) :: st
758 type(v_ks_t), intent(inout) :: ks
759 type(hamiltonian_elec_t), intent(inout) :: hm
760 integer, intent(in) :: iter
761 type(output_t), optional, intent(in) :: outp
762 integer, optional, intent(in) :: verbosity
763 type(restart_t), optional, intent(in) :: restart_dump
764
765 integer :: iqn, ib, ierr
766 class(convergence_criterion_t), pointer :: crit
767 type(criterion_iterator_t) :: iterator
768 logical :: is_crit_conv
769 real(real64) :: etime, itime
770
771 push_sub(scf_iter)
772
773 call profiling_in("SCF_CYCLE")
774 itime = loct_clock()
775
776 ! this initialization seems redundant but avoids improper optimization at -O3 by PGI 7 on chum,
777 ! which would cause a failure of testsuite/linear_response/04-vib_modes.03-vib_modes_fd.inp
778 scf%eigens%converged = 0
779
780 ! keep the information about the spectrum up to date, needed e.g. for Chebyshev expansion for imaginary time
781 call hm%update_span(gr%spacing(1:space%dim), minval(st%eigenval(:, :)), namespace)
782
783 !We update the quantities at the begining of the scf cycle
784 if (iter == 1) then
785 scf%evsum_in = states_elec_eigenvalues_sum(st)
786 end if
787 call iterator%start(scf%criterion_list)
788 do while (iterator%has_next())
789 crit => iterator%get_next()
790 call scf_update_initial_quantity(scf, hm, crit)
791 end do
792
793 if (scf%calc_force .or. scf%output_forces) then
794 !Used for computing the imperfect convegence contribution to the forces
795 scf%vhxc_old(1:gr%np, 1:st%d%nspin) = hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin)
796 end if
797
798 !We check if the system is coupled with a partner that requires self-consistency
799 ! if(hamiltonian_has_scf_partner(hm)) then
800 if (allocated(hm%vberry)) then
801 !In this case, v_Hxc is frozen and we do an internal SCF loop over the
802 ! partners that require SCF
803 ks%frozen_hxc = .true.
804 ! call perform_scf_partners()
805 call berry_perform_internal_scf(scf%berry, namespace, space, scf%eigens, gr, st, hm, iter, ks, ions, ext_partners)
806 !and we unfreeze the potential once finished
807 ks%frozen_hxc = .false.
808 else
809 scf%eigens%converged = 0
810 call scf%eigens%run(namespace, gr, st, hm, space, ext_partners, iter)
811 end if
812
813 scf%matvec = scf%matvec + scf%eigens%matvec
814
815 ! occupations
816 call states_elec_fermi(st, namespace, gr)
817 call lda_u_update_occ_matrices(hm%lda_u, namespace, gr, st, hm%hm_base, hm%phase, hm%energy)
818
819 ! compute output density, potential (if needed) and eigenvalues sum
820 call density_calc(st, gr, st%rho)
821
822 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoout)
823
824 select case (scf%mix_field)
825 case (option__mixfield__potential)
826 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
827 call mixfield_set_vout(scf%mixfield, hm%ks_pot%vhxc)
828 call vtau_mixer_set_vout(scf%vtau_mix, hm)
829 case (option__mixfield__density)
830 call mixfield_set_vout(scf%mixfield, scf%rhoout)
831 case(option__mixfield__states)
832 do iqn = st%d%kpt%start, st%d%kpt%end
833 do ib = st%group%block_start, st%group%block_end
834 call st%group%psib(ib, iqn)%copy_data_to(gr%np, scf%psioutb(ib, iqn))
835 end do
836 end do
837 end select
838
839 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none) then
840 call lda_u_mixer_set_vout(hm%lda_u, scf%lda_u_mix)
841 endif
842
843 ! recalculate total energy
844 call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit = -1)
845
846 if (present(outp)) then
847 ! compute forces only if requested
848 if (outp%duringscf .and. outp%what_now(option__output__forces, iter)) then
849 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
850 end if
851 end if
852
853 !We update the quantities at the end of the scf cycle
854 call iterator%start(scf%criterion_list)
855 do while (iterator%has_next())
856 crit => iterator%get_next()
857 call scf_update_diff_quantity(scf, hm, st, gr, scf%rhoout, scf%rhoin, crit)
858 end do
859
860 ! are we finished?
861 scf%converged_last = scf%converged_current
862
863 scf%converged_current = scf%check_conv .and. &
864 (.not. scf%conv_eigen_error .or. all(scf%eigens%converged >= st%nst_conv))
865 !Loop over the different criteria
866 call iterator%start(scf%criterion_list)
867 do while (iterator%has_next())
868 crit => iterator%get_next()
869 call crit%is_converged(is_crit_conv)
870 scf%converged_current = scf%converged_current .and. is_crit_conv
871 end do
872
873 ! only finish if the convergence criteria are fulfilled in two
874 ! consecutive iterations
875 scf%finish = scf%converged_last .and. scf%converged_current
876
877 etime = loct_clock() - itime
878 call scf_write_iter(namespace)
879
880 ! mixing
881 select case (scf%mix_field)
882 case (option__mixfield__density)
883 ! mix input and output densities and compute new potential
884 call mixing(namespace, scf%smix)
885 call mixfield_get_vnew(scf%mixfield, st%rho)
886 ! for spinors, having components 3 or 4 be negative is not unphysical
887 if (minval(st%rho(1:gr%np, 1:st%d%spin_channels)) < -1e-6_real64) then
888 write(message(1),*) 'Negative density after mixing. Minimum value = ', &
889 minval(st%rho(1:gr%np, 1:st%d%spin_channels))
890 call messages_warning(1, namespace=namespace)
891 end if
892 call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
893 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
894
895 case (option__mixfield__potential)
896 ! mix input and output potentials
897 call mixing(namespace, scf%smix)
898 call mixfield_get_vnew(scf%mixfield, hm%ks_pot%vhxc)
899 call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
900 call vtau_mixer_get_vnew(scf%vtau_mix, hm)
901 call hamiltonian_elec_update_pot(hm, gr)
902
903 case(option__mixfield__states)
904 do iqn = st%d%kpt%start, st%d%kpt%end
905 do ib = st%group%block_start, st%group%block_end
906 call batch_scal(gr%np, m_one - mix_coefficient(scf%smix), st%group%psib(ib, iqn))
907 call batch_axpy(gr%np, mix_coefficient(scf%smix), scf%psioutb(ib, iqn), st%group%psib(ib, iqn))
908 end do
909 end do
910 call density_calc(st, gr, st%rho)
911 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
912
913 case (option__mixfield__none)
914 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
915 end select
916
917 ! Are we asked to stop? (Whenever Fortran is ready for signals, this should go away)
918 scf%forced_finish = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
919
920 if (scf%finish .and. st%modelmbparticles%nparticle > 0) then
921 call modelmb_sym_all_states(space, gr, st)
922 end if
923
924 if (present(outp) .and. present(restart_dump)) then
925 ! save restart information
926
927 if ( (scf%finish .or. (modulo(iter, outp%restart_write_interval) == 0) &
928 .or. iter == scf%max_iter .or. scf%forced_finish) ) then
929
930 call states_elec_dump(scf%restart_dump, space, st, gr, hm%kpoints, ierr, iter=iter)
931 if (ierr /= 0) then
932 message(1) = 'Unable to write states wavefunctions.'
933 call messages_warning(1, namespace=namespace)
934 end if
935
936 call states_elec_dump_rho(scf%restart_dump, space, st, gr, ierr, iter=iter)
937 if (ierr /= 0) then
938 message(1) = 'Unable to write density.'
939 call messages_warning(1, namespace=namespace)
940 end if
941
942 if(hm%lda_u_level /= dft_u_none) then
943 call lda_u_dump(scf%restart_dump, namespace, hm%lda_u, st, gr, ierr)
944 if (ierr /= 0) then
945 message(1) = 'Unable to write DFT+U information.'
946 call messages_warning(1, namespace=namespace)
947 end if
948 end if
949
950 select case (scf%mix_field)
951 case (option__mixfield__density)
952 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
953 if (ierr /= 0) then
954 message(1) = 'Unable to write mixing information.'
955 call messages_warning(1, namespace=namespace)
956 end if
957 case (option__mixfield__potential)
958 call hm%ks_pot%dump(scf%restart_dump, space, gr, ierr)
959 if (ierr /= 0) then
960 message(1) = 'Unable to write Vhxc.'
961 call messages_warning(1, namespace=namespace)
962 end if
963
964 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
965 if (ierr /= 0) then
966 message(1) = 'Unable to write mixing information.'
967 call messages_warning(1, namespace=namespace)
968 end if
969 end select
970
971 end if
972 end if
973
974 call write_convergence_file(static_dir, "convergence")
975
976 call profiling_out("SCF_CYCLE")
977
978 pop_sub(scf_iter)
979 contains
980
981 ! ---------------------------------------------------------
982 subroutine scf_write_iter(namespace)
983 type(namespace_t), intent(in) :: namespace
984
985 character(len=50) :: str
986 real(real64) :: dipole(1:space%dim)
987
988 push_sub(scf_run.scf_write_iter)
989
990 if ( scf%verbosity_ == verb_full ) then
991
992 write(str, '(a,i5)') 'SCF CYCLE ITER #' ,iter
993 call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
994 write(message(1),'(a,es15.8,2(a,es9.2))') ' etot = ', units_from_atomic(units_out%energy, hm%energy%total), &
995 ' abs_ev = ', units_from_atomic(units_out%energy, scf%evsum_diff), &
996 ' rel_ev = ', scf%evsum_diff/(abs(scf%evsum_out)+1e-20)
997 write(message(2),'(a,es15.2,2(a,es9.2))') &
998 ' ediff = ', scf%energy_diff, ' abs_dens = ', scf%abs_dens_diff, &
999 ' rel_dens = ', scf%abs_dens_diff/st%qtot
1000 call messages_info(2, namespace=namespace)
1001
1002 write(message(1),'(a,i6)') 'Matrix vector products: ', scf%eigens%matvec
1003 write(message(2),'(a,i6)') 'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%nik))
1004 call messages_info(2, namespace=namespace)
1005 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, scf%eigens%diff, compact = .true., namespace=namespace)
1006
1007 if (allocated(hm%vberry)) then
1008 call calc_dipole(dipole, space, gr, st, ions)
1009 call write_dipole(st, hm, space, dipole, namespace=namespace)
1010 end if
1011
1012 if(st%d%ispin > unpolarized) then
1013 call write_magnetic_moments(gr, st, ions, gr%der%boundaries, scf%lmm_r, namespace=namespace)
1014 end if
1015
1016 if(hm%lda_u_level == dft_u_acbn0) then
1017 call lda_u_write_u(hm%lda_u, namespace=namespace)
1018 call lda_u_write_v(hm%lda_u, namespace=namespace)
1019 end if
1020
1021 write(message(1),'(a)') ''
1022 write(message(2),'(a,i5,a,f14.2)') 'Elapsed time for SCF step ', iter,':', etime
1023 call messages_info(2, namespace=namespace)
1024
1025 call scf_print_mem_use(namespace)
1026
1027 call messages_print_with_emphasis(namespace=namespace)
1028
1029 end if
1030
1031 if ( scf%verbosity_ == verb_compact ) then
1032 write(message(1),'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') &
1033 'iter ', iter, &
1034 ' : etot ', units_from_atomic(units_out%energy, hm%energy%total), &
1035 ' : abs_dens', scf%abs_dens_diff, &
1036 ' : etime ', etime, 's'
1037 call messages_info(1, namespace=namespace)
1038 end if
1039
1040 pop_sub(scf_run.scf_write_iter)
1041 end subroutine scf_write_iter
1042
1043
1044 ! -----------------------------------------------------
1045 subroutine write_convergence_file(dir, fname)
1046 character(len=*), intent(in) :: dir
1047 character(len=*), intent(in) :: fname
1048
1049 integer :: iunit
1050
1051 if(mpi_world%is_root()) then ! this the absolute master writes
1052 call io_mkdir(dir, namespace)
1053 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write', position='append')
1054 write(iunit, '(i5,es18.8)', advance = 'no') iter, units_from_atomic(units_out%energy, hm%energy%total)
1055 call iterator%start(scf%criterion_list)
1056 do while (iterator%has_next())
1057 crit => iterator%get_next()
1058 select type (crit)
1059 type is (energy_criterion_t)
1060 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, crit%val_abs)
1061 type is (density_criterion_t)
1062 write(iunit, '(2es13.5)', advance = 'no') crit%val_abs, crit%val_rel
1063 type is (eigenval_criterion_t)
1064 write(iunit, '(es13.5)', advance = 'no') units_from_atomic(units_out%energy, crit%val_abs)
1065 write(iunit, '(es13.5)', advance = 'no') crit%val_rel
1066 class default
1067 assert(.false.)
1068 end select
1069 end do
1070 if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
1071 .and. ks%theory_level /= generalized_kohn_sham_dft) then
1072 if (ks%oep%level == oep_level_full) then
1073 write(iunit, '(es13.5)', advance = 'no') ks%oep%norm2ss
1074 end if
1075 end if
1076 write(iunit,'(a)') ''
1077 call io_close(iunit)
1078 end if
1079 end subroutine write_convergence_file
1080
1081 end subroutine scf_iter
1082
1083 logical function scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, &
1084 verbosity, iters_done) result(completed)
1085 type(scf_t), intent(inout) :: scf
1086 type(namespace_t), intent(in) :: namespace
1087 type(electron_space_t), intent(in) :: space
1088 type(grid_t), intent(inout) :: gr
1089 type(ions_t), intent(inout) :: ions
1090 type(states_elec_t), intent(inout) :: st
1091 type(v_ks_t), intent(inout) :: ks
1092 type(hamiltonian_elec_t), intent(inout) :: hm
1093 integer, intent(in) :: iter
1094 type(output_t), optional, intent(in) :: outp
1095 integer, optional, intent(in) :: verbosity
1096 integer, optional, intent(out) :: iters_done
1097
1098 character(len=MAX_PATH_LEN) :: dirname
1099 integer(int64) :: what_i
1100
1101 push_sub(scf_iter_finish)
1102
1103 completed = .false.
1104
1105 if(scf%finish) then
1106 if(present(iters_done)) iters_done = iter
1107 if(scf%verbosity_ >= verb_compact) then
1108 write(message(1), '(a, i4, a)') 'Info: SCF converged in ', iter, ' iterations'
1109 write(message(2), '(a)') ''
1110 call messages_info(2, namespace=namespace)
1111 end if
1112 completed = .true.
1113 pop_sub(scf_iter_finish)
1114 return
1115 end if
1116 if (present(outp)) then
1117 if (any(outp%what) .and. outp%duringscf) then
1118 do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
1119 if (outp%what_now(what_i, iter)) then
1120 write(dirname,'(a,a,i4.4)') trim(outp%iter_dir),"scf.", iter
1121 call output_all(outp, namespace, space, dirname, gr, ions, iter, st, hm, ks)
1122 call output_modelmb(outp, namespace, space, dirname, gr, ions, iter, st)
1123 exit
1124 end if
1125 end do
1126 end if
1127 end if
1128
1129 ! save information for the next iteration
1130 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
1131
1132 ! restart mixing
1133 if (scf%mix_field /= option__mixfield__none) then
1134 if (scf%smix%ns_restart > 0) then
1135 if (mod(iter, scf%smix%ns_restart) == 0) then
1136 message(1) = "Info: restarting mixing."
1137 call messages_info(1, namespace=namespace)
1138 call scf_mix_clear(scf)
1139 end if
1140 end if
1141 end if
1142
1143 select case(scf%mix_field)
1144 case(option__mixfield__potential)
1145 call mixfield_set_vin(scf%mixfield, hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin))
1146 call vtau_mixer_set_vin(scf%vtau_mix, hm)
1147 case (option__mixfield__density)
1148 call mixfield_set_vin(scf%mixfield, scf%rhoin)
1149 end select
1150
1151 !If we use LDA+U, we also have do mix it
1152 if (scf%mix_field /= option__mixfield__states) then
1153 call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
1154 end if
1155
1156 ! check if debug mode should be enabled or disabled on the fly
1157 call io_debug_on_the_fly(namespace)
1158
1159 pop_sub(scf_iter_finish)
1160 end function scf_iter_finish
1161
1162 ! ---------------------------------------------------------
1163 subroutine scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
1164 verbosity, iters_done, restart_dump)
1165 type(scf_t), intent(inout) :: scf
1166 type(namespace_t), intent(in) :: namespace
1167 type(electron_space_t), intent(in) :: space
1168 type(multicomm_t), intent(in) :: mc
1169 type(grid_t), intent(inout) :: gr
1170 type(ions_t), intent(inout) :: ions
1171 type(partner_list_t), intent(in) :: ext_partners
1172 type(states_elec_t), intent(inout) :: st
1173 type(v_ks_t), intent(inout) :: ks
1174 type(hamiltonian_elec_t), intent(inout) :: hm
1175 integer, intent(in) :: iter
1176 type(output_t), optional, intent(in) :: outp
1177 integer, optional, intent(in) :: verbosity
1178 integer, optional, intent(out) :: iters_done
1179 type(restart_t), optional, intent(in) :: restart_dump
1180
1181 integer :: iqn, ib
1182 class(convergence_criterion_t), pointer :: crit
1183 type(criterion_iterator_t) :: iterator
1184
1185
1186 push_sub(scf_finish)
1187
1188 ! Compute the KS potential corresponding to the final density
1189 ! This is critical for getting consistent TD calculations
1190 if ((scf%max_iter > 0 .and. scf%mix_field == option__mixfield__potential) .or. scf%calc_current) then
1191 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
1192 calc_current=scf%calc_current)
1193 end if
1194
1195 select case(scf%mix_field)
1196 case(option__mixfield__states)
1197
1198 do iqn = st%d%kpt%start, st%d%kpt%end
1199 do ib = st%group%block_start, st%group%block_end
1200 call scf%psioutb(ib, iqn)%end()
1201 end do
1202 end do
1203
1204 ! There is a ICE with foss2022a-serial. I am changing to deallocate - NTD
1205 deallocate(scf%psioutb)
1206 end select
1207
1208 safe_deallocate_a(scf%rhoout)
1209 safe_deallocate_a(scf%rhoin)
1210
1211 if (scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst)) then
1212 write(message(1),'(a)') 'Some of the states are not fully converged!'
1213 if (all(scf%eigens%converged >= st%nst_conv)) then
1214 write(message(2),'(a)') 'But all requested states to converge are converged.'
1215 call messages_info(2, namespace=namespace)
1216 else
1217 call messages_warning(1, namespace=namespace)
1218 end if
1219 end if
1220
1221 if (.not.scf%finish) then
1222 write(message(1), '(a,i4,a)') 'SCF *not* converged after ', iter - 1, ' iterations.'
1223 call messages_warning(1, namespace=namespace)
1224 end if
1225
1226 write(message(1), '(a,i10)') 'Info: Number of matrix-vector products: ', scf%matvec
1227 call messages_info(1)
1228
1229 if (scf%calc_force) then
1230 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
1231 end if
1232
1233 if (scf%calc_stress) call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
1234
1235 ! Update the eigenvalues, to match the KS potential that just got recomputed
1236 if (scf%mix_field == option__mixfield__potential) then
1237 call energy_calc_eigenvalues(namespace, hm, gr%der, st)
1238 call states_elec_fermi(st, namespace, gr)
1239 end if
1240
1241 if(present(outp)) then
1242 ! output final information
1243 call scf_write_static(static_dir, "info")
1244 call output_all(outp, namespace, space, static_dir, gr, ions, -1, st, hm, ks)
1245 call output_modelmb(outp, namespace, space, static_dir, gr, ions, -1, st)
1246 end if
1247
1248 if (space%is_periodic() .and. st%nik > st%d%nspin) then
1249 if (bitand(hm%kpoints%method, kpoints_path) /= 0) then
1250 call states_elec_write_bandstructure(static_dir, namespace, st%nst, st, &
1251 ions, gr, hm%kpoints, hm%phase, vec_pot = hm%hm_base%uniform_vector_potential, &
1252 vec_pot_var = hm%hm_base%vector_potential)
1253 end if
1254 end if
1255
1256 if (ks%vdw%vdw_correction == option__vdwcorrection__vdw_ts) then
1257 call vdw_ts_write_c6ab(ks%vdw%vdw_ts, ions, static_dir, 'c6ab_eff', namespace)
1258 end if
1259
1260 safe_deallocate_a(scf%vhxc_old)
1261
1262 pop_sub(scf_finish)
1263
1264 contains
1265
1266 ! ---------------------------------------------------------
1267 subroutine scf_write_static(dir, fname)
1268 character(len=*), intent(in) :: dir, fname
1269
1270 integer :: iunit
1271 real(real64) :: dipole(1:space%dim)
1272 real(real64) :: ex_virial
1273
1274 push_sub(scf_run.scf_write_static)
1275
1276 if(mpi_world%is_root()) then ! this the absolute master writes
1277 call io_mkdir(dir, namespace)
1278 iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
1279
1280 call grid_write_info(gr, iunit=iunit)
1281
1282 call symmetries_write_info(gr%symm, space, iunit=iunit)
1283
1284 if (space%is_periodic()) then
1285 call hm%kpoints%write_info(iunit=iunit)
1286 write(iunit,'(1x)')
1287 end if
1288
1289 call v_ks_write_info(ks, iunit=iunit)
1290
1291 ! scf information
1292 if(scf%finish) then
1293 write(iunit, '(a, i4, a)')'SCF converged in ', iter, ' iterations'
1294 else
1295 write(iunit, '(a)') 'SCF *not* converged!'
1296 end if
1297 write(iunit, '(1x)')
1298
1299 if(any(scf%eigens%converged < st%nst)) then
1300 write(iunit,'(a)') 'Some of the states are not fully converged!'
1301 if (all(scf%eigens%converged >= st%nst_conv)) then
1302 write(iunit,'(a)') 'But all requested states to converge are converged.'
1303 end if
1304 end if
1305
1306 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, iunit=iunit)
1307 write(iunit, '(1x)')
1308
1309 if (space%is_periodic()) then
1310 call states_elec_write_gaps(iunit, st, space)
1311 write(iunit, '(1x)')
1312 end if
1313
1314 write(iunit, '(3a)') 'Energy [', trim(units_abbrev(units_out%energy)), ']:'
1315 else
1316 iunit = -1
1317 end if
1318
1319 call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full = .true.)
1320
1321 if(mpi_world%is_root()) write(iunit, '(1x)')
1322 if(st%d%ispin > unpolarized) then
1323 call write_magnetic_moments(gr, st, ions, gr%der%boundaries, scf%lmm_r, iunit=iunit)
1324 if (mpi_world%is_root()) write(iunit, '(1x)')
1325 end if
1326
1327 if(st%d%ispin == spinors .and. space%dim == 3 .and. &
1328 (ks%theory_level == kohn_sham_dft .or. ks%theory_level == generalized_kohn_sham_dft) ) then
1329 call write_total_xc_torque(iunit, gr, hm%ks_pot%vxc, st)
1330 if(mpi_world%is_root()) write(iunit, '(1x)')
1331 end if
1332
1333 if(hm%lda_u_level == dft_u_acbn0) then
1334 call lda_u_write_u(hm%lda_u, iunit=iunit)
1335 call lda_u_write_v(hm%lda_u, iunit=iunit)
1336 if(mpi_world%is_root()) write(iunit, '(1x)')
1337 end if
1338
1339 if(scf%calc_dipole) then
1340 call calc_dipole(dipole, space, gr, st, ions)
1341 call write_dipole(st, hm, space, dipole, iunit=iunit)
1342 end if
1343
1344 ! This only works when we do not have a correlation part
1345 if(ks%theory_level == kohn_sham_dft .and. &
1346 hm%xc%functional(func_c,1)%family == xc_family_none .and. st%d%ispin /= spinors) then
1347 call energy_calc_virial_ex(gr%der, hm%ks_pot%vxc, st, ex_virial)
1348
1349 if (mpi_world%is_root()) then
1350 write(iunit, '(3a)') 'Virial relation for exchange [', trim(units_abbrev(units_out%energy)), ']:'
1351 write(iunit,'(a,es14.6)') "Energy from the orbitals ", units_from_atomic(units_out%energy, hm%energy%exchange)
1352 write(iunit,'(a,es14.6)') "Energy from the potential (virial) ", units_from_atomic(units_out%energy, ex_virial)
1353 write(iunit, '(1x)')
1354 end if
1355 end if
1356
1357 if(mpi_world%is_root()) then
1358 if(scf%max_iter > 0) then
1359 write(iunit, '(a)') 'Convergence:'
1360 call iterator%start(scf%criterion_list)
1361 do while (iterator%has_next())
1362 crit => iterator%get_next()
1363 call crit%write_info(iunit)
1364 end do
1365 write(iunit,'(1x)')
1366 end if
1367 ! otherwise, these values are uninitialized, and unknown.
1368
1369 if (bitand(ks%xc_family, xc_family_oep) /= 0 .and. ks%theory_level /= hartree_fock &
1370 .and. ks%theory_level /= generalized_kohn_sham_dft) then
1371 if ((ks%oep_photon%level == oep_level_full) .or. (ks%oep_photon%level == oep_level_kli)) then
1372 write(iunit, '(a)') 'Photon observables:'
1373 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon number = ', ks%oep_photon%pt%number(1)
1374 write(iunit, '(6x, a, es15.8,a,es15.8,a)') 'Photon ex. = ', ks%oep_photon%pt%ex
1375 write(iunit,'(1x)')
1376 end if
1377 end if
1378
1379 if (scf%calc_force) call forces_write_info(iunit, ions, dir, namespace)
1380
1381 if (scf%calc_stress) then
1382 call output_stress(iunit, space%periodic_dim, st%stress_tensors, all_terms=.false.)
1383 call output_pressure(iunit, space%periodic_dim, st%stress_tensors%total)
1384 end if
1385
1386 end if
1387
1388 if(scf%calc_partial_charges) then
1389 call partial_charges_compute_and_print_charges(gr, st, ions, iunit)
1390 end if
1391
1392 if(mpi_world%is_root()) then
1393 call io_close(iunit)
1394 end if
1395
1396 pop_sub(scf_run.scf_write_static)
1397 end subroutine scf_write_static
1398
1399 end subroutine scf_finish
1400
1401 ! ---------------------------------------------------------
1402 subroutine scf_state_info(namespace, st)
1403 type(namespace_t), intent(in) :: namespace
1404 class(states_abst_t), intent(in) :: st
1405
1406 push_sub(scf_state_info)
1407
1408 if (states_are_real(st)) then
1409 call messages_write('Info: SCF using real wavefunctions.')
1410 else
1411 call messages_write('Info: SCF using complex wavefunctions.')
1412 end if
1413 call messages_info(namespace=namespace)
1414
1415 pop_sub(scf_state_info)
1416
1417 end subroutine scf_state_info
1418
1419 ! ---------------------------------------------------------
1420 subroutine scf_print_mem_use(namespace)
1421 type(namespace_t), intent(in) :: namespace
1422 real(real64) :: mem
1423 real(real64) :: mem_tmp
1424
1425 push_sub(scf_print_mem_use)
1426
1427 if(conf%report_memory) then
1428 mem = loct_get_memory_usage()/(1024.0_real64**2)
1429 call mpi_world%allreduce(mem, mem_tmp, 1, mpi_double_precision, mpi_sum)
1430 mem = mem_tmp
1431 write(message(1),'(a,f14.2)') 'Memory usage [Mbytes] :', mem
1432 call messages_info(1, namespace=namespace)
1433 end if
1434
1435 pop_sub(scf_print_mem_use)
1436 end subroutine scf_print_mem_use
1437
1438 ! --------------------------------------------------------
1440 subroutine scf_update_initial_quantity(scf, hm, criterion)
1441 type(scf_t), intent(inout) :: scf
1442 type(hamiltonian_elec_t), intent(in) :: hm
1443 class(convergence_criterion_t), intent(in) :: criterion
1444
1446
1447 select type (criterion)
1448 type is (energy_criterion_t)
1449 scf%energy_in = hm%energy%total
1450 type is (density_criterion_t)
1451 !Do nothing here
1452 type is (eigenval_criterion_t)
1453 !Setting of the value is done in the scf_update_diff_quantity routine
1454 class default
1455 assert(.false.)
1456 end select
1457
1459 end subroutine scf_update_initial_quantity
1460
1461 ! --------------------------------------------------------
1463 subroutine scf_update_diff_quantity(scf, hm, st, gr, rhoout, rhoin, criterion)
1464 type(scf_t), intent(inout) :: scf
1465 type(hamiltonian_elec_t), intent(in) :: hm
1466 type(states_elec_t), intent(in) :: st
1467 type(grid_t), intent(in) :: gr
1468 real(real64), intent(in) :: rhoout(:,:), rhoin(:,:)
1469 class(convergence_criterion_t), intent(in) :: criterion
1470
1471 integer :: is
1472 real(real64), allocatable :: tmp(:)
1473
1474 push_sub(scf_update_diff_quantity)
1475
1476 select type (criterion)
1477 type is (energy_criterion_t)
1478 scf%energy_diff = abs(hm%energy%total - scf%energy_in)
1479
1480 type is (density_criterion_t)
1481 scf%abs_dens_diff = m_zero
1482 safe_allocate(tmp(1:gr%np))
1483 do is = 1, st%d%nspin
1484 tmp(:) = abs(rhoin(1:gr%np, is) - rhoout(1:gr%np, is))
1485 scf%abs_dens_diff = scf%abs_dens_diff + dmf_integrate(gr, tmp)
1486 end do
1487 safe_deallocate_a(tmp)
1488
1489 type is (eigenval_criterion_t)
1490 scf%evsum_out = states_elec_eigenvalues_sum(st)
1491 scf%evsum_diff = abs(scf%evsum_out - scf%evsum_in)
1492 scf%evsum_in = scf%evsum_out
1493
1494 class default
1495 assert(.false.)
1496 end select
1499 end subroutine scf_update_diff_quantity
1500
1501 ! ---------------------------------------------------------
1502 subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
1503 type(states_elec_t), intent(in) :: st
1504 type(hamiltonian_elec_t), intent(in) :: hm
1505 type(electron_space_t), intent(in) :: space
1506 real(real64), intent(in) :: dipole(:)
1507 integer, optional, intent(in) :: iunit
1508 type(namespace_t), optional, intent(in) :: namespace
1509
1510 push_sub(write_dipole)
1511
1512 if(mpi_world%is_root()) then
1513 call output_dipole(dipole, space%dim, iunit=iunit, namespace=namespace)
1514
1515 if (space%is_periodic()) then
1516 message(1) = "Defined only up to quantum of polarization (e * lattice vector)."
1517 message(2) = "Single-point Berry's phase method only accurate for large supercells."
1518 call messages_info(2, iunit=iunit, namespace=namespace)
1519
1520 if (hm%kpoints%full%npoints > 1) then
1521 message(1) = &
1522 "WARNING: Single-point Berry's phase method for dipole should not be used when there is more than one k-point."
1523 message(2) = "Instead, finite differences on k-points (not yet implemented) are needed."
1524 call messages_info(2, iunit=iunit, namespace=namespace)
1525 end if
1526
1527 if(.not. smear_is_semiconducting(st%smear)) then
1528 message(1) = "Single-point Berry's phase dipole calculation not correct without integer occupations."
1529 call messages_info(1, iunit=iunit, namespace=namespace)
1530 end if
1531 end if
1532
1533 call messages_info(iunit=iunit, namespace=namespace)
1534 end if
1536 pop_sub(write_dipole)
1537 end subroutine write_dipole
1538
1539
1540end module scf_oct_m
1541
1542
1543!! Local Variables:
1544!! mode: f90
1545!! coding: utf-8
1546!! 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
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:596
subroutine, public forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old, t, dt)
Definition: forces.F90:342
real(real64), parameter, public m_zero
Definition: global.F90:190
integer, parameter, public hartree_fock
Definition: global.F90:236
integer, parameter, public independent_particles
Theory level.
Definition: global.F90:236
integer, parameter, public generalized_kohn_sham_dft
Definition: global.F90:236
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
Definition: global.F90:220
integer, parameter, public kohn_sham_dft
Definition: global.F90:236
type(conf_t), public conf
Global instance of Octopus configuration.
Definition: global.F90:180
character(len= *), parameter, public static_dir
Definition: global.F90:265
real(real64), parameter, public m_half
Definition: global.F90:196
real(real64), parameter, public m_one
Definition: global.F90:191
This module implements the underlying real-space grid.
Definition: grid.F90:119
subroutine, public grid_write_info(gr, iunit, namespace)
Definition: grid.F90:530
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
integer, parameter, public dft_u_acbn0
Definition: lda_u.F90:203
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, hm_base, phase, energy)
Definition: lda_u.F90:800
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:904
character(len=512), private msg
Definition: messages.F90:167
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:531
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1029
subroutine, public messages_new_line()
Definition: messages.F90:1118
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:416
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:697
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1069
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:600
real(real64) pure function, public mix_coefficient(this)
Definition: mix.F90:804
subroutine, public mixing(namespace, smix)
Main entry-point to SCF mixer.
Definition: mix.F90:830
subroutine, public mix_get_field(this, mixfield)
Definition: mix.F90:822
subroutine, public mix_dump(namespace, restart, smix, space, mesh, ierr)
Definition: mix.F90:579
subroutine, public mix_load(namespace, restart, smix, space, mesh, ierr)
Definition: mix.F90:679
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_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_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:572
logical function, public scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
Definition: scf.F90:1180
subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
Definition: scf.F90:1598
subroutine scf_update_initial_quantity(scf, hm, criterion)
Update the quantity at the begining of a SCF cycle.
Definition: scf.F90:1536
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:1559
subroutine, public scf_state_info(namespace, st)
Definition: scf.F90:1498
subroutine, public scf_print_mem_use(namespace)
Definition: scf.F90:1516
subroutine, public scf_mix_clear(scf)
Definition: scf.F90:555
integer, parameter, public verb_full
Definition: scf.F90:205
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:1260
subroutine, public scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, verbosity, restart_dump)
Definition: scf.F90:845
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:801
subroutine, public scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
Preparation of the SCF cycle.
Definition: scf.F90:657
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, 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:120
subroutine, public output_pressure(iunit, space_dim, total_stress_tensor)
Definition: stress.F90:1121
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:1059
subroutine, public symmetries_write_info(this, space, iunit, namespace)
Definition: symmetries.F90:606
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:1507
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:763
subroutine scf_write_iter(namespace)
Definition: scf.F90:1078
subroutine write_convergence_file(dir, fname)
Definition: scf.F90:1141
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)