41 use,
intrinsic :: iso_fortran_env
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
175 logical :: deactivate_oracle
198 if (
allocated(hm%vberry))
then
205 call iter%start(scf%criterion_list)
206 do while (iter%has_next())
207 crit => iter%get_next()
210 call crit%set_pointers(scf%energy_diff, scf%energy_in)
212 call crit%set_pointers(scf%abs_dens_diff, st%qtot)
214 call crit%set_pointers(scf%evsum_diff, scf%evsum_out)
221 if(.not. scf%check_conv .and. scf%max_iter < 0)
then
222 call messages_write(
"All convergence criteria are disabled. Octopus is cowardly refusing")
227 call messages_write(
"Please set one of the following variables to a positive value:")
250 call parse_variable(namespace,
'ConvEigenError', .false., scf%conv_eigen_error)
252 if(scf%max_iter < 0) scf%max_iter = huge(scf%max_iter)
258 call eigensolver_init(scf%eigens, namespace, gr, st, hm, mc, space, deactivate_oracle)
260 if(scf%eigens%es_type /=
rs_evo)
then
284 mixdefault = option__mixfield__potential
287 call parse_variable(namespace,
'MixField', mixdefault, scf%mix_field)
291 if (scf%mix_field == option__mixfield__potential .and. hm%theory_level ==
independent_particles)
then
292 call messages_write(
'Input: Cannot mix the potential for non-interacting particles.')
296 if (scf%mix_field == option__mixfield__potential .and. hm%pcm%run_pcm)
then
297 call messages_write(
'Input: You have selected to mix the potential.', new_line = .
true.)
298 call messages_write(
' This might produce convergence problems for solvated systems.', new_line = .
true.)
303 if(scf%mix_field == option__mixfield__density &
306 call messages_write(
'Input: You have selected to mix the density with OEP or MGGA XC functionals.', new_line = .
true.)
307 call messages_write(
' This might produce convergence problems. Mix the potential instead.')
311 if(scf%mix_field == option__mixfield__states)
then
316 select case(scf%mix_field)
317 case (option__mixfield__potential, option__mixfield__density)
319 case(option__mixfield__states)
326 if (scf%mix_field /= option__mixfield__none)
then
327 call mix_init(scf%smix, namespace, space, gr%der, scf%mixdim1, st%d%nspin, func_type_ = mix_type)
331 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none )
then
337 if(scf%mix_field == option__mixfield__potential)
then
343 scf%mix_field = option__mixfield__none
356 call parse_variable(namespace,
'SCFinLCAO', .false., scf%lcao_restricted)
357 if(scf%lcao_restricted)
then
359 message(1) =
'Info: SCF restricted to LCAO subspace.'
362 if(scf%conv_eigen_error)
then
363 message(1) =
"ConvEigenError cannot be used with SCFinLCAO, since error is unknown."
378 call parse_variable(namespace,
'SCFCalculateForces', .not. ions%only_user_def, scf%calc_force)
380 if(scf%calc_force .and. gr%der%boundaries%spiralBC)
then
381 message(1) =
'Forces cannot be calculated when using spiral boundary conditions.'
382 write(
message(2),
'(a)')
'Please use SCFCalculateForces = no.'
385 if(scf%calc_force)
then
386 if (
allocated(hm%ep%b_field) .or.
allocated(hm%ep%a_static))
then
387 write(
message(1),
'(a)')
'The forces are currently not properly calculated if static'
388 write(
message(2),
'(a)')
'magnetic fields or static vector potentials are present.'
389 write(
message(3),
'(a)')
'Please use SCFCalculateForces = no.'
402 call parse_variable(namespace,
'SCFCalculateStress', .false. , scf%calc_stress)
416 call parse_variable(namespace,
'SCFCalculateDipole', .not. space%is_periodic(), scf%calc_dipole)
417 if (
allocated(hm%vberry)) scf%calc_dipole = .
true.
427 call parse_variable(namespace,
'SCFCalculatePartialCharges', .false., scf%calc_partial_charges)
428 if (scf%calc_partial_charges)
call messages_experimental(
'SCFCalculatePartialCharges', namespace=namespace)
430 rmin = ions%min_distance()
446 scf%forced_finish = .false.
454 type(
scf_t),
intent(inout) :: scf
463 if(scf%mix_field /= option__mixfield__none)
call mix_end(scf%smix)
465 nullify(scf%mixfield)
467 if (scf%mix_field /= option__mixfield__states)
then
472 call iter%start(scf%criterion_list)
473 do while (iter%has_next())
474 crit => iter%get_next()
475 safe_deallocate_p(crit)
484 type(
scf_t),
intent(inout) :: scf
490 if (scf%mix_field /= option__mixfield__states)
then
500 subroutine scf_load(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, restart_load)
501 type(
scf_t),
intent(inout) :: scf
505 type(
grid_t),
intent(inout) :: gr
506 type(
ions_t),
intent(in) :: ions
509 type(
v_ks_t),
intent(inout) :: ks
511 type(
restart_t),
intent(in) :: restart_load
521 message(1) =
'Unable to read density. Density will be calculated from states.'
524 if (
bitand(ks%xc_family, xc_family_oep) == 0)
then
525 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
528 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
535 call hm%ks_pot%load(restart_load, space, gr, ierr)
537 message(1) =
'Unable to read Vhxc. Vhxc will be calculated from states.'
540 call hm%update(gr, namespace, space, ext_partners)
541 if (
bitand(ks%xc_family, xc_family_oep) /= 0)
then
543 do is = 1, st%d%nspin
544 ks%oep%vxc(1:gr%np, is) = hm%ks_pot%vhxc(1:gr%np, is) - hm%ks_pot%vhartree(1:gr%np)
546 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
553 if (scf%mix_field == option__mixfield__density .or. scf%mix_field == option__mixfield__potential)
then
554 call mix_load(namespace, restart_load, scf%smix, space, gr, ierr)
557 message(1) =
"Unable to read mixing information. Mixing will start from scratch."
563 call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr)
565 message(1) =
"Unable to read DFT+U information. DFT+U data will be calculated from states."
585 subroutine scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
586 type(
scf_t),
intent(inout) :: scf
589 type(
grid_t),
intent(inout) :: gr
590 type(
ions_t),
intent(inout) :: ions
592 type(
v_ks_t),
intent(inout) :: ks
594 type(
output_t),
optional,
intent(in) :: outp
595 integer,
optional,
intent(in) :: verbosity
601 if(scf%forced_finish)
then
602 message(1) =
"Previous clean stop, not doing SCF and quitting."
608 scf%output_during_scf = .false.
609 scf%output_forces = .false.
610 scf%calc_current = .false.
612 if (
present(outp))
then
615 if (outp%what(option__output__stress))
then
616 scf%calc_stress = .
true.
619 scf%output_during_scf = outp%duringscf
622 if (outp%duringscf .and. outp%what(option__output__forces))
then
623 scf%output_forces = .
true.
627 if(scf%lcao_restricted)
then
628 call lcao_init(scf%lcao, namespace, space, gr, ions, st, 1)
630 message(1) =
'LCAO is not available. Cannot do SCF in LCAO.'
635 safe_allocate(scf%rhoout(1:gr%np, 1:st%d%nspin))
636 safe_allocate(scf%rhoin (1:gr%np, 1:st%d%nspin))
638 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
641 if (scf%calc_force .or. scf%output_forces)
then
643 safe_allocate(scf%vhxc_old(1:gr%np, 1:st%d%nspin))
644 call lalg_copy(gr%np, st%d%nspin, hm%ks_pot%vhxc, scf%vhxc_old)
648 select case(scf%mix_field)
649 case(option__mixfield__potential)
652 case(option__mixfield__density)
655 case(option__mixfield__states)
658 allocate(
wfs_elec_t::scf%psioutb (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
660 do iqn = st%d%kpt%start, st%d%kpt%end
661 do ib = st%group%block_start, st%group%block_end
662 call st%group%psib(ib, iqn)%copy_to(scf%psioutb(ib, iqn))
670 if (scf%mix_field /= option__mixfield__states)
call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
674 if ( scf%verbosity_ /= verb_no )
then
675 if(scf%max_iter > 0)
then
676 write(
message(1),
'(a)')
'Info: Starting SCF iteration.'
678 write(
message(1),
'(a)')
'Info: No SCF iterations will be done.'
685 scf%converged_current = .false.
695 character(len=*),
intent(in) :: dir
696 character(len=*),
intent(in) :: fname
699 character(len=12) :: label
702 iunit =
io_open(trim(dir) //
"/" // trim(fname), namespace, action=
'write')
703 write(iunit,
'(a)', advance =
'no')
'#iter energy '
704 label =
'energy_diff'
705 write(iunit,
'(1x,a)', advance =
'no') label
707 write(iunit,
'(1x,a)', advance =
'no') label
709 write(iunit,
'(1x,a)', advance =
'no') label
711 write(iunit,
'(1x,a)', advance =
'no') label
713 write(iunit,
'(1x,a)', advance =
'no') label
717 label =
'OEP norm2ss'
718 write(iunit,
'(1x,a)', advance =
'no') label
721 write(iunit,
'(a)')
''
731 subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, &
732 verbosity, iters_done, restart_dump)
733 type(
scf_t),
intent(inout) :: scf
737 type(
grid_t),
intent(inout) :: gr
738 type(
ions_t),
intent(inout) :: ions
741 type(
v_ks_t),
intent(inout) :: ks
743 type(
output_t),
optional,
intent(in) :: outp
744 integer,
optional,
intent(in) :: verbosity
745 integer,
optional,
intent(out) :: iters_done
746 type(
restart_t),
optional,
intent(in) :: restart_dump
753 call scf_start(scf, namespace, space, gr, ions, st, ks, hm, outp, verbosity)
756 do iter = 1, scf%max_iter
758 call scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
759 verbosity, iters_done, restart_dump)
761 completed =
scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, verbosity, iters_done)
763 if(scf%forced_finish .or. completed)
then
768 call scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
769 verbosity, iters_done, restart_dump)
775 subroutine scf_iter(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
776 verbosity, iters_done, restart_dump)
777 type(
scf_t),
intent(inout) :: scf
781 type(
grid_t),
intent(inout) :: gr
782 type(
ions_t),
intent(inout) :: ions
785 type(
v_ks_t),
intent(inout) :: ks
787 integer,
intent(in) :: iter
788 type(
output_t),
optional,
intent(in) :: outp
789 integer,
optional,
intent(in) :: verbosity
790 integer,
optional,
intent(out) :: iters_done
791 type(
restart_t),
optional,
intent(in) :: restart_dump
793 integer :: iqn, ib, ierr
796 logical :: is_crit_conv
797 real(real64) :: etime, itime
807 scf%eigens%converged = 0
810 call hm%update_span(gr%spacing(1:space%dim), minval(st%eigenval(:, :)), namespace)
816 call iterator%start(scf%criterion_list)
817 do while (iterator%has_next())
818 crit => iterator%get_next()
822 if (scf%calc_force .or. scf%output_forces)
then
824 scf%vhxc_old(1:gr%np, 1:st%d%nspin) = hm%ks_pot%vhxc(1:gr%np, 1:st%d%nspin)
827 if(scf%lcao_restricted)
then
829 call lcao_wf(scf%lcao, st, gr, ions, hm, namespace)
834 if (
allocated(hm%vberry))
then
837 ks%frozen_hxc = .
true.
839 call berry_perform_internal_scf(scf%berry, namespace, space, scf%eigens, gr, st, hm, iter, ks, ions, ext_partners)
841 ks%frozen_hxc = .false.
843 scf%eigens%converged = 0
844 call scf%eigens%run(namespace, gr, st, hm, space, ext_partners, iter)
848 scf%matvec = scf%matvec + scf%eigens%matvec
857 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoout)
859 select case (scf%mix_field)
860 case (option__mixfield__potential)
861 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
864 case (option__mixfield__density)
866 case(option__mixfield__states)
868 do iqn = st%d%kpt%start, st%d%kpt%end
869 do ib = st%group%block_start, st%group%block_end
870 call st%group%psib(ib, iqn)%copy_data_to(gr%np, scf%psioutb(ib, iqn))
875 if (scf%mix_field /= option__mixfield__states .and. scf%mix_field /= option__mixfield__none)
then
882 if (
present(outp))
then
884 if (outp%duringscf .and. outp%what_now(option__output__forces, iter))
then
885 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
890 call iterator%start(scf%criterion_list)
891 do while (iterator%has_next())
892 crit => iterator%get_next()
897 scf%converged_last = scf%converged_current
899 scf%converged_current = scf%check_conv .and. &
900 (.not. scf%conv_eigen_error .or. all(scf%eigens%converged >= st%nst_conv))
902 call iterator%start(scf%criterion_list)
903 do while (iterator%has_next())
904 crit => iterator%get_next()
905 call crit%is_converged(is_crit_conv)
906 scf%converged_current = scf%converged_current .and. is_crit_conv
911 scf%finish = scf%converged_last .and. scf%converged_current
917 select case (scf%mix_field)
918 case (option__mixfield__density)
920 call mixing(namespace, scf%smix)
923 if (minval(st%rho(1:gr%np, 1:st%d%spin_channels)) < -1e-6_real64)
then
924 write(
message(1),*)
'Negative density after mixing. Minimum value = ', &
925 minval(st%rho(1:gr%np, 1:st%d%spin_channels))
929 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
930 case (option__mixfield__potential)
932 call mixing(namespace, scf%smix)
938 case(option__mixfield__states)
940 do iqn = st%d%kpt%start, st%d%kpt%end
941 do ib = st%group%block_start, st%group%block_end
948 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
950 case (option__mixfield__none)
951 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=scf%output_during_scf)
958 if (scf%finish .and. st%modelmbparticles%nparticle > 0)
then
962 if (
present(outp) .and.
present(restart_dump))
then
965 if ( (scf%finish .or. (modulo(iter, outp%restart_write_interval) == 0) &
966 .or. iter == scf%max_iter .or. scf%forced_finish) )
then
968 call states_elec_dump(scf%restart_dump, space, st, gr, hm%kpoints, ierr, iter=iter)
970 message(1) =
'Unable to write states wavefunctions.'
976 message(1) =
'Unable to write density.'
981 call lda_u_dump(scf%restart_dump, namespace, hm%lda_u, st, gr, ierr)
983 message(1) =
'Unable to write DFT+U information.'
988 select case (scf%mix_field)
989 case (option__mixfield__density)
990 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
992 message(1) =
'Unable to write mixing information.'
995 case (option__mixfield__potential)
996 call hm%ks_pot%dump(scf%restart_dump, space, gr, ierr)
998 message(1) =
'Unable to write Vhxc.'
1002 call mix_dump(namespace, scf%restart_dump, scf%smix, space, gr, ierr)
1004 message(1) =
'Unable to write mixing information.'
1022 character(len=50) :: str
1023 real(real64) :: dipole(1:space%dim)
1029 write(str,
'(a,i5)')
'SCF CYCLE ITER #' ,iter
1033 ' rel_ev = ', scf%evsum_diff/(abs(scf%evsum_out)+1e-20)
1034 write(
message(2),
'(a,es15.2,2(a,es9.2))') &
1035 ' ediff = ', scf%energy_diff,
' abs_dens = ', scf%abs_dens_diff, &
1036 ' rel_dens = ', scf%abs_dens_diff/st%qtot
1039 if(.not.scf%lcao_restricted)
then
1040 write(
message(1),
'(a,i6)')
'Matrix vector products: ', scf%eigens%matvec
1041 write(
message(2),
'(a,i6)')
'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%nik))
1048 if (
allocated(hm%vberry))
then
1050 call write_dipole(st, hm, space, dipole, namespace=namespace)
1063 write(
message(2),
'(a,i5,a,f14.2)')
'Elapsed time for SCF step ', iter,
':', etime
1073 write(
message(1),
'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') &
1076 ' : abs_dens', scf%abs_dens_diff, &
1077 ' : etime ', etime,
's'
1087 character(len=*),
intent(in) :: dir
1088 character(len=*),
intent(in) :: fname
1094 iunit =
io_open(trim(dir) //
"/" // trim(fname), namespace, action=
'write', position=
'append')
1096 call iterator%start(scf%criterion_list)
1097 do while (iterator%has_next())
1098 crit => iterator%get_next()
1103 write(iunit,
'(2es13.5)', advance =
'no') crit%val_abs, crit%val_rel
1106 write(iunit,
'(es13.5)', advance =
'no') crit%val_rel
1114 write(iunit,
'(es13.5)', advance =
'no') ks%oep%norm2ss
1117 write(iunit,
'(a)')
''
1124 logical function scf_iter_finish(scf, namespace, space, gr, ions, st, ks, hm, iter, outp, &
1125 verbosity, iters_done)
result(completed)
1126 type(
scf_t),
intent(inout) :: scf
1129 type(
grid_t),
intent(inout) :: gr
1130 type(
ions_t),
intent(inout) :: ions
1132 type(
v_ks_t),
intent(inout) :: ks
1134 integer,
intent(in) :: iter
1135 type(
output_t),
optional,
intent(in) :: outp
1136 integer,
optional,
intent(in) :: verbosity
1137 integer,
optional,
intent(out) :: iters_done
1139 character(len=MAX_PATH_LEN) :: dirname
1140 integer(int64) :: what_i
1147 if(
present(iters_done)) iters_done = iter
1149 write(
message(1),
'(a, i4, a)')
'Info: SCF converged in ', iter,
' iterations'
1157 if (
present(outp))
then
1158 if (any(outp%what) .and. outp%duringscf)
then
1159 do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
1160 if (outp%what_now(what_i, iter))
then
1161 write(dirname,
'(a,a,i4.4)') trim(outp%iter_dir),
"scf.", iter
1162 call output_all(outp, namespace, space, dirname, gr, ions, iter, st, hm, ks)
1163 call output_modelmb(outp, namespace, space, dirname, gr, ions, iter, st)
1171 call lalg_copy(gr%np, st%d%nspin, st%rho, scf%rhoin)
1174 if (scf%mix_field /= option__mixfield__none)
then
1175 if (scf%smix%ns_restart > 0)
then
1176 if (mod(iter, scf%smix%ns_restart) == 0)
then
1177 message(1) =
"Info: restarting mixing."
1184 select case(scf%mix_field)
1185 case(option__mixfield__potential)
1188 case (option__mixfield__density)
1193 if (scf%mix_field /= option__mixfield__states)
then
1204 subroutine scf_finish(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, iter, outp, &
1205 verbosity, iters_done, restart_dump)
1206 type(
scf_t),
intent(inout) :: scf
1210 type(
grid_t),
intent(inout) :: gr
1211 type(
ions_t),
intent(inout) :: ions
1214 type(
v_ks_t),
intent(inout) :: ks
1216 integer,
intent(in) :: iter
1218 integer,
optional,
intent(in) :: verbosity
1219 integer,
optional,
intent(out) :: iters_done
1220 type(
restart_t),
optional,
intent(in) :: restart_dump
1229 if(scf%lcao_restricted)
call lcao_end(scf%lcao)
1233 if ((scf%max_iter > 0 .and. scf%mix_field == option__mixfield__potential) .or. scf%calc_current)
then
1234 call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
1235 calc_current=scf%calc_current)
1238 select case(scf%mix_field)
1239 case(option__mixfield__states)
1241 do iqn = st%d%kpt%start, st%d%kpt%end
1242 do ib = st%group%block_start, st%group%block_end
1243 call scf%psioutb(ib, iqn)%end()
1248 deallocate(scf%psioutb)
1251 safe_deallocate_a(scf%rhoout)
1252 safe_deallocate_a(scf%rhoin)
1254 if (scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted)
then
1255 write(
message(1),
'(a)')
'Some of the states are not fully converged!'
1256 if (all(scf%eigens%converged >= st%nst_conv))
then
1257 write(
message(2),
'(a)')
'But all requested states to converge are converged.'
1264 if (.not.scf%finish)
then
1265 write(
message(1),
'(a,i4,a)')
'SCF *not* converged after ', iter - 1,
' iterations.'
1269 write(
message(1),
'(a,i10)')
'Info: Number of matrix-vector products: ', scf%matvec
1272 if (scf%calc_force)
then
1273 call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=scf%vhxc_old)
1276 if (scf%calc_stress)
call stress_calculate(namespace, gr, hm, st, ions, ks, ext_partners)
1279 if (scf%mix_field == option__mixfield__potential)
then
1285 if(
present(outp))
then
1292 if (space%is_periodic() .and. st%nik > st%d%nspin)
then
1295 ions, gr, hm%kpoints, hm%phase, vec_pot = hm%hm_base%uniform_vector_potential, &
1296 vec_pot_var = hm%hm_base%vector_potential)
1300 if (ks%vdw%vdw_correction == option__vdwcorrection__vdw_ts)
then
1304 safe_deallocate_a(scf%vhxc_old)
1312 character(len=*),
intent(in) :: dir, fname
1314 integer :: iunit, iatom
1315 real(real64),
allocatable :: hirshfeld_charges(:)
1316 real(real64) :: dipole(1:space%dim)
1317 real(real64) :: ex_virial
1323 iunit =
io_open(trim(dir) //
"/" // trim(fname), namespace, action=
'write')
1329 if (space%is_periodic())
then
1330 call hm%kpoints%write_info(iunit=iunit)
1338 write(iunit,
'(a, i4, a)')
'SCF converged in ', iter,
' iterations'
1340 write(iunit,
'(a)')
'SCF *not* converged!'
1342 write(iunit,
'(1x)')
1344 if(any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted)
then
1345 write(iunit,
'(a)')
'Some of the states are not fully converged!'
1346 if (all(scf%eigens%converged >= st%nst_conv))
then
1347 write(iunit,
'(a)')
'But all requested states to converge are converged.'
1352 write(iunit,
'(1x)')
1354 if (space%is_periodic())
then
1356 write(iunit,
'(1x)')
1372 if(st%d%ispin ==
spinors .and. space%dim == 3 .and. &
1384 if(scf%calc_dipole)
then
1391 hm%xc%functional(
func_c,1)%family == xc_family_none .and. st%d%ispin /=
spinors)
then
1398 write(iunit,
'(1x)')
1403 if(scf%max_iter > 0)
then
1404 write(iunit,
'(a)')
'Convergence:'
1405 call iterator%start(scf%criterion_list)
1406 do while (iterator%has_next())
1407 crit => iterator%get_next()
1408 call crit%write_info(iunit)
1417 write(iunit,
'(a)')
'Photon observables:'
1418 write(iunit,
'(6x, a, es15.8,a,es15.8,a)')
'Photon number = ', ks%oep_photon%pt%number(1)
1419 write(iunit,
'(6x, a, es15.8,a,es15.8,a)')
'Photon ex. = ', ks%oep_photon%pt%ex
1426 if (scf%calc_stress)
then
1427 call output_stress(iunit, space%periodic_dim, st%stress_tensors, all_terms=.false.)
1428 call output_pressure(iunit, space%periodic_dim, st%stress_tensors%total)
1433 if(scf%calc_partial_charges)
then
1434 safe_allocate(hirshfeld_charges(1:ions%natoms))
1440 write(iunit,
'(a)')
'Partial ionic charges'
1441 write(iunit,
'(a)')
' Ion Hirshfeld'
1443 do iatom = 1, ions%natoms
1444 write(iunit,
'(i4,a10,f16.3)') iatom, trim(ions%atom(iatom)%species%get_label()), hirshfeld_charges(iatom)
1450 safe_deallocate_a(hirshfeld_charges)
1485 real(real64) :: mem_tmp
1489 if(
conf%report_memory)
then
1491 call mpi_world%allreduce(mem, mem_tmp, 1, mpi_double_precision, mpi_sum)
1493 write(
message(1),
'(a,f14.2)')
'Memory usage [Mbytes] :', mem
1503 type(
scf_t),
intent(inout) :: scf
1509 select type (criterion)
1511 scf%energy_in = hm%energy%total
1526 type(
scf_t),
intent(inout) :: scf
1529 type(
grid_t),
intent(in) :: gr
1530 real(real64),
intent(in) :: rhoout(:,:), rhoin(:,:)
1534 real(real64),
allocatable :: tmp(:)
1538 select type (criterion)
1540 scf%energy_diff = abs(hm%energy%total - scf%energy_in)
1543 scf%abs_dens_diff =
m_zero
1544 safe_allocate(tmp(1:gr%np))
1545 do is = 1, st%d%nspin
1546 tmp(:) = abs(rhoin(1:gr%np, is) - rhoout(1:gr%np, is))
1547 scf%abs_dens_diff = scf%abs_dens_diff +
dmf_integrate(gr, tmp)
1549 safe_deallocate_a(tmp)
1553 scf%evsum_diff = abs(scf%evsum_out - scf%evsum_in)
1554 scf%evsum_in = scf%evsum_out
1564 subroutine write_dipole(st, hm, space, dipole, iunit, namespace)
1568 real(real64),
intent(in) :: dipole(:)
1569 integer,
optional,
intent(in) :: iunit
1570 type(
namespace_t),
optional,
intent(in) :: namespace
1577 if (space%is_periodic())
then
1578 message(1) =
"Defined only up to quantum of polarization (e * lattice vector)."
1579 message(2) =
"Single-point Berry's phase method only accurate for large supercells."
1582 if (hm%kpoints%full%npoints > 1)
then
1584 "WARNING: Single-point Berry's phase method for dipole should not be used when there is more than one k-point."
1585 message(2) =
"Instead, finite differences on k-points (not yet implemented) are needed."
1590 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.
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)
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, st_start)
logical function, public lcao_is_available(this)
Returns true if LCAO can be done.
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.
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_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.
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 xc_family_nc_mgga
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