28 USE dbcsr_api,
ONLY: dbcsr_p_type,&
69 calculate_subspace_eigenvalues
84 #include "./base/base_uses.f90"
90 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_scf_output'
110 INTEGER,
INTENT(IN) :: output_unit
111 TYPE(qs_environment_type),
POINTER :: qs_env
113 INTEGER :: nelectron_total
114 LOGICAL :: gapw, gapw_xc, qmmm
115 TYPE(dft_control_type),
POINTER :: dft_control
116 TYPE(qs_charges_type),
POINTER :: qs_charges
117 TYPE(qs_energy_type),
POINTER :: energy
118 TYPE(qs_rho_type),
POINTER :: rho
119 TYPE(qs_scf_env_type),
POINTER :: scf_env
121 NULLIFY (rho, energy, dft_control, scf_env, qs_charges)
122 CALL get_qs_env(qs_env=qs_env, rho=rho, energy=energy, dft_control=dft_control, &
123 scf_env=scf_env, qs_charges=qs_charges)
125 gapw = dft_control%qs_control%gapw
126 gapw_xc = dft_control%qs_control%gapw_xc
128 nelectron_total = scf_env%nelectron
130 CALL qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
131 dft_control, qmmm, qs_env, gapw, gapw_xc)
142 INTEGER :: output_unit
143 TYPE(mo_set_type),
DIMENSION(:),
INTENT(IN) :: mos
144 TYPE(dft_control_type),
POINTER :: dft_control
146 INTEGER :: homo, ispin, nao, nelectron_spin, nmo
148 IF (output_unit > 0)
THEN
149 DO ispin = 1, dft_control%nspins
152 nelectron=nelectron_spin, &
155 IF (dft_control%nspins > 1)
THEN
156 WRITE (unit=output_unit, fmt=
"(/,T2,A,I2)")
"Spin", ispin
158 WRITE (unit=output_unit, fmt=
"(/,(T2,A,T71,I10))") &
159 "Number of electrons:", nelectron_spin, &
160 "Number of occupied orbitals:", homo, &
161 "Number of molecular orbitals:", nmo
163 WRITE (unit=output_unit, fmt=
"(/,T2,A,T71,I10)") &
164 "Number of orbital functions:", nao
178 TYPE(qs_environment_type),
POINTER :: qs_env
179 TYPE(qs_scf_env_type),
POINTER :: scf_env
180 LOGICAL,
INTENT(IN) :: final_mos
182 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_write_mos'
184 CHARACTER(LEN=2) :: solver_method
185 CHARACTER(LEN=3*default_string_length) :: message
186 CHARACTER(LEN=5) :: spin
187 CHARACTER(LEN=default_string_length), &
188 DIMENSION(:),
POINTER :: tmpstringlist
189 INTEGER :: handle, homo, ikp, ispin, iw, kpoint, &
190 nao, nelectron, nkp, nmo, nspin, numo
191 INTEGER,
DIMENSION(2) :: nmos_occ
192 INTEGER,
DIMENSION(:),
POINTER :: mo_index_range
193 LOGICAL :: do_kpoints, print_eigvals, &
194 print_eigvecs, print_mo_info, &
195 print_occup, print_occup_stats
196 REAL(kind=
dp) :: flexible_electron_count, maxocc, n_el_f, &
197 occup_stats_occ_threshold
198 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_eigenvalues, umo_eigenvalues
199 TYPE(admm_type),
POINTER :: admm_env
200 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
201 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
202 TYPE(cp_fm_struct_type),
POINTER :: fm_struct_tmp
203 TYPE(cp_fm_type),
POINTER :: mo_coeff, umo_coeff
204 TYPE(cp_logger_type),
POINTER :: logger
205 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: ks, s
206 TYPE(dbcsr_type),
POINTER :: matrix_ks, matrix_s
207 TYPE(dft_control_type),
POINTER :: dft_control
208 TYPE(kpoint_type),
POINTER :: kpoints
209 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
210 TYPE(mo_set_type),
POINTER :: mo_set, umo_set
211 TYPE(mp_para_env_type),
POINTER :: para_env
212 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
213 TYPE(preconditioner_type),
POINTER :: local_preconditioner
214 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
215 TYPE(scf_control_type),
POINTER :: scf_control
216 TYPE(section_vals_type),
POINTER :: dft_section, input
218 CALL timeset(routinen, handle)
220 cpassert(
ASSOCIATED(qs_env))
224 atomic_kind_set=atomic_kind_set, &
225 blacs_env=blacs_env, &
226 dft_control=dft_control, &
227 do_kpoints=do_kpoints, &
229 qs_kind_set=qs_kind_set, &
231 particle_set=particle_set, &
232 scf_control=scf_control)
239 CALL section_vals_val_get(dft_section,
"PRINT%MO%OCCUPATION_NUMBERS_STATS", c_vals=tmpstringlist)
241 print_occup_stats = .false.
242 occup_stats_occ_threshold = 1e-6_dp
243 IF (
SIZE(tmpstringlist) > 0)
THEN
244 print_occup_stats = .true.
245 IF (len_trim(tmpstringlist(1)) > 0) &
246 READ (tmpstringlist(1), *) print_occup_stats
248 IF (
SIZE(tmpstringlist) > 1) &
249 READ (tmpstringlist(2), *) occup_stats_occ_threshold
254 IF ((.NOT. print_mo_info) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup .OR. print_occup_stats)))
THEN
255 CALL timestop(handle)
259 NULLIFY (fm_struct_tmp)
261 NULLIFY (mo_eigenvalues)
264 NULLIFY (umo_eigenvalues)
267 nspin = dft_control%nspins
273 nkp =
SIZE(kpoints%kp_env)
275 CALL get_qs_env(qs_env, matrix_ks=ks, matrix_s=s)
276 cpassert(
ASSOCIATED(ks))
277 cpassert(
ASSOCIATED(s))
284 mos => kpoints%kp_env(ikp)%kpoint_env%mos(1, :)
287 CALL get_qs_env(qs_env, matrix_ks=ks, mos=mos)
290 cpassert(
ASSOCIATED(mos))
301 cpabort(
"The OT method is not implemented for k points")
304 matrix_ks => ks(ispin)%matrix
305 matrix_s => s(1)%matrix
308 IF (dft_control%do_admm)
THEN
316 eigenvalues=mo_eigenvalues, &
319 nelectron=nelectron, &
323 flexible_electron_count=flexible_electron_count)
327 cpassert(
ASSOCIATED(mo_index_range))
328 numo = min(mo_index_range(2) - homo, nao - homo)
330 IF (.NOT. final_mos)
THEN
332 message =
"The MO information for unoccupied MOs is only calculated after "// &
333 "SCF convergence is achieved when the orbital transformation (OT) "// &
335 cpwarn(trim(message))
354 flexible_electron_count=flexible_electron_count)
356 fm_struct=fm_struct_tmp, &
357 name=
"Temporary MO set (unoccupied MOs only)for printout")
360 mo_coeff=umo_coeff, &
361 eigenvalues=umo_eigenvalues)
368 NULLIFY (local_preconditioner)
369 IF (
ASSOCIATED(scf_env%ot_preconditioner))
THEN
370 local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
372 NULLIFY (local_preconditioner)
377 matrix_c_fm=umo_coeff, &
378 matrix_orthogonal_space_fm=mo_coeff, &
379 eps_gradient=scf_control%eps_lumos, &
381 iter_max=scf_control%max_iter_lumos, &
382 size_ortho_space=nmo)
385 CALL calculate_subspace_eigenvalues(orbitals=umo_coeff, &
386 ks_matrix=matrix_ks, &
387 evals_arg=umo_eigenvalues, &
389 CALL set_mo_occupation(mo_set=umo_set)
392 IF (dft_control%do_admm)
THEN
418 cpabort(
"Invalid spin")
420 IF (
ASSOCIATED(umo_set))
THEN
422 dft_section, 4, kpoint, final_mos=final_mos, spin=trim(spin), &
423 solver_method=solver_method, umo_set=umo_set)
426 dft_section, 4, kpoint, final_mos=final_mos, spin=trim(spin), &
427 solver_method=solver_method)
430 IF (
ASSOCIATED(umo_set))
THEN
432 dft_section, 4, kpoint, final_mos=final_mos, &
433 solver_method=solver_method, umo_set=umo_set)
436 dft_section, 4, kpoint, final_mos=final_mos, &
437 solver_method=solver_method)
441 nmos_occ(ispin) = max(nmos_occ(ispin), count(mo_set%occupation_numbers > occup_stats_occ_threshold))
445 IF (
ASSOCIATED(umo_set))
THEN
458 IF (print_mo_info .AND. print_occup_stats)
THEN
460 ignore_should_output=print_mo_info, &
463 IF (
SIZE(mos) > 1)
THEN
464 WRITE (unit=iw, fmt=
"(A,I4)")
" MO| Total occupied (ALPHA):", nmos_occ(1)
465 WRITE (unit=iw, fmt=
"(A,I4)")
" MO| Total occupied (BETA): ", nmos_occ(2)
467 WRITE (unit=iw, fmt=
"(A,I4)")
" MO| Total occupied: ", nmos_occ(1)
469 WRITE (unit=iw, fmt=
"(A)")
""
472 ignore_should_output=print_mo_info)
475 CALL timestop(handle)
490 energy, total_steps, should_stop, outer_loop_converged)
491 INTEGER :: output_unit
492 TYPE(scf_control_type),
POINTER :: scf_control
493 TYPE(qs_scf_env_type),
POINTER :: scf_env
494 TYPE(qs_energy_type),
POINTER :: energy
495 INTEGER :: total_steps
496 LOGICAL,
INTENT(IN) :: should_stop, outer_loop_converged
498 REAL(kind=
dp) :: outer_loop_eps
500 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
501 IF (output_unit > 0)
WRITE (output_unit,
'(/,T3,A,I4,A,E10.2,A,F22.10)') &
502 "outer SCF iter = ", scf_env%outer_scf%iter_count, &
503 " RMS gradient = ", outer_loop_eps,
" energy =", energy%total
505 IF (outer_loop_converged)
THEN
506 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
507 "outer SCF loop converged in", scf_env%outer_scf%iter_count, &
508 " iterations or ", total_steps,
" steps"
509 ELSE IF (scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf &
510 .OR. should_stop)
THEN
511 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
512 "outer SCF loop FAILED to converge after ", &
513 scf_env%outer_scf%iter_count,
" iterations or ", total_steps,
" steps"
529 TYPE(qs_scf_env_type),
POINTER :: scf_env
530 INTEGER :: output_unit
531 LOGICAL :: just_energy
532 REAL(kind=
dp) :: t1, t2
533 TYPE(qs_energy_type),
POINTER :: energy
535 IF ((output_unit > 0) .AND. scf_env%print_iter_line)
THEN
536 IF (just_energy)
THEN
537 WRITE (unit=output_unit, &
538 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") &
539 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
540 t2 - t1, energy%total
542 IF ((abs(scf_env%iter_delta) < 1.0e-8_dp) .OR. &
543 (abs(scf_env%iter_delta) >= 1.0e5_dp))
THEN
544 WRITE (unit=output_unit, &
545 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") &
546 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
547 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
549 WRITE (unit=output_unit, &
550 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
551 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
552 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
577 SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
578 dft_control, qmmm, qs_env, gapw, gapw_xc)
579 INTEGER,
INTENT(IN) :: output_unit
580 TYPE(qs_rho_type),
POINTER :: rho
581 TYPE(qs_charges_type),
POINTER :: qs_charges
582 TYPE(qs_energy_type),
POINTER :: energy
583 INTEGER,
INTENT(IN) :: nelectron_total
584 TYPE(dft_control_type),
POINTER :: dft_control
585 LOGICAL,
INTENT(IN) :: qmmm
586 TYPE(qs_environment_type),
POINTER :: qs_env
587 LOGICAL,
INTENT(IN) :: gapw, gapw_xc
589 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_print_scf_summary'
591 INTEGER :: bc, handle, ispin, psolver
592 REAL(kind=
dp) :: exc1_energy, exc_energy, &
593 implicit_ps_ehartree, tot1_h, tot1_s
594 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tot_rho_r
595 TYPE(pw_env_type),
POINTER :: pw_env
597 NULLIFY (tot_rho_r, pw_env)
598 CALL timeset(routinen, handle)
601 psolver = pw_env%poisson_env%parameters%solver
603 IF (output_unit > 0)
THEN
605 IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
606 dft_control%qs_control%xtb .OR. &
607 dft_control%qs_control%dftb))
THEN
608 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T41,2F20.10))") &
609 "Electronic density on regular grids: ", &
610 accurate_sum(tot_rho_r), &
611 accurate_sum(tot_rho_r) + nelectron_total, &
612 "Core density on regular grids:", &
613 qs_charges%total_rho_core_rspace, &
614 qs_charges%total_rho_core_rspace - real(nelectron_total + dft_control%charge,
dp)
616 IF (dft_control%correct_surf_dip)
THEN
617 WRITE (unit=output_unit, fmt=
"((T3,A,/,T3,A,T41,F20.10))") &
618 "Total dipole moment perpendicular to ", &
619 "the slab [electrons-Angstroem]: ", &
620 qs_env%surface_dipole_moment
624 tot1_h = qs_charges%total_rho1_hard(1)
625 tot1_s = qs_charges%total_rho1_soft(1)
626 DO ispin = 2, dft_control%nspins
627 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
628 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
630 WRITE (unit=output_unit, fmt=
"((T3,A,T41,2F20.10))") &
631 "Hard and soft densities (Lebedev):", &
633 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
634 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
635 accurate_sum(tot_rho_r) + tot1_h - tot1_s, &
636 "Total charge density (r-space): ", &
637 accurate_sum(tot_rho_r) + tot1_h - tot1_s &
638 + qs_charges%total_rho_core_rspace, &
639 "Total Rho_soft + Rho0_soft (g-space):", &
640 qs_charges%total_rho_gspace
642 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
643 "Total charge density on r-space grids: ", &
644 accurate_sum(tot_rho_r) + &
645 qs_charges%total_rho_core_rspace, &
646 "Total charge density g-space grids: ", &
647 qs_charges%total_rho_gspace
650 IF (dft_control%qs_control%semi_empirical)
THEN
651 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
652 "Core-core repulsion energy [eV]: ", energy%core_overlap*
evolt, &
653 "Core Hamiltonian energy [eV]: ", energy%core*
evolt, &
654 "Two-electron integral energy [eV]: ", energy%hartree*
evolt, &
655 "Electronic energy [eV]: ", &
656 (energy%core + 0.5_dp*energy%hartree)*
evolt
657 IF (energy%dispersion /= 0.0_dp) &
658 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
659 "Dispersion energy [eV]: ", energy%dispersion*
evolt
660 ELSEIF (dft_control%qs_control%dftb)
THEN
661 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
662 "Core Hamiltonian energy: ", energy%core, &
663 "Repulsive potential energy: ", energy%repulsive, &
664 "Electronic energy: ", energy%hartree, &
665 "Dispersion energy: ", energy%dispersion
666 IF (energy%dftb3 /= 0.0_dp) &
667 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
668 "DFTB3 3rd order energy: ", energy%dftb3
669 IF (energy%efield /= 0.0_dp) &
670 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
671 "Electric field interaction energy: ", energy%efield
672 ELSEIF (dft_control%qs_control%xtb)
THEN
673 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
674 "Core Hamiltonian energy: ", energy%core, &
675 "Repulsive potential energy: ", energy%repulsive, &
676 "Electronic energy: ", energy%hartree, &
677 "DFTB3 3rd order energy: ", energy%dftb3, &
678 "Dispersion energy: ", energy%dispersion
679 IF (dft_control%qs_control%xtb_control%xb_interaction) &
680 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
681 "Correction for halogen bonding: ", energy%xtb_xb_inter
682 IF (dft_control%qs_control%xtb_control%do_nonbonded) &
683 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
684 "Correction for nonbonded interactions: ", energy%xtb_nonbonded
685 IF (energy%efield /= 0.0_dp) &
686 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
687 "Electric field interaction energy: ", energy%efield
689 IF (dft_control%do_admm)
THEN
690 exc_energy = energy%exc + energy%exc_aux_fit
691 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
693 exc_energy = energy%exc
694 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
698 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
699 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
702 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
703 "Overlap energy of the core charge distribution:", energy%core_overlap, &
704 "Self energy of the core charge distribution: ", energy%core_self, &
705 "Core Hamiltonian energy: ", energy%core, &
706 "Hartree energy: ", implicit_ps_ehartree, &
707 "Electric enthalpy: ", energy%hartree, &
708 "Exchange-correlation energy: ", exc_energy
710 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
711 "Overlap energy of the core charge distribution:", energy%core_overlap, &
712 "Self energy of the core charge distribution: ", energy%core_self, &
713 "Core Hamiltonian energy: ", energy%core, &
714 "Hartree energy: ", energy%hartree, &
715 "Exchange-correlation energy: ", exc_energy
718 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
719 "Overlap energy of the core charge distribution:", energy%core_overlap, &
720 "Self energy of the core charge distribution: ", energy%core_self, &
721 "Core Hamiltonian energy: ", energy%core, &
722 "Hartree energy: ", energy%hartree, &
723 "Exchange-correlation energy: ", exc_energy
725 IF (energy%e_hartree /= 0.0_dp) &
726 WRITE (unit=output_unit, fmt=
"(T3,A,/,T3,A,T56,F25.14)") &
727 "Coulomb Electron-Electron Interaction Energy ", &
728 "- Already included in the total Hartree term ", energy%e_hartree
729 IF (energy%ex /= 0.0_dp) &
730 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
731 "Hartree-Fock Exchange energy: ", energy%ex
732 IF (energy%dispersion /= 0.0_dp) &
733 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
734 "Dispersion energy: ", energy%dispersion
735 IF (energy%gcp /= 0.0_dp) &
736 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
737 "gCP energy: ", energy%gcp
738 IF (energy%efield /= 0.0_dp) &
739 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
740 "Electric field interaction energy: ", energy%efield
742 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
743 "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
744 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
747 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
748 "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
751 IF (dft_control%smear)
THEN
752 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
753 "Electronic entropic energy:", energy%kTS
754 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
755 "Fermi energy:", energy%efermi
757 IF (dft_control%dft_plus_u)
THEN
758 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
759 "DFT+U energy:", energy%dft_plus_u
761 IF (dft_control%do_sccs)
THEN
762 WRITE (unit=output_unit, fmt=
"(A)")
""
766 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
767 "QM/MM Electrostatic energy: ", energy%qmmm_el
768 IF (qs_env%qmmm_env_qm%image_charge)
THEN
769 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
770 "QM/MM image charge energy: ", energy%image_charge
773 IF (dft_control%qs_control%mulliken_restraint)
THEN
774 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
775 "Mulliken restraint energy: ", energy%mulliken
777 IF (dft_control%qs_control%semi_empirical)
THEN
778 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
779 "Total energy [eV]: ", energy%total*
evolt
780 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
781 "Atomic reference energy [eV]: ", energy%core_self*
evolt, &
782 "Heat of formation [kcal/mol]: ", &
783 (energy%total + energy%core_self)*
kcalmol
785 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
786 "Total energy: ", energy%total
789 IF (qs_env%qmmm_env_qm%image_charge)
THEN
796 CALL timestop(handle)
798 END SUBROUTINE qs_scf_print_scf_summary
809 TYPE(qs_environment_type),
POINTER :: qs_env
810 TYPE(qs_scf_env_type),
POINTER :: scf_env
811 TYPE(mp_para_env_type),
POINTER :: para_env
813 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_loop_print'
815 INTEGER :: after, handle, ic, ispin, iw
816 LOGICAL :: do_kpoints, omit_headers
817 REAL(kind=
dp) :: mo_mag_max, mo_mag_min, orthonormality
818 TYPE(cp_logger_type),
POINTER :: logger
819 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_p, matrix_s
820 TYPE(dft_control_type),
POINTER :: dft_control
821 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
822 TYPE(qs_rho_type),
POINTER :: rho
823 TYPE(section_vals_type),
POINTER :: dft_section, input, scf_section
826 CALL timeset(routinen, handle)
828 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
829 do_kpoints=do_kpoints)
835 DO ispin = 1, dft_control%nspins
838 dft_section,
"PRINT%AO_MATRICES/DENSITY"),
cp_p_file))
THEN
844 after = min(max(after, 1), 16)
845 DO ic = 1,
SIZE(matrix_p, 2)
847 output_unit=iw, omit_headers=omit_headers)
850 "PRINT%AO_MATRICES/DENSITY")
854 dft_section,
"PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"),
cp_p_file))
THEN
858 after = min(max(after, 1), 16)
859 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
860 DO ic = 1,
SIZE(matrix_ks, 2)
861 IF (dft_control%qs_control%semi_empirical)
THEN
863 scale=
evolt, output_unit=iw, omit_headers=omit_headers)
866 output_unit=iw, omit_headers=omit_headers)
870 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
876 scf_section,
"PRINT%MO_ORTHONORMALITY"),
cp_p_file))
THEN
881 WRITE (iw,
'(T8,A)') &
882 " K-points: Maximum deviation from MO S-orthonormality not determined"
885 "PRINT%MO_ORTHONORMALITY")
891 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
897 WRITE (iw,
'(T8,A,T61,E20.4)') &
898 " Maximum deviation from MO S-orthonormality", orthonormality
901 "PRINT%MO_ORTHONORMALITY")
905 scf_section,
"PRINT%MO_MAGNITUDE"),
cp_p_file))
THEN
910 WRITE (iw,
'(T8,A)') &
911 " K-points: Minimum/Maximum MO magnitude not determined"
914 "PRINT%MO_MAGNITUDE")
921 WRITE (iw,
'(T8,A,T41,2E20.4)') &
922 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
925 "PRINT%MO_MAGNITUDE")
929 CALL timestop(handle)
948 energy, total_steps, should_stop, outer_loop_converged, &
950 INTEGER :: output_unit
951 TYPE(scf_control_type),
POINTER :: scf_control
952 TYPE(qs_scf_env_type),
POINTER :: scf_env
953 TYPE(cdft_control_type),
POINTER :: cdft_control
954 TYPE(qs_energy_type),
POINTER :: energy
955 INTEGER :: total_steps
956 LOGICAL,
INTENT(IN) :: should_stop, outer_loop_converged, &
959 REAL(kind=
dp) :: outer_loop_eps
962 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
963 IF (output_unit > 0)
WRITE (output_unit,
'(/,T3,A,I4,A,E10.2,A,F22.10)') &
964 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
965 " RMS gradient = ", outer_loop_eps,
" energy =", energy%total
966 IF (outer_loop_converged)
THEN
967 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
968 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
969 " iterations or ", total_steps,
" steps"
971 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
972 .AND. .NOT. outer_loop_converged)
THEN
973 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
974 "CDFT SCF loop FAILED to converge after ", &
975 scf_env%outer_scf%iter_count,
" iterations or ", total_steps,
" steps"
990 INTEGER :: output_unit
991 TYPE(cdft_control_type),
POINTER :: cdft_control
993 IF (output_unit > 0)
THEN
994 WRITE (output_unit,
'(/,A)') &
995 " ---------------------------------- CDFT --------------------------------------"
996 WRITE (output_unit,
'(A)') &
997 " Optimizing a density constraint in an external SCF loop "
998 WRITE (output_unit,
'(A)')
" "
999 SELECT CASE (cdft_control%type)
1001 WRITE (output_unit,
'(A)')
" Type of constraint: Hirshfeld"
1003 WRITE (output_unit,
'(A)')
" Type of constraint: Becke"
1005 WRITE (output_unit,
'(A,I8)')
" Number of constraints: ",
SIZE(cdft_control%group)
1006 WRITE (output_unit,
'(A,L8)')
" Using fragment densities:", cdft_control%fragment_density
1007 WRITE (output_unit,
'(A)')
" "
1008 IF (cdft_control%atomic_charges)
WRITE (output_unit,
'(A,/)')
" Calculating atomic CDFT charges"
1009 SELECT CASE (cdft_control%constraint_control%optimizer)
1011 WRITE (output_unit,
'(A)') &
1012 " Minimizer : SD : steepest descent"
1014 WRITE (output_unit,
'(A)') &
1015 " Minimizer : DIIS : direct inversion"
1016 WRITE (output_unit,
'(A)') &
1017 " in the iterative subspace"
1018 WRITE (output_unit,
'(A,I3,A)') &
1020 cdft_control%constraint_control%diis_buffer_length,
" DIIS vectors"
1022 WRITE (output_unit,
'(A)') &
1023 " Minimizer : BISECT : gradient bisection"
1024 WRITE (output_unit,
'(A,I3)') &
1025 " using a trust count of", &
1026 cdft_control%constraint_control%bisect_trust_count
1030 cdft_control%constraint_control%optimizer, output_unit)
1032 WRITE (output_unit,
'(A)')
" Minimizer : Secant"
1036 WRITE (output_unit,
'(/,A,L7)') &
1037 " Reusing OT preconditioner: ", cdft_control%reuse_precond
1038 IF (cdft_control%reuse_precond)
THEN
1039 WRITE (output_unit,
'(A,I3,A,I3,A)') &
1040 " using old preconditioner for up to ", &
1041 cdft_control%max_reuse,
" subsequent CDFT SCF"
1042 WRITE (output_unit,
'(A,I3,A,I3,A)') &
1043 " iterations if the relevant loop converged in less than ", &
1044 cdft_control%precond_freq,
" steps"
1046 SELECT CASE (cdft_control%type)
1048 WRITE (output_unit,
'(/,A)')
" Hirshfeld constraint settings"
1049 WRITE (output_unit,
'(A)')
" "
1050 SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1052 WRITE (output_unit,
'(A, A8)') &
1053 " Shape function type: ",
"Gaussian"
1054 WRITE (output_unit,
'(A)', advance=
'NO') &
1055 " Type of Gaussian: "
1056 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1058 WRITE (output_unit,
'(A13)')
"Default"
1060 WRITE (output_unit,
'(A13)')
"Covalent"
1062 WRITE (output_unit,
'(A13)')
"Fixed radius"
1064 WRITE (output_unit,
'(A13)')
"Van der Waals"
1066 WRITE (output_unit,
'(A13)')
"User-defined"
1070 WRITE (output_unit,
'(A, A8)') &
1071 " Shape function type: ",
"Density"
1074 WRITE (output_unit,
'(/, A)')
" Becke constraint settings"
1075 WRITE (output_unit,
'(A)')
" "
1076 SELECT CASE (cdft_control%becke_control%cutoff_type)
1078 WRITE (output_unit,
'(A,F8.3,A)') &
1079 " Cutoff for partitioning :",
cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1080 "angstrom"),
" angstrom"
1082 WRITE (output_unit,
'(A)') &
1083 " Using element specific cutoffs for partitioning"
1085 WRITE (output_unit,
'(A,L7)') &
1086 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1087 WRITE (output_unit,
'(A,L7)') &
1088 " Precompute gradients : ", cdft_control%becke_control%in_memory
1089 WRITE (output_unit,
'(A)')
" "
1090 IF (cdft_control%becke_control%adjust) &
1091 WRITE (output_unit,
'(A)') &
1092 " Using atomic radii to generate a heteronuclear charge partitioning"
1093 WRITE (output_unit,
'(A)')
" "
1094 IF (.NOT. cdft_control%becke_control%cavity_confine)
THEN
1095 WRITE (output_unit,
'(A)') &
1096 " No confinement is active"
1098 WRITE (output_unit,
'(A)')
" Confinement using a Gaussian shaped cavity is active"
1099 SELECT CASE (cdft_control%becke_control%cavity_shape)
1101 WRITE (output_unit,
'(A,F8.4, A)') &
1102 " Type of Gaussian : Fixed radius: ", &
1105 WRITE (output_unit,
'(A)') &
1106 " Type of Gaussian : Covalent radius "
1108 WRITE (output_unit,
'(A)') &
1109 " Type of Gaussian : vdW radius "
1111 WRITE (output_unit,
'(A)') &
1112 " Type of Gaussian : User radius "
1114 WRITE (output_unit,
'(A,ES12.4)') &
1115 " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1118 WRITE (output_unit,
'(/,A)') &
1119 " ---------------------------------- CDFT --------------------------------------"
1132 INTEGER :: output_unit
1133 TYPE(cdft_control_type),
POINTER :: cdft_control
1137 IF (output_unit > 0)
THEN
1138 SELECT CASE (cdft_control%type)
1140 WRITE (output_unit,
'(/,T3,A,T60)') &
1141 '------------------- Hirshfeld constraint information -------------------'
1143 WRITE (output_unit,
'(/,T3,A,T60)') &
1144 '--------------------- Becke constraint information ---------------------'
1146 cpabort(
"Unknown CDFT constraint.")
1148 DO igroup = 1,
SIZE(cdft_control%target)
1149 IF (igroup > 1)
WRITE (output_unit,
'(T3,A)')
' '
1150 WRITE (output_unit,
'(T3,A,T54,(3X,I18))') &
1151 'Atomic group :', igroup
1152 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1154 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1155 WRITE (output_unit,
'(T3,A,T42,A)') &
1156 'Type of constraint :', adjustr(
'Charge density constraint (frag.)')
1158 WRITE (output_unit,
'(T3,A,T50,A)') &
1159 'Type of constraint :', adjustr(
'Charge density constraint')
1162 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1163 WRITE (output_unit,
'(T3,A,T35,A)') &
1164 'Type of constraint :', adjustr(
'Magnetization density constraint (frag.)')
1166 WRITE (output_unit,
'(T3,A,T43,A)') &
1167 'Type of constraint :', adjustr(
'Magnetization density constraint')
1170 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1171 WRITE (output_unit,
'(T3,A,T38,A)') &
1172 'Type of constraint :', adjustr(
'Alpha spin density constraint (frag.)')
1174 WRITE (output_unit,
'(T3,A,T46,A)') &
1175 'Type of constraint :', adjustr(
'Alpha spin density constraint')
1178 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1179 WRITE (output_unit,
'(T3,A,T39,A)') &
1180 'Type of constraint :', adjustr(
'Beta spin density constraint (frag.)')
1182 WRITE (output_unit,
'(T3,A,T47,A)') &
1183 'Type of constraint :', adjustr(
'Beta spin density constraint')
1186 cpabort(
"Unknown constraint type.")
1188 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1189 'Target value of constraint :', cdft_control%target(igroup)
1190 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1191 'Current value of constraint :', cdft_control%value(igroup)
1192 WRITE (output_unit,
'(T3,A,T59,(3X,ES13.3))') &
1193 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1194 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1195 'Strength of constraint :', cdft_control%strength(igroup)
1197 WRITE (output_unit,
'(T3,A)') &
1198 '------------------------------------------------------------------------'
Types and set/get functions for auxiliary density matrix methods.
Contains methods used in the context of density fitting.
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Define the atomic kind types and their sub types.
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Types and basic routines needed for a kpoint calculation.
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Interface to the message passing library MPI.
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public kcalmol
real(kind=dp), parameter, public evolt
computes preconditioners, and implements methods to apply them currently used in qs_ot
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
integer, parameter, public mixed_bc
integer, parameter, public mixed_periodic_bc
integer, parameter, public periodic_bc
container for various plainwaves related things
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_implicit
Routines for image charge calculation within QM/MM.
subroutine, public print_image_coefficients(image_coeff, qs_env)
Print image coefficients.
Control parameters for optimizers that work with CDFT constraints.
subroutine, public cdft_opt_type_write(cdft_opt_control, optimizer, output_unit)
writes information about the CDFT optimizer object
Defines CDFT control structures.
container for information about total charges on the grids
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition and initialisation of the mo data type.
subroutine, public write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, dft_section, before, kpoint, final_mos, spin, solver_method, rtp, cpart, sim_step, umo_set)
Write MO information to output file (eigenvalues, occupation numbers, coefficients)
collects routines that perform operations directly related to MOs
subroutine, public calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
...
subroutine, public calculate_orthonormality(orthonormality, mo_array, matrix_s)
...
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Self-consistent continuum solvation (SCCS) model implementation.
subroutine, public print_sccs_results(energy, sccs_control, output_unit)
Print SCCS results.
subroutine, public qs_scf_print_summary(output_unit, qs_env)
writes a summary of information after scf
subroutine, public qs_scf_cdft_constraint_info(output_unit, cdft_control)
writes CDFT constraint information
subroutine, public qs_scf_loop_print(qs_env, scf_env, para_env)
collects the 'heavy duty' printing tasks out of the SCF loop
subroutine, public qs_scf_initial_info(output_unit, mos, dft_control)
writes basic information at the beginning of an scf run
subroutine, public qs_scf_outer_loop_info(output_unit, scf_control, scf_env, energy, total_steps, should_stop, outer_loop_converged)
writes basic information obtained in a scf outer loop step
subroutine, public qs_scf_write_mos(qs_env, scf_env, final_mos)
Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit.
subroutine, public qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
writes basic information obtained in a scf step
subroutine, public qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, energy, total_steps, should_stop, outer_loop_converged, cdft_loop)
writes CDFT constraint information and optionally CDFT scf loop info
subroutine, public qs_scf_cdft_initial_info(output_unit, cdft_control)
writes information about the CDFT env
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
integer, parameter, public special_diag_method_nr
parameters that control an scf iteration