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