84#include "./base/base_uses.f90"
90 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_scf_output'
110 INTEGER,
INTENT(IN) :: output_unit
113 INTEGER :: nelectron_total
114 LOGICAL :: gapw, gapw_xc, qmmm
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
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
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
203 TYPE(
cp_fm_type),
POINTER :: mo_coeff, umo_coeff
206 TYPE(
dbcsr_type),
POINTER :: matrix_ks, matrix_s
214 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
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)
386 ks_matrix=matrix_ks, &
387 evals_arg=umo_eigenvalues, &
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
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"
530 INTEGER :: output_unit
531 LOGICAL :: just_energy
532 REAL(kind=
dp) :: t1, t2
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,A,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") &
539 " -", trim(scf_env%iter_method), scf_env%iter_param, t2 - t1, energy%total
541 IF ((abs(scf_env%iter_delta) < 1.0e-8_dp) .OR. &
542 (abs(scf_env%iter_delta) >= 1.0e5_dp))
THEN
543 WRITE (unit=output_unit, &
544 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") &
545 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
546 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
548 WRITE (unit=output_unit, &
549 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
550 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
551 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
576 SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
577 dft_control, qmmm, qs_env, gapw, gapw_xc)
578 INTEGER,
INTENT(IN) :: output_unit
582 INTEGER,
INTENT(IN) :: nelectron_total
584 LOGICAL,
INTENT(IN) :: qmmm
586 LOGICAL,
INTENT(IN) :: gapw, gapw_xc
588 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_print_scf_summary'
590 INTEGER :: bc, handle, ispin, psolver
591 REAL(kind=
dp) :: exc1_energy, exc_energy, &
592 implicit_ps_ehartree, tot1_h, tot1_s
593 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tot_rho_r
596 NULLIFY (tot_rho_r, pw_env)
597 CALL timeset(routinen, handle)
600 psolver = pw_env%poisson_env%parameters%solver
602 IF (output_unit > 0)
THEN
604 IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
605 dft_control%qs_control%xtb .OR. &
606 dft_control%qs_control%dftb))
THEN
607 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T41,2F20.10))") &
608 "Electronic density on regular grids: ", &
611 "Core density on regular grids:", &
612 qs_charges%total_rho_core_rspace, &
613 qs_charges%total_rho_core_rspace - real(nelectron_total + dft_control%charge,
dp)
615 IF (dft_control%correct_surf_dip)
THEN
616 WRITE (unit=output_unit, fmt=
"((T3,A,/,T3,A,T41,F20.10))") &
617 "Total dipole moment perpendicular to ", &
618 "the slab [electrons-Angstroem]: ", &
619 qs_env%surface_dipole_moment
623 tot1_h = qs_charges%total_rho1_hard(1)
624 tot1_s = qs_charges%total_rho1_soft(1)
625 DO ispin = 2, dft_control%nspins
626 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
627 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
629 WRITE (unit=output_unit, fmt=
"((T3,A,T41,2F20.10))") &
630 "Hard and soft densities (Lebedev):", &
632 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
633 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
635 "Total charge density (r-space): ", &
637 + qs_charges%total_rho_core_rspace, &
638 "Total Rho_soft + Rho0_soft (g-space):", &
639 qs_charges%total_rho_gspace
641 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
642 "Total charge density on r-space grids: ", &
644 qs_charges%total_rho_core_rspace, &
645 "Total charge density g-space grids: ", &
646 qs_charges%total_rho_gspace
649 IF (dft_control%qs_control%semi_empirical)
THEN
650 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
651 "Core-core repulsion energy [eV]: ", energy%core_overlap*
evolt, &
652 "Core Hamiltonian energy [eV]: ", energy%core*
evolt, &
653 "Two-electron integral energy [eV]: ", energy%hartree*
evolt, &
654 "Electronic energy [eV]: ", &
655 (energy%core + 0.5_dp*energy%hartree)*
evolt
656 IF (energy%dispersion /= 0.0_dp) &
657 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
658 "Dispersion energy [eV]: ", energy%dispersion*
evolt
659 ELSEIF (dft_control%qs_control%dftb)
THEN
660 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
661 "Core Hamiltonian energy: ", energy%core, &
662 "Repulsive potential energy: ", energy%repulsive, &
663 "Electronic energy: ", energy%hartree, &
664 "Dispersion energy: ", energy%dispersion
665 IF (energy%dftb3 /= 0.0_dp) &
666 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
667 "DFTB3 3rd order energy: ", energy%dftb3
668 IF (energy%efield /= 0.0_dp) &
669 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
670 "Electric field interaction energy: ", energy%efield
671 ELSEIF (dft_control%qs_control%xtb)
THEN
672 IF (dft_control%qs_control%xtb_control%do_tblite)
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 "Electrostatic energy: ", energy%el_stat, &
677 "Non-self consistent dispersion energy: ", energy%dispersion, &
678 "Self-consistent dispersion energy: ", energy%dispersion_sc, &
679 "Correction for halogen bonding: ", energy%xtb_xb_inter
681 IF (dft_control%qs_control%xtb_control%gfn_type == 0)
THEN
682 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
683 "Core Hamiltonian energy: ", energy%core, &
684 "Repulsive potential energy: ", energy%repulsive, &
685 "SRB Correction energy: ", energy%srb, &
686 "Charge equilibration energy: ", energy%eeq, &
687 "Dispersion energy: ", energy%dispersion
688 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 1)
THEN
689 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
690 "Core Hamiltonian energy: ", energy%core, &
691 "Repulsive potential energy: ", energy%repulsive, &
692 "Electronic energy: ", energy%hartree, &
693 "DFTB3 3rd order energy: ", energy%dftb3, &
694 "Dispersion energy: ", energy%dispersion
695 IF (dft_control%qs_control%xtb_control%xb_interaction) &
696 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
697 "Correction for halogen bonding: ", energy%xtb_xb_inter
698 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 2)
THEN
699 cpabort(
"gfn_typ 2 NYA")
701 cpabort(
"invalid gfn_typ")
704 IF (dft_control%qs_control%xtb_control%do_nonbonded) &
705 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
706 "Correction for nonbonded interactions: ", energy%xtb_nonbonded
707 IF (energy%efield /= 0.0_dp) &
708 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
709 "Electric field interaction energy: ", energy%efield
711 IF (dft_control%do_admm)
THEN
712 exc_energy = energy%exc + energy%exc_aux_fit
713 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
715 exc_energy = energy%exc
716 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
720 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
721 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
724 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
725 "Overlap energy of the core charge distribution:", energy%core_overlap, &
726 "Self energy of the core charge distribution: ", energy%core_self, &
727 "Core Hamiltonian energy: ", energy%core, &
728 "Hartree energy: ", implicit_ps_ehartree, &
729 "Electric enthalpy: ", energy%hartree, &
730 "Exchange-correlation energy: ", exc_energy
732 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
733 "Overlap energy of the core charge distribution:", energy%core_overlap, &
734 "Self energy of the core charge distribution: ", energy%core_self, &
735 "Core Hamiltonian energy: ", energy%core, &
736 "Hartree energy: ", energy%hartree, &
737 "Exchange-correlation energy: ", exc_energy
740 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
741 "Overlap energy of the core charge distribution:", energy%core_overlap, &
742 "Self energy of the core charge distribution: ", energy%core_self, &
743 "Core Hamiltonian energy: ", energy%core, &
744 "Hartree energy: ", energy%hartree, &
745 "Exchange-correlation energy: ", exc_energy
747 IF (energy%e_hartree /= 0.0_dp) &
748 WRITE (unit=output_unit, fmt=
"(T3,A,/,T3,A,T56,F25.14)") &
749 "Coulomb Electron-Electron Interaction Energy ", &
750 "- Already included in the total Hartree term ", energy%e_hartree
751 IF (energy%ex /= 0.0_dp) &
752 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
753 "Hartree-Fock Exchange energy: ", energy%ex
754 IF (energy%dispersion /= 0.0_dp) &
755 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
756 "Dispersion energy: ", energy%dispersion
757 IF (energy%gcp /= 0.0_dp) &
758 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
759 "gCP energy: ", energy%gcp
760 IF (energy%efield /= 0.0_dp) &
761 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
762 "Electric field interaction energy: ", energy%efield
764 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
765 "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
766 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
769 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
770 "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
773 IF (dft_control%hairy_probes .EQV. .true.)
THEN
774 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
775 "Electronic entropic energy:", energy%kTS
776 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
777 "Fermi energy:", energy%efermi
779 IF (dft_control%smear)
THEN
780 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
781 "Electronic entropic energy:", energy%kTS
782 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
783 "Fermi energy:", energy%efermi
785 IF (dft_control%dft_plus_u)
THEN
786 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
787 "DFT+U energy:", energy%dft_plus_u
789 IF (dft_control%do_sccs)
THEN
790 WRITE (unit=output_unit, fmt=
"(A)")
""
794 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
795 "QM/MM Electrostatic energy: ", energy%qmmm_el
796 IF (qs_env%qmmm_env_qm%image_charge)
THEN
797 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
798 "QM/MM image charge energy: ", energy%image_charge
801 IF (dft_control%qs_control%mulliken_restraint)
THEN
802 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
803 "Mulliken restraint energy: ", energy%mulliken
805 IF (dft_control%qs_control%semi_empirical)
THEN
806 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
807 "Total energy [eV]: ", energy%total*
evolt
808 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
809 "Atomic reference energy [eV]: ", energy%core_self*
evolt, &
810 "Heat of formation [kcal/mol]: ", &
811 (energy%total + energy%core_self)*
kcalmol
813 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
814 "Total energy: ", energy%total
817 IF (qs_env%qmmm_env_qm%image_charge)
THEN
824 CALL timestop(handle)
826 END SUBROUTINE qs_scf_print_scf_summary
841 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_loop_print'
843 INTEGER :: after, handle, ic, ispin, iw
844 LOGICAL :: do_kpoints, omit_headers
845 REAL(kind=
dp) :: mo_mag_max, mo_mag_min, orthonormality
847 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_p, matrix_s
854 CALL timeset(routinen, handle)
856 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
857 do_kpoints=do_kpoints)
863 DO ispin = 1, dft_control%nspins
866 dft_section,
"PRINT%AO_MATRICES/DENSITY"),
cp_p_file))
THEN
872 after = min(max(after, 1), 16)
873 DO ic = 1,
SIZE(matrix_p, 2)
875 output_unit=iw, omit_headers=omit_headers)
878 "PRINT%AO_MATRICES/DENSITY")
882 dft_section,
"PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"),
cp_p_file))
THEN
886 after = min(max(after, 1), 16)
887 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
888 DO ic = 1,
SIZE(matrix_ks, 2)
889 IF (dft_control%qs_control%semi_empirical)
THEN
891 scale=
evolt, output_unit=iw, omit_headers=omit_headers)
894 output_unit=iw, omit_headers=omit_headers)
898 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
904 scf_section,
"PRINT%MO_ORTHONORMALITY"),
cp_p_file))
THEN
909 WRITE (iw,
'(T8,A)') &
910 " K-points: Maximum deviation from MO S-orthonormality not determined"
913 "PRINT%MO_ORTHONORMALITY")
919 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
925 WRITE (iw,
'(T8,A,T61,E20.4)') &
926 " Maximum deviation from MO S-orthonormality", orthonormality
929 "PRINT%MO_ORTHONORMALITY")
933 scf_section,
"PRINT%MO_MAGNITUDE"),
cp_p_file))
THEN
938 WRITE (iw,
'(T8,A)') &
939 " K-points: Minimum/Maximum MO magnitude not determined"
942 "PRINT%MO_MAGNITUDE")
949 WRITE (iw,
'(T8,A,T41,2E20.4)') &
950 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
953 "PRINT%MO_MAGNITUDE")
957 CALL timestop(handle)
976 energy, total_steps, should_stop, outer_loop_converged, &
978 INTEGER :: output_unit
983 INTEGER :: total_steps
984 LOGICAL,
INTENT(IN) :: should_stop, outer_loop_converged, &
987 REAL(kind=
dp) :: outer_loop_eps
990 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
991 IF (output_unit > 0)
WRITE (output_unit,
'(/,T3,A,I4,A,E10.2,A,F22.10)') &
992 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
993 " RMS gradient = ", outer_loop_eps,
" energy =", energy%total
994 IF (outer_loop_converged)
THEN
995 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
996 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
997 " iterations or ", total_steps,
" steps"
999 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
1000 .AND. .NOT. outer_loop_converged)
THEN
1001 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
1002 "CDFT SCF loop FAILED to converge after ", &
1003 scf_env%outer_scf%iter_count,
" iterations or ", total_steps,
" steps"
1018 INTEGER :: output_unit
1021 IF (output_unit > 0)
THEN
1022 WRITE (output_unit,
'(/,A)') &
1023 " ---------------------------------- CDFT --------------------------------------"
1024 WRITE (output_unit,
'(A)') &
1025 " Optimizing a density constraint in an external SCF loop "
1026 WRITE (output_unit,
'(A)')
" "
1027 SELECT CASE (cdft_control%type)
1029 WRITE (output_unit,
'(A)')
" Type of constraint: Hirshfeld"
1031 WRITE (output_unit,
'(A)')
" Type of constraint: Becke"
1033 WRITE (output_unit,
'(A,I8)')
" Number of constraints: ",
SIZE(cdft_control%group)
1034 WRITE (output_unit,
'(A,L8)')
" Using fragment densities:", cdft_control%fragment_density
1035 WRITE (output_unit,
'(A)')
" "
1036 IF (cdft_control%atomic_charges)
WRITE (output_unit,
'(A,/)')
" Calculating atomic CDFT charges"
1037 SELECT CASE (cdft_control%constraint_control%optimizer)
1039 WRITE (output_unit,
'(A)') &
1040 " Minimizer : SD : steepest descent"
1042 WRITE (output_unit,
'(A)') &
1043 " Minimizer : DIIS : direct inversion"
1044 WRITE (output_unit,
'(A)') &
1045 " in the iterative subspace"
1046 WRITE (output_unit,
'(A,I3,A)') &
1048 cdft_control%constraint_control%diis_buffer_length,
" DIIS vectors"
1050 WRITE (output_unit,
'(A)') &
1051 " Minimizer : BISECT : gradient bisection"
1052 WRITE (output_unit,
'(A,I3)') &
1053 " using a trust count of", &
1054 cdft_control%constraint_control%bisect_trust_count
1058 cdft_control%constraint_control%optimizer, output_unit)
1060 WRITE (output_unit,
'(A)')
" Minimizer : Secant"
1064 WRITE (output_unit,
'(/,A,L7)') &
1065 " Reusing OT preconditioner: ", cdft_control%reuse_precond
1066 IF (cdft_control%reuse_precond)
THEN
1067 WRITE (output_unit,
'(A,I3,A,I3,A)') &
1068 " using old preconditioner for up to ", &
1069 cdft_control%max_reuse,
" subsequent CDFT SCF"
1070 WRITE (output_unit,
'(A,I3,A,I3,A)') &
1071 " iterations if the relevant loop converged in less than ", &
1072 cdft_control%precond_freq,
" steps"
1074 SELECT CASE (cdft_control%type)
1076 WRITE (output_unit,
'(/,A)')
" Hirshfeld constraint settings"
1077 WRITE (output_unit,
'(A)')
" "
1078 SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1080 WRITE (output_unit,
'(A, A8)') &
1081 " Shape function type: ",
"Gaussian"
1082 WRITE (output_unit,
'(A)', advance=
'NO') &
1083 " Type of Gaussian: "
1084 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1086 WRITE (output_unit,
'(A13)')
"Default"
1088 WRITE (output_unit,
'(A13)')
"Covalent"
1090 WRITE (output_unit,
'(A13)')
"Fixed radius"
1092 WRITE (output_unit,
'(A13)')
"Van der Waals"
1094 WRITE (output_unit,
'(A13)')
"User-defined"
1098 WRITE (output_unit,
'(A, A8)') &
1099 " Shape function type: ",
"Density"
1102 WRITE (output_unit,
'(/, A)')
" Becke constraint settings"
1103 WRITE (output_unit,
'(A)')
" "
1104 SELECT CASE (cdft_control%becke_control%cutoff_type)
1106 WRITE (output_unit,
'(A,F8.3,A)') &
1107 " Cutoff for partitioning :",
cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1108 "angstrom"),
" angstrom"
1110 WRITE (output_unit,
'(A)') &
1111 " Using element specific cutoffs for partitioning"
1113 WRITE (output_unit,
'(A,L7)') &
1114 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1115 WRITE (output_unit,
'(A,L7)') &
1116 " Precompute gradients : ", cdft_control%becke_control%in_memory
1117 WRITE (output_unit,
'(A)')
" "
1118 IF (cdft_control%becke_control%adjust) &
1119 WRITE (output_unit,
'(A)') &
1120 " Using atomic radii to generate a heteronuclear charge partitioning"
1121 WRITE (output_unit,
'(A)')
" "
1122 IF (.NOT. cdft_control%becke_control%cavity_confine)
THEN
1123 WRITE (output_unit,
'(A)') &
1124 " No confinement is active"
1126 WRITE (output_unit,
'(A)')
" Confinement using a Gaussian shaped cavity is active"
1127 SELECT CASE (cdft_control%becke_control%cavity_shape)
1129 WRITE (output_unit,
'(A,F8.4, A)') &
1130 " Type of Gaussian : Fixed radius: ", &
1133 WRITE (output_unit,
'(A)') &
1134 " Type of Gaussian : Covalent radius "
1136 WRITE (output_unit,
'(A)') &
1137 " Type of Gaussian : vdW radius "
1139 WRITE (output_unit,
'(A)') &
1140 " Type of Gaussian : User radius "
1142 WRITE (output_unit,
'(A,ES12.4)') &
1143 " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1146 WRITE (output_unit,
'(/,A)') &
1147 " ---------------------------------- CDFT --------------------------------------"
1160 INTEGER :: output_unit
1165 IF (output_unit > 0)
THEN
1166 SELECT CASE (cdft_control%type)
1168 WRITE (output_unit,
'(/,T3,A,T60)') &
1169 '------------------- Hirshfeld constraint information -------------------'
1171 WRITE (output_unit,
'(/,T3,A,T60)') &
1172 '--------------------- Becke constraint information ---------------------'
1174 cpabort(
"Unknown CDFT constraint.")
1176 DO igroup = 1,
SIZE(cdft_control%target)
1177 IF (igroup > 1)
WRITE (output_unit,
'(T3,A)')
' '
1178 WRITE (output_unit,
'(T3,A,T54,(3X,I18))') &
1179 'Atomic group :', igroup
1180 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1182 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1183 WRITE (output_unit,
'(T3,A,T42,A)') &
1184 'Type of constraint :', adjustr(
'Charge density constraint (frag.)')
1186 WRITE (output_unit,
'(T3,A,T50,A)') &
1187 'Type of constraint :', adjustr(
'Charge density constraint')
1190 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1191 WRITE (output_unit,
'(T3,A,T35,A)') &
1192 'Type of constraint :', adjustr(
'Magnetization density constraint (frag.)')
1194 WRITE (output_unit,
'(T3,A,T43,A)') &
1195 'Type of constraint :', adjustr(
'Magnetization density constraint')
1198 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1199 WRITE (output_unit,
'(T3,A,T38,A)') &
1200 'Type of constraint :', adjustr(
'Alpha spin density constraint (frag.)')
1202 WRITE (output_unit,
'(T3,A,T46,A)') &
1203 'Type of constraint :', adjustr(
'Alpha spin density constraint')
1206 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1207 WRITE (output_unit,
'(T3,A,T39,A)') &
1208 'Type of constraint :', adjustr(
'Beta spin density constraint (frag.)')
1210 WRITE (output_unit,
'(T3,A,T47,A)') &
1211 'Type of constraint :', adjustr(
'Beta spin density constraint')
1214 cpabort(
"Unknown constraint type.")
1216 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1217 'Target value of constraint :', cdft_control%target(igroup)
1218 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1219 'Current value of constraint :', cdft_control%value(igroup)
1220 WRITE (output_unit,
'(T3,A,T59,(3X,ES13.3))') &
1221 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1222 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1223 'Strength of constraint :', cdft_control%strength(igroup)
1225 WRITE (output_unit,
'(T3,A)') &
1226 '------------------------------------------------------------------------'
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_pp, 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, harris_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, eeq, rhs, tb_tblite)
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
stores some data used in wavefunction fitting
Provides all information about an atomic kind.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
stores all the informations relevant to an mpi environment
contained for different pw related things
Container for information about total charges on the grids.
Provides all information about a quickstep kind.
keeps the density in various representations, keeping track of which ones are valid.