(git:419edc0)
Loading...
Searching...
No Matches
qs_scf_output.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
9 USE admm_types, ONLY: admm_type
15 USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
25 USE cp_output_handling, ONLY: cp_p_file,&
30 USE input_constants, ONLY: &
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
49 USE physcon, ONLY: evolt,&
52 USE ps_implicit_types, ONLY: mixed_bc,&
56 USE pw_env_types, ONLY: pw_env_type
71 USE qs_mo_types, ONLY: allocate_mo_set,&
77 USE qs_rho_types, ONLY: qs_rho_get,&
80 USE qs_scf_types, ONLY: ot_method_nr,&
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, &
101
102CONTAINS
103
104! **************************************************************************************************
105!> \brief writes a summary of information after scf
106!> \param output_unit ...
107!> \param qs_env ...
108! **************************************************************************************************
109 SUBROUTINE qs_scf_print_summary(output_unit, qs_env)
110 INTEGER, INTENT(IN) :: output_unit
111 TYPE(qs_environment_type), POINTER :: qs_env
112
113 INTEGER :: nelectron_total
114 LOGICAL :: gapw, gapw_xc, qmmm
115 TYPE(dft_control_type), POINTER :: dft_control
116 TYPE(qs_charges_type), POINTER :: qs_charges
117 TYPE(qs_energy_type), POINTER :: energy
118 TYPE(qs_rho_type), POINTER :: rho
119 TYPE(qs_scf_env_type), POINTER :: scf_env
120
121 NULLIFY (rho, energy, dft_control, scf_env, qs_charges)
122 CALL get_qs_env(qs_env=qs_env, rho=rho, energy=energy, dft_control=dft_control, &
123 scf_env=scf_env, qs_charges=qs_charges)
124
125 gapw = dft_control%qs_control%gapw
126 gapw_xc = dft_control%qs_control%gapw_xc
127 qmmm = qs_env%qmmm
128 nelectron_total = scf_env%nelectron
129
130 CALL qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
131 dft_control, qmmm, qs_env, gapw, gapw_xc)
132
133 END SUBROUTINE qs_scf_print_summary
134
135! **************************************************************************************************
136!> \brief writes basic information at the beginning of an scf run
137!> \param output_unit ...
138!> \param mos ...
139!> \param dft_control ...
140! **************************************************************************************************
141 SUBROUTINE qs_scf_initial_info(output_unit, mos, dft_control)
142 INTEGER :: output_unit
143 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
144 TYPE(dft_control_type), POINTER :: dft_control
145
146 INTEGER :: homo, ispin, nao, nelectron_spin, nmo
147
148 IF (output_unit > 0) THEN
149 DO ispin = 1, dft_control%nspins
150 CALL get_mo_set(mo_set=mos(ispin), &
151 homo=homo, &
152 nelectron=nelectron_spin, &
153 nao=nao, &
154 nmo=nmo)
155 IF (dft_control%nspins > 1) THEN
156 WRITE (unit=output_unit, fmt="(/,T2,A,I2)") "Spin", ispin
157 END IF
158 WRITE (unit=output_unit, fmt="(/,(T2,A,T71,I10))") &
159 "Number of electrons:", nelectron_spin, &
160 "Number of occupied orbitals:", homo, &
161 "Number of molecular orbitals:", nmo
162 END DO
163 WRITE (unit=output_unit, fmt="(/,T2,A,T71,I10)") &
164 "Number of orbital functions:", nao
165 END IF
166
167 END SUBROUTINE qs_scf_initial_info
168
169! **************************************************************************************************
170!> \brief Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit
171!> \param qs_env ...
172!> \param scf_env ...
173!> \param final_mos ...
174!> \par History
175!> - Revise MO printout to enable eigenvalues with OT (05.05.2021, MK)
176! **************************************************************************************************
177 SUBROUTINE qs_scf_write_mos(qs_env, scf_env, final_mos)
178 TYPE(qs_environment_type), POINTER :: qs_env
179 TYPE(qs_scf_env_type), POINTER :: scf_env
180 LOGICAL, INTENT(IN) :: final_mos
181
182 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_scf_write_mos'
183
184 CHARACTER(LEN=2) :: solver_method
185 CHARACTER(LEN=3*default_string_length) :: message
186 CHARACTER(LEN=5) :: spin
187 CHARACTER(LEN=default_string_length), &
188 DIMENSION(:), POINTER :: tmpstringlist
189 INTEGER :: handle, homo, ikp, ispin, iw, kpoint, &
190 nao, nelectron, nkp, nmo, nspin, numo
191 INTEGER, DIMENSION(2) :: nmos_occ
192 INTEGER, DIMENSION(:), POINTER :: mo_index_range
193 LOGICAL :: do_kpoints, print_eigvals, &
194 print_eigvecs, print_mo_info, &
195 print_occup, print_occup_stats
196 REAL(kind=dp) :: flexible_electron_count, maxocc, n_el_f, &
197 occup_stats_occ_threshold
198 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues, umo_eigenvalues
199 TYPE(admm_type), POINTER :: admm_env
200 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
201 TYPE(cp_blacs_env_type), POINTER :: blacs_env
202 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
203 TYPE(cp_fm_type), POINTER :: mo_coeff, umo_coeff
204 TYPE(cp_logger_type), POINTER :: logger
205 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks, s
206 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_s
207 TYPE(dft_control_type), POINTER :: dft_control
208 TYPE(kpoint_type), POINTER :: kpoints
209 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
210 TYPE(mo_set_type), POINTER :: mo_set, umo_set
211 TYPE(mp_para_env_type), POINTER :: para_env
212 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
213 TYPE(preconditioner_type), POINTER :: local_preconditioner
214 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
215 TYPE(scf_control_type), POINTER :: scf_control
216 TYPE(section_vals_type), POINTER :: dft_section, input
217
218 CALL timeset(routinen, handle)
219
220 cpassert(ASSOCIATED(qs_env))
221
222 ! Retrieve the required information for the requested print output
223 CALL get_qs_env(qs_env, &
224 atomic_kind_set=atomic_kind_set, &
225 blacs_env=blacs_env, &
226 dft_control=dft_control, &
227 do_kpoints=do_kpoints, &
228 input=input, &
229 qs_kind_set=qs_kind_set, &
230 para_env=para_env, &
231 particle_set=particle_set, &
232 scf_control=scf_control)
233
234 ! Quick return, if no printout of MO information is requested
235 dft_section => section_vals_get_subs_vals(input, "DFT")
236 CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVALUES", l_val=print_eigvals)
237 CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
238 CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occup)
239 CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS_STATS", c_vals=tmpstringlist)
240
241 print_occup_stats = .false.
242 occup_stats_occ_threshold = 1e-6_dp
243 IF (SIZE(tmpstringlist) > 0) THEN ! the lone_keyword_c_vals doesn't work as advertised, handle it manually
244 print_occup_stats = .true.
245 IF (len_trim(tmpstringlist(1)) > 0) &
246 READ (tmpstringlist(1), *) print_occup_stats
247 END IF
248 IF (SIZE(tmpstringlist) > 1) &
249 READ (tmpstringlist(2), *) occup_stats_occ_threshold
250
251 logger => cp_get_default_logger()
252 print_mo_info = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO") /= 0) .OR. final_mos
253
254 IF ((.NOT. print_mo_info) .OR. (.NOT. (print_eigvals .OR. print_eigvecs .OR. print_occup .OR. print_occup_stats))) THEN
255 CALL timestop(handle)
256 RETURN
257 END IF
258
259 NULLIFY (fm_struct_tmp)
260 NULLIFY (mo_coeff)
261 NULLIFY (mo_eigenvalues)
262 NULLIFY (mo_set)
263 NULLIFY (umo_coeff)
264 NULLIFY (umo_eigenvalues)
265 NULLIFY (umo_set)
266
267 nspin = dft_control%nspins
268 nmos_occ = 0
269
270 ! Check, if we have k points
271 IF (do_kpoints) THEN
272 CALL get_qs_env(qs_env, kpoints=kpoints)
273 nkp = SIZE(kpoints%kp_env)
274 ELSE
275 CALL get_qs_env(qs_env, matrix_ks=ks, matrix_s=s)
276 cpassert(ASSOCIATED(ks))
277 cpassert(ASSOCIATED(s))
278 nkp = 1
279 END IF
280
281 DO ikp = 1, nkp
282
283 IF (do_kpoints) THEN
284 mos => kpoints%kp_env(ikp)%kpoint_env%mos(1, :)
285 kpoint = ikp
286 ELSE
287 CALL get_qs_env(qs_env, matrix_ks=ks, mos=mos)
288 kpoint = 0 ! Gamma point only
289 END IF
290 cpassert(ASSOCIATED(mos))
291
292 ! Prepare MO information for printout
293 DO ispin = 1, nspin
294
295 ! Calculate MO eigenvalues and eigenvector when OT is used
296 IF (scf_env%method == ot_method_nr) THEN
297
298 solver_method = "OT"
299
300 IF (do_kpoints) THEN
301 cpabort("The OT method is not implemented for k points")
302 END IF
303
304 matrix_ks => ks(ispin)%matrix
305 matrix_s => s(1)%matrix
306
307 ! With ADMM, we have to modify the Kohn-Sham matrix
308 IF (dft_control%do_admm) THEN
309 CALL get_qs_env(qs_env, admm_env=admm_env)
310 CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks)
311 END IF
312
313 mo_set => mos(ispin)
314 CALL get_mo_set(mo_set=mo_set, &
315 mo_coeff=mo_coeff, &
316 eigenvalues=mo_eigenvalues, &
317 homo=homo, &
318 maxocc=maxocc, &
319 nelectron=nelectron, &
320 n_el_f=n_el_f, &
321 nao=nao, &
322 nmo=nmo, &
323 flexible_electron_count=flexible_electron_count)
324
325 ! Retrieve the index of the last MO for which a printout is requested
326 mo_index_range => section_get_ivals(dft_section, "PRINT%MO%MO_INDEX_RANGE")
327 cpassert(ASSOCIATED(mo_index_range))
328 numo = min(mo_index_range(2) - homo, nao - homo)
329
330 IF (.NOT. final_mos) THEN
331 numo = 0
332 message = "The MO information for unoccupied MOs is only calculated after "// &
333 "SCF convergence is achieved when the orbital transformation (OT) "// &
334 "method is used"
335 cpwarn(trim(message))
336 END IF
337
338 ! Calculate the unoccupied MO set (umo_set) with OT if needed
339 IF (numo > 0) THEN
340
341 ! Create temporary virtual MO set for printout
342 CALL cp_fm_struct_create(fm_struct_tmp, &
343 context=blacs_env, &
344 para_env=para_env, &
345 nrow_global=nao, &
346 ncol_global=numo)
347 ALLOCATE (umo_set)
348 CALL allocate_mo_set(mo_set=umo_set, &
349 nao=nao, &
350 nmo=numo, &
351 nelectron=0, &
352 n_el_f=n_el_f, &
353 maxocc=maxocc, &
354 flexible_electron_count=flexible_electron_count)
355 CALL init_mo_set(mo_set=umo_set, &
356 fm_struct=fm_struct_tmp, &
357 name="Temporary MO set (unoccupied MOs only)for printout")
358 CALL cp_fm_struct_release(fm_struct_tmp)
359 CALL get_mo_set(mo_set=umo_set, &
360 mo_coeff=umo_coeff, &
361 eigenvalues=umo_eigenvalues)
362
363 ! Calculate the MO information for the request MO index range
364 IF (final_mos) THEN
365 ! Prepare printout of the additional unoccupied MOs when OT is being employed
366 CALL cp_fm_init_random(umo_coeff)
367 ! The FULL_ALL preconditioner makes not much sense for the unoccupied orbitals
368 NULLIFY (local_preconditioner)
369 IF (ASSOCIATED(scf_env%ot_preconditioner)) THEN
370 local_preconditioner => scf_env%ot_preconditioner(1)%preconditioner
371 IF (local_preconditioner%in_use == ot_precond_full_all) THEN
372 NULLIFY (local_preconditioner)
373 END IF
374 END IF
375 CALL ot_eigensolver(matrix_h=matrix_ks, &
376 matrix_s=matrix_s, &
377 matrix_c_fm=umo_coeff, &
378 matrix_orthogonal_space_fm=mo_coeff, &
379 eps_gradient=scf_control%eps_lumos, &
380 preconditioner=local_preconditioner, &
381 iter_max=scf_control%max_iter_lumos, &
382 size_ortho_space=nmo)
383 END IF
384
385 CALL calculate_subspace_eigenvalues(orbitals=umo_coeff, &
386 ks_matrix=matrix_ks, &
387 evals_arg=umo_eigenvalues, &
388 do_rotation=.true.)
389 CALL set_mo_occupation(mo_set=umo_set)
390
391 ! With ADMM, we have to undo the modification of the Kohn-Sham matrix
392 IF (dft_control%do_admm) THEN
393 CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks)
394 END IF
395
396 ELSE
397
398 NULLIFY (umo_set)
399
400 END IF ! numo > 0
401
402 ELSE
403
404 solver_method = "TD"
405 mo_set => mos(ispin)
406 NULLIFY (umo_set)
407
408 END IF ! OT is used
409
410 ! Print MO information
411 IF (nspin > 1) THEN
412 SELECT CASE (ispin)
413 CASE (1)
414 spin = "ALPHA"
415 CASE (2)
416 spin = "BETA"
417 CASE DEFAULT
418 cpabort("Invalid spin")
419 END SELECT
420 IF (ASSOCIATED(umo_set)) THEN
421 CALL write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
422 dft_section, 4, kpoint, final_mos=final_mos, spin=trim(spin), &
423 solver_method=solver_method, umo_set=umo_set)
424 ELSE
425 CALL write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
426 dft_section, 4, kpoint, final_mos=final_mos, spin=trim(spin), &
427 solver_method=solver_method)
428 END IF
429 ELSE
430 IF (ASSOCIATED(umo_set)) THEN
431 CALL write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
432 dft_section, 4, kpoint, final_mos=final_mos, &
433 solver_method=solver_method, umo_set=umo_set)
434 ELSE
435 CALL write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, &
436 dft_section, 4, kpoint, final_mos=final_mos, &
437 solver_method=solver_method)
438 END IF
439 END IF
440
441 nmos_occ(ispin) = max(nmos_occ(ispin), count(mo_set%occupation_numbers > occup_stats_occ_threshold))
442
443 ! Deallocate temporary objects needed for OT
444 IF (scf_env%method == ot_method_nr) THEN
445 IF (ASSOCIATED(umo_set)) THEN
446 CALL deallocate_mo_set(umo_set)
447 DEALLOCATE (umo_set)
448 END IF
449 NULLIFY (matrix_ks)
450 NULLIFY (matrix_s)
451 END IF
452 NULLIFY (mo_set)
453
454 END DO ! ispin
455
456 END DO ! k point loop
457
458 IF (print_mo_info .AND. print_occup_stats) THEN
459 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO", &
460 ignore_should_output=print_mo_info, &
461 extension=".MOLog")
462 IF (iw > 0) THEN
463 IF (SIZE(mos) > 1) THEN
464 WRITE (unit=iw, fmt="(A,I4)") " MO| Total occupied (ALPHA):", nmos_occ(1)
465 WRITE (unit=iw, fmt="(A,I4)") " MO| Total occupied (BETA): ", nmos_occ(2)
466 ELSE
467 WRITE (unit=iw, fmt="(A,I4)") " MO| Total occupied: ", nmos_occ(1)
468 END IF
469 WRITE (unit=iw, fmt="(A)") ""
470 END IF
471 CALL cp_print_key_finished_output(iw, logger, dft_section, "PRINT%MO", &
472 ignore_should_output=print_mo_info)
473 END IF
474
475 CALL timestop(handle)
476
477 END SUBROUTINE qs_scf_write_mos
478
479! **************************************************************************************************
480!> \brief writes basic information obtained in a scf outer loop step
481!> \param output_unit ...
482!> \param scf_control ...
483!> \param scf_env ...
484!> \param energy ...
485!> \param total_steps ...
486!> \param should_stop ...
487!> \param outer_loop_converged ...
488! **************************************************************************************************
489 SUBROUTINE qs_scf_outer_loop_info(output_unit, scf_control, scf_env, &
490 energy, total_steps, should_stop, outer_loop_converged)
491 INTEGER :: output_unit
492 TYPE(scf_control_type), POINTER :: scf_control
493 TYPE(qs_scf_env_type), POINTER :: scf_env
494 TYPE(qs_energy_type), POINTER :: energy
495 INTEGER :: total_steps
496 LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged
497
498 REAL(kind=dp) :: outer_loop_eps
499
500 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
501 IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') &
502 "outer SCF iter = ", scf_env%outer_scf%iter_count, &
503 " RMS gradient = ", outer_loop_eps, " energy =", energy%total
504
505 IF (outer_loop_converged) THEN
506 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
507 "outer SCF loop converged in", scf_env%outer_scf%iter_count, &
508 " iterations or ", total_steps, " steps"
509 ELSE IF (scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf &
510 .OR. should_stop) THEN
511 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
512 "outer SCF loop FAILED to converge after ", &
513 scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps"
514 END IF
515
516 END SUBROUTINE qs_scf_outer_loop_info
517
518! **************************************************************************************************
519!> \brief writes basic information obtained in a scf step
520!> \param scf_env ...
521!> \param output_unit ...
522!> \param just_energy ...
523!> \param t1 ...
524!> \param t2 ...
525!> \param energy ...
526! **************************************************************************************************
527 SUBROUTINE qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
528
529 TYPE(qs_scf_env_type), POINTER :: scf_env
530 INTEGER :: output_unit
531 LOGICAL :: just_energy
532 REAL(kind=dp) :: t1, t2
533 TYPE(qs_energy_type), POINTER :: energy
534
535 IF ((output_unit > 0) .AND. scf_env%print_iter_line) THEN
536 IF (just_energy) THEN
537 WRITE (unit=output_unit, &
538 fmt="(T2,A,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") &
539 " -", trim(scf_env%iter_method), scf_env%iter_param, t2 - t1, energy%total
540 ELSE
541 IF ((abs(scf_env%iter_delta) < 1.0e-8_dp) .OR. &
542 (abs(scf_env%iter_delta) >= 1.0e5_dp)) THEN
543 WRITE (unit=output_unit, &
544 fmt="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") &
545 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
546 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
547 ELSE
548 WRITE (unit=output_unit, &
549 fmt="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
550 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
551 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
552 END IF
553 END IF
554 END IF
555
556 END SUBROUTINE qs_scf_loop_info
557
558! **************************************************************************************************
559!> \brief writes rather detailed summary of densities and energies
560!> after the SCF
561!> \param output_unit ...
562!> \param rho ...
563!> \param qs_charges ...
564!> \param energy ...
565!> \param nelectron_total ...
566!> \param dft_control ...
567!> \param qmmm ...
568!> \param qs_env ...
569!> \param gapw ...
570!> \param gapw_xc ...
571!> \par History
572!> 03.2006 created [Joost VandeVondele]
573!> 10.2019 print dipole moment [SGh]
574!> 11.2022 print SCCS results [MK]
575! **************************************************************************************************
576 SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
577 dft_control, qmmm, qs_env, gapw, gapw_xc)
578 INTEGER, INTENT(IN) :: output_unit
579 TYPE(qs_rho_type), POINTER :: rho
580 TYPE(qs_charges_type), POINTER :: qs_charges
581 TYPE(qs_energy_type), POINTER :: energy
582 INTEGER, INTENT(IN) :: nelectron_total
583 TYPE(dft_control_type), POINTER :: dft_control
584 LOGICAL, INTENT(IN) :: qmmm
585 TYPE(qs_environment_type), POINTER :: qs_env
586 LOGICAL, INTENT(IN) :: gapw, gapw_xc
587
588 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_scf_print_scf_summary'
589
590 INTEGER :: bc, handle, ispin, psolver
591 REAL(kind=dp) :: exc1_energy, exc_energy, &
592 implicit_ps_ehartree, tot1_h, tot1_s
593 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
594 TYPE(pw_env_type), POINTER :: pw_env
595
596 NULLIFY (tot_rho_r, pw_env)
597 CALL timeset(routinen, handle)
598
599 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
600 psolver = pw_env%poisson_env%parameters%solver
601
602 IF (output_unit > 0) THEN
603 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r)
604 IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
605 dft_control%qs_control%xtb .OR. &
606 dft_control%qs_control%dftb)) THEN
607 WRITE (unit=output_unit, fmt="(/,(T3,A,T41,2F20.10))") &
608 "Electronic density on regular grids: ", &
609 accurate_sum(tot_rho_r), &
610 accurate_sum(tot_rho_r) + nelectron_total, &
611 "Core density on regular grids:", &
612 qs_charges%total_rho_core_rspace, &
613 qs_charges%total_rho_core_rspace - real(nelectron_total + dft_control%charge, dp)
614
615 IF (dft_control%correct_surf_dip) THEN
616 WRITE (unit=output_unit, fmt="((T3,A,/,T3,A,T41,F20.10))") &
617 "Total dipole moment perpendicular to ", &
618 "the slab [electrons-Angstroem]: ", &
619 qs_env%surface_dipole_moment
620 END IF
621
622 IF (gapw) THEN
623 tot1_h = qs_charges%total_rho1_hard(1)
624 tot1_s = qs_charges%total_rho1_soft(1)
625 DO ispin = 2, dft_control%nspins
626 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
627 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
628 END DO
629 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
630 "Hard and soft densities (Lebedev):", &
631 tot1_h, tot1_s
632 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
633 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
634 accurate_sum(tot_rho_r) + tot1_h - tot1_s, &
635 "Total charge density (r-space): ", &
636 accurate_sum(tot_rho_r) + tot1_h - tot1_s &
637 + qs_charges%total_rho_core_rspace, &
638 "Total Rho_soft + Rho0_soft (g-space):", &
639 qs_charges%total_rho_gspace
640 ELSE
641 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
642 "Total charge density on r-space grids: ", &
643 accurate_sum(tot_rho_r) + &
644 qs_charges%total_rho_core_rspace, &
645 "Total charge density g-space grids: ", &
646 qs_charges%total_rho_gspace
647 END IF
648 END IF
649 IF (dft_control%qs_control%semi_empirical) THEN
650 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
651 "Core-core repulsion energy [eV]: ", energy%core_overlap*evolt, &
652 "Core Hamiltonian energy [eV]: ", energy%core*evolt, &
653 "Two-electron integral energy [eV]: ", energy%hartree*evolt, &
654 "Electronic energy [eV]: ", &
655 (energy%core + 0.5_dp*energy%hartree)*evolt
656 IF (energy%dispersion /= 0.0_dp) &
657 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
658 "Dispersion energy [eV]: ", energy%dispersion*evolt
659 ELSEIF (dft_control%qs_control%dftb) THEN
660 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
661 "Core Hamiltonian energy: ", energy%core, &
662 "Repulsive potential energy: ", energy%repulsive, &
663 "Electronic energy: ", energy%hartree, &
664 "Dispersion energy: ", energy%dispersion
665 IF (energy%dftb3 /= 0.0_dp) &
666 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
667 "DFTB3 3rd order energy: ", energy%dftb3
668 IF (energy%efield /= 0.0_dp) &
669 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
670 "Electric field interaction energy: ", energy%efield
671 ELSEIF (dft_control%qs_control%xtb) THEN
672 IF (dft_control%qs_control%xtb_control%gfn_type == 0) THEN
673 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
674 "Core Hamiltonian energy: ", energy%core, &
675 "Repulsive potential energy: ", energy%repulsive, &
676 "SRB Correction energy: ", energy%srb, &
677 "Charge equilibration energy: ", energy%eeq, &
678 "Dispersion energy: ", energy%dispersion
679 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 1) THEN
680 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
681 "Core Hamiltonian energy: ", energy%core, &
682 "Repulsive potential energy: ", energy%repulsive, &
683 "Electronic energy: ", energy%hartree, &
684 "DFTB3 3rd order energy: ", energy%dftb3, &
685 "Dispersion energy: ", energy%dispersion
686 IF (dft_control%qs_control%xtb_control%xb_interaction) &
687 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
688 "Correction for halogen bonding: ", energy%xtb_xb_inter
689 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 2) THEN
690 cpabort("gfn_typ 2 NYA")
691 ELSE
692 cpabort("invalid gfn_typ")
693 END IF
694 IF (dft_control%qs_control%xtb_control%do_nonbonded) &
695 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
696 "Correction for nonbonded interactions: ", energy%xtb_nonbonded
697 IF (energy%efield /= 0.0_dp) &
698 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
699 "Electric field interaction energy: ", energy%efield
700 ELSE
701 IF (dft_control%do_admm) THEN
702 exc_energy = energy%exc + energy%exc_aux_fit
703 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
704 ELSE
705 exc_energy = energy%exc
706 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
707 END IF
708
709 IF (psolver .EQ. pw_poisson_implicit) THEN
710 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
711 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
712 SELECT CASE (bc)
714 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
715 "Overlap energy of the core charge distribution:", energy%core_overlap, &
716 "Self energy of the core charge distribution: ", energy%core_self, &
717 "Core Hamiltonian energy: ", energy%core, &
718 "Hartree energy: ", implicit_ps_ehartree, &
719 "Electric enthalpy: ", energy%hartree, &
720 "Exchange-correlation energy: ", exc_energy
721 CASE (periodic_bc, neumann_bc)
722 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
723 "Overlap energy of the core charge distribution:", energy%core_overlap, &
724 "Self energy of the core charge distribution: ", energy%core_self, &
725 "Core Hamiltonian energy: ", energy%core, &
726 "Hartree energy: ", energy%hartree, &
727 "Exchange-correlation energy: ", exc_energy
728 END SELECT
729 ELSE
730 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
731 "Overlap energy of the core charge distribution:", energy%core_overlap, &
732 "Self energy of the core charge distribution: ", energy%core_self, &
733 "Core Hamiltonian energy: ", energy%core, &
734 "Hartree energy: ", energy%hartree, &
735 "Exchange-correlation energy: ", exc_energy
736 END IF
737 IF (energy%e_hartree /= 0.0_dp) &
738 WRITE (unit=output_unit, fmt="(T3,A,/,T3,A,T56,F25.14)") &
739 "Coulomb Electron-Electron Interaction Energy ", &
740 "- Already included in the total Hartree term ", energy%e_hartree
741 IF (energy%ex /= 0.0_dp) &
742 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
743 "Hartree-Fock Exchange energy: ", energy%ex
744 IF (energy%dispersion /= 0.0_dp) &
745 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
746 "Dispersion energy: ", energy%dispersion
747 IF (energy%gcp /= 0.0_dp) &
748 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
749 "gCP energy: ", energy%gcp
750 IF (energy%efield /= 0.0_dp) &
751 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
752 "Electric field interaction energy: ", energy%efield
753 IF (gapw) THEN
754 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
755 "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
756 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
757 END IF
758 IF (gapw_xc) THEN
759 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
760 "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
761 END IF
762 END IF
763 IF (dft_control%smear) THEN
764 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
765 "Electronic entropic energy:", energy%kTS
766 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
767 "Fermi energy:", energy%efermi
768 END IF
769 IF (dft_control%dft_plus_u) THEN
770 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
771 "DFT+U energy:", energy%dft_plus_u
772 END IF
773 IF (dft_control%do_sccs) THEN
774 WRITE (unit=output_unit, fmt="(A)") ""
775 CALL print_sccs_results(energy, dft_control%sccs_control, output_unit)
776 END IF
777 IF (qmmm) THEN
778 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
779 "QM/MM Electrostatic energy: ", energy%qmmm_el
780 IF (qs_env%qmmm_env_qm%image_charge) THEN
781 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
782 "QM/MM image charge energy: ", energy%image_charge
783 END IF
784 END IF
785 IF (dft_control%qs_control%mulliken_restraint) THEN
786 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
787 "Mulliken restraint energy: ", energy%mulliken
788 END IF
789 IF (dft_control%qs_control%semi_empirical) THEN
790 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
791 "Total energy [eV]: ", energy%total*evolt
792 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
793 "Atomic reference energy [eV]: ", energy%core_self*evolt, &
794 "Heat of formation [kcal/mol]: ", &
795 (energy%total + energy%core_self)*kcalmol
796 ELSE
797 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
798 "Total energy: ", energy%total
799 END IF
800 IF (qmmm) THEN
801 IF (qs_env%qmmm_env_qm%image_charge) THEN
802 CALL print_image_coefficients(qs_env%image_coeff, qs_env)
803 END IF
804 END IF
805 CALL m_flush(output_unit)
806 END IF
807
808 CALL timestop(handle)
809
810 END SUBROUTINE qs_scf_print_scf_summary
811
812! **************************************************************************************************
813!> \brief collects the 'heavy duty' printing tasks out of the SCF loop
814!> \param qs_env ...
815!> \param scf_env ...
816!> \param para_env ...
817!> \par History
818!> 03.2006 created [Joost VandeVondele]
819! **************************************************************************************************
820 SUBROUTINE qs_scf_loop_print(qs_env, scf_env, para_env)
821 TYPE(qs_environment_type), POINTER :: qs_env
822 TYPE(qs_scf_env_type), POINTER :: scf_env
823 TYPE(mp_para_env_type), POINTER :: para_env
824
825 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_scf_loop_print'
826
827 INTEGER :: after, handle, ic, ispin, iw
828 LOGICAL :: do_kpoints, omit_headers
829 REAL(kind=dp) :: mo_mag_max, mo_mag_min, orthonormality
830 TYPE(cp_logger_type), POINTER :: logger
831 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
832 TYPE(dft_control_type), POINTER :: dft_control
833 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
834 TYPE(qs_rho_type), POINTER :: rho
835 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
836
837 logger => cp_get_default_logger()
838 CALL timeset(routinen, handle)
839
840 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
841 do_kpoints=do_kpoints)
842
843 dft_section => section_vals_get_subs_vals(input, "DFT")
844 scf_section => section_vals_get_subs_vals(dft_section, "SCF")
845
846 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
847 DO ispin = 1, dft_control%nspins
848
849 IF (btest(cp_print_key_should_output(logger%iter_info, &
850 dft_section, "PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
851 CALL get_qs_env(qs_env, rho=rho)
852 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
853 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/DENSITY", &
854 extension=".Log")
855 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
856 after = min(max(after, 1), 16)
857 DO ic = 1, SIZE(matrix_p, 2)
858 CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, ic)%matrix, 4, after, qs_env, para_env, &
859 output_unit=iw, omit_headers=omit_headers)
860 END DO
861 CALL cp_print_key_finished_output(iw, logger, dft_section, &
862 "PRINT%AO_MATRICES/DENSITY")
863 END IF
864
865 IF (btest(cp_print_key_should_output(logger%iter_info, &
866 dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN
867 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
868 extension=".Log")
869 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
870 after = min(max(after, 1), 16)
871 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
872 DO ic = 1, SIZE(matrix_ks, 2)
873 IF (dft_control%qs_control%semi_empirical) THEN
874 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
875 scale=evolt, output_unit=iw, omit_headers=omit_headers)
876 ELSE
877 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
878 output_unit=iw, omit_headers=omit_headers)
879 END IF
880 END DO
881 CALL cp_print_key_finished_output(iw, logger, dft_section, &
882 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
883 END IF
884
885 END DO
886
887 IF (btest(cp_print_key_should_output(logger%iter_info, &
888 scf_section, "PRINT%MO_ORTHONORMALITY"), cp_p_file)) THEN
889 IF (do_kpoints) THEN
890 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
891 extension=".scfLog")
892 IF (iw > 0) THEN
893 WRITE (iw, '(T8,A)') &
894 " K-points: Maximum deviation from MO S-orthonormality not determined"
895 END IF
896 CALL cp_print_key_finished_output(iw, logger, scf_section, &
897 "PRINT%MO_ORTHONORMALITY")
898 ELSE
899 CALL get_qs_env(qs_env, mos=mos)
900 IF (scf_env%method == special_diag_method_nr) THEN
901 CALL calculate_orthonormality(orthonormality, mos)
902 ELSE
903 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
904 CALL calculate_orthonormality(orthonormality, mos, matrix_s(1, 1)%matrix)
905 END IF
906 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
907 extension=".scfLog")
908 IF (iw > 0) THEN
909 WRITE (iw, '(T8,A,T61,E20.4)') &
910 " Maximum deviation from MO S-orthonormality", orthonormality
911 END IF
912 CALL cp_print_key_finished_output(iw, logger, scf_section, &
913 "PRINT%MO_ORTHONORMALITY")
914 END IF
915 END IF
916 IF (btest(cp_print_key_should_output(logger%iter_info, &
917 scf_section, "PRINT%MO_MAGNITUDE"), cp_p_file)) THEN
918 IF (do_kpoints) THEN
919 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
920 extension=".scfLog")
921 IF (iw > 0) THEN
922 WRITE (iw, '(T8,A)') &
923 " K-points: Minimum/Maximum MO magnitude not determined"
924 END IF
925 CALL cp_print_key_finished_output(iw, logger, scf_section, &
926 "PRINT%MO_MAGNITUDE")
927 ELSE
928 CALL get_qs_env(qs_env, mos=mos)
929 CALL calculate_magnitude(mos, mo_mag_min, mo_mag_max)
930 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
931 extension=".scfLog")
932 IF (iw > 0) THEN
933 WRITE (iw, '(T8,A,T41,2E20.4)') &
934 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
935 END IF
936 CALL cp_print_key_finished_output(iw, logger, scf_section, &
937 "PRINT%MO_MAGNITUDE")
938 END IF
939 END IF
940
941 CALL timestop(handle)
942
943 END SUBROUTINE qs_scf_loop_print
944
945! **************************************************************************************************
946!> \brief writes CDFT constraint information and optionally CDFT scf loop info
947!> \param output_unit where to write the information
948!> \param scf_control settings of the SCF loop
949!> \param scf_env the env which holds convergence data
950!> \param cdft_control the env which holds information about the constraint
951!> \param energy the total energy
952!> \param total_steps the total number of performed SCF iterations
953!> \param should_stop if the calculation should stop
954!> \param outer_loop_converged logical which determines if the CDFT SCF loop converged
955!> \param cdft_loop logical which determines a CDFT SCF loop is active
956!> \par History
957!> 12.2015 created [Nico Holmberg]
958! **************************************************************************************************
959 SUBROUTINE qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
960 energy, total_steps, should_stop, outer_loop_converged, &
961 cdft_loop)
962 INTEGER :: output_unit
963 TYPE(scf_control_type), POINTER :: scf_control
964 TYPE(qs_scf_env_type), POINTER :: scf_env
965 TYPE(cdft_control_type), POINTER :: cdft_control
966 TYPE(qs_energy_type), POINTER :: energy
967 INTEGER :: total_steps
968 LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged, &
969 cdft_loop
970
971 REAL(kind=dp) :: outer_loop_eps
972
973 IF (cdft_loop) THEN
974 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
975 IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') &
976 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
977 " RMS gradient = ", outer_loop_eps, " energy =", energy%total
978 IF (outer_loop_converged) THEN
979 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
980 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
981 " iterations or ", total_steps, " steps"
982 END IF
983 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
984 .AND. .NOT. outer_loop_converged) THEN
985 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
986 "CDFT SCF loop FAILED to converge after ", &
987 scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps"
988 END IF
989 END IF
990 CALL qs_scf_cdft_constraint_info(output_unit, cdft_control)
991
992 END SUBROUTINE qs_scf_cdft_info
993
994! **************************************************************************************************
995!> \brief writes information about the CDFT env
996!> \param output_unit where to write the information
997!> \param cdft_control the CDFT env that stores information about the constraint calculation
998!> \par History
999!> 12.2015 created [Nico Holmberg]
1000! **************************************************************************************************
1001 SUBROUTINE qs_scf_cdft_initial_info(output_unit, cdft_control)
1002 INTEGER :: output_unit
1003 TYPE(cdft_control_type), POINTER :: cdft_control
1004
1005 IF (output_unit > 0) THEN
1006 WRITE (output_unit, '(/,A)') &
1007 " ---------------------------------- CDFT --------------------------------------"
1008 WRITE (output_unit, '(A)') &
1009 " Optimizing a density constraint in an external SCF loop "
1010 WRITE (output_unit, '(A)') " "
1011 SELECT CASE (cdft_control%type)
1013 WRITE (output_unit, '(A)') " Type of constraint: Hirshfeld"
1015 WRITE (output_unit, '(A)') " Type of constraint: Becke"
1016 END SELECT
1017 WRITE (output_unit, '(A,I8)') " Number of constraints: ", SIZE(cdft_control%group)
1018 WRITE (output_unit, '(A,L8)') " Using fragment densities:", cdft_control%fragment_density
1019 WRITE (output_unit, '(A)') " "
1020 IF (cdft_control%atomic_charges) WRITE (output_unit, '(A,/)') " Calculating atomic CDFT charges"
1021 SELECT CASE (cdft_control%constraint_control%optimizer)
1023 WRITE (output_unit, '(A)') &
1024 " Minimizer : SD : steepest descent"
1026 WRITE (output_unit, '(A)') &
1027 " Minimizer : DIIS : direct inversion"
1028 WRITE (output_unit, '(A)') &
1029 " in the iterative subspace"
1030 WRITE (output_unit, '(A,I3,A)') &
1031 " using ", &
1032 cdft_control%constraint_control%diis_buffer_length, " DIIS vectors"
1034 WRITE (output_unit, '(A)') &
1035 " Minimizer : BISECT : gradient bisection"
1036 WRITE (output_unit, '(A,I3)') &
1037 " using a trust count of", &
1038 cdft_control%constraint_control%bisect_trust_count
1041 CALL cdft_opt_type_write(cdft_control%constraint_control%cdft_opt_control, &
1042 cdft_control%constraint_control%optimizer, output_unit)
1044 WRITE (output_unit, '(A)') " Minimizer : Secant"
1045 CASE DEFAULT
1046 cpabort("")
1047 END SELECT
1048 WRITE (output_unit, '(/,A,L7)') &
1049 " Reusing OT preconditioner: ", cdft_control%reuse_precond
1050 IF (cdft_control%reuse_precond) THEN
1051 WRITE (output_unit, '(A,I3,A,I3,A)') &
1052 " using old preconditioner for up to ", &
1053 cdft_control%max_reuse, " subsequent CDFT SCF"
1054 WRITE (output_unit, '(A,I3,A,I3,A)') &
1055 " iterations if the relevant loop converged in less than ", &
1056 cdft_control%precond_freq, " steps"
1057 END IF
1058 SELECT CASE (cdft_control%type)
1060 WRITE (output_unit, '(/,A)') " Hirshfeld constraint settings"
1061 WRITE (output_unit, '(A)') " "
1062 SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1064 WRITE (output_unit, '(A, A8)') &
1065 " Shape function type: ", "Gaussian"
1066 WRITE (output_unit, '(A)', advance='NO') &
1067 " Type of Gaussian: "
1068 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1069 CASE (radius_default)
1070 WRITE (output_unit, '(A13)') "Default"
1071 CASE (radius_covalent)
1072 WRITE (output_unit, '(A13)') "Covalent"
1073 CASE (radius_single)
1074 WRITE (output_unit, '(A13)') "Fixed radius"
1075 CASE (radius_vdw)
1076 WRITE (output_unit, '(A13)') "Van der Waals"
1077 CASE (radius_user)
1078 WRITE (output_unit, '(A13)') "User-defined"
1079
1080 END SELECT
1082 WRITE (output_unit, '(A, A8)') &
1083 " Shape function type: ", "Density"
1084 END SELECT
1086 WRITE (output_unit, '(/, A)') " Becke constraint settings"
1087 WRITE (output_unit, '(A)') " "
1088 SELECT CASE (cdft_control%becke_control%cutoff_type)
1089 CASE (becke_cutoff_global)
1090 WRITE (output_unit, '(A,F8.3,A)') &
1091 " Cutoff for partitioning :", cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1092 "angstrom"), " angstrom"
1094 WRITE (output_unit, '(A)') &
1095 " Using element specific cutoffs for partitioning"
1096 END SELECT
1097 WRITE (output_unit, '(A,L7)') &
1098 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1099 WRITE (output_unit, '(A,L7)') &
1100 " Precompute gradients : ", cdft_control%becke_control%in_memory
1101 WRITE (output_unit, '(A)') " "
1102 IF (cdft_control%becke_control%adjust) &
1103 WRITE (output_unit, '(A)') &
1104 " Using atomic radii to generate a heteronuclear charge partitioning"
1105 WRITE (output_unit, '(A)') " "
1106 IF (.NOT. cdft_control%becke_control%cavity_confine) THEN
1107 WRITE (output_unit, '(A)') &
1108 " No confinement is active"
1109 ELSE
1110 WRITE (output_unit, '(A)') " Confinement using a Gaussian shaped cavity is active"
1111 SELECT CASE (cdft_control%becke_control%cavity_shape)
1112 CASE (radius_single)
1113 WRITE (output_unit, '(A,F8.4, A)') &
1114 " Type of Gaussian : Fixed radius: ", &
1115 cp_unit_from_cp2k(cdft_control%becke_control%rcavity, "angstrom"), " angstrom"
1116 CASE (radius_covalent)
1117 WRITE (output_unit, '(A)') &
1118 " Type of Gaussian : Covalent radius "
1119 CASE (radius_vdw)
1120 WRITE (output_unit, '(A)') &
1121 " Type of Gaussian : vdW radius "
1122 CASE (radius_user)
1123 WRITE (output_unit, '(A)') &
1124 " Type of Gaussian : User radius "
1125 END SELECT
1126 WRITE (output_unit, '(A,ES12.4)') &
1127 " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1128 END IF
1129 END SELECT
1130 WRITE (output_unit, '(/,A)') &
1131 " ---------------------------------- CDFT --------------------------------------"
1132 END IF
1133
1134 END SUBROUTINE qs_scf_cdft_initial_info
1135
1136! **************************************************************************************************
1137!> \brief writes CDFT constraint information
1138!> \param output_unit where to write the information
1139!> \param cdft_control the env which holds information about the constraint
1140!> \par History
1141!> 08.2018 separated from qs_scf_cdft_info to make code callable elsewhere [Nico Holmberg]
1142! **************************************************************************************************
1143 SUBROUTINE qs_scf_cdft_constraint_info(output_unit, cdft_control)
1144 INTEGER :: output_unit
1145 TYPE(cdft_control_type), POINTER :: cdft_control
1146
1147 INTEGER :: igroup
1148
1149 IF (output_unit > 0) THEN
1150 SELECT CASE (cdft_control%type)
1152 WRITE (output_unit, '(/,T3,A,T60)') &
1153 '------------------- Hirshfeld constraint information -------------------'
1155 WRITE (output_unit, '(/,T3,A,T60)') &
1156 '--------------------- Becke constraint information ---------------------'
1157 CASE DEFAULT
1158 cpabort("Unknown CDFT constraint.")
1159 END SELECT
1160 DO igroup = 1, SIZE(cdft_control%target)
1161 IF (igroup > 1) WRITE (output_unit, '(T3,A)') ' '
1162 WRITE (output_unit, '(T3,A,T54,(3X,I18))') &
1163 'Atomic group :', igroup
1164 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1166 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1167 WRITE (output_unit, '(T3,A,T42,A)') &
1168 'Type of constraint :', adjustr('Charge density constraint (frag.)')
1169 ELSE
1170 WRITE (output_unit, '(T3,A,T50,A)') &
1171 'Type of constraint :', adjustr('Charge density constraint')
1172 END IF
1174 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1175 WRITE (output_unit, '(T3,A,T35,A)') &
1176 'Type of constraint :', adjustr('Magnetization density constraint (frag.)')
1177 ELSE
1178 WRITE (output_unit, '(T3,A,T43,A)') &
1179 'Type of constraint :', adjustr('Magnetization density constraint')
1180 END IF
1182 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1183 WRITE (output_unit, '(T3,A,T38,A)') &
1184 'Type of constraint :', adjustr('Alpha spin density constraint (frag.)')
1185 ELSE
1186 WRITE (output_unit, '(T3,A,T46,A)') &
1187 'Type of constraint :', adjustr('Alpha spin density constraint')
1188 END IF
1190 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1191 WRITE (output_unit, '(T3,A,T39,A)') &
1192 'Type of constraint :', adjustr('Beta spin density constraint (frag.)')
1193 ELSE
1194 WRITE (output_unit, '(T3,A,T47,A)') &
1195 'Type of constraint :', adjustr('Beta spin density constraint')
1196 END IF
1197 CASE DEFAULT
1198 cpabort("Unknown constraint type.")
1199 END SELECT
1200 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1201 'Target value of constraint :', cdft_control%target(igroup)
1202 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1203 'Current value of constraint :', cdft_control%value(igroup)
1204 WRITE (output_unit, '(T3,A,T59,(3X,ES13.3))') &
1205 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1206 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1207 'Strength of constraint :', cdft_control%strength(igroup)
1208 END DO
1209 WRITE (output_unit, '(T3,A)') &
1210 '------------------------------------------------------------------------'
1211 END IF
1212
1213 END SUBROUTINE qs_scf_cdft_constraint_info
1214
1215END MODULE qs_scf_output
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
Contains methods used in the context of density fitting.
Definition admm_utils.F:15
subroutine, public admm_uncorrect_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:127
subroutine, public admm_correct_for_eigenvalues(ispin, admm_env, ks_matrix)
...
Definition admm_utils.F:53
Define the atomic kind types and their sub types.
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR output in CP2K.
subroutine, public cp_dbcsr_write_sparse_matrix(sparse_matrix, before, after, qs_env, para_env, first_row, last_row, first_col, last_col, scale, output_unit, omit_headers)
...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1179
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public radius_vdw
integer, parameter, public outer_scf_optimizer_sd
integer, parameter, public cdft_beta_constraint
integer, parameter, public cdft_magnetization_constraint
integer, parameter, public outer_scf_optimizer_bisect
integer, parameter, public outer_scf_optimizer_secant
integer, parameter, public becke_cutoff_element
integer, parameter, public radius_default
integer, parameter, public outer_scf_optimizer_broyden
integer, parameter, public cdft_charge_constraint
integer, parameter, public outer_scf_becke_constraint
integer, parameter, public radius_user
integer, parameter, public shape_function_density
integer, parameter, public outer_scf_hirshfeld_constraint
integer, parameter, public radius_covalent
integer, parameter, public shape_function_gaussian
integer, parameter, public radius_single
integer, parameter, public outer_scf_optimizer_newton_ls
integer, parameter, public outer_scf_optimizer_newton
integer, parameter, public becke_cutoff_global
integer, parameter, public cdft_alpha_constraint
integer, parameter, public outer_scf_optimizer_diis
integer, parameter, public ot_precond_full_all
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Types and basic routines needed for a kpoint calculation.
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
Interface to the message passing library MPI.
Define the data structure for the particle information.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public kcalmol
Definition physcon.F:171
real(kind=dp), parameter, public evolt
Definition physcon.F:183
types of preconditioners
computes preconditioners, and implements methods to apply them currently used in qs_ot
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
integer, parameter, public mixed_bc
integer, parameter, public mixed_periodic_bc
integer, parameter, public periodic_bc
container for various plainwaves related things
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_implicit
Routines for image charge calculation within QM/MM.
subroutine, public print_image_coefficients(image_coeff, qs_env)
Print image coefficients.
Control parameters for optimizers that work with CDFT constraints.
subroutine, public cdft_opt_type_write(cdft_opt_control, optimizer, output_unit)
writes information about the CDFT optimizer object
Defines CDFT control structures.
container for information about total charges on the grids
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition and initialisation of the mo data type.
Definition qs_mo_io.F:21
subroutine, public write_mo_set_to_output_unit(mo_set, atomic_kind_set, qs_kind_set, particle_set, dft_section, before, kpoint, final_mos, spin, solver_method, rtp, cpart, sim_step, umo_set)
Write MO information to output file (eigenvalues, occupation numbers, coefficients)
Definition qs_mo_io.F:995
collects routines that perform operations directly related to MOs
subroutine, public calculate_magnitude(mo_array, mo_mag_min, mo_mag_max)
...
subroutine, public calculate_orthonormality(orthonormality, mo_array, matrix_s)
...
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Self-consistent continuum solvation (SCCS) model implementation.
Definition qs_sccs.F:29
subroutine, public print_sccs_results(energy, sccs_control, output_unit)
Print SCCS results.
Definition qs_sccs.F:1212
subroutine, public qs_scf_print_summary(output_unit, qs_env)
writes a summary of information after scf
subroutine, public qs_scf_cdft_constraint_info(output_unit, cdft_control)
writes CDFT constraint information
subroutine, public qs_scf_loop_print(qs_env, scf_env, para_env)
collects the 'heavy duty' printing tasks out of the SCF loop
subroutine, public qs_scf_initial_info(output_unit, mos, dft_control)
writes basic information at the beginning of an scf run
subroutine, public qs_scf_outer_loop_info(output_unit, scf_control, scf_env, energy, total_steps, should_stop, outer_loop_converged)
writes basic information obtained in a scf outer loop step
subroutine, public qs_scf_write_mos(qs_env, scf_env, final_mos)
Write the MO eigenvector, eigenvalues, and occupation numbers to the output unit.
subroutine, public qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy)
writes basic information obtained in a scf step
subroutine, public qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, energy, total_steps, should_stop, outer_loop_converged, cdft_loop)
writes CDFT constraint information and optionally CDFT scf loop info
subroutine, public qs_scf_cdft_initial_info(output_unit, cdft_control)
writes information about the CDFT env
module that contains the definitions of the scf types
integer, parameter, public ot_method_nr
integer, parameter, public special_diag_method_nr
parameters that control an scf iteration
stores some data used in wavefunction fitting
Definition admm_types.F:120
Provides all information about an atomic kind.
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Contains information about kpoints.
stores all the informations relevant to an mpi environment
contained for different pw related things
Container for information about total charges on the grids.
Provides all information about a quickstep kind.
keeps the density in various representations, keeping track of which ones are valid.