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