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