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