41 use,
intrinsic :: iso_fortran_env
88 use xc_functional_oct_m
110 integer,
public,
parameter :: &
118 integer,
public :: max_iter
120 real(real64),
public :: lmm_r
123 logical :: conv_eigen_error
124 logical :: check_conv
127 logical :: lcao_restricted
128 logical :: calc_force
129 logical,
public :: calc_stress
130 logical :: calc_dipole
131 logical :: calc_partial_charges
133 type(mixfield_t),
pointer :: mixfield
134 type(eigensolver_t) :: eigens
136 logical :: forced_finish
137 type(lda_u_mixer_t) :: lda_u_mix
138 type(vtau_mixer_t) :: vtau_mix
139 type(berry_t) :: berry
142 type(restart_t),
public :: restart_load, restart_dump
144 type(criterion_list_t),
public :: criterion_list
145 real(real64) :: energy_in, energy_diff, abs_dens_diff, evsum_in, evsum_out, evsum_diff
148 logical :: converged_current, converged_last
149 integer :: verbosity_
151 real(real64),
allocatable :: rhoout(:,:), rhoin(:,:)
152 real(real64),
allocatable :: vhxc_old(:,:)
153 class(wfs_elec_t),
allocatable :: psioutb(:, :)
154 logical :: output_forces, calc_current, output_during_scf, finish
160 subroutine scf_init(scf, namespace, gr, ions, st, mc, hm, space)
161 type(scf_t),
intent(inout) :: scf
162 type(grid_t),
intent(in) :: gr
163 type(namespace_t),
intent(in) :: namespace
164 type(ions_t),
intent(in) :: ions
165 type(states_elec_t),
intent(in) :: st
166 type(multicomm_t),
intent(in) :: mc
167 type(hamiltonian_elec_t),
intent(inout) :: hm
168 class(space_t),
intent(in) :: space
171 integer :: mixdefault
172 type(type_t) :: mix_type
173 class(convergence_criterion_t),
pointer :: crit
174 type(criterion_iterator_t) :: iter
197 if (
allocated(hm%vberry))
then
204 call iter%start(scf%criterion_list)
205 do while (iter%has_next())
206 crit => iter%get_next()
209 call crit%set_pointers(scf%energy_diff, scf%energy_in)
211 call crit%set_pointers(scf%abs_dens_diff, st%qtot)
213 call crit%set_pointers(scf%evsum_diff, scf%evsum_out)
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")
226 call messages_write(
"Please set one of the following variables to a positive value:")
249 if(scf%max_iter < 0) scf%max_iter = huge(scf%max_iter)
256 if(scf%eigens%es_type /=
rs_evo)
then
280 mixdefault = option__mixfield__potential
283 call parse_variable(namespace,
'MixField', mixdefault, scf%mix_field)
287 if (scf%mix_field == option__mixfield__potential .and. hm%theory_level ==
independent_particles)
then
288 call messages_write(
'Input: Cannot mix the potential for non-interacting particles.')
292 if (scf%mix_field == option__mixfield__potential .and. hm%pcm%run_pcm)
then
293 call messages_write(
'Input: You have selected to mix the potential.', new_line = .
true.)
294 call messages_write(
' This might produce convergence problems for solvated systems.', new_line = .
true.)
299 if(scf%mix_field == option__mixfield__density &
300 .and.
bitand(hm%xc%family, xc_family_oep + xc_family_mgga + xc_family_hyb_mgga + xc_family_nc_mgga) /= 0)
then
302 call messages_write(
'Input: You have selected to mix the density with OEP or MGGA XC functionals.', new_line = .
true.)
303 call messages_write(
' This might produce convergence problems. Mix the potential instead.')
307 if(scf%mix_field == option__mixfield__states)
then
312 select case(scf%mix_field)
313 case (option__mixfield__potential, option__mixfield__density)
315 case(option__mixfield__states)
322 if (scf%mix_field /= option__mixfield__none)
then
323 call mix_init(scf%smix, namespace, space, gr%der, scf%mixdim1, st%d%nspin, func_type_ = mix_type)
327 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none )
then
333 if(scf%mix_field == option__mixfield__potential)
then
339 scf%mix_field = option__mixfield__none
352 call parse_variable(namespace,
'SCFinLCAO', .false., scf%lcao_restricted)
353 if(scf%lcao_restricted)
then
355 message(1) =
'Info: SCF restricted to LCAO subspace.'
358 if(scf%conv_eigen_error)
then
359 message(1) =
"ConvEigenError cannot be used with SCFinLCAO, since error is unknown."
374 call parse_variable(namespace,
'SCFCalculateForces', .not. ions%only_user_def, scf%calc_force)
376 if(scf%calc_force .and. gr%der%boundaries%spiralBC)
then
377 message(1) =
'Forces cannot be calculated when using spiral boundary conditions.'
378 write(
message(2),
'(a)')
'Please use SCFCalculateForces = no.'
381 if(scf%calc_force)
then
382 if (
allocated(hm%ep%b_field) .or.
allocated(hm%ep%a_static))
then
383 write(
message(1),
'(a)')
'The forces are currently not properly calculated if static'
384 write(
message(2),
'(a)')
'magnetic fields or static vector potentials are present.'
385 write(
message(3),
'(a)')
'Please use SCFCalculateForces = no.'
398 call parse_variable(namespace,
'SCFCalculateStress', .false. , scf%calc_stress)
412 call parse_variable(namespace,
'SCFCalculateDipole', .not. space%is_periodic(), scf%calc_dipole)
413 if (
allocated(hm%vberry)) scf%calc_dipole = .
true.
423 call parse_variable(namespace,
'SCFCalculatePartialCharges', .false., scf%calc_partial_charges)
424 if (scf%calc_partial_charges)
call messages_experimental(
'SCFCalculatePartialCharges', namespace=namespace)
426 rmin = ions%min_distance()
442 scf%forced_finish = .false.
450 type(
scf_t),
intent(inout) :: scf
459 if(scf%mix_field /= option__mixfield__none)
call mix_end(scf%smix)
461 nullify(scf%mixfield)
463 if (scf%mix_field /= option__mixfield__states)
then
468 call iter%start(scf%criterion_list)
469 do while (iter%has_next())
470 crit => iter%get_next()
471 safe_deallocate_p(crit)
480 type(
scf_t),
intent(inout) :: scf
486 if (scf%mix_field /= option__mixfield__states)
then
496 subroutine scf_load(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, restart_load)
497 type(
scf_t),
intent(inout) :: scf
501 type(
grid_t),
intent(inout) :: gr
502 type(
ions_t),
intent(in) :: ions
505 type(
v_ks_t),
intent(inout) :: ks
507 type(
restart_t),
intent(in) :: restart_load
517 message(1) =
'Unable to read density. Density will be calculated from states.'
520 if (
bitand(ks%xc_family, xc_family_oep) == 0)
then
521 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
524 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
531 call hm%ks_pot%load(restart_load, space, gr, ierr)
533 message(1) =
'Unable to read Vhxc. Vhxc will be calculated from states.'
536 call hm%update(gr, namespace, space, ext_partners)
537 if (
bitand(ks%xc_family, xc_family_oep) /= 0)
then
539 do is = 1, st%d%nspin
540 ks%oep%vxc(1:gr%np, is) = hm%ks_pot%vhxc(1:gr%np, is) - hm%ks_pot%vhartree(1:gr%np)
542 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
549 if (scf%mix_field == option__mixfield__density .or. scf%mix_field == option__mixfield__potential)
then
550 call mix_load(namespace, restart_load, scf%smix, space, gr, ierr)
553 message(1) =
"Unable to read mixing information. Mixing will start from scratch."
559 call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr)
561 message(1) =
"Unable to read DFT+U information. DFT+U data will be calculated from states."
578 subroutine scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
579 type(
scf_t),
intent(inout) :: scf
582 type(
grid_t),
intent(inout) :: gr
583 type(
ions_t),
intent(inout) :: ions
585 type(
v_ks_t),
intent(inout) :: ks
587 type(
output_t),
optional,
intent(in) :: outp
588 integer,
optional,
intent(in) :: verbosity
594 if(scf%forced_finish)
then
595 message(1) =
"Previous clean stop, not doing SCF and quitting."
601 scf%output_during_scf = .false.
602 scf%output_forces = .false.
603 scf%calc_current = .false.
605 if (
present(outp))
then
608 if (outp%what(option__output__stress))
then
609 scf%calc_stress = .
true.
612 scf%output_during_scf = outp%duringscf
615 if (outp%duringscf .and. outp%what(option__output__forces))
then
616 scf%output_forces = .
true.
620 if(scf%lcao_restricted)
then
621 call lcao_init(scf%lcao, namespace, space, gr, ions, st)
623 message(1) =
'LCAO is not available. Cannot do SCF in LCAO.'
628 safe_allocate(scf%rhoout(1:gr%np, 1:st%d%nspin))
629 safe_allocate(scf%rhoin (1:gr%np, 1:st%d%nspin))
631 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
634 if (scf%calc_force .or. scf%output_forces)
then
636 safe_allocate(scf%vhxc_old(1:gr%np, 1:st%d%nspin))
637 call lalg_copy(gr%np, st%d%nspin, hm%ks_pot%vhxc, scf%vhxc_old)
641 select case(scf%mix_field)
642 case(option__mixfield__potential)
645 case(option__mixfield__density)
648 case(option__mixfield__states)
651 allocate(
wfs_elec_t::scf%psioutb (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
653 do iqn = st%d%kpt%start, st%d%kpt%end
654 do ib = st%group%block_start, st%group%block_end
655 call st%group%psib(ib, iqn)%copy_to(scf%psioutb(ib, iqn))
663 if (scf%mix_field /= option__mixfield__states)
call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
667 if ( scf%verbosity_ /= verb_no )
then
668 if(scf%max_iter > 0)
then
669 write(
message(1),
'(a)')
'Info: Starting SCF iteration.'
671 write(
message(1),
'(a)')
'Info: No SCF iterations will be done.'
678 scf%converged_current = .false.
688 character(len=*),
intent(in) :: dir
689 character(len=*),
intent(in) :: fname
692 character(len=12) :: label
695 iunit =
io_open(trim(dir) //
"/" // trim(fname), namespace, action=
'write')
696 write(iunit,
'(a)', advance =
'no')
'#iter energy '
697 label =
'energy_diff'
698 write(iunit,
'(1x,a)', advance =
'no') label
700 write(iunit,
'(1x,a)', advance =
'no') label
702 write(iunit,
'(1x,a)', advance =
'no') label
704 write(iunit,
'(1x,a)', advance =
'no') label
706 write(iunit,
'(1x,a)', advance =
'no') label
710 label =
'OEP norm2ss'
711 write(iunit,
'(1x,a)', advance =
'no') label
714 write(iunit,
'(a)')
''
724 subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
725 verbosity, iters_done, restart_dump)
726 type(
scf_t),
intent(inout) :: scf
730 type(
grid_t),
intent(inout) :: gr
731 type(
ions_t),
intent(inout) :: ions
734 type(
v_ks_t),
intent(inout) :: ks
736 type(
output_t),
optional,
intent(in) :: outp
737 integer,
optional,
intent(in) :: verbosity
738 integer,
optional,
intent(out) :: iters_done
739 type(
restart_t),
optional,
intent(in) :: restart_dump
746 call scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
749 do iter = 1, scf%max_iter
751 call scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
752 verbosity, iters_done, restart_dump)
754 completed =
scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
756 if(scf%forced_finish .or. completed)
then
761 call scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
762 verbosity, iters_done, restart_dump)
768 subroutine scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
769 verbosity, iters_done, restart_dump)
770 type(
scf_t),
intent(inout) :: scf
774 type(
grid_t),
intent(inout) :: gr
775 type(
ions_t),
intent(inout) :: ions
778 type(
v_ks_t),
intent(inout) :: ks
780 integer,
intent(in) :: iter
781 type(
output_t),
optional,
intent(in) :: outp
782 integer,
optional,
intent(in) :: verbosity
783 integer,
optional,
intent(out) :: iters_done
784 type(
restart_t),
optional,
intent(in) :: restart_dump
786 integer :: iqn, ib, ierr
789 logical :: is_crit_conv
790 real(real64) :: etime, itime
800 scf%eigens%converged = 0
803 call hm%update_span(gr%spacing(1:space%dim), minval(st%eigenval(:, :)), namespace)
809 call iterator%start(scf%criterion_list)
810 do while (iterator%has_next())
811 crit => iterator%get_next()
815 if (scf%calc_force .or. scf%output_forces)
then
817 scf%vhxc_old(1:gr%np, 1:st%d%nspin) = hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin)
820 if(scf%lcao_restricted)
then
822 call lcao_wf(scf%lcao, st, gr, ions, hm, namespace)
827 if (
allocated(hm%vberry))
then
830 ks%frozen_hxc = .
true.
832 call berry_perform_internal_scf(scf%berry, namespace, space, scf%eigens, gr, st, hm, iter, ks, ions, ext_partners)
834 ks%frozen_hxc = .false.
836 scf%eigens%converged = 0
837 call scf%eigens%run(namespace, gr, st, hm, iter)
841 scf%matvec = scf%matvec + scf%eigens%matvec
850 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoout)
852 select case (scf%mix_field)
853 case (option__mixfield__potential)
854 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
857 case (option__mixfield__density)
859 case(option__mixfield__states)
861 do iqn = st%d%kpt%start, st%d%kpt%end
862 do ib = st%group%block_start, st%group%block_end
863 call st%group%psib(ib, iqn)%copy_data_to(gr%np, scf%psioutb(ib, iqn))
868 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none)
then
875 if (
present(outp))
then
877 if (outp%duringscf .and. outp%what_now(option__output__forces, iter))
then
878 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
883 call iterator%start(scf%criterion_list)
884 do while (iterator%has_next())
885 crit => iterator%get_next()
890 scf%converged_last = scf%converged_current
892 scf%converged_current = scf%check_conv .and. &
893 (.not. scf%conv_eigen_error .or. all(scf%eigens%converged == st%nst))
895 call iterator%start(scf%criterion_list)
896 do while (iterator%has_next())
897 crit => iterator%get_next()
898 call crit%is_converged(is_crit_conv)
899 scf%converged_current = scf%converged_current .and. is_crit_conv
904 scf%finish = scf%converged_last .and. scf%converged_current
910 select case (scf%mix_field)
911 case (option__mixfield__density)
913 call mixing(namespace, scf%smix)
916 if (minval(st%rho(1:gr%np, 1:st%d%spin_channels)) < -1e-6_real64)
then
917 write(
message(1),*)
'Negative density after mixing. Minimum value = ', &
918 minval(st%rho(1:gr%np, 1:st%d%spin_channels))
922 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
923 case (option__mixfield__potential)
925 call mixing(namespace, scf%smix)
931 case(option__mixfield__states)
933 do iqn = st%d%kpt%start, st%d%kpt%end
934 do ib = st%group%block_start, st%group%block_end
941 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
943 case (option__mixfield__none)
944 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
951 if (scf%finish .and. st%modelmbparticles%nparticle > 0)
then
955 if (
present(outp) .and.
present(restart_dump))
then
958 if ( (scf%finish .or. (modulo(iter, outp%restart_write_interval) == 0) &
959 .or. iter == scf%max_iter .or. scf%forced_finish) )
then
961 call states_elec_dump(scf%restart_dump, space, st, gr, hm%kpoints, ierr, iter=iter)
963 message(1) =
'Unable to write states wavefunctions.'
969 message(1) =
'Unable to write density.'
974 call lda_u_dump(scf%restart_dump, namespace, hm%lda_u, st, gr, ierr)
976 message(1) =
'Unable to write DFT+U information.'
981 select case (scf%mix_field)
982 case (option__mixfield__density)
983 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
985 message(1) =
'Unable to write mixing information.'
988 case (option__mixfield__potential)
989 call hm%ks_pot%dump(scf%restart_dump, space, gr, ierr)
991 message(1) =
'Unable to write Vhxc.'
995 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
997 message(1) =
'Unable to write mixing information.'
1015 character(len=50) :: str
1016 real(real64) :: dipole(1:space%dim)
1022 write(str,
'(a,i5)')
'SCF CYCLE ITER #' ,iter
1026 ' rel_ev = ', scf%evsum_diff/(abs(scf%evsum_out)+1e-20)
1027 write(
message(2),
'(a,es15.2,2(a,es9.2))') &
1028 ' ediff = ', scf%energy_diff,
' abs_dens = ', scf%abs_dens_diff, &
1029 ' rel_dens = ', scf%abs_dens_diff/st%qtot
1032 if(.not.scf%lcao_restricted)
then
1033 write(
message(1),
'(a,i6)')
'Matrix vector products: ', scf%eigens%matvec
1034 write(
message(2),
'(a,i6)')
'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%nik))
1041 if (
allocated(hm%vberry))
then
1043 call write_dipole(st, hm, space, dipole, namespace=namespace)
1056 write(
message(2),
'(a,i5,a,f14.2)')
'Elapsed time for SCF step ', iter,
':', etime
1066 write(
message(1),
'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') &
1069 ' : abs_dens', scf%abs_dens_diff, &
1070 ' : etime ', etime,
's'
1080 character(len=*),
intent(in) :: dir
1081 character(len=*),
intent(in) :: fname
1087 iunit =
io_open(trim(dir) //
"/" // trim(fname), namespace, action=
'write', position=
'append')
1089 call iterator%start(scf%criterion_list)
1090 do while (iterator%has_next())
1091 crit => iterator%get_next()
1096 write(iunit,
'(2es13.5)', advance =
'no') crit%val_abs, crit%val_rel
1099 write(iunit,
'(es13.5)', advance =
'no') crit%val_rel
1107 write(iunit,
'(es13.5)', advance =
'no') ks%oep%norm2ss
1110 write(iunit,
'(a)')
''
1117 logical function scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, &
1118 verbosity, iters_done)
result(completed)
1119 type(
scf_t),
intent(inout) :: scf
1122 type(
grid_t),
intent(inout) :: gr
1123 type(
ions_t),
intent(inout) :: ions
1125 type(
v_ks_t),
intent(inout) :: ks
1127 integer,
intent(in) :: iter
1128 type(
output_t),
optional,
intent(in) :: outp
1129 integer,
optional,
intent(in) :: verbosity
1130 integer,
optional,
intent(out) :: iters_done
1132 character(len=MAX_PATH_LEN) :: dirname
1133 integer(int64) :: what_i
1140 if(
present(iters_done)) iters_done = iter
1142 write(
message(1),
'(a, i4, a)')
'Info: SCF converged in ', iter,
' iterations'
1150 if (
present(outp))
then
1151 if (any(outp%what) .and. outp%duringscf)
then
1152 do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
1153 if (outp%what_now(what_i, iter))
then
1154 write(dirname,
'(a,a,i4.4)') trim(outp%iter_dir),
"scf.", iter
1155 call output_all(outp, namespace, space, dirname, gr, ions, iter, st, hm, ks)
1156 call output_modelmb(outp, namespace, space, dirname, gr, ions, iter, st)
1164 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
1167 if (scf%mix_field /= option__mixfield__none)
then
1168 if (scf%smix%ns_restart > 0)
then
1169 if (mod(iter, scf%smix%ns_restart) == 0)
then
1170 message(1) =
"Info: restarting mixing."
1177 select case(scf%mix_field)
1178 case(option__mixfield__potential)
1181 case (option__mixfield__density)
1186 if (scf%mix_field /= option__mixfield__states)
then
1197 subroutine scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
1198 verbosity, iters_done, restart_dump)
1199 type(
scf_t),
intent(inout) :: scf
1203 type(
grid_t),
intent(inout) :: gr
1204 type(
ions_t),
intent(inout) :: ions
1207 type(
v_ks_t),
intent(inout) :: ks
1209 integer,
intent(in) :: iter
1211 integer,
optional,
intent(in) :: verbosity
1212 integer,
optional,
intent(out) :: iters_done
1213 type(
restart_t),
optional,
intent(in) :: restart_dump
1222 if(scf%lcao_restricted)
call lcao_end(scf%lcao)
1224 if ((scf%max_iter > 0 .and. scf%mix_field == option__mixfield__potential) .or. scf%calc_current)
then
1225 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
1226 calc_current=scf%calc_current)
1229 select case(scf%mix_field)
1230 case(option__mixfield__states)
1232 do iqn = st%d%kpt%start, st%d%kpt%end
1233 do ib = st%group%block_start, st%group%block_end
1234 call scf%psioutb(ib, iqn)%end()
1239 deallocate(scf%psioutb)
1242 safe_deallocate_a(scf%rhoout)
1243 safe_deallocate_a(scf%rhoin)
1245 if (scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted)
then
1246 write(
message(1),
'(a)')
'Some of the states are not fully converged!'
1250 if (.not.scf%finish)
then
1251 write(
message(1),
'(a,i4,a)')
'SCF *not* converged after ', iter - 1,
' iterations.'
1255 write(
message(1),
'(a,i10)')
'Info: Number of matrix-vector products: ', scf%matvec
1258 if (scf%calc_force)
then
1259 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
1262 if (scf%calc_stress)
call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
1264 if(scf%max_iter == 0)
then
1270 if(
present(outp))
then
1277 if (space%is_periodic() .and. st%nik > st%d%nspin)
then
1280 ions, gr, hm%kpoints, hm%phase, vec_pot = hm%hm_base%uniform_vector_potential, &
1281 vec_pot_var = hm%hm_base%vector_potential)
1285 if (ks%vdw%vdw_correction == option__vdwcorrection__vdw_ts)
then
1289 safe_deallocate_a(scf%vhxc_old)
1297 character(len=*),
intent(in) :: dir, fname
1299 integer :: iunit, iatom
1300 real(real64),
allocatable :: hirshfeld_charges(:)
1301 real(real64) :: dipole(1:space%dim)
1302 real(real64) :: ex_virial
1308 iunit =
io_open(trim(dir) //
"/" // trim(fname), namespace, action=
'write')
1314 if (space%is_periodic())
then
1315 call hm%kpoints%write_info(iunit=iunit)
1323 write(iunit,
'(a, i4, a)')
'SCF converged in ', iter,
' iterations'
1325 write(iunit,
'(a)')
'SCF *not* converged!'
1327 write(iunit,
'(1x)')
1329 if(any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted)
then
1330 write(iunit,
'(a)')
'Some of the states are not fully converged!'
1334 write(iunit,
'(1x)')
1336 if (space%is_periodic())
then
1338 write(iunit,
'(1x)')
1354 if(st%d%ispin ==
spinors .and. space%dim == 3 .and. &
1366 if(scf%calc_dipole)
then
1373 hm%xc%functional(
func_c,1)%family == xc_family_none .and. st%d%ispin /=
spinors)
then
1380 write(iunit,
'(1x)')
1385 if(scf%max_iter > 0)
then
1386 write(iunit,
'(a)')
'Convergence:'
1387 call iterator%start(scf%criterion_list)
1388 do while (iterator%has_next())
1389 crit => iterator%get_next()
1390 call crit%write_info(iunit)
1399 write(iunit,
'(a)')
'Photon observables:'
1400 write(iunit,
'(6x, a, es15.8,a,es15.8,a)')
'Photon number = ', ks%oep_photon%pt%number(1)
1401 write(iunit,
'(6x, a, es15.8,a,es15.8,a)')
'Photon ex. = ', ks%oep_photon%pt%ex
1408 if (scf%calc_stress)
then
1409 call output_stress(iunit, space%periodic_dim, st%stress_tensors, all_terms=.false.)
1410 call output_pressure(iunit, space%periodic_dim, st%stress_tensors%total)
1415 if(scf%calc_partial_charges)
then
1416 safe_allocate(hirshfeld_charges(1:ions%natoms))
1422 write(iunit,
'(a)')
'Partial ionic charges'
1423 write(iunit,
'(a)')
' Ion Hirshfeld'
1425 do iatom = 1, ions%natoms
1426 write(iunit,
'(i4,a10,f16.3)') iatom, trim(ions%atom(iatom)%species%get_label()), hirshfeld_charges(iatom)
1432 safe_deallocate_a(hirshfeld_charges)
1467 real(real64) :: mem_tmp
1471 if(
conf%report_memory)
then
1473 call mpi_world%allreduce(mem, mem_tmp, 1, mpi_double_precision, mpi_sum)
1475 write(
message(1),
'(a,f14.2)')
'Memory usage [Mbytes] :', mem
1485 type(
scf_t),
intent(inout) :: scf
1491 select type (criterion)
1493 scf%energy_in = hm%energy%total
1508 type(
scf_t),
intent(inout) :: scf
1511 type(
grid_t),
intent(in) :: gr
1512 real(real64),
intent(in) :: rhoout(:,:), rhoin(:,:)
1516 real(real64),
allocatable :: tmp(:)
1520 select type (criterion)
1522 scf%energy_diff = abs(hm%energy%total - scf%energy_in)
1525 scf%abs_dens_diff =
m_zero
1526 safe_allocate(tmp(1:gr%np))
1527 do is = 1, st%d%nspin
1528 tmp = abs(rhoin(1:gr%np, is) - rhoout(1:gr%np, is))
1529 scf%abs_dens_diff = scf%abs_dens_diff +
dmf_integrate(gr, tmp)
1531 safe_deallocate_a(tmp)
1535 scf%evsum_diff = abs(scf%evsum_out - scf%evsum_in)
1536 scf%evsum_in = scf%evsum_out
1546 subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
1550 real(real64),
intent(in) :: dipole(:)
1551 integer,
optional,
intent(in) :: iunit
1552 type(
namespace_t),
optional,
intent(in) :: namespace
1559 if (space%is_periodic())
then
1560 message(1) =
"Defined only up to quantum of polarization (e * lattice vector)."
1561 message(2) =
"Single-point Berry's phase method only accurate for large supercells."
1564 if (hm%kpoints%full%npoints > 1)
then
1566 "WARNING: Single-point Berry's phase method for dipole should not be used when there is more than one k-point."
1567 message(2) =
"Instead, finite differences on k-points (not yet implemented) are needed."
1572 message(1) =
"Single-point Berry's phase dipole calculation not correct without integer occupations."
batchified version of the BLAS axpy routine:
scale a batch by a constant or vector
Copies a vector x, to a vector y.
This module implements common operations on batches of mesh functions.
subroutine, public berry_perform_internal_scf(this, namespace, space, eigensolver, gr, st, hm, iter, ks, ions, ext_partners)
subroutine, public berry_init(this, namespace)
subroutine, public calc_dipole(dipole, space, mesh, st, ions)
subroutine, public criteria_factory_init(list, namespace, check_conv)
This module implements a calculator for the density and defines related functions.
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
subroutine, public eigensolver_init(eigens, namespace, gr, st, mc, space)
integer, parameter, public rs_evo
subroutine, public eigensolver_end(eigens)
integer, parameter, public unpolarized
Parameters...
integer, parameter, public spinors
subroutine, public energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit, full)
This subroutine calculates the total energy of the system. Basically, it adds up the KS eigenvalues,...
subroutine, public energy_calc_virial_ex(der, vxc, st, ex)
subroutine, public energy_calc_eigenvalues(namespace, hm, der, st)
subroutine, public forces_write_info(iunit, ions, dir, namespace)
subroutine, public forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old, t, dt)
real(real64), parameter, public m_zero
real(real64), parameter, public lmm_r_single_atom
Default local magnetic moments sphere radius for an isolated system.
type(conf_t), public conf
Global instance of Octopus configuration.
character(len= *), parameter, public static_dir
real(real64), parameter, public m_half
real(real64), parameter, public m_one
This module implements the underlying real-space grid.
subroutine, public grid_write_info(gr, iunit, namespace)
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.
subroutine, public io_close(iunit, grp)
subroutine, public io_debug_on_the_fly(namespace)
check if debug mode should be enabled or disabled on the fly
subroutine, public io_mkdir(fname, namespace, parents)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
integer, parameter, public kpoints_path
A module to handle KS potential, without the external potential.
integer, parameter, public hartree_fock
integer, parameter, public independent_particles
integer, parameter, public generalized_kohn_sham_dft
integer, parameter, public kohn_sham_dft
subroutine, public lcao_init_orbitals(this, namespace, st, gr, ions, start)
subroutine, public lcao_wf(this, st, gr, ions, hm, namespace, start)
subroutine, public lcao_end(this)
subroutine, public lcao_init(this, namespace, space, gr, ions, st)
logical function, public lcao_is_available(this)
subroutine, public lda_u_dump(restart, namespace, this, st, mesh, ierr)
subroutine, public lda_u_write_u(this, iunit, namespace)
subroutine, public lda_u_load(restart, this, st, dftu_energy, ierr, occ_only, u_only)
subroutine, public lda_u_write_v(this, iunit, namespace)
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
integer, parameter, public dft_u_acbn0
subroutine, public lda_u_update_occ_matrices(this, namespace, mesh, st, hm_base, phase, energy)
subroutine, public write_magnetic_moments(mesh, st, ions, boundaries, lmm_r, iunit, namespace)
subroutine, public write_total_xc_torque(iunit, mesh, vxc, st)
This module is intended to contain "only mathematical" functions and procedures.
This module defines various routines, operating on mesh functions.
real(real64) function, public dmf_integrate(mesh, ff, mask, reduce)
Integrate a function on the mesh.
This module defines the meshes, which are used in Octopus.
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
character(len=512), private msg
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
subroutine, public messages_new_line()
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
real(real64) pure function, public mix_coefficient(this)
subroutine, public mixing(namespace, smix)
Main entry-point to SCF mixer.
subroutine, public mix_get_field(this, mixfield)
subroutine, public mix_dump(namespace, restart, smix, space, mesh, ierr)
subroutine, public mix_load(namespace, restart, smix, space, mesh, ierr)
subroutine, public mix_init(smix, namespace, space, der, d1, d2, def_, func_type_, prefix_)
Initialise mix_t instance.
subroutine, public mix_end(smix)
subroutine, public mix_clear(smix)
subroutine, public modelmb_sym_all_states(space, mesh, st)
logical function mpi_grp_is_root(grp)
Is the current MPI process of grpcomm, root.
type(mpi_grp_t), public mpi_world
This module handles the communicators for the various parallelization strategies.
this module contains the low-level part of the output system
subroutine, public output_modelmb(outp, namespace, space, dir, gr, ions, iter, st)
this module contains the output system
logical function, public output_needs_current(outp, states_are_real)
subroutine, public output_all(outp, namespace, space, dir, gr, ions, iter, st, hm, ks)
subroutine, public partial_charges_calculate(mesh, st, ions, hirshfeld_charges)
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
logical function, public clean_stop(comm)
returns true if a file named stop exists
integer, parameter, public restart_flag_mix
integer, parameter, public restart_flag_rho
integer, parameter, public restart_flag_vhxc
logical pure function, public restart_has_flag(restart, flag)
Returns true if...
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.
logical function, public scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
subroutine scf_update_initial_quantity(scf, hm, criterion)
Update the quantity at the begining of a SCF cycle.
subroutine scf_update_diff_quantity(scf, hm, st, gr, rhoout, rhoin, criterion)
Update the quantity at the begining of a SCF cycle.
subroutine, public scf_state_info(namespace, st)
subroutine, public scf_print_mem_use(namespace)
subroutine, public scf_mix_clear(scf)
integer, parameter, public verb_full
subroutine, public scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, verbosity, iters_done, restart_dump)
integer, parameter, public verb_compact
subroutine, public scf_init(scf, namespace, gr, ions, st, mc, hm, space)
subroutine, public scf_end(scf)
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.
subroutine, public scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, verbosity, iters_done, restart_dump)
subroutine, public scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
Preparation of the SCF cycle.
logical pure function, public smear_is_semiconducting(this)
pure logical function, public states_are_real(st)
This module defines routines to write information about states.
subroutine, public states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
write the eigenvalues for some states to a file.
subroutine, public states_elec_write_gaps(iunit, st, space)
calculate gaps and write to a file.
subroutine, public states_elec_write_bandstructure(dir, namespace, nst, st, ions, mesh, kpoints, phase, vec_pot, vec_pot_var)
calculate and write the bandstructure
subroutine, public states_elec_fermi(st, namespace, mesh, compute_spin)
calculate the Fermi level for the states in this object
real(real64) function, public states_elec_eigenvalues_sum(st, alt_eig)
function to calculate the eigenvalues sum using occupations as weights
This module handles reading and writing restart information for the states_elec_t.
subroutine, public states_elec_dump(restart, space, st, mesh, kpoints, ierr, iter, lr, st_start_writing, verbose)
subroutine, public states_elec_load_rho(restart, space, st, mesh, ierr)
subroutine, public states_elec_dump_rho(restart, space, st, mesh, ierr, iter)
This module implements the calculation of the stress tensor.
subroutine, public output_pressure(iunit, space_dim, total_stress_tensor)
subroutine, public stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
This computes the total stress on the lattice.
subroutine, public output_stress(iunit, space_dim, stress_tensors, all_terms)
subroutine, public symmetries_write_info(this, space, iunit, namespace)
type(type_t), public type_float
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
character(len=20) pure function, public units_abbrev(this)
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.
subroutine, public output_dipole(dipole, ndim, iunit, namespace)
subroutine, public v_ks_write_info(ks, iunit, namespace)
subroutine, public v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_eigenval, time, calc_energy, calc_current, force_semilocal)
Tkatchenko-Scheffler pairwise method for van der Waals (vdW, dispersion) interactions.
subroutine, public vdw_ts_write_c6ab(this, ions, dir, fname, namespace)
subroutine, public vtau_mixer_end(mixer, smix)
subroutine, public vtau_mixer_init_auxmixer(namespace, mixer, smix, hm, np, nspin)
subroutine, public vtau_mixer_set_vout(mixer, hm)
subroutine, public vtau_mixer_get_vnew(mixer, hm)
subroutine, public vtau_mixer_clear(mixer, smix)
subroutine, public vtau_mixer_set_vin(mixer, hm)
This module provices a simple timer class which can be used to trigger the writing of a restart file ...
logical function, public walltimer_alarm(comm, print)
indicate whether time is up
integer, parameter, public func_c
integer, parameter, public oep_level_full
integer, parameter, public oep_level_kli
subroutine scf_write_static(dir, fname)
subroutine create_convergence_file(dir, fname)
subroutine scf_write_iter(namespace)
subroutine write_convergence_file(dir, fname)
Extension of space that contains the knowledge of the spin dimension.
Description of the grid, containing information on derivatives, stencil, and symmetries.
Stores all communicators and groups.
some variables used for the SCF cycle
abstract class for states
The states_elec_t class contains all electronic wave functions.
batches of electronic states