The Self-consistency cycle
Definition of scf_t
type scf_t
private
integer, public :: max_iter !< maximum number of SCF iterations
FLOAT, public :: lmm_r
! several convergence criteria
logical :: conv_eigen_error
logical :: check_conv
integer :: mix_field
logical :: lcao_restricted
logical :: calc_force
logical :: calc_stress
logical :: calc_dipole
logical :: calc_partial_charges
type(mix_t) :: smix
type(mixfield_t), pointer :: mixfield
type(eigensolver_t) :: eigens
integer :: mixdim1
logical :: forced_finish !< remember if 'touch stop' was triggered earlier.
type(lda_u_mixer_t) :: lda_u_mix
type(berry_t) :: berry
integer :: matvec !< number matrix-vector products
type(criterion_list_t), public :: criterion_list
FLOAT :: energy_in, energy_diff, abs_dens_diff, evsum_in, evsum_out, evsum_diff
end type scf_t
Definition of scf_run()
subroutine scf_run(scf, namespace, space, mc, gr, ions, ext_partners, st, ks, hm, outp, gs_run, &
verbosity, iters_done, restart_load, restart_dump)
type(scf_t), intent(inout) :: scf !< self consistent cycle
type(namespace_t), intent(in) :: namespace
type(space_t), intent(in) :: space
type(multicomm_t), intent(in) :: mc
type(grid_t), intent(inout) :: gr !< grid
type(ions_t), intent(inout) :: ions !< geometry
type(partner_list_t), intent(in) :: ext_partners
type(states_elec_t), intent(inout) :: st !< States
type(v_ks_t), intent(inout) :: ks !< Kohn-Sham
type(hamiltonian_elec_t), intent(inout) :: hm !< Hamiltonian
type(output_t), intent(in) :: outp
logical, optional, intent(in) :: gs_run
integer, optional, intent(in) :: verbosity
integer, optional, intent(out) :: iters_done
type(restart_t), optional, intent(in) :: restart_load
type(restart_t), optional, intent(in) :: restart_dump
logical :: finish, converged_current, converged_last, gs_run_
integer :: iter, is, nspin, ierr, verbosity_, ib, iqn
integer(i8) :: what_i
FLOAT :: etime, itime
character(len=MAX_PATH_LEN) :: dirname
type(lcao_t) :: lcao !< Linear combination of atomic orbitals
type(profile_t), save :: prof
FLOAT, allocatable :: rhoout(:,:,:), rhoin(:,:,:)
FLOAT, allocatable :: vhxc_old(:,:)
class(wfs_elec_t), allocatable :: psioutb(:, :)
class(convergence_criterion_t), pointer :: crit
type(criterion_iterator_t) :: iterator
logical :: is_crit_conv
PUSH_SUB(scf_run)
if(scf%forced_finish) then
message(1) = "Previous clean stop, not doing SCF and quitting."
call messages_fatal(1, only_root_writes = .true., namespace=namespace)
end if
gs_run_ = .true.
if(present(gs_run)) gs_run_ = gs_run
verbosity_ = VERB_FULL
if(present(verbosity)) verbosity_ = verbosity
if(scf%lcao_restricted) then
call lcao_init(lcao, namespace, space, gr, ions, st)
if(.not. lcao_is_available(lcao)) then
message(1) = 'LCAO is not available. Cannot do SCF in LCAO.'
call messages_fatal(1, namespace=namespace)
end if
end if
nspin = st%d%nspin
if (present(restart_load)) then
if (restart_has_flag(restart_load, RESTART_FLAG_RHO)) then
! Load density and used it to recalculated the KS potential.
call states_elec_load_rho(restart_load, space, st, gr%mesh, ierr)
if (ierr /= 0) then
message(1) = 'Unable to read density. Density will be calculated from states.'
call messages_warning(1, namespace=namespace)
else
if (bitand(ks%xc_family, XC_FAMILY_OEP) == 0) then
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
else
if (.not. restart_has_flag(restart_load, RESTART_FLAG_VHXC) .and. ks%oep%level /= XC_OEP_FULL) then
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
end if
end if
end if
end if
if (restart_has_flag(restart_load, RESTART_FLAG_VHXC)) then
call hamiltonian_elec_load_vhxc(restart_load, hm, space, gr%mesh, ierr)
if (ierr /= 0) then
message(1) = 'Unable to read Vhxc. Vhxc will be calculated from states.'
call messages_warning(1, namespace=namespace)
else
call hamiltonian_elec_update(hm, gr%mesh, namespace, space, ext_partners)
if (bitand(ks%xc_family, XC_FAMILY_OEP) /= 0) then
if (ks%oep%level == XC_OEP_FULL) then
do is = 1, st%d%nspin
ks%oep%vxc(1:gr%mesh%np, is) = hm%vhxc(1:gr%mesh%np, is) - hm%vhartree(1:gr%mesh%np)
end do
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners)
end if
end if
end if
end if
if (restart_has_flag(restart_load, RESTART_FLAG_MIX)) then
if (scf%mix_field == OPTION__MIXFIELD__DENSITY .or. scf%mix_field == OPTION__MIXFIELD__POTENTIAL) then
call mix_load(namespace, restart_load, scf%smix, space, gr%mesh, ierr)
end if
if (ierr /= 0) then
message(1) = "Unable to read mixing information. Mixing will start from scratch."
call messages_warning(1, namespace=namespace)
end if
end if
if(hm%lda_u_level /= DFT_U_NONE) then
call lda_u_load(restart_load, hm%lda_u, st, hm%energy%dft_u, ierr)
if (ierr /= 0) then
message(1) = "Unable to read LDA+U information. LDA+U data will be calculated from states."
call messages_warning(1, namespace=namespace)
end if
end if
end if
SAFE_ALLOCATE(rhoout(1:gr%mesh%np, 1, 1:nspin))
SAFE_ALLOCATE(rhoin (1:gr%mesh%np, 1, 1:nspin))
rhoin(1:gr%mesh%np, 1, 1:nspin) = st%rho(1:gr%mesh%np, 1:nspin)
rhoout = M_ZERO
!We store the Hxc potential for the contribution to the forces
if (scf%calc_force .or. (outp%duringscf .and. outp%what(OPTION__OUTPUT__FORCES))) then
SAFE_ALLOCATE(vhxc_old(1:gr%mesh%np, 1:nspin))
vhxc_old(1:gr%mesh%np, 1:nspin) = hm%vhxc(1:gr%mesh%np, 1:nspin)
end if
select case(scf%mix_field)
case(OPTION__MIXFIELD__POTENTIAL)
call mixfield_set_vin(scf%mixfield, hm%vhxc)
case(OPTION__MIXFIELD__DENSITY)
call mixfield_set_vin(scf%mixfield, rhoin)
case(OPTION__MIXFIELD__STATES)
SAFE_ALLOCATE_TYPE_ARRAY(wfs_elec_t, psioutb, (st%group%block_start:st%group%block_end, st%d%kpt%start:st%d%kpt%end))
do iqn = st%d%kpt%start, st%d%kpt%end
do ib = st%group%block_start, st%group%block_end
call st%group%psib(ib, iqn)%copy_to(psioutb(ib, iqn))
end do
end do
end select
call lda_u_update_occ_matrices(hm%lda_u, namespace, gr%mesh, st, hm%hm_base, hm%energy)
!If we use LDA+U, we also have do mix it
if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
call create_convergence_file(STATIC_DIR, "convergence")
if ( verbosity_ /= VERB_NO ) then
if(scf%max_iter > 0) then
write(message(1),'(a)') 'Info: Starting SCF iteration.'
else
write(message(1),'(a)') 'Info: No SCF iterations will be done.'
! we cannot tell whether it is converged.
finish = .false.
end if
call messages_info(1, namespace=namespace)
end if
converged_current = .false.
scf%matvec = 0
! SCF cycle
itime = loct_clock()
do iter = 1, scf%max_iter
call profiling_in(prof, "SCF_CYCLE")
! this initialization seems redundant but avoids improper optimization at -O3 by PGI 7 on chum,
! which would cause a failure of testsuite/linear_response/04-vib_modes.03-vib_modes_fd.inp
scf%eigens%converged = 0
!We update the quantities at the begining of the scf cycle
if (iter == 1) then
scf%evsum_in = states_elec_eigenvalues_sum(st)
end if
call iterator%start(scf%criterion_list)
do while (iterator%has_next())
crit => iterator%get_next()
call scf_update_initial_quantity(scf, hm, crit)
end do
!Used for computing the imperfect convegence contribution to the forces
if (scf%calc_force .or. (outp%duringscf .and. outp%what(OPTION__OUTPUT__FORCES))) then
vhxc_old(1:gr%mesh%np, 1:nspin) = hm%vhxc(1:gr%mesh%np, 1:nspin)
end if
if(scf%lcao_restricted) then
call lcao_init_orbitals(lcao, namespace, st, gr, ions)
call lcao_wf(lcao, st, gr, ions, hm, namespace)
else
!We check if the system is coupled with a partner that requires self-consistency
! if(hamiltonian_has_scf_partner(hm)) then
if (allocated(hm%vberry)) then
!In this case, v_Hxc is frozen and we do an internal SCF loop over the
! partners that require SCF
ks%frozen_hxc = .true.
! call perform_scf_partners()
call berry_perform_internal_scf(scf%berry, namespace, space, scf%eigens, gr, st, hm, iter, ks, ions, ext_partners)
!and we unfreeze the potential once finished
ks%frozen_hxc = .false.
else
scf%eigens%converged = 0
call eigensolver_run(scf%eigens, namespace, gr, st, hm, iter)
end if
end if
scf%matvec = scf%matvec + scf%eigens%matvec
! occupations
call states_elec_fermi(st, namespace, gr%mesh)
call lda_u_update_occ_matrices(hm%lda_u, namespace, gr%mesh, st, hm%hm_base, hm%energy)
! compute output density, potential (if needed) and eigenvalues sum
call density_calc(st, gr, st%rho)
rhoout(1:gr%mesh%np, 1, 1:nspin) = st%rho(1:gr%mesh%np, 1:nspin)
select case (scf%mix_field)
case (OPTION__MIXFIELD__POTENTIAL)
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=outp%duringscf)
call mixfield_set_vout(scf%mixfield, hm%vhxc)
case (OPTION__MIXFIELD__DENSITY)
call mixfield_set_vout(scf%mixfield, rhoout)
case(OPTION__MIXFIELD__STATES)
do iqn = st%d%kpt%start, st%d%kpt%end
do ib = st%group%block_start, st%group%block_end
call st%group%psib(ib, iqn)%copy_data_to(gr%mesh%np, psioutb(ib, iqn))
end do
end do
end select
if (scf%mix_field /= OPTION__MIXFIELD__STATES .and. scf%mix_field /= OPTION__MIXFIELD__NONE) then
call lda_u_mixer_set_vout(hm%lda_u, scf%lda_u_mix)
endif
! recalculate total energy
call energy_calc_total(namespace, space, hm, gr, st, ext_partners, iunit = 0)
! compute forces only if requested
if (outp%duringscf .and. outp%what_now(OPTION__OUTPUT__FORCES, iter) &
.and. gs_run_) then
call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=vhxc_old)
end if
!We update the quantities at the end of the scf cycle
call iterator%start(scf%criterion_list)
do while (iterator%has_next())
crit => iterator%get_next()
call scf_update_diff_quantity(scf, hm, st, gr, rhoout, rhoin, crit)
end do
! are we finished?
converged_last = converged_current
converged_current = scf%check_conv .and. &
(.not. scf%conv_eigen_error .or. all(scf%eigens%converged == st%nst))
!Loop over the different criteria
call iterator%start(scf%criterion_list)
do while (iterator%has_next())
crit => iterator%get_next()
call crit%is_converged(is_crit_conv)
converged_current = converged_current .and. is_crit_conv
end do
! only finish if the convergence criteria are fulfilled in two
! consecutive iterations
finish = converged_last .and. converged_current
etime = loct_clock() - itime
itime = etime + itime
call scf_write_iter(namespace)
! mixing
select case (scf%mix_field)
case (OPTION__MIXFIELD__DENSITY)
! mix input and output densities and compute new potential
call mixing(namespace, scf%smix)
call mixfield_get_vnew(scf%mixfield, st%rho)
! for spinors, having components 3 or 4 be negative is not unphysical
if (minval(st%rho(1:gr%mesh%np, 1:st%d%spin_channels)) < -CNST(1e-6)) then
write(message(1),*) 'Negative density after mixing. Minimum value = ', &
minval(st%rho(1:gr%mesh%np, 1:st%d%spin_channels))
call messages_warning(1, namespace=namespace)
end if
call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=outp%duringscf)
case (OPTION__MIXFIELD__POTENTIAL)
! mix input and output potentials
call mixing(namespace, scf%smix)
call mixfield_get_vnew(scf%mixfield, hm%vhxc)
call lda_u_mixer_get_vnew(hm%lda_u, scf%lda_u_mix, st)
call hamiltonian_elec_update_pot(hm, gr%mesh, accel_copy=.true.)
case(OPTION__MIXFIELD__STATES)
do iqn = st%d%kpt%start, st%d%kpt%end
do ib = st%group%block_start, st%group%block_end
call batch_scal(gr%mesh%np, M_ONE - mix_coefficient(scf%smix), st%group%psib(ib, iqn))
call batch_axpy(gr%mesh%np, mix_coefficient(scf%smix), psioutb(ib, iqn), st%group%psib(ib, iqn))
end do
end do
call density_calc(st, gr, st%rho)
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=outp%duringscf)
case (OPTION__MIXFIELD__NONE)
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, calc_current=outp%duringscf)
end select
! Are we asked to stop? (Whenever Fortran is ready for signals, this should go away)
scf%forced_finish = clean_stop(mc%master_comm) .or. walltimer_alarm(mc%master_comm)
if (finish .and. st%modelmbparticles%nparticle > 0) then
call modelmb_sym_all_states(space, gr%mesh, st)
end if
if (gs_run_ .and. present(restart_dump)) then
! save restart information
if ( (finish .or. (modulo(iter, outp%restart_write_interval) == 0) &
.or. iter == scf%max_iter .or. scf%forced_finish) ) then
call states_elec_dump(restart_dump, space, st, gr%mesh, hm%kpoints, ierr, iter=iter)
if (ierr /= 0) then
message(1) = 'Unable to write states wavefunctions.'
call messages_warning(1, namespace=namespace)
end if
call states_elec_dump_rho(restart_dump, space, st, gr%mesh, ierr, iter=iter)
if (ierr /= 0) then
message(1) = 'Unable to write density.'
call messages_warning(1, namespace=namespace)
end if
if(hm%lda_u_level /= DFT_U_NONE) then
call lda_u_dump(restart_dump, hm%lda_u, st, ierr)
if (ierr /= 0) then
message(1) = 'Unable to write LDA+U information.'
call messages_warning(1, namespace=namespace)
end if
end if
select case (scf%mix_field)
case (OPTION__MIXFIELD__DENSITY)
call mix_dump(namespace, restart_dump, scf%smix, space, gr%mesh, ierr)
if (ierr /= 0) then
message(1) = 'Unable to write mixing information.'
call messages_warning(1, namespace=namespace)
end if
case (OPTION__MIXFIELD__POTENTIAL)
call hamiltonian_elec_dump_vhxc(restart_dump, hm, space, gr%mesh, ierr)
if (ierr /= 0) then
message(1) = 'Unable to write Vhxc.'
call messages_warning(1, namespace=namespace)
end if
call mix_dump(namespace, restart_dump, scf%smix, space, gr%mesh, ierr)
if (ierr /= 0) then
message(1) = 'Unable to write mixing information.'
call messages_warning(1, namespace=namespace)
end if
end select
end if
end if
call write_convergence_file(STATIC_DIR, "convergence")
if(finish) then
if(present(iters_done)) iters_done = iter
if(verbosity_ >= VERB_COMPACT) then
write(message(1), '(a, i4, a)') 'Info: SCF converged in ', iter, ' iterations'
write(message(2), '(a)') ''
call messages_info(2, namespace=namespace)
end if
call profiling_out(prof)
exit
end if
if (any(outp%what) .and. outp%duringscf .and. gs_run_) then
do what_i = lbound(outp%what, 1), ubound(outp%what, 1)
if (outp%what_now(what_i, iter)) then
write(dirname,'(a,a,i4.4)') trim(outp%iter_dir),"scf.",iter
call output_all(outp, namespace, space, dirname, gr, ions, iter, st, hm, ks)
exit
end if
end do
end if
! save information for the next iteration
rhoin(1:gr%mesh%np, 1, 1:nspin) = st%rho(1:gr%mesh%np, 1:nspin)
! restart mixing
if (scf%mix_field /= OPTION__MIXFIELD__NONE) then
if (scf%smix%ns_restart > 0) then
if (mod(iter, scf%smix%ns_restart) == 0) then
message(1) = "Info: restarting mixing."
call messages_info(1, namespace=namespace)
call scf_mix_clear(scf)
end if
end if
end if
select case(scf%mix_field)
case(OPTION__MIXFIELD__POTENTIAL)
call mixfield_set_vin(scf%mixfield, hm%vhxc(1:gr%mesh%np, 1:nspin))
case (OPTION__MIXFIELD__DENSITY)
call mixfield_set_vin(scf%mixfield, rhoin)
end select
if(scf%mix_field /= OPTION__MIXFIELD__STATES) call lda_u_mixer_set_vin(hm%lda_u, scf%lda_u_mix)
if(scf%forced_finish) then
call profiling_out(prof)
exit
end if
! check if debug mode should be enabled or disabled on the fly
call io_debug_on_the_fly(namespace)
call profiling_out(prof)
end do !iter
if(scf%lcao_restricted) call lcao_end(lcao)
if ((scf%max_iter > 0 .and. scf%mix_field == OPTION__MIXFIELD__POTENTIAL) .or. output_needs_current(outp, states_are_real(st))) then
call v_ks_calc(ks, namespace, space, hm, st, ions, ext_partners, &
calc_current=output_needs_current(outp, states_are_real(st)))
end if
select case(scf%mix_field)
case(OPTION__MIXFIELD__STATES)
do iqn = st%d%kpt%start, st%d%kpt%end
do ib = st%group%block_start, st%group%block_end
call psioutb(ib, iqn)%end()
end do
end do
SAFE_DEALLOCATE_A(psioutb)
end select
SAFE_DEALLOCATE_A(rhoout)
SAFE_DEALLOCATE_A(rhoin)
if(scf%max_iter > 0 .and. any(scf%eigens%converged < st%nst) .and. .not. scf%lcao_restricted) then
write(message(1),'(a)') 'Some of the states are not fully converged!'
call messages_warning(1, namespace=namespace)
end if
if(.not.finish) then
write(message(1), '(a,i4,a)') 'SCF *not* converged after ', iter - 1, ' iterations.'
call messages_warning(1, namespace=namespace)
end if
write(message(1), '(a,i10)') 'Info: Number of matrix-vector products: ', scf%matvec
call messages_info(1)
! calculate forces
if (scf%calc_force) then
call forces_calculate(gr, namespace, ions, hm, ext_partners, st, ks, vhxc_old=vhxc_old)
end if
! calculate stress
if (scf%calc_stress) call stress_calculate(namespace, gr, hm, st, ions, ks)
if(scf%max_iter == 0) then
call energy_calc_eigenvalues(namespace, hm, gr%der, st)
call states_elec_fermi(st, namespace, gr%mesh)
call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, namespace=namespace)
end if
if(gs_run_) then
! output final information
call scf_write_static(STATIC_DIR, "info")
call output_all(outp, namespace, space, STATIC_DIR, gr, ions, -1, st, hm, ks)
end if
if (space%is_periodic() .and. st%d%nik > st%d%nspin) then
if (bitand(hm%kpoints%method, KPOINTS_PATH) /= 0) then
call states_elec_write_bandstructure(STATIC_DIR, namespace, st%nst, st, gr%box, &
ions, gr%mesh, hm%kpoints, hm%hm_base%phase, vec_pot = hm%hm_base%uniform_vector_potential, &
vec_pot_var = hm%hm_base%vector_potential)
end if
end if
if( ks%vdw_correction == OPTION__VDWCORRECTION__VDW_TS) then
call vdw_ts_write_c6ab(ks%vdw_ts, ions, STATIC_DIR, 'c6ab_eff', namespace)
end if
SAFE_DEALLOCATE_A(vhxc_old)
POP_SUB(scf_run)
contains
! ---------------------------------------------------------
subroutine scf_write_iter(namespace)
type(namespace_t), intent(in) :: namespace
character(len=50) :: str
FLOAT :: dipole(1:space%dim)
PUSH_SUB(scf_run.scf_write_iter)
if ( verbosity_ == VERB_FULL ) then
write(str, '(a,i5)') 'SCF CYCLE ITER #' ,iter
call messages_print_stress(msg=trim(str), namespace=namespace)
write(message(1),'(a,es15.8,2(a,es9.2))') ' etot = ', units_from_atomic(units_out%energy, hm%energy%total), &
' abs_ev = ', units_from_atomic(units_out%energy, scf%evsum_diff), &
' rel_ev = ', scf%evsum_diff/abs(scf%evsum_out)
write(message(2),'(a,es15.2,2(a,es9.2))') &
' ediff = ', scf%energy_diff, ' abs_dens = ', scf%abs_dens_diff, &
' rel_dens = ', scf%abs_dens_diff/st%qtot
call messages_info(2, namespace=namespace)
if(.not.scf%lcao_restricted) then
write(message(1),'(a,i6)') 'Matrix vector products: ', scf%eigens%matvec
write(message(2),'(a,i6)') 'Converged eigenvectors: ', sum(scf%eigens%converged(1:st%d%nik))
call messages_info(2, namespace=namespace)
call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, scf%eigens%diff, compact = .true., namespace=namespace)
else
call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, compact = .true., namespace=namespace)
end if
if (allocated(hm%vberry)) then
call calc_dipole(dipole, space, gr%mesh, st, ions)
call write_dipole(dipole, namespace=namespace)
end if
if(st%d%ispin > UNPOLARIZED) then
call write_magnetic_moments(gr%mesh, st, ions, gr%der%boundaries, scf%lmm_r, namespace=namespace)
end if
if(hm%lda_u_level == DFT_U_ACBN0) then
call lda_u_write_U(hm%lda_u, namespace=namespace)
call lda_u_write_V(hm%lda_u, namespace=namespace)
end if
write(message(1),'(a)') ''
write(message(2),'(a,i5,a,f14.2)') 'Elapsed time for SCF step ', iter,':', etime
call messages_info(2, namespace=namespace)
call scf_print_mem_use(namespace)
call messages_print_stress(namespace=namespace)
end if
if ( verbosity_ == VERB_COMPACT ) then
write(message(1),'(a,i4,a,es15.8, a,es9.2, a, f7.1, a)') &
'iter ', iter, &
' : etot ', units_from_atomic(units_out%energy, hm%energy%total), &
' : abs_dens', scf%abs_dens_diff, &
' : etime ', etime, 's'
call messages_info(1, namespace=namespace)
end if
POP_SUB(scf_run.scf_write_iter)
end subroutine scf_write_iter