Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE qs_scf_output
9 : USE admm_types, ONLY: admm_type
10 : USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
11 : admm_uncorrect_for_eigenvalues
12 : USE cp_blacs_env, ONLY: cp_blacs_env_type
13 : USE cp_control_types, ONLY: dft_control_type
14 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
15 : dbcsr_type
16 : USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix
17 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
18 : cp_fm_struct_release,&
19 : cp_fm_struct_type
20 : USE cp_fm_types, ONLY: cp_fm_init_random,&
21 : cp_fm_type
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_type
24 : USE cp_output_handling, ONLY: cp_p_file,&
25 : cp_print_key_finished_output,&
26 : cp_print_key_should_output,&
27 : cp_print_key_unit_nr
28 : USE cp_units, ONLY: cp_unit_from_cp2k
29 : USE input_constants, ONLY: &
30 : becke_cutoff_element, becke_cutoff_global, cdft_alpha_constraint, cdft_beta_constraint, &
31 : cdft_charge_constraint, cdft_magnetization_constraint, ot_precond_full_all, &
32 : outer_scf_becke_constraint, outer_scf_hirshfeld_constraint, outer_scf_optimizer_bisect, &
33 : outer_scf_optimizer_broyden, outer_scf_optimizer_diis, outer_scf_optimizer_newton, &
34 : outer_scf_optimizer_newton_ls, outer_scf_optimizer_sd, outer_scf_optimizer_secant, &
35 : radius_covalent, radius_default, radius_single, radius_user, radius_vdw, &
36 : shape_function_density, shape_function_gaussian, smear_fermi_dirac, smear_gaussian, &
37 : smear_mp, smear_mv
38 : USE input_section_types, ONLY: section_get_ivals,&
39 : section_vals_get_subs_vals,&
40 : section_vals_type,&
41 : section_vals_val_get
42 : USE kahan_sum, ONLY: accurate_sum
43 : USE kinds, ONLY: default_string_length,&
44 : dp
45 : USE kpoint_types, ONLY: kpoint_type
46 : USE machine, ONLY: m_flush
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE particle_types, ONLY: particle_type
49 : USE physcon, ONLY: evolt,&
50 : kcalmol
51 : USE preconditioner_types, ONLY: preconditioner_type
52 : USE ps_implicit_types, ONLY: MIXED_BC,&
53 : MIXED_PERIODIC_BC,&
54 : NEUMANN_BC,&
55 : PERIODIC_BC
56 : USE pw_env_types, ONLY: pw_env_type
57 : USE pw_poisson_types, ONLY: pw_poisson_implicit
58 : USE qmmm_image_charge, ONLY: print_image_coefficients
59 : USE qs_cdft_opt_types, ONLY: cdft_opt_type_write
60 : USE qs_cdft_types, ONLY: cdft_control_type
61 : USE qs_charges_types, ONLY: qs_charges_type
62 : USE qs_energy_types, ONLY: qs_energy_type
63 : USE qs_environment_types, ONLY: get_qs_env,&
64 : qs_environment_type
65 : USE qs_kind_types, ONLY: qs_kind_type
66 : USE qs_mo_io, ONLY: write_mo_set_to_output_unit
67 : USE qs_mo_methods, ONLY: calculate_magnitude,&
68 : calculate_orthonormality,&
69 : calculate_subspace_eigenvalues
70 : USE qs_mo_occupation, ONLY: set_mo_occupation
71 : USE qs_mo_types, ONLY: allocate_mo_set,&
72 : deallocate_mo_set,&
73 : get_mo_set,&
74 : init_mo_set,&
75 : mo_set_type
76 : USE qs_ot_eigensolver, ONLY: ot_eigensolver
77 : USE qs_rho_types, ONLY: qs_rho_get,&
78 : qs_rho_type
79 : USE qs_sccs, ONLY: print_sccs_results
80 : USE qs_scf_types, ONLY: ot_method_nr,&
81 : qs_scf_env_type,&
82 : special_diag_method_nr
83 : USE scf_control_types, ONLY: scf_control_type
84 : #include "./base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_output'
91 :
92 : PUBLIC :: qs_scf_loop_info, &
93 : qs_scf_print_summary, &
94 : qs_scf_loop_print, &
95 : qs_scf_outer_loop_info, &
96 : qs_scf_initial_info, &
97 : qs_scf_write_mos, &
98 : qs_scf_cdft_info, &
99 : qs_scf_cdft_initial_info, &
100 : qs_scf_cdft_constraint_info
101 :
102 : CONTAINS
103 :
104 : ! **************************************************************************************************
105 : !> \brief writes a summary of information after scf
106 : !> \param output_unit ...
107 : !> \param qs_env ...
108 : ! **************************************************************************************************
109 20265 : SUBROUTINE qs_scf_print_summary(output_unit, qs_env)
110 : INTEGER, INTENT(IN) :: output_unit
111 : TYPE(qs_environment_type), POINTER :: qs_env
112 :
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
120 :
121 20265 : 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 20265 : scf_env=scf_env, qs_charges=qs_charges)
124 :
125 20265 : gapw = dft_control%qs_control%gapw
126 20265 : gapw_xc = dft_control%qs_control%gapw_xc
127 20265 : qmmm = qs_env%qmmm
128 20265 : nelectron_total = scf_env%nelectron
129 :
130 : CALL qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
131 20265 : dft_control, qmmm, qs_env, gapw, gapw_xc)
132 :
133 20265 : END SUBROUTINE qs_scf_print_summary
134 :
135 : ! **************************************************************************************************
136 : !> \brief writes basic information at the beginning of an scf run
137 : !> \param output_unit ...
138 : !> \param mos ...
139 : !> \param dft_control ...
140 : !> \param ndep ...
141 : ! **************************************************************************************************
142 21251 : SUBROUTINE qs_scf_initial_info(output_unit, mos, dft_control, ndep)
143 : INTEGER :: output_unit
144 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
145 : TYPE(dft_control_type), POINTER :: dft_control
146 : INTEGER, INTENT(IN) :: ndep
147 :
148 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_initial_info'
149 :
150 : INTEGER :: handle, homo, ispin, nao, &
151 : nelectron_spin, nmo
152 :
153 21251 : CALL timeset(routineN, handle)
154 :
155 21251 : IF (output_unit > 0) THEN
156 22922 : DO ispin = 1, dft_control%nspins
157 : CALL get_mo_set(mo_set=mos(ispin), &
158 : homo=homo, &
159 : nelectron=nelectron_spin, &
160 : nao=nao, &
161 12114 : nmo=nmo)
162 12114 : IF (dft_control%nspins > 1) THEN
163 2612 : WRITE (UNIT=output_unit, FMT="(/,T2,A,I2)") "Spin", ispin
164 : END IF
165 : WRITE (UNIT=output_unit, FMT="(/,(T2,A,T71,I10))") &
166 12114 : "Number of electrons:", nelectron_spin, &
167 12114 : "Number of occupied orbitals:", homo, &
168 47150 : "Number of molecular orbitals:", nmo
169 : END DO
170 : WRITE (UNIT=output_unit, FMT="(/,(T2,A,T71,I10))") &
171 10808 : "Number of orbital functions:", nao, &
172 21616 : "Number of independent orbital functions:", nao - ndep
173 : END IF
174 :
175 21251 : CALL timestop(handle)
176 :
177 21251 : END SUBROUTINE qs_scf_initial_info
178 :
179 : ! **************************************************************************************************
180 : !> \brief Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit
181 : !> \param qs_env ...
182 : !> \param scf_env ...
183 : !> \param final_mos ...
184 : !> \par History
185 : !> - Revise MO printout to enable eigenvalues with OT (05.05.2021, MK)
186 : ! **************************************************************************************************
187 762000 : SUBROUTINE qs_scf_write_mos(qs_env, scf_env, final_mos)
188 : TYPE(qs_environment_type), POINTER :: qs_env
189 : TYPE(qs_scf_env_type), POINTER :: scf_env
190 : LOGICAL, INTENT(IN) :: final_mos
191 :
192 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_write_mos'
193 :
194 : CHARACTER(LEN=2) :: solver_method
195 : CHARACTER(LEN=3*default_string_length) :: message
196 : CHARACTER(LEN=5) :: spin
197 : CHARACTER(LEN=default_string_length), &
198 190500 : DIMENSION(:), POINTER :: tmpstringlist
199 : INTEGER :: handle, homo, ikp, ispin, iw, kpoint, &
200 : nao, nelectron, nkp, nmo, nspin, numo
201 : INTEGER, DIMENSION(2) :: nmos_occ
202 190500 : INTEGER, DIMENSION(:), POINTER :: mo_index_range
203 : LOGICAL :: do_kpoints, do_printout, print_eigvals, &
204 : print_eigvecs, print_mo_info, &
205 : print_occup, print_occup_stats
206 : REAL(KIND=dp) :: flexible_electron_count, maxocc, n_el_f, &
207 : occup_stats_occ_threshold
208 190500 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, umo_eigenvalues
209 : TYPE(admm_type), POINTER :: admm_env
210 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
211 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
212 : TYPE(cp_fm_type), POINTER :: mo_coeff, umo_coeff
213 : TYPE(cp_logger_type), POINTER :: logger
214 190500 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks, s
215 : TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s, mo_coeff_deriv
216 : TYPE(dft_control_type), POINTER :: dft_control
217 : TYPE(kpoint_type), POINTER :: kpoints
218 190500 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
219 : TYPE(mo_set_type), POINTER :: mo_set, umo_set
220 : TYPE(mp_para_env_type), POINTER :: para_env
221 190500 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
222 : TYPE(preconditioner_type), POINTER :: local_preconditioner
223 190500 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
224 : TYPE(scf_control_type), POINTER :: scf_control
225 : TYPE(section_vals_type), POINTER :: dft_section, input
226 :
227 190500 : CALL timeset(routineN, handle)
228 :
229 190500 : CPASSERT(ASSOCIATED(qs_env))
230 :
231 : ! Retrieve the required information for the requested print output
232 : CALL get_qs_env(qs_env, &
233 : blacs_env=blacs_env, &
234 : dft_control=dft_control, &
235 : do_kpoints=do_kpoints, &
236 : input=input, &
237 : qs_kind_set=qs_kind_set, &
238 : para_env=para_env, &
239 : particle_set=particle_set, &
240 190500 : scf_control=scf_control)
241 :
242 : ! Quick return, if no printout of MO information is requested
243 190500 : dft_section => section_vals_get_subs_vals(input, "DFT")
244 190500 : CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVALUES", l_val=print_eigvals)
245 190500 : CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
246 190500 : CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occup)
247 190500 : CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS_STATS", c_vals=tmpstringlist)
248 :
249 190500 : print_occup_stats = .FALSE.
250 190500 : occup_stats_occ_threshold = 1e-6_dp
251 190500 : IF (SIZE(tmpstringlist) > 0) THEN ! the lone_keyword_c_vals doesn't work as advertised, handle it manually
252 190500 : print_occup_stats = .TRUE.
253 190500 : IF (LEN_TRIM(tmpstringlist(1)) > 0) &
254 190500 : READ (tmpstringlist(1), *) print_occup_stats
255 : END IF
256 190500 : IF (SIZE(tmpstringlist) > 1) &
257 190500 : READ (tmpstringlist(2), *) occup_stats_occ_threshold
258 :
259 190500 : logger => cp_get_default_logger()
260 190500 : print_mo_info = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO") /= 0)
261 :
262 190500 : IF ((.NOT. print_mo_info) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup .OR. print_occup_stats))) THEN
263 183078 : CALL timestop(handle)
264 183078 : RETURN
265 : END IF
266 :
267 7422 : NULLIFY (fm_struct_tmp)
268 7422 : NULLIFY (mo_coeff)
269 7422 : NULLIFY (mo_coeff_deriv)
270 7422 : NULLIFY (mo_eigenvalues)
271 7422 : NULLIFY (mo_set)
272 7422 : NULLIFY (umo_coeff)
273 7422 : NULLIFY (umo_eigenvalues)
274 7422 : NULLIFY (umo_set)
275 :
276 7422 : do_printout = .TRUE.
277 7422 : nspin = dft_control%nspins
278 7422 : nmos_occ = 0
279 :
280 : ! Check, if we have k points
281 7422 : IF (do_kpoints) THEN
282 22 : CALL get_qs_env(qs_env, kpoints=kpoints)
283 22 : nkp = SIZE(kpoints%kp_env)
284 : ELSE
285 7400 : CALL get_qs_env(qs_env, matrix_ks=ks, matrix_s=s)
286 7400 : CPASSERT(ASSOCIATED(ks))
287 7400 : CPASSERT(ASSOCIATED(s))
288 : nkp = 1
289 : END IF
290 :
291 12436 : kp_loop: DO ikp = 1, nkp
292 :
293 7986 : IF (do_kpoints) THEN
294 586 : mos => kpoints%kp_env(ikp)%kpoint_env%mos(1, :)
295 586 : kpoint = ikp
296 : ELSE
297 7400 : CALL get_qs_env(qs_env, matrix_ks=ks, mos=mos)
298 7400 : kpoint = 0 ! Gamma point only
299 : END IF
300 7986 : CPASSERT(ASSOCIATED(mos))
301 :
302 : ! Prepare MO information for printout
303 17744 : DO ispin = 1, nspin
304 :
305 : ! Calculate MO eigenvalues and eigenvector when OT is used
306 8280 : IF (scf_env%method == ot_method_nr) THEN
307 :
308 3194 : solver_method = "OT"
309 :
310 3194 : IF (do_kpoints) THEN
311 0 : CPABORT("The OT method is not implemented for k points")
312 : END IF
313 :
314 3194 : IF (final_mos) THEN
315 :
316 222 : matrix_ks => ks(ispin)%matrix
317 222 : matrix_s => s(1)%matrix
318 :
319 : ! With ADMM, we have to modify the Kohn-Sham matrix
320 222 : IF (dft_control%do_admm) THEN
321 0 : CALL get_qs_env(qs_env, admm_env=admm_env)
322 0 : CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks)
323 : END IF
324 :
325 222 : mo_set => mos(ispin)
326 : CALL get_mo_set(mo_set=mo_set, &
327 : mo_coeff=mo_coeff, &
328 : eigenvalues=mo_eigenvalues, &
329 : homo=homo, &
330 : maxocc=maxocc, &
331 : nelectron=nelectron, &
332 : n_el_f=n_el_f, &
333 : nao=nao, &
334 : nmo=nmo, &
335 222 : flexible_electron_count=flexible_electron_count)
336 :
337 222 : IF (ASSOCIATED(qs_env%mo_derivs)) THEN
338 222 : mo_coeff_deriv => qs_env%mo_derivs(ispin)%matrix
339 : ELSE
340 0 : mo_coeff_deriv => NULL()
341 : END IF
342 :
343 : ! Update the eigenvalues of the occupied orbitals
344 : CALL calculate_subspace_eigenvalues(orbitals=mo_coeff, &
345 : ks_matrix=matrix_ks, &
346 : evals_arg=mo_eigenvalues, &
347 222 : co_rotate_dbcsr=mo_coeff_deriv)
348 222 : CALL set_mo_occupation(mo_set=mo_set)
349 :
350 : ! Retrieve the index of the last MO for which a printout is requested
351 222 : mo_index_range => section_get_ivals(dft_section, "PRINT%MO%MO_INDEX_RANGE")
352 222 : CPASSERT(ASSOCIATED(mo_index_range))
353 222 : IF (mo_index_range(2) < 0) THEN
354 0 : numo = nao - homo
355 : ELSE
356 222 : numo = MIN(mo_index_range(2) - homo, nao - homo)
357 : END IF
358 :
359 : ! Calculate the unoccupied MO set (umo_set) with OT if needed
360 222 : IF (numo > 0) THEN
361 :
362 : ! Create temporary virtual MO set for printout
363 : CALL cp_fm_struct_create(fm_struct_tmp, &
364 : context=blacs_env, &
365 : para_env=para_env, &
366 : nrow_global=nao, &
367 20 : ncol_global=numo)
368 20 : ALLOCATE (umo_set)
369 : CALL allocate_mo_set(mo_set=umo_set, &
370 : nao=nao, &
371 : nmo=numo, &
372 : nelectron=0, &
373 : n_el_f=n_el_f, &
374 : maxocc=maxocc, &
375 20 : flexible_electron_count=flexible_electron_count)
376 : CALL init_mo_set(mo_set=umo_set, &
377 : fm_struct=fm_struct_tmp, &
378 20 : name="Temporary MO set (unoccupied MOs only) for printout")
379 20 : CALL cp_fm_struct_release(fm_struct_tmp)
380 : CALL get_mo_set(mo_set=umo_set, &
381 : mo_coeff=umo_coeff, &
382 20 : eigenvalues=umo_eigenvalues)
383 :
384 : ! Prepare printout of the additional unoccupied MOs when OT is being employed
385 20 : CALL cp_fm_init_random(umo_coeff)
386 :
387 : ! The FULL_ALL preconditioner makes not much sense for the unoccupied orbitals
388 20 : NULLIFY (local_preconditioner)
389 20 : IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
390 20 : local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
391 20 : IF (local_preconditioner%in_use == ot_precond_full_all) THEN
392 0 : NULLIFY (local_preconditioner)
393 : END IF
394 : END IF
395 :
396 : ! Calculate the MO information for the request MO index range
397 : CALL ot_eigensolver(matrix_h=matrix_ks, &
398 : matrix_s=matrix_s, &
399 : matrix_c_fm=umo_coeff, &
400 : matrix_orthogonal_space_fm=mo_coeff, &
401 : eps_gradient=scf_control%eps_lumos, &
402 : preconditioner=local_preconditioner, &
403 : iter_max=scf_control%max_iter_lumos, &
404 20 : size_ortho_space=nmo)
405 :
406 : CALL calculate_subspace_eigenvalues(orbitals=umo_coeff, &
407 : ks_matrix=matrix_ks, &
408 20 : evals_arg=umo_eigenvalues)
409 20 : CALL set_mo_occupation(mo_set=umo_set)
410 :
411 : END IF ! numo > 0
412 :
413 : ! With ADMM, we have to undo the modification of the Kohn-Sham matrix
414 222 : IF (dft_control%do_admm) THEN
415 0 : CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks)
416 : END IF
417 :
418 : ELSE
419 :
420 : message = "The MO information is only calculated after SCF convergence "// &
421 2972 : "is achieved when the orbital transformation (OT) method is used"
422 2972 : CPWARN(TRIM(message))
423 2972 : do_printout = .FALSE.
424 2972 : EXIT kp_loop
425 :
426 : END IF ! final MOs
427 :
428 : ELSE
429 :
430 5086 : solver_method = "TD"
431 5086 : mo_set => mos(ispin)
432 5086 : NULLIFY (umo_set)
433 :
434 : END IF ! OT is used
435 :
436 : ! Print MO information
437 5308 : IF (nspin > 1) THEN
438 294 : SELECT CASE (ispin)
439 : CASE (1)
440 294 : spin = "ALPHA"
441 : CASE (2)
442 294 : spin = "BETA"
443 : CASE DEFAULT
444 588 : CPABORT("Invalid spin")
445 : END SELECT
446 588 : IF (ASSOCIATED(umo_set)) THEN
447 : CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, &
448 : final_mos=final_mos, spin=TRIM(spin), solver_method=solver_method, &
449 12 : umo_set=umo_set)
450 : ELSE
451 : CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, &
452 576 : final_mos=final_mos, spin=TRIM(spin), solver_method=solver_method)
453 : END IF
454 : ELSE
455 4720 : IF (ASSOCIATED(umo_set)) THEN
456 : CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, &
457 : final_mos=final_mos, solver_method=solver_method, &
458 8 : umo_set=umo_set)
459 : ELSE
460 : CALL write_mo_set_to_output_unit(mo_set, qs_kind_set, particle_set, dft_section, 4, kpoint, &
461 4712 : final_mos=final_mos, solver_method=solver_method)
462 : END IF
463 : END IF
464 :
465 46350 : nmos_occ(ispin) = MAX(nmos_occ(ispin), COUNT(mo_set%occupation_numbers > occup_stats_occ_threshold))
466 :
467 : ! Deallocate temporary objects needed for OT
468 5308 : IF (scf_env%method == ot_method_nr) THEN
469 222 : IF (ASSOCIATED(umo_set)) THEN
470 20 : CALL deallocate_mo_set(umo_set)
471 20 : DEALLOCATE (umo_set)
472 : END IF
473 222 : NULLIFY (matrix_ks)
474 222 : NULLIFY (matrix_s)
475 : END IF
476 10322 : NULLIFY (mo_set)
477 :
478 : END DO ! ispin
479 :
480 : END DO kp_loop
481 :
482 7422 : IF (do_printout .AND. print_mo_info .AND. print_occup_stats) THEN
483 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO", &
484 : ignore_should_output=print_mo_info, &
485 0 : extension=".MOLog")
486 0 : IF (iw > 0) THEN
487 0 : IF (SIZE(mos) > 1) THEN
488 0 : WRITE (UNIT=iw, FMT="(A,I4)") " MO| Total occupied (ALPHA):", nmos_occ(1)
489 0 : WRITE (UNIT=iw, FMT="(A,I4)") " MO| Total occupied (BETA): ", nmos_occ(2)
490 : ELSE
491 0 : WRITE (UNIT=iw, FMT="(A,I4)") " MO| Total occupied: ", nmos_occ(1)
492 : END IF
493 0 : WRITE (UNIT=iw, FMT="(A)") ""
494 : END IF
495 : CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%MO", &
496 0 : ignore_should_output=print_mo_info)
497 : END IF
498 :
499 7422 : CALL timestop(handle)
500 :
501 190500 : END SUBROUTINE qs_scf_write_mos
502 :
503 : ! **************************************************************************************************
504 : !> \brief writes basic information obtained in a scf outer loop step
505 : !> \param output_unit ...
506 : !> \param scf_control ...
507 : !> \param scf_env ...
508 : !> \param energy ...
509 : !> \param total_steps ...
510 : !> \param should_stop ...
511 : !> \param outer_loop_converged ...
512 : ! **************************************************************************************************
513 5291 : SUBROUTINE qs_scf_outer_loop_info(output_unit, scf_control, scf_env, &
514 : energy, total_steps, should_stop, outer_loop_converged)
515 : INTEGER :: output_unit
516 : TYPE(scf_control_type), POINTER :: scf_control
517 : TYPE(qs_scf_env_type), POINTER :: scf_env
518 : TYPE(qs_energy_type), POINTER :: energy
519 : INTEGER :: total_steps
520 : LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged
521 :
522 : REAL(KIND=dp) :: outer_loop_eps
523 :
524 15873 : outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
525 5291 : IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') &
526 2762 : "outer SCF iter = ", scf_env%outer_scf%iter_count, &
527 5524 : " RMS gradient = ", outer_loop_eps, " energy =", energy%total
528 :
529 5291 : IF (outer_loop_converged) THEN
530 4327 : IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
531 2273 : "outer SCF loop converged in", scf_env%outer_scf%iter_count, &
532 4546 : " iterations or ", total_steps, " steps"
533 : ELSE IF (scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf &
534 964 : .OR. should_stop) THEN
535 108 : IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
536 54 : "outer SCF loop FAILED to converge after ", &
537 108 : scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps"
538 : END IF
539 :
540 5291 : END SUBROUTINE qs_scf_outer_loop_info
541 :
542 : ! **************************************************************************************************
543 : !> \brief writes basic information obtained in a scf step
544 : !> \param scf_env ...
545 : !> \param output_unit ...
546 : !> \param just_energy ...
547 : !> \param t1 ...
548 : !> \param t2 ...
549 : !> \param energy ...
550 : ! **************************************************************************************************
551 173987 : SUBROUTINE qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
552 :
553 : TYPE(qs_scf_env_type), POINTER :: scf_env
554 : INTEGER :: output_unit
555 : LOGICAL :: just_energy
556 : REAL(KIND=dp) :: t1, t2
557 : TYPE(qs_energy_type), POINTER :: energy
558 :
559 173987 : IF ((output_unit > 0) .AND. scf_env%print_iter_line) THEN
560 88252 : IF (just_energy) THEN
561 : WRITE (UNIT=output_unit, &
562 : FMT="(T2,A,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") &
563 7226 : " -", TRIM(scf_env%iter_method), scf_env%iter_param, t2 - t1, energy%total
564 : ELSE
565 78121 : IF ((ABS(scf_env%iter_delta) < 1.0E-8_dp) .OR. &
566 81026 : (ABS(scf_env%iter_delta) >= 1.0E5_dp)) THEN
567 : WRITE (UNIT=output_unit, &
568 : FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") &
569 2905 : scf_env%iter_count, TRIM(scf_env%iter_method), scf_env%iter_param, &
570 5810 : t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
571 : ELSE
572 : WRITE (UNIT=output_unit, &
573 : FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
574 78121 : scf_env%iter_count, TRIM(scf_env%iter_method), scf_env%iter_param, &
575 156242 : t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
576 : END IF
577 : END IF
578 : END IF
579 :
580 173987 : END SUBROUTINE qs_scf_loop_info
581 :
582 : ! **************************************************************************************************
583 : !> \brief writes rather detailed summary of densities and energies
584 : !> after the SCF
585 : !> \param output_unit ...
586 : !> \param rho ...
587 : !> \param qs_charges ...
588 : !> \param energy ...
589 : !> \param nelectron_total ...
590 : !> \param dft_control ...
591 : !> \param qmmm ...
592 : !> \param qs_env ...
593 : !> \param gapw ...
594 : !> \param gapw_xc ...
595 : !> \par History
596 : !> 03.2006 created [Joost VandeVondele]
597 : !> 10.2019 print dipole moment [SGh]
598 : !> 11.2022 print SCCS results [MK]
599 : ! **************************************************************************************************
600 20265 : SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
601 : dft_control, qmmm, qs_env, gapw, gapw_xc)
602 : INTEGER, INTENT(IN) :: output_unit
603 : TYPE(qs_rho_type), POINTER :: rho
604 : TYPE(qs_charges_type), POINTER :: qs_charges
605 : TYPE(qs_energy_type), POINTER :: energy
606 : INTEGER, INTENT(IN) :: nelectron_total
607 : TYPE(dft_control_type), POINTER :: dft_control
608 : LOGICAL, INTENT(IN) :: qmmm
609 : TYPE(qs_environment_type), POINTER :: qs_env
610 : LOGICAL, INTENT(IN) :: gapw, gapw_xc
611 :
612 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_print_scf_summary'
613 :
614 : INTEGER :: bc, handle, ispin, psolver
615 : REAL(kind=dp) :: e_extrapolated, exc1_energy, exc_energy, &
616 : implicit_ps_ehartree, tot1_h, tot1_s
617 20265 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
618 : TYPE(pw_env_type), POINTER :: pw_env
619 : TYPE(scf_control_type), POINTER :: scf_control
620 :
621 20265 : NULLIFY (tot_rho_r, pw_env)
622 20265 : CALL timeset(routineN, handle)
623 :
624 20265 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, scf_control=scf_control)
625 20265 : psolver = pw_env%poisson_env%parameters%solver
626 :
627 20265 : IF (output_unit > 0) THEN
628 10322 : CALL qs_rho_get(rho, tot_rho_r=tot_rho_r)
629 10322 : IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
630 : dft_control%qs_control%xtb .OR. &
631 : dft_control%qs_control%dftb)) THEN
632 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T41,2F20.10))") &
633 5716 : "Electronic density on regular grids: ", &
634 5716 : accurate_sum(tot_rho_r), &
635 5716 : accurate_sum(tot_rho_r) + nelectron_total, &
636 5716 : "Core density on regular grids:", &
637 5716 : qs_charges%total_rho_core_rspace, &
638 : qs_charges%total_rho_core_rspace + &
639 : qs_charges%total_rho1_hard_nuc - &
640 11432 : REAL(nelectron_total + dft_control%charge, dp)
641 :
642 5716 : IF (dft_control%correct_surf_dip) THEN
643 : WRITE (UNIT=output_unit, FMT="((T3,A,/,T3,A,T41,F20.10))") &
644 5 : "Total dipole moment perpendicular to ", &
645 5 : "the slab [electrons-Angstroem]: ", &
646 10 : qs_env%surface_dipole_moment
647 : END IF
648 :
649 5716 : IF (gapw) THEN
650 1013 : tot1_h = qs_charges%total_rho1_hard(1)
651 1013 : tot1_s = qs_charges%total_rho1_soft(1)
652 1219 : DO ispin = 2, dft_control%nspins
653 206 : tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
654 1219 : tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
655 : END DO
656 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
657 1013 : "Hard and soft densities (Lebedev):", &
658 2026 : tot1_h, tot1_s
659 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
660 1013 : "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
661 1013 : accurate_sum(tot_rho_r) + tot1_h - tot1_s, &
662 1013 : "Total charge density (r-space): ", &
663 : accurate_sum(tot_rho_r) + tot1_h - tot1_s &
664 : + qs_charges%total_rho_core_rspace &
665 2026 : + qs_charges%total_rho1_hard_nuc
666 1013 : IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp) THEN
667 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
668 4 : "Total CNEO nuc. char. den. (Lebedev): ", &
669 4 : qs_charges%total_rho1_hard_nuc, &
670 4 : "Total CNEO soft char. den. (Lebedev): ", &
671 4 : qs_charges%total_rho1_soft_nuc_lebedev, &
672 4 : "Total CNEO soft char. den. (r-space): ", &
673 4 : qs_charges%total_rho1_soft_nuc_rspace, &
674 4 : "Total soft Rho_e+n+0 (g-space):", &
675 8 : qs_charges%total_rho_gspace
676 : ELSE
677 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
678 1009 : "Total Rho_soft + Rho0_soft (g-space):", &
679 2018 : qs_charges%total_rho_gspace
680 : END IF
681 : ! only add total_rho1_hard_nuc for gapw as cneo requires gapw
682 : ELSE
683 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
684 4703 : "Total charge density on r-space grids: ", &
685 : accurate_sum(tot_rho_r) + &
686 4703 : qs_charges%total_rho_core_rspace, &
687 4703 : "Total charge density g-space grids: ", &
688 9406 : qs_charges%total_rho_gspace
689 : END IF
690 : END IF
691 10322 : IF (dft_control%qs_control%semi_empirical) THEN
692 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
693 1917 : "Core-core repulsion energy [eV]: ", energy%core_overlap*evolt, &
694 1917 : "Core Hamiltonian energy [eV]: ", energy%core*evolt, &
695 1917 : "Two-electron integral energy [eV]: ", energy%hartree*evolt, &
696 1917 : "Electronic energy [eV]: ", &
697 3834 : (energy%core + 0.5_dp*energy%hartree)*evolt
698 1917 : IF (energy%dispersion /= 0.0_dp) &
699 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
700 8 : "Dispersion energy [eV]: ", energy%dispersion*evolt
701 8405 : ELSEIF (dft_control%qs_control%dftb) THEN
702 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
703 681 : "Core Hamiltonian energy: ", energy%core, &
704 681 : "Repulsive potential energy: ", energy%repulsive, &
705 681 : "Electronic energy: ", energy%hartree, &
706 1362 : "Dispersion energy: ", energy%dispersion
707 681 : IF (energy%dftb3 /= 0.0_dp) &
708 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
709 135 : "DFTB3 3rd order energy: ", energy%dftb3
710 681 : IF (energy%efield /= 0.0_dp) &
711 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
712 16 : "Electric field interaction energy: ", energy%efield
713 7724 : ELSEIF (dft_control%qs_control%xtb) THEN
714 2008 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
715 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
716 447 : "Core Hamiltonian energy: ", energy%core, &
717 447 : "Repulsive potential energy: ", energy%repulsive, &
718 447 : "Electrostatic energy: ", energy%el_stat, &
719 447 : "Self-consistent dispersion energy: ", energy%dispersion_sc, &
720 447 : "Non-self consistent dispersion energy: ", energy%dispersion, &
721 894 : "Correction for halogen bonding: ", energy%xtb_xb_inter
722 : ELSE
723 1561 : IF (dft_control%qs_control%xtb_control%gfn_type == 0) THEN
724 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
725 0 : "Core Hamiltonian energy: ", energy%core, &
726 0 : "Repulsive potential energy: ", energy%repulsive, &
727 0 : "SRB Correction energy: ", energy%srb, &
728 0 : "Charge equilibration energy: ", energy%eeq, &
729 0 : "Dispersion energy: ", energy%dispersion
730 1561 : ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 1) THEN
731 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
732 1561 : "Core Hamiltonian energy: ", energy%core, &
733 1561 : "Repulsive potential energy: ", energy%repulsive, &
734 1561 : "Electronic energy: ", energy%hartree, &
735 1561 : "DFTB3 3rd order energy: ", energy%dftb3, &
736 3122 : "Dispersion energy: ", energy%dispersion
737 1561 : IF (dft_control%qs_control%xtb_control%xb_interaction) &
738 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
739 1561 : "Correction for halogen bonding: ", energy%xtb_xb_inter
740 0 : ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 2) THEN
741 0 : CPABORT("gfn_typ 2 NYA")
742 : ELSE
743 0 : CPABORT("invalid gfn_typ")
744 : END IF
745 : END IF
746 2008 : IF (dft_control%qs_control%xtb_control%do_nonbonded) &
747 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
748 12 : "Correction for nonbonded interactions: ", energy%xtb_nonbonded
749 2008 : IF (energy%efield /= 0.0_dp) &
750 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
751 394 : "Electric field interaction energy: ", energy%efield
752 : ELSE
753 5716 : IF (dft_control%do_admm) THEN
754 504 : exc_energy = energy%exc + energy%exc_aux_fit
755 504 : IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
756 : ELSE
757 5212 : exc_energy = energy%exc
758 5212 : IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
759 : END IF
760 :
761 5716 : IF (psolver == pw_poisson_implicit) THEN
762 60 : implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
763 60 : bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
764 41 : SELECT CASE (bc)
765 : CASE (MIXED_PERIODIC_BC, MIXED_BC)
766 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
767 41 : "Overlap energy of the core charge distribution:", energy%core_overlap, &
768 41 : "Self energy of the core charge distribution: ", energy%core_self, &
769 41 : "Core Hamiltonian energy: ", energy%core, &
770 41 : "Hartree energy: ", implicit_ps_ehartree, &
771 41 : "Electric enthalpy: ", energy%hartree, &
772 82 : "Exchange-correlation energy: ", exc_energy
773 : CASE (PERIODIC_BC, NEUMANN_BC)
774 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
775 19 : "Overlap energy of the core charge distribution:", energy%core_overlap, &
776 19 : "Self energy of the core charge distribution: ", energy%core_self, &
777 19 : "Core Hamiltonian energy: ", energy%core, &
778 19 : "Hartree energy: ", energy%hartree, &
779 79 : "Exchange-correlation energy: ", exc_energy
780 : END SELECT
781 : ELSE
782 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
783 5656 : "Overlap energy of the core charge distribution:", energy%core_overlap, &
784 5656 : "Self energy of the core charge distribution: ", energy%core_self, &
785 5656 : "Core Hamiltonian energy: ", energy%core, &
786 5656 : "Hartree energy: ", energy%hartree, &
787 11312 : "Exchange-correlation energy: ", exc_energy
788 : END IF
789 5716 : IF (energy%e_hartree /= 0.0_dp) &
790 : WRITE (UNIT=output_unit, FMT="(T3,A,/,T3,A,T56,F25.14)") &
791 44 : "Coulomb Electron-Electron Interaction Energy ", &
792 88 : "- Already included in the total Hartree term ", energy%e_hartree
793 5716 : IF (energy%ex /= 0.0_dp) &
794 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
795 1141 : "Hartree-Fock Exchange energy: ", energy%ex
796 5716 : IF (energy%dispersion /= 0.0_dp) &
797 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
798 160 : "Dispersion energy: ", energy%dispersion
799 5716 : IF (energy%gcp /= 0.0_dp) &
800 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
801 2 : "gCP energy: ", energy%gcp
802 5716 : IF (energy%efield /= 0.0_dp) &
803 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
804 447 : "Electric field interaction energy: ", energy%efield
805 5716 : IF (gapw) THEN
806 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
807 1013 : "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
808 2026 : "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
809 : END IF
810 5716 : IF (gapw_xc) THEN
811 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
812 187 : "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
813 : END IF
814 5716 : IF (energy%core_cneo /= 0.0_dp) THEN
815 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
816 4 : "CNEO| quantum nuclear core energy: ", energy%core_cneo
817 : END IF
818 : END IF
819 10322 : IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
820 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
821 2 : "Electronic entropic energy:", energy%kTS
822 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
823 2 : "Fermi energy:", energy%efermi
824 : END IF
825 10322 : IF (dft_control%smear) THEN
826 757 : SELECT CASE (scf_control%smear%method)
827 : CASE (smear_gaussian, smear_mp, smear_mv)
828 : ! kTS does not have physical meaning in these smearing methods
829 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
830 61 : "Smearing free energy correction:", energy%kTS
831 : CASE DEFAULT
832 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
833 696 : "Electronic entropic energy:", energy%kTS
834 : END SELECT
835 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
836 696 : "Fermi energy:", energy%efermi
837 : END IF
838 10322 : IF (dft_control%dft_plus_u) THEN
839 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
840 50 : "DFT+U energy:", energy%dft_plus_u
841 : END IF
842 10322 : IF (dft_control%do_sccs) THEN
843 6 : WRITE (UNIT=output_unit, FMT="(A)") ""
844 6 : CALL print_sccs_results(energy, dft_control%sccs_control, output_unit)
845 : END IF
846 10322 : IF (qmmm) THEN
847 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
848 1746 : "QM/MM Electrostatic energy: ", energy%qmmm_el
849 1746 : IF (qs_env%qmmm_env_qm%image_charge) THEN
850 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
851 10 : "QM/MM image charge energy: ", energy%image_charge
852 : END IF
853 : END IF
854 10322 : IF (dft_control%qs_control%mulliken_restraint) THEN
855 : WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
856 3 : "Mulliken restraint energy: ", energy%mulliken
857 : END IF
858 10322 : IF (dft_control%qs_control%semi_empirical) THEN
859 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
860 1917 : "Total energy [eV]: ", energy%total*evolt
861 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
862 1917 : "Atomic reference energy [eV]: ", energy%core_self*evolt, &
863 1917 : "Heat of formation [kcal/mol]: ", &
864 3834 : (energy%total + energy%core_self)*kcalmol
865 : ELSE
866 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
867 8405 : "Total energy: ", energy%total
868 8405 : IF (dft_control%smear) THEN
869 1324 : SELECT CASE (scf_control%smear%method)
870 : CASE (smear_fermi_dirac)
871 628 : e_extrapolated = energy%total - 0.5_dp*energy%kTS
872 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
873 628 : "Total energy (extrapolated to T->0): ", e_extrapolated
874 : CASE (smear_gaussian)
875 59 : e_extrapolated = energy%total - 0.5_dp*energy%kTS
876 : WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &
877 696 : "Total energy (extrapolated to sigma->0): ", e_extrapolated
878 : CASE (smear_mp, smear_mv)
879 : ! Sigma->0 extrapolation does not apply to MP or MV method.
880 : END SELECT
881 : END IF
882 : END IF
883 10322 : IF (qmmm) THEN
884 1746 : IF (qs_env%qmmm_env_qm%image_charge) THEN
885 10 : CALL print_image_coefficients(qs_env%image_coeff, qs_env)
886 : END IF
887 : END IF
888 10322 : CALL m_flush(output_unit)
889 : END IF
890 :
891 20265 : CALL timestop(handle)
892 :
893 20265 : END SUBROUTINE qs_scf_print_scf_summary
894 :
895 : ! **************************************************************************************************
896 : !> \brief collects the 'heavy duty' printing tasks out of the SCF loop
897 : !> \param qs_env ...
898 : !> \param scf_env ...
899 : !> \param para_env ...
900 : !> \par History
901 : !> 03.2006 created [Joost VandeVondele]
902 : ! **************************************************************************************************
903 528381 : SUBROUTINE qs_scf_loop_print(qs_env, scf_env, para_env)
904 : TYPE(qs_environment_type), POINTER :: qs_env
905 : TYPE(qs_scf_env_type), POINTER :: scf_env
906 : TYPE(mp_para_env_type), POINTER :: para_env
907 :
908 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_print'
909 :
910 : INTEGER :: after, handle, ic, ispin, iw
911 : LOGICAL :: do_kpoints, omit_headers
912 : REAL(KIND=dp) :: mo_mag_max, mo_mag_min, orthonormality
913 : TYPE(cp_logger_type), POINTER :: logger
914 176127 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
915 : TYPE(dft_control_type), POINTER :: dft_control
916 176127 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
917 : TYPE(qs_rho_type), POINTER :: rho
918 : TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
919 :
920 352254 : logger => cp_get_default_logger()
921 176127 : CALL timeset(routineN, handle)
922 :
923 : CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
924 176127 : do_kpoints=do_kpoints)
925 :
926 176127 : dft_section => section_vals_get_subs_vals(input, "DFT")
927 176127 : scf_section => section_vals_get_subs_vals(dft_section, "SCF")
928 :
929 176127 : CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
930 376775 : DO ispin = 1, dft_control%nspins
931 :
932 200648 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
933 : dft_section, "PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
934 6892 : CALL get_qs_env(qs_env, rho=rho)
935 6892 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
936 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/DENSITY", &
937 6892 : extension=".Log")
938 6892 : CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
939 6892 : after = MIN(MAX(after, 1), 16)
940 13784 : DO ic = 1, SIZE(matrix_p, 2)
941 : CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, ic)%matrix, 4, after, qs_env, para_env, &
942 13784 : output_unit=iw, omit_headers=omit_headers)
943 : END DO
944 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
945 6892 : "PRINT%AO_MATRICES/DENSITY")
946 : END IF
947 :
948 200648 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
949 176127 : dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN
950 : iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
951 5766 : extension=".Log")
952 5766 : CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
953 5766 : after = MIN(MAX(after, 1), 16)
954 5766 : CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
955 11532 : DO ic = 1, SIZE(matrix_ks, 2)
956 11532 : IF (dft_control%qs_control%semi_empirical) THEN
957 : CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
958 5762 : scale=evolt, output_unit=iw, omit_headers=omit_headers)
959 : ELSE
960 : CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
961 4 : output_unit=iw, omit_headers=omit_headers)
962 : END IF
963 : END DO
964 : CALL cp_print_key_finished_output(iw, logger, dft_section, &
965 5766 : "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
966 : END IF
967 :
968 : END DO
969 :
970 176127 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
971 : scf_section, "PRINT%MO_ORTHONORMALITY"), cp_p_file)) THEN
972 1366 : IF (do_kpoints) THEN
973 : iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
974 10 : extension=".scfLog")
975 10 : IF (iw > 0) THEN
976 : WRITE (iw, '(T8,A)') &
977 5 : " K-points: Maximum deviation from MO S-orthonormality not determined"
978 : END IF
979 : CALL cp_print_key_finished_output(iw, logger, scf_section, &
980 10 : "PRINT%MO_ORTHONORMALITY")
981 : ELSE
982 1356 : CALL get_qs_env(qs_env, mos=mos)
983 1356 : IF (scf_env%method == special_diag_method_nr) THEN
984 58 : CALL calculate_orthonormality(orthonormality, mos)
985 : ELSE
986 1298 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
987 1298 : CALL calculate_orthonormality(orthonormality, mos, matrix_s(1, 1)%matrix)
988 : END IF
989 : iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
990 1356 : extension=".scfLog")
991 1356 : IF (iw > 0) THEN
992 : WRITE (iw, '(T8,A,T61,E20.4)') &
993 678 : " Maximum deviation from MO S-orthonormality", orthonormality
994 : END IF
995 : CALL cp_print_key_finished_output(iw, logger, scf_section, &
996 1356 : "PRINT%MO_ORTHONORMALITY")
997 : END IF
998 : END IF
999 176127 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
1000 : scf_section, "PRINT%MO_MAGNITUDE"), cp_p_file)) THEN
1001 1366 : IF (do_kpoints) THEN
1002 : iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
1003 10 : extension=".scfLog")
1004 10 : IF (iw > 0) THEN
1005 : WRITE (iw, '(T8,A)') &
1006 5 : " K-points: Minimum/Maximum MO magnitude not determined"
1007 : END IF
1008 : CALL cp_print_key_finished_output(iw, logger, scf_section, &
1009 10 : "PRINT%MO_MAGNITUDE")
1010 : ELSE
1011 1356 : CALL get_qs_env(qs_env, mos=mos)
1012 1356 : CALL calculate_magnitude(mos, mo_mag_min, mo_mag_max)
1013 : iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
1014 1356 : extension=".scfLog")
1015 1356 : IF (iw > 0) THEN
1016 : WRITE (iw, '(T8,A,T41,2E20.4)') &
1017 678 : " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
1018 : END IF
1019 : CALL cp_print_key_finished_output(iw, logger, scf_section, &
1020 1356 : "PRINT%MO_MAGNITUDE")
1021 : END IF
1022 : END IF
1023 :
1024 176127 : CALL timestop(handle)
1025 :
1026 176127 : END SUBROUTINE qs_scf_loop_print
1027 :
1028 : ! **************************************************************************************************
1029 : !> \brief writes CDFT constraint information and optionally CDFT scf loop info
1030 : !> \param output_unit where to write the information
1031 : !> \param scf_control settings of the SCF loop
1032 : !> \param scf_env the env which holds convergence data
1033 : !> \param cdft_control the env which holds information about the constraint
1034 : !> \param energy the total energy
1035 : !> \param total_steps the total number of performed SCF iterations
1036 : !> \param should_stop if the calculation should stop
1037 : !> \param outer_loop_converged logical which determines if the CDFT SCF loop converged
1038 : !> \param cdft_loop logical which determines a CDFT SCF loop is active
1039 : !> \par History
1040 : !> 12.2015 created [Nico Holmberg]
1041 : ! **************************************************************************************************
1042 626 : SUBROUTINE qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
1043 : energy, total_steps, should_stop, outer_loop_converged, &
1044 : cdft_loop)
1045 : INTEGER :: output_unit
1046 : TYPE(scf_control_type), POINTER :: scf_control
1047 : TYPE(qs_scf_env_type), POINTER :: scf_env
1048 : TYPE(cdft_control_type), POINTER :: cdft_control
1049 : TYPE(qs_energy_type), POINTER :: energy
1050 : INTEGER :: total_steps
1051 : LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged, &
1052 : cdft_loop
1053 :
1054 : REAL(KIND=dp) :: outer_loop_eps
1055 :
1056 626 : IF (cdft_loop) THEN
1057 1622 : outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
1058 512 : IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') &
1059 274 : "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
1060 548 : " RMS gradient = ", outer_loop_eps, " energy =", energy%total
1061 512 : IF (outer_loop_converged) THEN
1062 270 : IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
1063 153 : "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
1064 306 : " iterations or ", total_steps, " steps"
1065 : END IF
1066 : IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
1067 512 : .AND. .NOT. outer_loop_converged) THEN
1068 56 : IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
1069 28 : "CDFT SCF loop FAILED to converge after ", &
1070 56 : scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps"
1071 : END IF
1072 : END IF
1073 626 : CALL qs_scf_cdft_constraint_info(output_unit, cdft_control)
1074 :
1075 626 : END SUBROUTINE qs_scf_cdft_info
1076 :
1077 : ! **************************************************************************************************
1078 : !> \brief writes information about the CDFT env
1079 : !> \param output_unit where to write the information
1080 : !> \param cdft_control the CDFT env that stores information about the constraint calculation
1081 : !> \par History
1082 : !> 12.2015 created [Nico Holmberg]
1083 : ! **************************************************************************************************
1084 181 : SUBROUTINE qs_scf_cdft_initial_info(output_unit, cdft_control)
1085 : INTEGER :: output_unit
1086 : TYPE(cdft_control_type), POINTER :: cdft_control
1087 :
1088 181 : IF (output_unit > 0) THEN
1089 : WRITE (output_unit, '(/,A)') &
1090 181 : " ---------------------------------- CDFT --------------------------------------"
1091 : WRITE (output_unit, '(A)') &
1092 181 : " Optimizing a density constraint in an external SCF loop "
1093 181 : WRITE (output_unit, '(A)') " "
1094 196 : SELECT CASE (cdft_control%type)
1095 : CASE (outer_scf_hirshfeld_constraint)
1096 15 : WRITE (output_unit, '(A)') " Type of constraint: Hirshfeld"
1097 : CASE (outer_scf_becke_constraint)
1098 181 : WRITE (output_unit, '(A)') " Type of constraint: Becke"
1099 : END SELECT
1100 181 : WRITE (output_unit, '(A,I8)') " Number of constraints: ", SIZE(cdft_control%group)
1101 181 : WRITE (output_unit, '(A,L8)') " Using fragment densities:", cdft_control%fragment_density
1102 181 : WRITE (output_unit, '(A)') " "
1103 181 : IF (cdft_control%atomic_charges) WRITE (output_unit, '(A,/)') " Calculating atomic CDFT charges"
1104 181 : SELECT CASE (cdft_control%constraint_control%optimizer)
1105 : CASE (outer_scf_optimizer_sd)
1106 : WRITE (output_unit, '(A)') &
1107 0 : " Minimizer : SD : steepest descent"
1108 : CASE (outer_scf_optimizer_diis)
1109 : WRITE (output_unit, '(A)') &
1110 5 : " Minimizer : DIIS : direct inversion"
1111 : WRITE (output_unit, '(A)') &
1112 5 : " in the iterative subspace"
1113 : WRITE (output_unit, '(A,I3,A)') &
1114 5 : " using ", &
1115 10 : cdft_control%constraint_control%diis_buffer_length, " DIIS vectors"
1116 : CASE (outer_scf_optimizer_bisect)
1117 : WRITE (output_unit, '(A)') &
1118 115 : " Minimizer : BISECT : gradient bisection"
1119 : WRITE (output_unit, '(A,I3)') &
1120 115 : " using a trust count of", &
1121 230 : cdft_control%constraint_control%bisect_trust_count
1122 : CASE (outer_scf_optimizer_broyden, outer_scf_optimizer_newton, &
1123 : outer_scf_optimizer_newton_ls)
1124 : CALL cdft_opt_type_write(cdft_control%constraint_control%cdft_opt_control, &
1125 60 : cdft_control%constraint_control%optimizer, output_unit)
1126 : CASE (outer_scf_optimizer_secant)
1127 1 : WRITE (output_unit, '(A)') " Minimizer : Secant"
1128 : CASE DEFAULT
1129 181 : CPABORT("")
1130 : END SELECT
1131 : WRITE (output_unit, '(/,A,L7)') &
1132 181 : " Reusing OT preconditioner: ", cdft_control%reuse_precond
1133 181 : IF (cdft_control%reuse_precond) THEN
1134 : WRITE (output_unit, '(A,I3,A,I3,A)') &
1135 0 : " using old preconditioner for up to ", &
1136 0 : cdft_control%max_reuse, " subsequent CDFT SCF"
1137 : WRITE (output_unit, '(A,I3,A,I3,A)') &
1138 0 : " iterations if the relevant loop converged in less than ", &
1139 0 : cdft_control%precond_freq, " steps"
1140 : END IF
1141 196 : SELECT CASE (cdft_control%type)
1142 : CASE (outer_scf_hirshfeld_constraint)
1143 15 : WRITE (output_unit, '(/,A)') " Hirshfeld constraint settings"
1144 15 : WRITE (output_unit, '(A)') " "
1145 194 : SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1146 : CASE (shape_function_gaussian)
1147 : WRITE (output_unit, '(A, A8)') &
1148 13 : " Shape function type: ", "Gaussian"
1149 : WRITE (output_unit, '(A)', ADVANCE='NO') &
1150 13 : " Type of Gaussian: "
1151 17 : SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1152 : CASE (radius_default)
1153 2 : WRITE (output_unit, '(A13)') "Default"
1154 : CASE (radius_covalent)
1155 11 : WRITE (output_unit, '(A13)') "Covalent"
1156 : CASE (radius_single)
1157 0 : WRITE (output_unit, '(A13)') "Fixed radius"
1158 : CASE (radius_vdw)
1159 0 : WRITE (output_unit, '(A13)') "Van der Waals"
1160 : CASE (radius_user)
1161 13 : WRITE (output_unit, '(A13)') "User-defined"
1162 :
1163 : END SELECT
1164 : CASE (shape_function_density)
1165 : WRITE (output_unit, '(A, A8)') &
1166 15 : " Shape function type: ", "Density"
1167 : END SELECT
1168 : CASE (outer_scf_becke_constraint)
1169 166 : WRITE (output_unit, '(/, A)') " Becke constraint settings"
1170 166 : WRITE (output_unit, '(A)') " "
1171 263 : SELECT CASE (cdft_control%becke_control%cutoff_type)
1172 : CASE (becke_cutoff_global)
1173 : WRITE (output_unit, '(A,F8.3,A)') &
1174 97 : " Cutoff for partitioning :", cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1175 194 : "angstrom"), " angstrom"
1176 : CASE (becke_cutoff_element)
1177 : WRITE (output_unit, '(A)') &
1178 166 : " Using element specific cutoffs for partitioning"
1179 : END SELECT
1180 : WRITE (output_unit, '(A,L7)') &
1181 166 : " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1182 : WRITE (output_unit, '(A,L7)') &
1183 166 : " Precompute gradients : ", cdft_control%becke_control%in_memory
1184 166 : WRITE (output_unit, '(A)') " "
1185 166 : IF (cdft_control%becke_control%adjust) &
1186 : WRITE (output_unit, '(A)') &
1187 110 : " Using atomic radii to generate a heteronuclear charge partitioning"
1188 166 : WRITE (output_unit, '(A)') " "
1189 347 : IF (.NOT. cdft_control%becke_control%cavity_confine) THEN
1190 : WRITE (output_unit, '(A)') &
1191 9 : " No confinement is active"
1192 : ELSE
1193 157 : WRITE (output_unit, '(A)') " Confinement using a Gaussian shaped cavity is active"
1194 158 : SELECT CASE (cdft_control%becke_control%cavity_shape)
1195 : CASE (radius_single)
1196 : WRITE (output_unit, '(A,F8.4, A)') &
1197 1 : " Type of Gaussian : Fixed radius: ", &
1198 2 : cp_unit_from_cp2k(cdft_control%becke_control%rcavity, "angstrom"), " angstrom"
1199 : CASE (radius_covalent)
1200 : WRITE (output_unit, '(A)') &
1201 1 : " Type of Gaussian : Covalent radius "
1202 : CASE (radius_vdw)
1203 : WRITE (output_unit, '(A)') &
1204 154 : " Type of Gaussian : vdW radius "
1205 : CASE (radius_user)
1206 : WRITE (output_unit, '(A)') &
1207 157 : " Type of Gaussian : User radius "
1208 : END SELECT
1209 : WRITE (output_unit, '(A,ES12.4)') &
1210 157 : " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1211 : END IF
1212 : END SELECT
1213 : WRITE (output_unit, '(/,A)') &
1214 181 : " ---------------------------------- CDFT --------------------------------------"
1215 : END IF
1216 :
1217 181 : END SUBROUTINE qs_scf_cdft_initial_info
1218 :
1219 : ! **************************************************************************************************
1220 : !> \brief writes CDFT constraint information
1221 : !> \param output_unit where to write the information
1222 : !> \param cdft_control the env which holds information about the constraint
1223 : !> \par History
1224 : !> 08.2018 separated from qs_scf_cdft_info to make code callable elsewhere [Nico Holmberg]
1225 : ! **************************************************************************************************
1226 3660 : SUBROUTINE qs_scf_cdft_constraint_info(output_unit, cdft_control)
1227 : INTEGER :: output_unit
1228 : TYPE(cdft_control_type), POINTER :: cdft_control
1229 :
1230 : INTEGER :: igroup
1231 :
1232 3660 : IF (output_unit > 0) THEN
1233 1955 : SELECT CASE (cdft_control%type)
1234 : CASE (outer_scf_hirshfeld_constraint)
1235 : WRITE (output_unit, '(/,T3,A,T60)') &
1236 61 : '------------------- Hirshfeld constraint information -------------------'
1237 : CASE (outer_scf_becke_constraint)
1238 : WRITE (output_unit, '(/,T3,A,T60)') &
1239 1833 : '--------------------- Becke constraint information ---------------------'
1240 : CASE DEFAULT
1241 1894 : CPABORT("Unknown CDFT constraint.")
1242 : END SELECT
1243 4343 : DO igroup = 1, SIZE(cdft_control%target)
1244 2449 : IF (igroup > 1) WRITE (output_unit, '(T3,A)') ' '
1245 : WRITE (output_unit, '(T3,A,T54,(3X,I18))') &
1246 2449 : 'Atomic group :', igroup
1247 3788 : SELECT CASE (cdft_control%group(igroup)%constraint_type)
1248 : CASE (cdft_charge_constraint)
1249 1339 : IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1250 : WRITE (output_unit, '(T3,A,T42,A)') &
1251 6 : 'Type of constraint :', ADJUSTR('Charge density constraint (frag.)')
1252 : ELSE
1253 : WRITE (output_unit, '(T3,A,T50,A)') &
1254 1333 : 'Type of constraint :', ADJUSTR('Charge density constraint')
1255 : END IF
1256 : CASE (cdft_magnetization_constraint)
1257 8 : IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1258 : WRITE (output_unit, '(T3,A,T35,A)') &
1259 6 : 'Type of constraint :', ADJUSTR('Magnetization density constraint (frag.)')
1260 : ELSE
1261 : WRITE (output_unit, '(T3,A,T43,A)') &
1262 2 : 'Type of constraint :', ADJUSTR('Magnetization density constraint')
1263 : END IF
1264 : CASE (cdft_alpha_constraint)
1265 551 : IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1266 : WRITE (output_unit, '(T3,A,T38,A)') &
1267 0 : 'Type of constraint :', ADJUSTR('Alpha spin density constraint (frag.)')
1268 : ELSE
1269 : WRITE (output_unit, '(T3,A,T46,A)') &
1270 551 : 'Type of constraint :', ADJUSTR('Alpha spin density constraint')
1271 : END IF
1272 : CASE (cdft_beta_constraint)
1273 551 : IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1274 : WRITE (output_unit, '(T3,A,T39,A)') &
1275 0 : 'Type of constraint :', ADJUSTR('Beta spin density constraint (frag.)')
1276 : ELSE
1277 : WRITE (output_unit, '(T3,A,T47,A)') &
1278 551 : 'Type of constraint :', ADJUSTR('Beta spin density constraint')
1279 : END IF
1280 : CASE DEFAULT
1281 2449 : CPABORT("Unknown constraint type.")
1282 : END SELECT
1283 : WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1284 2449 : 'Target value of constraint :', cdft_control%target(igroup)
1285 : WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1286 2449 : 'Current value of constraint :', cdft_control%value(igroup)
1287 : WRITE (output_unit, '(T3,A,T59,(3X,ES13.3))') &
1288 2449 : 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1289 : WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1290 4343 : 'Strength of constraint :', cdft_control%strength(igroup)
1291 : END DO
1292 : WRITE (output_unit, '(T3,A)') &
1293 1894 : '------------------------------------------------------------------------'
1294 : END IF
1295 :
1296 3660 : END SUBROUTINE qs_scf_cdft_constraint_info
1297 :
1298 : END MODULE qs_scf_output
|
