(git:9754b87)
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,I5,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") &
539 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
540 t2 - t1, energy%total
541 ELSE
542 IF ((abs(scf_env%iter_delta) < 1.0e-8_dp) .OR. &
543 (abs(scf_env%iter_delta) >= 1.0e5_dp)) THEN
544 WRITE (unit=output_unit, &
545 fmt="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") &
546 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
547 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
548 ELSE
549 WRITE (unit=output_unit, &
550 fmt="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") &
551 scf_env%iter_count, trim(scf_env%iter_method), scf_env%iter_param, &
552 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old
553 END IF
554 END IF
555 END IF
556
557 END SUBROUTINE qs_scf_loop_info
558
559! **************************************************************************************************
560!> \brief writes rather detailed summary of densities and energies
561!> after the SCF
562!> \param output_unit ...
563!> \param rho ...
564!> \param qs_charges ...
565!> \param energy ...
566!> \param nelectron_total ...
567!> \param dft_control ...
568!> \param qmmm ...
569!> \param qs_env ...
570!> \param gapw ...
571!> \param gapw_xc ...
572!> \par History
573!> 03.2006 created [Joost VandeVondele]
574!> 10.2019 print dipole moment [SGh]
575!> 11.2022 print SCCS results [MK]
576! **************************************************************************************************
577 SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, &
578 dft_control, qmmm, qs_env, gapw, gapw_xc)
579 INTEGER, INTENT(IN) :: output_unit
580 TYPE(qs_rho_type), POINTER :: rho
581 TYPE(qs_charges_type), POINTER :: qs_charges
582 TYPE(qs_energy_type), POINTER :: energy
583 INTEGER, INTENT(IN) :: nelectron_total
584 TYPE(dft_control_type), POINTER :: dft_control
585 LOGICAL, INTENT(IN) :: qmmm
586 TYPE(qs_environment_type), POINTER :: qs_env
587 LOGICAL, INTENT(IN) :: gapw, gapw_xc
588
589 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_scf_print_scf_summary'
590
591 INTEGER :: bc, handle, ispin, psolver
592 REAL(kind=dp) :: exc1_energy, exc_energy, &
593 implicit_ps_ehartree, tot1_h, tot1_s
594 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
595 TYPE(pw_env_type), POINTER :: pw_env
596
597 NULLIFY (tot_rho_r, pw_env)
598 CALL timeset(routinen, handle)
599
600 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
601 psolver = pw_env%poisson_env%parameters%solver
602
603 IF (output_unit > 0) THEN
604 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r)
605 IF (.NOT. (dft_control%qs_control%semi_empirical .OR. &
606 dft_control%qs_control%xtb .OR. &
607 dft_control%qs_control%dftb)) THEN
608 WRITE (unit=output_unit, fmt="(/,(T3,A,T41,2F20.10))") &
609 "Electronic density on regular grids: ", &
610 accurate_sum(tot_rho_r), &
611 accurate_sum(tot_rho_r) + nelectron_total, &
612 "Core density on regular grids:", &
613 qs_charges%total_rho_core_rspace, &
614 qs_charges%total_rho_core_rspace - real(nelectron_total + dft_control%charge, dp)
615
616 IF (dft_control%correct_surf_dip) THEN
617 WRITE (unit=output_unit, fmt="((T3,A,/,T3,A,T41,F20.10))") &
618 "Total dipole moment perpendicular to ", &
619 "the slab [electrons-Angstroem]: ", &
620 qs_env%surface_dipole_moment
621 END IF
622
623 IF (gapw) THEN
624 tot1_h = qs_charges%total_rho1_hard(1)
625 tot1_s = qs_charges%total_rho1_soft(1)
626 DO ispin = 2, dft_control%nspins
627 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
628 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
629 END DO
630 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
631 "Hard and soft densities (Lebedev):", &
632 tot1_h, tot1_s
633 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
634 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
635 accurate_sum(tot_rho_r) + tot1_h - tot1_s, &
636 "Total charge density (r-space): ", &
637 accurate_sum(tot_rho_r) + tot1_h - tot1_s &
638 + qs_charges%total_rho_core_rspace, &
639 "Total Rho_soft + Rho0_soft (g-space):", &
640 qs_charges%total_rho_gspace
641 ELSE
642 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
643 "Total charge density on r-space grids: ", &
644 accurate_sum(tot_rho_r) + &
645 qs_charges%total_rho_core_rspace, &
646 "Total charge density g-space grids: ", &
647 qs_charges%total_rho_gspace
648 END IF
649 END IF
650 IF (dft_control%qs_control%semi_empirical) THEN
651 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
652 "Core-core repulsion energy [eV]: ", energy%core_overlap*evolt, &
653 "Core Hamiltonian energy [eV]: ", energy%core*evolt, &
654 "Two-electron integral energy [eV]: ", energy%hartree*evolt, &
655 "Electronic energy [eV]: ", &
656 (energy%core + 0.5_dp*energy%hartree)*evolt
657 IF (energy%dispersion /= 0.0_dp) &
658 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
659 "Dispersion energy [eV]: ", energy%dispersion*evolt
660 ELSEIF (dft_control%qs_control%dftb) THEN
661 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
662 "Core Hamiltonian energy: ", energy%core, &
663 "Repulsive potential energy: ", energy%repulsive, &
664 "Electronic energy: ", energy%hartree, &
665 "Dispersion energy: ", energy%dispersion
666 IF (energy%dftb3 /= 0.0_dp) &
667 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
668 "DFTB3 3rd order energy: ", energy%dftb3
669 IF (energy%efield /= 0.0_dp) &
670 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
671 "Electric field interaction energy: ", energy%efield
672 ELSEIF (dft_control%qs_control%xtb) THEN
673 IF (dft_control%qs_control%xtb_control%gfn_type == 0) THEN
674 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
675 "Core Hamiltonian energy: ", energy%core, &
676 "Repulsive potential energy: ", energy%repulsive, &
677 "SRB Correction energy: ", energy%srb, &
678 "Charge equilibration energy: ", energy%eeq, &
679 "Dispersion energy: ", energy%dispersion
680 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 1) THEN
681 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
682 "Core Hamiltonian energy: ", energy%core, &
683 "Repulsive potential energy: ", energy%repulsive, &
684 "Electronic energy: ", energy%hartree, &
685 "DFTB3 3rd order energy: ", energy%dftb3, &
686 "Dispersion energy: ", energy%dispersion
687 IF (dft_control%qs_control%xtb_control%xb_interaction) &
688 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
689 "Correction for halogen bonding: ", energy%xtb_xb_inter
690 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 2) THEN
691 cpabort("gfn_typ 2 NYA")
692 ELSE
693 cpabort("invalid gfn_typ")
694 END IF
695 IF (dft_control%qs_control%xtb_control%do_nonbonded) &
696 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
697 "Correction for nonbonded interactions: ", energy%xtb_nonbonded
698 IF (energy%efield /= 0.0_dp) &
699 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
700 "Electric field interaction energy: ", energy%efield
701 ELSE
702 IF (dft_control%do_admm) THEN
703 exc_energy = energy%exc + energy%exc_aux_fit
704 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
705 ELSE
706 exc_energy = energy%exc
707 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
708 END IF
709
710 IF (psolver .EQ. pw_poisson_implicit) THEN
711 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
712 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
713 SELECT CASE (bc)
715 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
716 "Overlap energy of the core charge distribution:", energy%core_overlap, &
717 "Self energy of the core charge distribution: ", energy%core_self, &
718 "Core Hamiltonian energy: ", energy%core, &
719 "Hartree energy: ", implicit_ps_ehartree, &
720 "Electric enthalpy: ", energy%hartree, &
721 "Exchange-correlation energy: ", exc_energy
722 CASE (periodic_bc, neumann_bc)
723 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
724 "Overlap energy of the core charge distribution:", energy%core_overlap, &
725 "Self energy of the core charge distribution: ", energy%core_self, &
726 "Core Hamiltonian energy: ", energy%core, &
727 "Hartree energy: ", energy%hartree, &
728 "Exchange-correlation energy: ", exc_energy
729 END SELECT
730 ELSE
731 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
732 "Overlap energy of the core charge distribution:", energy%core_overlap, &
733 "Self energy of the core charge distribution: ", energy%core_self, &
734 "Core Hamiltonian energy: ", energy%core, &
735 "Hartree energy: ", energy%hartree, &
736 "Exchange-correlation energy: ", exc_energy
737 END IF
738 IF (energy%e_hartree /= 0.0_dp) &
739 WRITE (unit=output_unit, fmt="(T3,A,/,T3,A,T56,F25.14)") &
740 "Coulomb Electron-Electron Interaction Energy ", &
741 "- Already included in the total Hartree term ", energy%e_hartree
742 IF (energy%ex /= 0.0_dp) &
743 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
744 "Hartree-Fock Exchange energy: ", energy%ex
745 IF (energy%dispersion /= 0.0_dp) &
746 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
747 "Dispersion energy: ", energy%dispersion
748 IF (energy%gcp /= 0.0_dp) &
749 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
750 "gCP energy: ", energy%gcp
751 IF (energy%efield /= 0.0_dp) &
752 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
753 "Electric field interaction energy: ", energy%efield
754 IF (gapw) THEN
755 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
756 "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
757 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
758 END IF
759 IF (gapw_xc) THEN
760 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
761 "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
762 END IF
763 END IF
764 IF (dft_control%smear) THEN
765 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
766 "Electronic entropic energy:", energy%kTS
767 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
768 "Fermi energy:", energy%efermi
769 END IF
770 IF (dft_control%dft_plus_u) THEN
771 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
772 "DFT+U energy:", energy%dft_plus_u
773 END IF
774 IF (dft_control%do_sccs) THEN
775 WRITE (unit=output_unit, fmt="(A)") ""
776 CALL print_sccs_results(energy, dft_control%sccs_control, output_unit)
777 END IF
778 IF (qmmm) THEN
779 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
780 "QM/MM Electrostatic energy: ", energy%qmmm_el
781 IF (qs_env%qmmm_env_qm%image_charge) THEN
782 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
783 "QM/MM image charge energy: ", energy%image_charge
784 END IF
785 END IF
786 IF (dft_control%qs_control%mulliken_restraint) THEN
787 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
788 "Mulliken restraint energy: ", energy%mulliken
789 END IF
790 IF (dft_control%qs_control%semi_empirical) THEN
791 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
792 "Total energy [eV]: ", energy%total*evolt
793 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
794 "Atomic reference energy [eV]: ", energy%core_self*evolt, &
795 "Heat of formation [kcal/mol]: ", &
796 (energy%total + energy%core_self)*kcalmol
797 ELSE
798 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
799 "Total energy: ", energy%total
800 END IF
801 IF (qmmm) THEN
802 IF (qs_env%qmmm_env_qm%image_charge) THEN
803 CALL print_image_coefficients(qs_env%image_coeff, qs_env)
804 END IF
805 END IF
806 CALL m_flush(output_unit)
807 END IF
808
809 CALL timestop(handle)
810
811 END SUBROUTINE qs_scf_print_scf_summary
812
813! **************************************************************************************************
814!> \brief collects the 'heavy duty' printing tasks out of the SCF loop
815!> \param qs_env ...
816!> \param scf_env ...
817!> \param para_env ...
818!> \par History
819!> 03.2006 created [Joost VandeVondele]
820! **************************************************************************************************
821 SUBROUTINE qs_scf_loop_print(qs_env, scf_env, para_env)
822 TYPE(qs_environment_type), POINTER :: qs_env
823 TYPE(qs_scf_env_type), POINTER :: scf_env
824 TYPE(mp_para_env_type), POINTER :: para_env
825
826 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_scf_loop_print'
827
828 INTEGER :: after, handle, ic, ispin, iw
829 LOGICAL :: do_kpoints, omit_headers
830 REAL(kind=dp) :: mo_mag_max, mo_mag_min, orthonormality
831 TYPE(cp_logger_type), POINTER :: logger
832 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
833 TYPE(dft_control_type), POINTER :: dft_control
834 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
835 TYPE(qs_rho_type), POINTER :: rho
836 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
837
838 logger => cp_get_default_logger()
839 CALL timeset(routinen, handle)
840
841 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
842 do_kpoints=do_kpoints)
843
844 dft_section => section_vals_get_subs_vals(input, "DFT")
845 scf_section => section_vals_get_subs_vals(dft_section, "SCF")
846
847 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
848 DO ispin = 1, dft_control%nspins
849
850 IF (btest(cp_print_key_should_output(logger%iter_info, &
851 dft_section, "PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
852 CALL get_qs_env(qs_env, rho=rho)
853 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
854 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/DENSITY", &
855 extension=".Log")
856 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
857 after = min(max(after, 1), 16)
858 DO ic = 1, SIZE(matrix_p, 2)
859 CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, ic)%matrix, 4, after, qs_env, para_env, &
860 output_unit=iw, omit_headers=omit_headers)
861 END DO
862 CALL cp_print_key_finished_output(iw, logger, dft_section, &
863 "PRINT%AO_MATRICES/DENSITY")
864 END IF
865
866 IF (btest(cp_print_key_should_output(logger%iter_info, &
867 dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN
868 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
869 extension=".Log")
870 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
871 after = min(max(after, 1), 16)
872 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
873 DO ic = 1, SIZE(matrix_ks, 2)
874 IF (dft_control%qs_control%semi_empirical) THEN
875 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
876 scale=evolt, output_unit=iw, omit_headers=omit_headers)
877 ELSE
878 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
879 output_unit=iw, omit_headers=omit_headers)
880 END IF
881 END DO
882 CALL cp_print_key_finished_output(iw, logger, dft_section, &
883 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
884 END IF
885
886 END DO
887
888 IF (btest(cp_print_key_should_output(logger%iter_info, &
889 scf_section, "PRINT%MO_ORTHONORMALITY"), cp_p_file)) THEN
890 IF (do_kpoints) THEN
891 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
892 extension=".scfLog")
893 IF (iw > 0) THEN
894 WRITE (iw, '(T8,A)') &
895 " K-points: Maximum deviation from MO S-orthonormality not determined"
896 END IF
897 CALL cp_print_key_finished_output(iw, logger, scf_section, &
898 "PRINT%MO_ORTHONORMALITY")
899 ELSE
900 CALL get_qs_env(qs_env, mos=mos)
901 IF (scf_env%method == special_diag_method_nr) THEN
902 CALL calculate_orthonormality(orthonormality, mos)
903 ELSE
904 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
905 CALL calculate_orthonormality(orthonormality, mos, matrix_s(1, 1)%matrix)
906 END IF
907 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
908 extension=".scfLog")
909 IF (iw > 0) THEN
910 WRITE (iw, '(T8,A,T61,E20.4)') &
911 " Maximum deviation from MO S-orthonormality", orthonormality
912 END IF
913 CALL cp_print_key_finished_output(iw, logger, scf_section, &
914 "PRINT%MO_ORTHONORMALITY")
915 END IF
916 END IF
917 IF (btest(cp_print_key_should_output(logger%iter_info, &
918 scf_section, "PRINT%MO_MAGNITUDE"), cp_p_file)) THEN
919 IF (do_kpoints) THEN
920 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
921 extension=".scfLog")
922 IF (iw > 0) THEN
923 WRITE (iw, '(T8,A)') &
924 " K-points: Minimum/Maximum MO magnitude not determined"
925 END IF
926 CALL cp_print_key_finished_output(iw, logger, scf_section, &
927 "PRINT%MO_MAGNITUDE")
928 ELSE
929 CALL get_qs_env(qs_env, mos=mos)
930 CALL calculate_magnitude(mos, mo_mag_min, mo_mag_max)
931 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
932 extension=".scfLog")
933 IF (iw > 0) THEN
934 WRITE (iw, '(T8,A,T41,2E20.4)') &
935 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
936 END IF
937 CALL cp_print_key_finished_output(iw, logger, scf_section, &
938 "PRINT%MO_MAGNITUDE")
939 END IF
940 END IF
941
942 CALL timestop(handle)
943
944 END SUBROUTINE qs_scf_loop_print
945
946! **************************************************************************************************
947!> \brief writes CDFT constraint information and optionally CDFT scf loop info
948!> \param output_unit where to write the information
949!> \param scf_control settings of the SCF loop
950!> \param scf_env the env which holds convergence data
951!> \param cdft_control the env which holds information about the constraint
952!> \param energy the total energy
953!> \param total_steps the total number of performed SCF iterations
954!> \param should_stop if the calculation should stop
955!> \param outer_loop_converged logical which determines if the CDFT SCF loop converged
956!> \param cdft_loop logical which determines a CDFT SCF loop is active
957!> \par History
958!> 12.2015 created [Nico Holmberg]
959! **************************************************************************************************
960 SUBROUTINE qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
961 energy, total_steps, should_stop, outer_loop_converged, &
962 cdft_loop)
963 INTEGER :: output_unit
964 TYPE(scf_control_type), POINTER :: scf_control
965 TYPE(qs_scf_env_type), POINTER :: scf_env
966 TYPE(cdft_control_type), POINTER :: cdft_control
967 TYPE(qs_energy_type), POINTER :: energy
968 INTEGER :: total_steps
969 LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged, &
970 cdft_loop
971
972 REAL(kind=dp) :: outer_loop_eps
973
974 IF (cdft_loop) THEN
975 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
976 IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') &
977 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
978 " RMS gradient = ", outer_loop_eps, " energy =", energy%total
979 IF (outer_loop_converged) THEN
980 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
981 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
982 " iterations or ", total_steps, " steps"
983 END IF
984 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
985 .AND. .NOT. outer_loop_converged) THEN
986 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
987 "CDFT SCF loop FAILED to converge after ", &
988 scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps"
989 END IF
990 END IF
991 CALL qs_scf_cdft_constraint_info(output_unit, cdft_control)
992
993 END SUBROUTINE qs_scf_cdft_info
994
995! **************************************************************************************************
996!> \brief writes information about the CDFT env
997!> \param output_unit where to write the information
998!> \param cdft_control the CDFT env that stores information about the constraint calculation
999!> \par History
1000!> 12.2015 created [Nico Holmberg]
1001! **************************************************************************************************
1002 SUBROUTINE qs_scf_cdft_initial_info(output_unit, cdft_control)
1003 INTEGER :: output_unit
1004 TYPE(cdft_control_type), POINTER :: cdft_control
1005
1006 IF (output_unit > 0) THEN
1007 WRITE (output_unit, '(/,A)') &
1008 " ---------------------------------- CDFT --------------------------------------"
1009 WRITE (output_unit, '(A)') &
1010 " Optimizing a density constraint in an external SCF loop "
1011 WRITE (output_unit, '(A)') " "
1012 SELECT CASE (cdft_control%type)
1014 WRITE (output_unit, '(A)') " Type of constraint: Hirshfeld"
1016 WRITE (output_unit, '(A)') " Type of constraint: Becke"
1017 END SELECT
1018 WRITE (output_unit, '(A,I8)') " Number of constraints: ", SIZE(cdft_control%group)
1019 WRITE (output_unit, '(A,L8)') " Using fragment densities:", cdft_control%fragment_density
1020 WRITE (output_unit, '(A)') " "
1021 IF (cdft_control%atomic_charges) WRITE (output_unit, '(A,/)') " Calculating atomic CDFT charges"
1022 SELECT CASE (cdft_control%constraint_control%optimizer)
1024 WRITE (output_unit, '(A)') &
1025 " Minimizer : SD : steepest descent"
1027 WRITE (output_unit, '(A)') &
1028 " Minimizer : DIIS : direct inversion"
1029 WRITE (output_unit, '(A)') &
1030 " in the iterative subspace"
1031 WRITE (output_unit, '(A,I3,A)') &
1032 " using ", &
1033 cdft_control%constraint_control%diis_buffer_length, " DIIS vectors"
1035 WRITE (output_unit, '(A)') &
1036 " Minimizer : BISECT : gradient bisection"
1037 WRITE (output_unit, '(A,I3)') &
1038 " using a trust count of", &
1039 cdft_control%constraint_control%bisect_trust_count
1042 CALL cdft_opt_type_write(cdft_control%constraint_control%cdft_opt_control, &
1043 cdft_control%constraint_control%optimizer, output_unit)
1045 WRITE (output_unit, '(A)') " Minimizer : Secant"
1046 CASE DEFAULT
1047 cpabort("")
1048 END SELECT
1049 WRITE (output_unit, '(/,A,L7)') &
1050 " Reusing OT preconditioner: ", cdft_control%reuse_precond
1051 IF (cdft_control%reuse_precond) THEN
1052 WRITE (output_unit, '(A,I3,A,I3,A)') &
1053 " using old preconditioner for up to ", &
1054 cdft_control%max_reuse, " subsequent CDFT SCF"
1055 WRITE (output_unit, '(A,I3,A,I3,A)') &
1056 " iterations if the relevant loop converged in less than ", &
1057 cdft_control%precond_freq, " steps"
1058 END IF
1059 SELECT CASE (cdft_control%type)
1061 WRITE (output_unit, '(/,A)') " Hirshfeld constraint settings"
1062 WRITE (output_unit, '(A)') " "
1063 SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1065 WRITE (output_unit, '(A, A8)') &
1066 " Shape function type: ", "Gaussian"
1067 WRITE (output_unit, '(A)', advance='NO') &
1068 " Type of Gaussian: "
1069 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1070 CASE (radius_default)
1071 WRITE (output_unit, '(A13)') "Default"
1072 CASE (radius_covalent)
1073 WRITE (output_unit, '(A13)') "Covalent"
1074 CASE (radius_single)
1075 WRITE (output_unit, '(A13)') "Fixed radius"
1076 CASE (radius_vdw)
1077 WRITE (output_unit, '(A13)') "Van der Waals"
1078 CASE (radius_user)
1079 WRITE (output_unit, '(A13)') "User-defined"
1080
1081 END SELECT
1083 WRITE (output_unit, '(A, A8)') &
1084 " Shape function type: ", "Density"
1085 END SELECT
1087 WRITE (output_unit, '(/, A)') " Becke constraint settings"
1088 WRITE (output_unit, '(A)') " "
1089 SELECT CASE (cdft_control%becke_control%cutoff_type)
1090 CASE (becke_cutoff_global)
1091 WRITE (output_unit, '(A,F8.3,A)') &
1092 " Cutoff for partitioning :", cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1093 "angstrom"), " angstrom"
1095 WRITE (output_unit, '(A)') &
1096 " Using element specific cutoffs for partitioning"
1097 END SELECT
1098 WRITE (output_unit, '(A,L7)') &
1099 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1100 WRITE (output_unit, '(A,L7)') &
1101 " Precompute gradients : ", cdft_control%becke_control%in_memory
1102 WRITE (output_unit, '(A)') " "
1103 IF (cdft_control%becke_control%adjust) &
1104 WRITE (output_unit, '(A)') &
1105 " Using atomic radii to generate a heteronuclear charge partitioning"
1106 WRITE (output_unit, '(A)') " "
1107 IF (.NOT. cdft_control%becke_control%cavity_confine) THEN
1108 WRITE (output_unit, '(A)') &
1109 " No confinement is active"
1110 ELSE
1111 WRITE (output_unit, '(A)') " Confinement using a Gaussian shaped cavity is active"
1112 SELECT CASE (cdft_control%becke_control%cavity_shape)
1113 CASE (radius_single)
1114 WRITE (output_unit, '(A,F8.4, A)') &
1115 " Type of Gaussian : Fixed radius: ", &
1116 cp_unit_from_cp2k(cdft_control%becke_control%rcavity, "angstrom"), " angstrom"
1117 CASE (radius_covalent)
1118 WRITE (output_unit, '(A)') &
1119 " Type of Gaussian : Covalent radius "
1120 CASE (radius_vdw)
1121 WRITE (output_unit, '(A)') &
1122 " Type of Gaussian : vdW radius "
1123 CASE (radius_user)
1124 WRITE (output_unit, '(A)') &
1125 " Type of Gaussian : User radius "
1126 END SELECT
1127 WRITE (output_unit, '(A,ES12.4)') &
1128 " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1129 END IF
1130 END SELECT
1131 WRITE (output_unit, '(/,A)') &
1132 " ---------------------------------- CDFT --------------------------------------"
1133 END IF
1134
1135 END SUBROUTINE qs_scf_cdft_initial_info
1136
1137! **************************************************************************************************
1138!> \brief writes CDFT constraint information
1139!> \param output_unit where to write the information
1140!> \param cdft_control the env which holds information about the constraint
1141!> \par History
1142!> 08.2018 separated from qs_scf_cdft_info to make code callable elsewhere [Nico Holmberg]
1143! **************************************************************************************************
1144 SUBROUTINE qs_scf_cdft_constraint_info(output_unit, cdft_control)
1145 INTEGER :: output_unit
1146 TYPE(cdft_control_type), POINTER :: cdft_control
1147
1148 INTEGER :: igroup
1149
1150 IF (output_unit > 0) THEN
1151 SELECT CASE (cdft_control%type)
1153 WRITE (output_unit, '(/,T3,A,T60)') &
1154 '------------------- Hirshfeld constraint information -------------------'
1156 WRITE (output_unit, '(/,T3,A,T60)') &
1157 '--------------------- Becke constraint information ---------------------'
1158 CASE DEFAULT
1159 cpabort("Unknown CDFT constraint.")
1160 END SELECT
1161 DO igroup = 1, SIZE(cdft_control%target)
1162 IF (igroup > 1) WRITE (output_unit, '(T3,A)') ' '
1163 WRITE (output_unit, '(T3,A,T54,(3X,I18))') &
1164 'Atomic group :', igroup
1165 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1167 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1168 WRITE (output_unit, '(T3,A,T42,A)') &
1169 'Type of constraint :', adjustr('Charge density constraint (frag.)')
1170 ELSE
1171 WRITE (output_unit, '(T3,A,T50,A)') &
1172 'Type of constraint :', adjustr('Charge density constraint')
1173 END IF
1175 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1176 WRITE (output_unit, '(T3,A,T35,A)') &
1177 'Type of constraint :', adjustr('Magnetization density constraint (frag.)')
1178 ELSE
1179 WRITE (output_unit, '(T3,A,T43,A)') &
1180 'Type of constraint :', adjustr('Magnetization density constraint')
1181 END IF
1183 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1184 WRITE (output_unit, '(T3,A,T38,A)') &
1185 'Type of constraint :', adjustr('Alpha spin density constraint (frag.)')
1186 ELSE
1187 WRITE (output_unit, '(T3,A,T46,A)') &
1188 'Type of constraint :', adjustr('Alpha spin density constraint')
1189 END IF
1191 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1192 WRITE (output_unit, '(T3,A,T39,A)') &
1193 'Type of constraint :', adjustr('Beta spin density constraint (frag.)')
1194 ELSE
1195 WRITE (output_unit, '(T3,A,T47,A)') &
1196 'Type of constraint :', adjustr('Beta spin density constraint')
1197 END IF
1198 CASE DEFAULT
1199 cpabort("Unknown constraint type.")
1200 END SELECT
1201 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1202 'Target value of constraint :', cdft_control%target(igroup)
1203 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1204 'Current value of constraint :', cdft_control%value(igroup)
1205 WRITE (output_unit, '(T3,A,T59,(3X,ES13.3))') &
1206 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1207 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1208 'Strength of constraint :', cdft_control%strength(igroup)
1209 END DO
1210 WRITE (output_unit, '(T3,A)') &
1211 '------------------------------------------------------------------------'
1212 END IF
1213
1214 END SUBROUTINE qs_scf_cdft_constraint_info
1215
1216END 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.