83#include "./base/base_uses.f90"
89 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_scf_output'
109 INTEGER,
INTENT(IN) :: output_unit
112 INTEGER :: nelectron_total
113 LOGICAL :: gapw, gapw_xc, qmmm
120 NULLIFY (rho, energy, dft_control, scf_env, qs_charges)
121 CALL get_qs_env(qs_env=qs_env, rho=rho, energy=energy, dft_control=dft_control, &
122 scf_env=scf_env, qs_charges=qs_charges)
124 gapw = dft_control%qs_control%gapw
125 gapw_xc = dft_control%qs_control%gapw_xc
127 nelectron_total = scf_env%nelectron
129 CALL qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
130 dft_control, qmmm, qs_env, gapw, gapw_xc)
142 INTEGER :: output_unit
145 INTEGER,
INTENT(IN) :: ndep
147 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_initial_info'
149 INTEGER :: handle, homo, ispin, nao, &
152 CALL timeset(routinen, handle)
154 IF (output_unit > 0)
THEN
155 DO ispin = 1, dft_control%nspins
158 nelectron=nelectron_spin, &
161 IF (dft_control%nspins > 1)
THEN
162 WRITE (unit=output_unit, fmt=
"(/,T2,A,I2)")
"Spin", ispin
164 WRITE (unit=output_unit, fmt=
"(/,(T2,A,T71,I10))") &
165 "Number of electrons:", nelectron_spin, &
166 "Number of occupied orbitals:", homo, &
167 "Number of molecular orbitals:", nmo
169 WRITE (unit=output_unit, fmt=
"(/,(T2,A,T71,I10))") &
170 "Number of orbital functions:", nao, &
171 "Number of independent orbital functions:", nao - ndep
174 CALL timestop(handle)
189 LOGICAL,
INTENT(IN) :: final_mos
191 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_write_mos'
193 CHARACTER(LEN=2) :: solver_method
194 CHARACTER(LEN=3*default_string_length) :: message
195 CHARACTER(LEN=5) :: spin
196 CHARACTER(LEN=default_string_length), &
197 DIMENSION(:),
POINTER :: tmpstringlist
198 INTEGER :: handle, homo, ikp, ispin, iw, kpoint, &
199 nao, nelectron, nkp, nmo, nspin, numo
200 INTEGER,
DIMENSION(2) :: nmos_occ
201 INTEGER,
DIMENSION(:),
POINTER :: mo_index_range
202 LOGICAL :: do_kpoints, do_printout, print_eigvals, &
203 print_eigvecs, print_mo_info, &
204 print_occup, print_occup_stats
205 REAL(kind=
dp) :: flexible_electron_count, maxocc, n_el_f, &
206 occup_stats_occ_threshold
207 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_eigenvalues, umo_eigenvalues
211 TYPE(
cp_fm_type),
POINTER :: mo_coeff, umo_coeff
214 TYPE(
dbcsr_type),
POINTER :: matrix_ks, matrix_s
222 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
226 CALL timeset(routinen, handle)
228 cpassert(
ASSOCIATED(qs_env))
232 blacs_env=blacs_env, &
233 dft_control=dft_control, &
234 do_kpoints=do_kpoints, &
236 qs_kind_set=qs_kind_set, &
238 particle_set=particle_set, &
239 scf_control=scf_control)
246 CALL section_vals_val_get(dft_section,
"PRINT%MO%OCCUPATION_NUMBERS_STATS", c_vals=tmpstringlist)
248 print_occup_stats = .false.
249 occup_stats_occ_threshold = 1e-6_dp
250 IF (
SIZE(tmpstringlist) > 0)
THEN
251 print_occup_stats = .true.
252 IF (len_trim(tmpstringlist(1)) > 0) &
253 READ (tmpstringlist(1), *) print_occup_stats
255 IF (
SIZE(tmpstringlist) > 1) &
256 READ (tmpstringlist(2), *) occup_stats_occ_threshold
261 IF ((.NOT. print_mo_info) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup .OR. print_occup_stats)))
THEN
262 CALL timestop(handle)
266 NULLIFY (fm_struct_tmp)
268 NULLIFY (mo_eigenvalues)
271 NULLIFY (umo_eigenvalues)
275 nspin = dft_control%nspins
281 nkp =
SIZE(kpoints%kp_env)
283 CALL get_qs_env(qs_env, matrix_ks=ks, matrix_s=s)
284 cpassert(
ASSOCIATED(ks))
285 cpassert(
ASSOCIATED(s))
289 kp_loop:
DO ikp = 1, nkp
292 mos => kpoints%kp_env(ikp)%kpoint_env%mos(1, :)
295 CALL get_qs_env(qs_env, matrix_ks=ks, mos=mos)
298 cpassert(
ASSOCIATED(mos))
309 cpabort(
"The OT method is not implemented for k points")
314 matrix_ks => ks(ispin)%matrix
315 matrix_s => s(1)%matrix
318 IF (dft_control%do_admm)
THEN
326 eigenvalues=mo_eigenvalues, &
329 nelectron=nelectron, &
333 flexible_electron_count=flexible_electron_count)
337 cpassert(
ASSOCIATED(mo_index_range))
338 IF (mo_index_range(2) == -1)
THEN
341 numo = min(mo_index_range(2) - homo, nao - homo)
360 flexible_electron_count=flexible_electron_count)
362 fm_struct=fm_struct_tmp, &
363 name=
"Temporary MO set (unoccupied MOs only) for printout")
366 mo_coeff=umo_coeff, &
367 eigenvalues=umo_eigenvalues)
373 NULLIFY (local_preconditioner)
374 IF (
ASSOCIATED(scf_env%ot_preconditioner))
THEN
375 local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
377 NULLIFY (local_preconditioner)
384 matrix_c_fm=umo_coeff, &
385 matrix_orthogonal_space_fm=mo_coeff, &
386 eps_gradient=scf_control%eps_lumos, &
388 iter_max=scf_control%max_iter_lumos, &
389 size_ortho_space=nmo)
392 ks_matrix=matrix_ks, &
393 evals_arg=umo_eigenvalues, &
398 IF (dft_control%do_admm)
THEN
406 message =
"The MO information is only calculated after SCF convergence "// &
407 "is achieved when the orbital transformation (OT) method is used"
408 cpwarn(trim(message))
409 do_printout = .false.
430 cpabort(
"Invalid spin")
432 IF (
ASSOCIATED(umo_set))
THEN
434 final_mos=final_mos, spin=trim(spin), solver_method=solver_method, &
438 final_mos=final_mos, spin=trim(spin), solver_method=solver_method)
441 IF (
ASSOCIATED(umo_set))
THEN
443 final_mos=final_mos, solver_method=solver_method, &
447 final_mos=final_mos, solver_method=solver_method)
451 nmos_occ(ispin) = max(nmos_occ(ispin), count(mo_set%occupation_numbers > occup_stats_occ_threshold))
455 IF (
ASSOCIATED(umo_set))
THEN
468 IF (do_printout .AND. print_mo_info .AND. print_occup_stats)
THEN
470 ignore_should_output=print_mo_info, &
473 IF (
SIZE(mos) > 1)
THEN
474 WRITE (unit=iw, fmt=
"(A,I4)")
" MO| Total occupied (ALPHA):", nmos_occ(1)
475 WRITE (unit=iw, fmt=
"(A,I4)")
" MO| Total occupied (BETA): ", nmos_occ(2)
477 WRITE (unit=iw, fmt=
"(A,I4)")
" MO| Total occupied: ", nmos_occ(1)
479 WRITE (unit=iw, fmt=
"(A)")
""
482 ignore_should_output=print_mo_info)
485 CALL timestop(handle)
500 energy, total_steps, should_stop, outer_loop_converged)
501 INTEGER :: output_unit
505 INTEGER :: total_steps
506 LOGICAL,
INTENT(IN) :: should_stop, outer_loop_converged
508 REAL(kind=
dp) :: outer_loop_eps
510 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
511 IF (output_unit > 0)
WRITE (output_unit,
'(/,T3,A,I4,A,E10.2,A,F22.10)') &
512 "outer SCF iter = ", scf_env%outer_scf%iter_count, &
513 " RMS gradient = ", outer_loop_eps,
" energy =", energy%total
515 IF (outer_loop_converged)
THEN
516 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
517 "outer SCF loop converged in", scf_env%outer_scf%iter_count, &
518 " iterations or ", total_steps,
" steps"
519 ELSE IF (scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf &
520 .OR. should_stop)
THEN
521 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
522 "outer SCF loop FAILED to converge after ", &
523 scf_env%outer_scf%iter_count,
" iterations or ", total_steps,
" steps"
540 INTEGER :: output_unit
541 LOGICAL :: just_energy
542 REAL(kind=
dp) :: t1, t2
545 IF ((output_unit > 0) .AND. scf_env%print_iter_line)
THEN
546 IF (just_energy)
THEN
547 WRITE (unit=output_unit, &
548 fmt=
"(T2,A,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") &
549 " -", trim(scf_env%iter_method), scf_env%iter_param, t2 - t1, energy%total
551 IF ((abs(scf_env%iter_delta) < 1.0e-8_dp) .OR. &
552 (abs(scf_env%iter_delta) >= 1.0e5_dp))
THEN
553 WRITE (unit=output_unit, &
554 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") &
555 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
556 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
558 WRITE (unit=output_unit, &
559 fmt=
"(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
560 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
561 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
586 SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
587 dft_control, qmmm, qs_env, gapw, gapw_xc)
588 INTEGER,
INTENT(IN) :: output_unit
592 INTEGER,
INTENT(IN) :: nelectron_total
594 LOGICAL,
INTENT(IN) :: qmmm
596 LOGICAL,
INTENT(IN) :: gapw, gapw_xc
598 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_print_scf_summary'
600 INTEGER :: bc, handle, ispin, psolver
601 REAL(kind=
dp) :: exc1_energy, exc_energy, &
602 implicit_ps_ehartree, tot1_h, tot1_s
603 REAL(kind=
dp),
DIMENSION(:),
POINTER :: tot_rho_r
606 NULLIFY (tot_rho_r, pw_env)
607 CALL timeset(routinen, handle)
610 psolver = pw_env%poisson_env%parameters%solver
612 IF (output_unit > 0)
THEN
614 IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
615 dft_control%qs_control%xtb .OR. &
616 dft_control%qs_control%dftb))
THEN
617 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T41,2F20.10))") &
618 "Electronic density on regular grids: ", &
621 "Core density on regular grids:", &
622 qs_charges%total_rho_core_rspace, &
623 qs_charges%total_rho_core_rspace - real(nelectron_total + dft_control%charge,
dp)
625 IF (dft_control%correct_surf_dip)
THEN
626 WRITE (unit=output_unit, fmt=
"((T3,A,/,T3,A,T41,F20.10))") &
627 "Total dipole moment perpendicular to ", &
628 "the slab [electrons-Angstroem]: ", &
629 qs_env%surface_dipole_moment
633 tot1_h = qs_charges%total_rho1_hard(1)
634 tot1_s = qs_charges%total_rho1_soft(1)
635 DO ispin = 2, dft_control%nspins
636 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
637 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
639 WRITE (unit=output_unit, fmt=
"((T3,A,T41,2F20.10))") &
640 "Hard and soft densities (Lebedev):", &
642 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
643 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
645 "Total charge density (r-space): ", &
647 + qs_charges%total_rho_core_rspace, &
648 "Total Rho_soft + Rho0_soft (g-space):", &
649 qs_charges%total_rho_gspace
651 WRITE (unit=output_unit, fmt=
"(T3,A,T41,F20.10)") &
652 "Total charge density on r-space grids: ", &
654 qs_charges%total_rho_core_rspace, &
655 "Total charge density g-space grids: ", &
656 qs_charges%total_rho_gspace
659 IF (dft_control%qs_control%semi_empirical)
THEN
660 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
661 "Core-core repulsion energy [eV]: ", energy%core_overlap*
evolt, &
662 "Core Hamiltonian energy [eV]: ", energy%core*
evolt, &
663 "Two-electron integral energy [eV]: ", energy%hartree*
evolt, &
664 "Electronic energy [eV]: ", &
665 (energy%core + 0.5_dp*energy%hartree)*
evolt
666 IF (energy%dispersion /= 0.0_dp) &
667 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
668 "Dispersion energy [eV]: ", energy%dispersion*
evolt
669 ELSEIF (dft_control%qs_control%dftb)
THEN
670 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
671 "Core Hamiltonian energy: ", energy%core, &
672 "Repulsive potential energy: ", energy%repulsive, &
673 "Electronic energy: ", energy%hartree, &
674 "Dispersion energy: ", energy%dispersion
675 IF (energy%dftb3 /= 0.0_dp) &
676 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
677 "DFTB3 3rd order energy: ", energy%dftb3
678 IF (energy%efield /= 0.0_dp) &
679 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
680 "Electric field interaction energy: ", energy%efield
681 ELSEIF (dft_control%qs_control%xtb)
THEN
682 IF (dft_control%qs_control%xtb_control%do_tblite)
THEN
683 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
684 "Core Hamiltonian energy: ", energy%core, &
685 "Repulsive potential energy: ", energy%repulsive, &
686 "Electrostatic energy: ", energy%el_stat, &
687 "Self-consistent dispersion energy: ", energy%dispersion_sc, &
688 "Non-self consistent dispersion energy: ", energy%dispersion, &
689 "Correction for halogen bonding: ", energy%xtb_xb_inter
691 IF (dft_control%qs_control%xtb_control%gfn_type == 0)
THEN
692 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
693 "Core Hamiltonian energy: ", energy%core, &
694 "Repulsive potential energy: ", energy%repulsive, &
695 "SRB Correction energy: ", energy%srb, &
696 "Charge equilibration energy: ", energy%eeq, &
697 "Dispersion energy: ", energy%dispersion
698 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 1)
THEN
699 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
700 "Core Hamiltonian energy: ", energy%core, &
701 "Repulsive potential energy: ", energy%repulsive, &
702 "Electronic energy: ", energy%hartree, &
703 "DFTB3 3rd order energy: ", energy%dftb3, &
704 "Dispersion energy: ", energy%dispersion
705 IF (dft_control%qs_control%xtb_control%xb_interaction) &
706 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
707 "Correction for halogen bonding: ", energy%xtb_xb_inter
708 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 2)
THEN
709 cpabort(
"gfn_typ 2 NYA")
711 cpabort(
"invalid gfn_typ")
714 IF (dft_control%qs_control%xtb_control%do_nonbonded) &
715 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
716 "Correction for nonbonded interactions: ", energy%xtb_nonbonded
717 IF (energy%efield /= 0.0_dp) &
718 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
719 "Electric field interaction energy: ", energy%efield
721 IF (dft_control%do_admm)
THEN
722 exc_energy = energy%exc + energy%exc_aux_fit
723 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
725 exc_energy = energy%exc
726 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
730 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
731 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
734 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
735 "Overlap energy of the core charge distribution:", energy%core_overlap, &
736 "Self energy of the core charge distribution: ", energy%core_self, &
737 "Core Hamiltonian energy: ", energy%core, &
738 "Hartree energy: ", implicit_ps_ehartree, &
739 "Electric enthalpy: ", energy%hartree, &
740 "Exchange-correlation energy: ", exc_energy
742 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
743 "Overlap energy of the core charge distribution:", energy%core_overlap, &
744 "Self energy of the core charge distribution: ", energy%core_self, &
745 "Core Hamiltonian energy: ", energy%core, &
746 "Hartree energy: ", energy%hartree, &
747 "Exchange-correlation energy: ", exc_energy
750 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
751 "Overlap energy of the core charge distribution:", energy%core_overlap, &
752 "Self energy of the core charge distribution: ", energy%core_self, &
753 "Core Hamiltonian energy: ", energy%core, &
754 "Hartree energy: ", energy%hartree, &
755 "Exchange-correlation energy: ", exc_energy
757 IF (energy%e_hartree /= 0.0_dp) &
758 WRITE (unit=output_unit, fmt=
"(T3,A,/,T3,A,T56,F25.14)") &
759 "Coulomb Electron-Electron Interaction Energy ", &
760 "- Already included in the total Hartree term ", energy%e_hartree
761 IF (energy%ex /= 0.0_dp) &
762 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
763 "Hartree-Fock Exchange energy: ", energy%ex
764 IF (energy%dispersion /= 0.0_dp) &
765 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
766 "Dispersion energy: ", energy%dispersion
767 IF (energy%gcp /= 0.0_dp) &
768 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
769 "gCP energy: ", energy%gcp
770 IF (energy%efield /= 0.0_dp) &
771 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
772 "Electric field interaction energy: ", energy%efield
774 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
775 "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
776 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
779 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
780 "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
783 IF (dft_control%hairy_probes .EQV. .true.)
THEN
784 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
785 "Electronic entropic energy:", energy%kTS
786 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
787 "Fermi energy:", energy%efermi
789 IF (dft_control%smear)
THEN
790 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
791 "Electronic entropic energy:", energy%kTS
792 WRITE (unit=output_unit, fmt=
"((T3,A,T56,F25.14))") &
793 "Fermi energy:", energy%efermi
795 IF (dft_control%dft_plus_u)
THEN
796 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
797 "DFT+U energy:", energy%dft_plus_u
799 IF (dft_control%do_sccs)
THEN
800 WRITE (unit=output_unit, fmt=
"(A)")
""
804 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
805 "QM/MM Electrostatic energy: ", energy%qmmm_el
806 IF (qs_env%qmmm_env_qm%image_charge)
THEN
807 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
808 "QM/MM image charge energy: ", energy%image_charge
811 IF (dft_control%qs_control%mulliken_restraint)
THEN
812 WRITE (unit=output_unit, fmt=
"(T3,A,T56,F25.14)") &
813 "Mulliken restraint energy: ", energy%mulliken
815 IF (dft_control%qs_control%semi_empirical)
THEN
816 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
817 "Total energy [eV]: ", energy%total*
evolt
818 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
819 "Atomic reference energy [eV]: ", energy%core_self*
evolt, &
820 "Heat of formation [kcal/mol]: ", &
821 (energy%total + energy%core_self)*
kcalmol
823 WRITE (unit=output_unit, fmt=
"(/,(T3,A,T56,F25.14))") &
824 "Total energy: ", energy%total
827 IF (qs_env%qmmm_env_qm%image_charge)
THEN
834 CALL timestop(handle)
836 END SUBROUTINE qs_scf_print_scf_summary
851 CHARACTER(LEN=*),
PARAMETER :: routinen =
'qs_scf_loop_print'
853 INTEGER :: after, handle, ic, ispin, iw
854 LOGICAL :: do_kpoints, omit_headers
855 REAL(kind=
dp) :: mo_mag_max, mo_mag_min, orthonormality
857 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, matrix_p, matrix_s
864 CALL timeset(routinen, handle)
866 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
867 do_kpoints=do_kpoints)
873 DO ispin = 1, dft_control%nspins
876 dft_section,
"PRINT%AO_MATRICES/DENSITY"),
cp_p_file))
THEN
882 after = min(max(after, 1), 16)
883 DO ic = 1,
SIZE(matrix_p, 2)
885 output_unit=iw, omit_headers=omit_headers)
888 "PRINT%AO_MATRICES/DENSITY")
892 dft_section,
"PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"),
cp_p_file))
THEN
896 after = min(max(after, 1), 16)
897 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
898 DO ic = 1,
SIZE(matrix_ks, 2)
899 IF (dft_control%qs_control%semi_empirical)
THEN
901 scale=
evolt, output_unit=iw, omit_headers=omit_headers)
904 output_unit=iw, omit_headers=omit_headers)
908 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
914 scf_section,
"PRINT%MO_ORTHONORMALITY"),
cp_p_file))
THEN
919 WRITE (iw,
'(T8,A)') &
920 " K-points: Maximum deviation from MO S-orthonormality not determined"
923 "PRINT%MO_ORTHONORMALITY")
929 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
935 WRITE (iw,
'(T8,A,T61,E20.4)') &
936 " Maximum deviation from MO S-orthonormality", orthonormality
939 "PRINT%MO_ORTHONORMALITY")
943 scf_section,
"PRINT%MO_MAGNITUDE"),
cp_p_file))
THEN
948 WRITE (iw,
'(T8,A)') &
949 " K-points: Minimum/Maximum MO magnitude not determined"
952 "PRINT%MO_MAGNITUDE")
959 WRITE (iw,
'(T8,A,T41,2E20.4)') &
960 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
963 "PRINT%MO_MAGNITUDE")
967 CALL timestop(handle)
986 energy, total_steps, should_stop, outer_loop_converged, &
988 INTEGER :: output_unit
993 INTEGER :: total_steps
994 LOGICAL,
INTENT(IN) :: should_stop, outer_loop_converged, &
997 REAL(kind=
dp) :: outer_loop_eps
1000 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
1001 IF (output_unit > 0)
WRITE (output_unit,
'(/,T3,A,I4,A,E10.2,A,F22.10)') &
1002 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
1003 " RMS gradient = ", outer_loop_eps,
" energy =", energy%total
1004 IF (outer_loop_converged)
THEN
1005 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
1006 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
1007 " iterations or ", total_steps,
" steps"
1009 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
1010 .AND. .NOT. outer_loop_converged)
THEN
1011 IF (output_unit > 0)
WRITE (output_unit,
'(T3,A,I4,A,I4,A,/)') &
1012 "CDFT SCF loop FAILED to converge after ", &
1013 scf_env%outer_scf%iter_count,
" iterations or ", total_steps,
" steps"
1028 INTEGER :: output_unit
1031 IF (output_unit > 0)
THEN
1032 WRITE (output_unit,
'(/,A)') &
1033 " ---------------------------------- CDFT --------------------------------------"
1034 WRITE (output_unit,
'(A)') &
1035 " Optimizing a density constraint in an external SCF loop "
1036 WRITE (output_unit,
'(A)')
" "
1037 SELECT CASE (cdft_control%type)
1039 WRITE (output_unit,
'(A)')
" Type of constraint: Hirshfeld"
1041 WRITE (output_unit,
'(A)')
" Type of constraint: Becke"
1043 WRITE (output_unit,
'(A,I8)')
" Number of constraints: ",
SIZE(cdft_control%group)
1044 WRITE (output_unit,
'(A,L8)')
" Using fragment densities:", cdft_control%fragment_density
1045 WRITE (output_unit,
'(A)')
" "
1046 IF (cdft_control%atomic_charges)
WRITE (output_unit,
'(A,/)')
" Calculating atomic CDFT charges"
1047 SELECT CASE (cdft_control%constraint_control%optimizer)
1049 WRITE (output_unit,
'(A)') &
1050 " Minimizer : SD : steepest descent"
1052 WRITE (output_unit,
'(A)') &
1053 " Minimizer : DIIS : direct inversion"
1054 WRITE (output_unit,
'(A)') &
1055 " in the iterative subspace"
1056 WRITE (output_unit,
'(A,I3,A)') &
1058 cdft_control%constraint_control%diis_buffer_length,
" DIIS vectors"
1060 WRITE (output_unit,
'(A)') &
1061 " Minimizer : BISECT : gradient bisection"
1062 WRITE (output_unit,
'(A,I3)') &
1063 " using a trust count of", &
1064 cdft_control%constraint_control%bisect_trust_count
1068 cdft_control%constraint_control%optimizer, output_unit)
1070 WRITE (output_unit,
'(A)')
" Minimizer : Secant"
1074 WRITE (output_unit,
'(/,A,L7)') &
1075 " Reusing OT preconditioner: ", cdft_control%reuse_precond
1076 IF (cdft_control%reuse_precond)
THEN
1077 WRITE (output_unit,
'(A,I3,A,I3,A)') &
1078 " using old preconditioner for up to ", &
1079 cdft_control%max_reuse,
" subsequent CDFT SCF"
1080 WRITE (output_unit,
'(A,I3,A,I3,A)') &
1081 " iterations if the relevant loop converged in less than ", &
1082 cdft_control%precond_freq,
" steps"
1084 SELECT CASE (cdft_control%type)
1086 WRITE (output_unit,
'(/,A)')
" Hirshfeld constraint settings"
1087 WRITE (output_unit,
'(A)')
" "
1088 SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1090 WRITE (output_unit,
'(A, A8)') &
1091 " Shape function type: ",
"Gaussian"
1092 WRITE (output_unit,
'(A)', advance=
'NO') &
1093 " Type of Gaussian: "
1094 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1096 WRITE (output_unit,
'(A13)')
"Default"
1098 WRITE (output_unit,
'(A13)')
"Covalent"
1100 WRITE (output_unit,
'(A13)')
"Fixed radius"
1102 WRITE (output_unit,
'(A13)')
"Van der Waals"
1104 WRITE (output_unit,
'(A13)')
"User-defined"
1108 WRITE (output_unit,
'(A, A8)') &
1109 " Shape function type: ",
"Density"
1112 WRITE (output_unit,
'(/, A)')
" Becke constraint settings"
1113 WRITE (output_unit,
'(A)')
" "
1114 SELECT CASE (cdft_control%becke_control%cutoff_type)
1116 WRITE (output_unit,
'(A,F8.3,A)') &
1117 " Cutoff for partitioning :",
cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1118 "angstrom"),
" angstrom"
1120 WRITE (output_unit,
'(A)') &
1121 " Using element specific cutoffs for partitioning"
1123 WRITE (output_unit,
'(A,L7)') &
1124 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1125 WRITE (output_unit,
'(A,L7)') &
1126 " Precompute gradients : ", cdft_control%becke_control%in_memory
1127 WRITE (output_unit,
'(A)')
" "
1128 IF (cdft_control%becke_control%adjust) &
1129 WRITE (output_unit,
'(A)') &
1130 " Using atomic radii to generate a heteronuclear charge partitioning"
1131 WRITE (output_unit,
'(A)')
" "
1132 IF (.NOT. cdft_control%becke_control%cavity_confine)
THEN
1133 WRITE (output_unit,
'(A)') &
1134 " No confinement is active"
1136 WRITE (output_unit,
'(A)')
" Confinement using a Gaussian shaped cavity is active"
1137 SELECT CASE (cdft_control%becke_control%cavity_shape)
1139 WRITE (output_unit,
'(A,F8.4, A)') &
1140 " Type of Gaussian : Fixed radius: ", &
1143 WRITE (output_unit,
'(A)') &
1144 " Type of Gaussian : Covalent radius "
1146 WRITE (output_unit,
'(A)') &
1147 " Type of Gaussian : vdW radius "
1149 WRITE (output_unit,
'(A)') &
1150 " Type of Gaussian : User radius "
1152 WRITE (output_unit,
'(A,ES12.4)') &
1153 " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1156 WRITE (output_unit,
'(/,A)') &
1157 " ---------------------------------- CDFT --------------------------------------"
1170 INTEGER :: output_unit
1175 IF (output_unit > 0)
THEN
1176 SELECT CASE (cdft_control%type)
1178 WRITE (output_unit,
'(/,T3,A,T60)') &
1179 '------------------- Hirshfeld constraint information -------------------'
1181 WRITE (output_unit,
'(/,T3,A,T60)') &
1182 '--------------------- Becke constraint information ---------------------'
1184 cpabort(
"Unknown CDFT constraint.")
1186 DO igroup = 1,
SIZE(cdft_control%target)
1187 IF (igroup > 1)
WRITE (output_unit,
'(T3,A)')
' '
1188 WRITE (output_unit,
'(T3,A,T54,(3X,I18))') &
1189 'Atomic group :', igroup
1190 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1192 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1193 WRITE (output_unit,
'(T3,A,T42,A)') &
1194 'Type of constraint :', adjustr(
'Charge density constraint (frag.)')
1196 WRITE (output_unit,
'(T3,A,T50,A)') &
1197 'Type of constraint :', adjustr(
'Charge density constraint')
1200 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1201 WRITE (output_unit,
'(T3,A,T35,A)') &
1202 'Type of constraint :', adjustr(
'Magnetization density constraint (frag.)')
1204 WRITE (output_unit,
'(T3,A,T43,A)') &
1205 'Type of constraint :', adjustr(
'Magnetization density constraint')
1208 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1209 WRITE (output_unit,
'(T3,A,T38,A)') &
1210 'Type of constraint :', adjustr(
'Alpha spin density constraint (frag.)')
1212 WRITE (output_unit,
'(T3,A,T46,A)') &
1213 'Type of constraint :', adjustr(
'Alpha spin density constraint')
1216 IF (cdft_control%group(igroup)%is_fragment_constraint)
THEN
1217 WRITE (output_unit,
'(T3,A,T39,A)') &
1218 'Type of constraint :', adjustr(
'Beta spin density constraint (frag.)')
1220 WRITE (output_unit,
'(T3,A,T47,A)') &
1221 'Type of constraint :', adjustr(
'Beta spin density constraint')
1224 cpabort(
"Unknown constraint type.")
1226 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1227 'Target value of constraint :', cdft_control%target(igroup)
1228 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1229 'Current value of constraint :', cdft_control%value(igroup)
1230 WRITE (output_unit,
'(T3,A,T59,(3X,ES13.3))') &
1231 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1232 WRITE (output_unit,
'(T3,A,T54,(3X,F18.12))') &
1233 'Strength of constraint :', cdft_control%strength(igroup)
1235 WRITE (output_unit,
'(T3,A)') &
1236 '------------------------------------------------------------------------'
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)
...
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, do_rixs, 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, 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_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_initial_info(output_unit, mos, dft_control, ndep)
writes basic information at the beginning of an scf run
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
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.