(git:e5b1968)
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%do_tblite) 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 "Electrostatic energy: ", energy%el_stat, &
677 "Non-self consistent dispersion energy: ", energy%dispersion, &
678 "Self-consistent dispersion energy: ", energy%dispersion_sc, &
679 "Correction for halogen bonding: ", energy%xtb_xb_inter
680 ELSE
681 IF (dft_control%qs_control%xtb_control%gfn_type == 0) THEN
682 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
683 "Core Hamiltonian energy: ", energy%core, &
684 "Repulsive potential energy: ", energy%repulsive, &
685 "SRB Correction energy: ", energy%srb, &
686 "Charge equilibration energy: ", energy%eeq, &
687 "Dispersion energy: ", energy%dispersion
688 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 1) THEN
689 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
690 "Core Hamiltonian energy: ", energy%core, &
691 "Repulsive potential energy: ", energy%repulsive, &
692 "Electronic energy: ", energy%hartree, &
693 "DFTB3 3rd order energy: ", energy%dftb3, &
694 "Dispersion energy: ", energy%dispersion
695 IF (dft_control%qs_control%xtb_control%xb_interaction) &
696 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
697 "Correction for halogen bonding: ", energy%xtb_xb_inter
698 ELSEIF (dft_control%qs_control%xtb_control%gfn_type == 2) THEN
699 cpabort("gfn_typ 2 NYA")
700 ELSE
701 cpabort("invalid gfn_typ")
702 END IF
703 END IF
704 IF (dft_control%qs_control%xtb_control%do_nonbonded) &
705 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
706 "Correction for nonbonded interactions: ", energy%xtb_nonbonded
707 IF (energy%efield /= 0.0_dp) &
708 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
709 "Electric field interaction energy: ", energy%efield
710 ELSE
711 IF (dft_control%do_admm) THEN
712 exc_energy = energy%exc + energy%exc_aux_fit
713 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1 + energy%exc1_aux_fit
714 ELSE
715 exc_energy = energy%exc
716 IF (gapw .OR. gapw_xc) exc1_energy = energy%exc1
717 END IF
718
719 IF (psolver .EQ. pw_poisson_implicit) THEN
720 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
721 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
722 SELECT CASE (bc)
724 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
725 "Overlap energy of the core charge distribution:", energy%core_overlap, &
726 "Self energy of the core charge distribution: ", energy%core_self, &
727 "Core Hamiltonian energy: ", energy%core, &
728 "Hartree energy: ", implicit_ps_ehartree, &
729 "Electric enthalpy: ", energy%hartree, &
730 "Exchange-correlation energy: ", exc_energy
731 CASE (periodic_bc, neumann_bc)
732 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
733 "Overlap energy of the core charge distribution:", energy%core_overlap, &
734 "Self energy of the core charge distribution: ", energy%core_self, &
735 "Core Hamiltonian energy: ", energy%core, &
736 "Hartree energy: ", energy%hartree, &
737 "Exchange-correlation energy: ", exc_energy
738 END SELECT
739 ELSE
740 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
741 "Overlap energy of the core charge distribution:", energy%core_overlap, &
742 "Self energy of the core charge distribution: ", energy%core_self, &
743 "Core Hamiltonian energy: ", energy%core, &
744 "Hartree energy: ", energy%hartree, &
745 "Exchange-correlation energy: ", exc_energy
746 END IF
747 IF (energy%e_hartree /= 0.0_dp) &
748 WRITE (unit=output_unit, fmt="(T3,A,/,T3,A,T56,F25.14)") &
749 "Coulomb Electron-Electron Interaction Energy ", &
750 "- Already included in the total Hartree term ", energy%e_hartree
751 IF (energy%ex /= 0.0_dp) &
752 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
753 "Hartree-Fock Exchange energy: ", energy%ex
754 IF (energy%dispersion /= 0.0_dp) &
755 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
756 "Dispersion energy: ", energy%dispersion
757 IF (energy%gcp /= 0.0_dp) &
758 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
759 "gCP energy: ", energy%gcp
760 IF (energy%efield /= 0.0_dp) &
761 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
762 "Electric field interaction energy: ", energy%efield
763 IF (gapw) THEN
764 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
765 "GAPW| Exc from hard and soft atomic rho1: ", exc1_energy, &
766 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
767 END IF
768 IF (gapw_xc) THEN
769 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
770 "GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
771 END IF
772 END IF
773 IF (dft_control%hairy_probes .EQV. .true.) THEN
774 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
775 "Electronic entropic energy:", energy%kTS
776 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
777 "Fermi energy:", energy%efermi
778 END IF
779 IF (dft_control%smear) THEN
780 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
781 "Electronic entropic energy:", energy%kTS
782 WRITE (unit=output_unit, fmt="((T3,A,T56,F25.14))") &
783 "Fermi energy:", energy%efermi
784 END IF
785 IF (dft_control%dft_plus_u) THEN
786 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
787 "DFT+U energy:", energy%dft_plus_u
788 END IF
789 IF (dft_control%do_sccs) THEN
790 WRITE (unit=output_unit, fmt="(A)") ""
791 CALL print_sccs_results(energy, dft_control%sccs_control, output_unit)
792 END IF
793 IF (qmmm) THEN
794 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
795 "QM/MM Electrostatic energy: ", energy%qmmm_el
796 IF (qs_env%qmmm_env_qm%image_charge) THEN
797 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
798 "QM/MM image charge energy: ", energy%image_charge
799 END IF
800 END IF
801 IF (dft_control%qs_control%mulliken_restraint) THEN
802 WRITE (unit=output_unit, fmt="(T3,A,T56,F25.14)") &
803 "Mulliken restraint energy: ", energy%mulliken
804 END IF
805 IF (dft_control%qs_control%semi_empirical) THEN
806 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
807 "Total energy [eV]: ", energy%total*evolt
808 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
809 "Atomic reference energy [eV]: ", energy%core_self*evolt, &
810 "Heat of formation [kcal/mol]: ", &
811 (energy%total + energy%core_self)*kcalmol
812 ELSE
813 WRITE (unit=output_unit, fmt="(/,(T3,A,T56,F25.14))") &
814 "Total energy: ", energy%total
815 END IF
816 IF (qmmm) THEN
817 IF (qs_env%qmmm_env_qm%image_charge) THEN
818 CALL print_image_coefficients(qs_env%image_coeff, qs_env)
819 END IF
820 END IF
821 CALL m_flush(output_unit)
822 END IF
823
824 CALL timestop(handle)
825
826 END SUBROUTINE qs_scf_print_scf_summary
827
828! **************************************************************************************************
829!> \brief collects the 'heavy duty' printing tasks out of the SCF loop
830!> \param qs_env ...
831!> \param scf_env ...
832!> \param para_env ...
833!> \par History
834!> 03.2006 created [Joost VandeVondele]
835! **************************************************************************************************
836 SUBROUTINE qs_scf_loop_print(qs_env, scf_env, para_env)
837 TYPE(qs_environment_type), POINTER :: qs_env
838 TYPE(qs_scf_env_type), POINTER :: scf_env
839 TYPE(mp_para_env_type), POINTER :: para_env
840
841 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_scf_loop_print'
842
843 INTEGER :: after, handle, ic, ispin, iw
844 LOGICAL :: do_kpoints, omit_headers
845 REAL(kind=dp) :: mo_mag_max, mo_mag_min, orthonormality
846 TYPE(cp_logger_type), POINTER :: logger
847 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s
848 TYPE(dft_control_type), POINTER :: dft_control
849 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
850 TYPE(qs_rho_type), POINTER :: rho
851 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section
852
853 logger => cp_get_default_logger()
854 CALL timeset(routinen, handle)
855
856 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, &
857 do_kpoints=do_kpoints)
858
859 dft_section => section_vals_get_subs_vals(input, "DFT")
860 scf_section => section_vals_get_subs_vals(dft_section, "SCF")
861
862 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
863 DO ispin = 1, dft_control%nspins
864
865 IF (btest(cp_print_key_should_output(logger%iter_info, &
866 dft_section, "PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN
867 CALL get_qs_env(qs_env, rho=rho)
868 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
869 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/DENSITY", &
870 extension=".Log")
871 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
872 after = min(max(after, 1), 16)
873 DO ic = 1, SIZE(matrix_p, 2)
874 CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, ic)%matrix, 4, after, qs_env, para_env, &
875 output_unit=iw, omit_headers=omit_headers)
876 END DO
877 CALL cp_print_key_finished_output(iw, logger, dft_section, &
878 "PRINT%AO_MATRICES/DENSITY")
879 END IF
880
881 IF (btest(cp_print_key_should_output(logger%iter_info, &
882 dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN
883 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", &
884 extension=".Log")
885 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after)
886 after = min(max(after, 1), 16)
887 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks)
888 DO ic = 1, SIZE(matrix_ks, 2)
889 IF (dft_control%qs_control%semi_empirical) THEN
890 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
891 scale=evolt, output_unit=iw, omit_headers=omit_headers)
892 ELSE
893 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, &
894 output_unit=iw, omit_headers=omit_headers)
895 END IF
896 END DO
897 CALL cp_print_key_finished_output(iw, logger, dft_section, &
898 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX")
899 END IF
900
901 END DO
902
903 IF (btest(cp_print_key_should_output(logger%iter_info, &
904 scf_section, "PRINT%MO_ORTHONORMALITY"), cp_p_file)) THEN
905 IF (do_kpoints) THEN
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)') &
910 " K-points: Maximum deviation from MO S-orthonormality not determined"
911 END IF
912 CALL cp_print_key_finished_output(iw, logger, scf_section, &
913 "PRINT%MO_ORTHONORMALITY")
914 ELSE
915 CALL get_qs_env(qs_env, mos=mos)
916 IF (scf_env%method == special_diag_method_nr) THEN
917 CALL calculate_orthonormality(orthonormality, mos)
918 ELSE
919 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
920 CALL calculate_orthonormality(orthonormality, mos, matrix_s(1, 1)%matrix)
921 END IF
922 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", &
923 extension=".scfLog")
924 IF (iw > 0) THEN
925 WRITE (iw, '(T8,A,T61,E20.4)') &
926 " Maximum deviation from MO S-orthonormality", orthonormality
927 END IF
928 CALL cp_print_key_finished_output(iw, logger, scf_section, &
929 "PRINT%MO_ORTHONORMALITY")
930 END IF
931 END IF
932 IF (btest(cp_print_key_should_output(logger%iter_info, &
933 scf_section, "PRINT%MO_MAGNITUDE"), cp_p_file)) THEN
934 IF (do_kpoints) THEN
935 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
936 extension=".scfLog")
937 IF (iw > 0) THEN
938 WRITE (iw, '(T8,A)') &
939 " K-points: Minimum/Maximum MO magnitude not determined"
940 END IF
941 CALL cp_print_key_finished_output(iw, logger, scf_section, &
942 "PRINT%MO_MAGNITUDE")
943 ELSE
944 CALL get_qs_env(qs_env, mos=mos)
945 CALL calculate_magnitude(mos, mo_mag_min, mo_mag_max)
946 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", &
947 extension=".scfLog")
948 IF (iw > 0) THEN
949 WRITE (iw, '(T8,A,T41,2E20.4)') &
950 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max
951 END IF
952 CALL cp_print_key_finished_output(iw, logger, scf_section, &
953 "PRINT%MO_MAGNITUDE")
954 END IF
955 END IF
956
957 CALL timestop(handle)
958
959 END SUBROUTINE qs_scf_loop_print
960
961! **************************************************************************************************
962!> \brief writes CDFT constraint information and optionally CDFT scf loop info
963!> \param output_unit where to write the information
964!> \param scf_control settings of the SCF loop
965!> \param scf_env the env which holds convergence data
966!> \param cdft_control the env which holds information about the constraint
967!> \param energy the total energy
968!> \param total_steps the total number of performed SCF iterations
969!> \param should_stop if the calculation should stop
970!> \param outer_loop_converged logical which determines if the CDFT SCF loop converged
971!> \param cdft_loop logical which determines a CDFT SCF loop is active
972!> \par History
973!> 12.2015 created [Nico Holmberg]
974! **************************************************************************************************
975 SUBROUTINE qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, &
976 energy, total_steps, should_stop, outer_loop_converged, &
977 cdft_loop)
978 INTEGER :: output_unit
979 TYPE(scf_control_type), POINTER :: scf_control
980 TYPE(qs_scf_env_type), POINTER :: scf_env
981 TYPE(cdft_control_type), POINTER :: cdft_control
982 TYPE(qs_energy_type), POINTER :: energy
983 INTEGER :: total_steps
984 LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged, &
985 cdft_loop
986
987 REAL(kind=dp) :: outer_loop_eps
988
989 IF (cdft_loop) THEN
990 outer_loop_eps = sqrt(maxval(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2))
991 IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') &
992 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, &
993 " RMS gradient = ", outer_loop_eps, " energy =", energy%total
994 IF (outer_loop_converged) THEN
995 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
996 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, &
997 " iterations or ", total_steps, " steps"
998 END IF
999 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) &
1000 .AND. .NOT. outer_loop_converged) THEN
1001 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') &
1002 "CDFT SCF loop FAILED to converge after ", &
1003 scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps"
1004 END IF
1005 END IF
1006 CALL qs_scf_cdft_constraint_info(output_unit, cdft_control)
1007
1008 END SUBROUTINE qs_scf_cdft_info
1009
1010! **************************************************************************************************
1011!> \brief writes information about the CDFT env
1012!> \param output_unit where to write the information
1013!> \param cdft_control the CDFT env that stores information about the constraint calculation
1014!> \par History
1015!> 12.2015 created [Nico Holmberg]
1016! **************************************************************************************************
1017 SUBROUTINE qs_scf_cdft_initial_info(output_unit, cdft_control)
1018 INTEGER :: output_unit
1019 TYPE(cdft_control_type), POINTER :: cdft_control
1020
1021 IF (output_unit > 0) THEN
1022 WRITE (output_unit, '(/,A)') &
1023 " ---------------------------------- CDFT --------------------------------------"
1024 WRITE (output_unit, '(A)') &
1025 " Optimizing a density constraint in an external SCF loop "
1026 WRITE (output_unit, '(A)') " "
1027 SELECT CASE (cdft_control%type)
1029 WRITE (output_unit, '(A)') " Type of constraint: Hirshfeld"
1031 WRITE (output_unit, '(A)') " Type of constraint: Becke"
1032 END SELECT
1033 WRITE (output_unit, '(A,I8)') " Number of constraints: ", SIZE(cdft_control%group)
1034 WRITE (output_unit, '(A,L8)') " Using fragment densities:", cdft_control%fragment_density
1035 WRITE (output_unit, '(A)') " "
1036 IF (cdft_control%atomic_charges) WRITE (output_unit, '(A,/)') " Calculating atomic CDFT charges"
1037 SELECT CASE (cdft_control%constraint_control%optimizer)
1039 WRITE (output_unit, '(A)') &
1040 " Minimizer : SD : steepest descent"
1042 WRITE (output_unit, '(A)') &
1043 " Minimizer : DIIS : direct inversion"
1044 WRITE (output_unit, '(A)') &
1045 " in the iterative subspace"
1046 WRITE (output_unit, '(A,I3,A)') &
1047 " using ", &
1048 cdft_control%constraint_control%diis_buffer_length, " DIIS vectors"
1050 WRITE (output_unit, '(A)') &
1051 " Minimizer : BISECT : gradient bisection"
1052 WRITE (output_unit, '(A,I3)') &
1053 " using a trust count of", &
1054 cdft_control%constraint_control%bisect_trust_count
1057 CALL cdft_opt_type_write(cdft_control%constraint_control%cdft_opt_control, &
1058 cdft_control%constraint_control%optimizer, output_unit)
1060 WRITE (output_unit, '(A)') " Minimizer : Secant"
1061 CASE DEFAULT
1062 cpabort("")
1063 END SELECT
1064 WRITE (output_unit, '(/,A,L7)') &
1065 " Reusing OT preconditioner: ", cdft_control%reuse_precond
1066 IF (cdft_control%reuse_precond) THEN
1067 WRITE (output_unit, '(A,I3,A,I3,A)') &
1068 " using old preconditioner for up to ", &
1069 cdft_control%max_reuse, " subsequent CDFT SCF"
1070 WRITE (output_unit, '(A,I3,A,I3,A)') &
1071 " iterations if the relevant loop converged in less than ", &
1072 cdft_control%precond_freq, " steps"
1073 END IF
1074 SELECT CASE (cdft_control%type)
1076 WRITE (output_unit, '(/,A)') " Hirshfeld constraint settings"
1077 WRITE (output_unit, '(A)') " "
1078 SELECT CASE (cdft_control%hirshfeld_control%shape_function)
1080 WRITE (output_unit, '(A, A8)') &
1081 " Shape function type: ", "Gaussian"
1082 WRITE (output_unit, '(A)', advance='NO') &
1083 " Type of Gaussian: "
1084 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape)
1085 CASE (radius_default)
1086 WRITE (output_unit, '(A13)') "Default"
1087 CASE (radius_covalent)
1088 WRITE (output_unit, '(A13)') "Covalent"
1089 CASE (radius_single)
1090 WRITE (output_unit, '(A13)') "Fixed radius"
1091 CASE (radius_vdw)
1092 WRITE (output_unit, '(A13)') "Van der Waals"
1093 CASE (radius_user)
1094 WRITE (output_unit, '(A13)') "User-defined"
1095
1096 END SELECT
1098 WRITE (output_unit, '(A, A8)') &
1099 " Shape function type: ", "Density"
1100 END SELECT
1102 WRITE (output_unit, '(/, A)') " Becke constraint settings"
1103 WRITE (output_unit, '(A)') " "
1104 SELECT CASE (cdft_control%becke_control%cutoff_type)
1105 CASE (becke_cutoff_global)
1106 WRITE (output_unit, '(A,F8.3,A)') &
1107 " Cutoff for partitioning :", cp_unit_from_cp2k(cdft_control%becke_control%rglobal, &
1108 "angstrom"), " angstrom"
1110 WRITE (output_unit, '(A)') &
1111 " Using element specific cutoffs for partitioning"
1112 END SELECT
1113 WRITE (output_unit, '(A,L7)') &
1114 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip
1115 WRITE (output_unit, '(A,L7)') &
1116 " Precompute gradients : ", cdft_control%becke_control%in_memory
1117 WRITE (output_unit, '(A)') " "
1118 IF (cdft_control%becke_control%adjust) &
1119 WRITE (output_unit, '(A)') &
1120 " Using atomic radii to generate a heteronuclear charge partitioning"
1121 WRITE (output_unit, '(A)') " "
1122 IF (.NOT. cdft_control%becke_control%cavity_confine) THEN
1123 WRITE (output_unit, '(A)') &
1124 " No confinement is active"
1125 ELSE
1126 WRITE (output_unit, '(A)') " Confinement using a Gaussian shaped cavity is active"
1127 SELECT CASE (cdft_control%becke_control%cavity_shape)
1128 CASE (radius_single)
1129 WRITE (output_unit, '(A,F8.4, A)') &
1130 " Type of Gaussian : Fixed radius: ", &
1131 cp_unit_from_cp2k(cdft_control%becke_control%rcavity, "angstrom"), " angstrom"
1132 CASE (radius_covalent)
1133 WRITE (output_unit, '(A)') &
1134 " Type of Gaussian : Covalent radius "
1135 CASE (radius_vdw)
1136 WRITE (output_unit, '(A)') &
1137 " Type of Gaussian : vdW radius "
1138 CASE (radius_user)
1139 WRITE (output_unit, '(A)') &
1140 " Type of Gaussian : User radius "
1141 END SELECT
1142 WRITE (output_unit, '(A,ES12.4)') &
1143 " Cavity threshold : ", cdft_control%becke_control%eps_cavity
1144 END IF
1145 END SELECT
1146 WRITE (output_unit, '(/,A)') &
1147 " ---------------------------------- CDFT --------------------------------------"
1148 END IF
1149
1150 END SUBROUTINE qs_scf_cdft_initial_info
1151
1152! **************************************************************************************************
1153!> \brief writes CDFT constraint information
1154!> \param output_unit where to write the information
1155!> \param cdft_control the env which holds information about the constraint
1156!> \par History
1157!> 08.2018 separated from qs_scf_cdft_info to make code callable elsewhere [Nico Holmberg]
1158! **************************************************************************************************
1159 SUBROUTINE qs_scf_cdft_constraint_info(output_unit, cdft_control)
1160 INTEGER :: output_unit
1161 TYPE(cdft_control_type), POINTER :: cdft_control
1162
1163 INTEGER :: igroup
1164
1165 IF (output_unit > 0) THEN
1166 SELECT CASE (cdft_control%type)
1168 WRITE (output_unit, '(/,T3,A,T60)') &
1169 '------------------- Hirshfeld constraint information -------------------'
1171 WRITE (output_unit, '(/,T3,A,T60)') &
1172 '--------------------- Becke constraint information ---------------------'
1173 CASE DEFAULT
1174 cpabort("Unknown CDFT constraint.")
1175 END SELECT
1176 DO igroup = 1, SIZE(cdft_control%target)
1177 IF (igroup > 1) WRITE (output_unit, '(T3,A)') ' '
1178 WRITE (output_unit, '(T3,A,T54,(3X,I18))') &
1179 'Atomic group :', igroup
1180 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1182 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1183 WRITE (output_unit, '(T3,A,T42,A)') &
1184 'Type of constraint :', adjustr('Charge density constraint (frag.)')
1185 ELSE
1186 WRITE (output_unit, '(T3,A,T50,A)') &
1187 'Type of constraint :', adjustr('Charge density constraint')
1188 END IF
1190 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1191 WRITE (output_unit, '(T3,A,T35,A)') &
1192 'Type of constraint :', adjustr('Magnetization density constraint (frag.)')
1193 ELSE
1194 WRITE (output_unit, '(T3,A,T43,A)') &
1195 'Type of constraint :', adjustr('Magnetization density constraint')
1196 END IF
1198 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1199 WRITE (output_unit, '(T3,A,T38,A)') &
1200 'Type of constraint :', adjustr('Alpha spin density constraint (frag.)')
1201 ELSE
1202 WRITE (output_unit, '(T3,A,T46,A)') &
1203 'Type of constraint :', adjustr('Alpha spin density constraint')
1204 END IF
1206 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN
1207 WRITE (output_unit, '(T3,A,T39,A)') &
1208 'Type of constraint :', adjustr('Beta spin density constraint (frag.)')
1209 ELSE
1210 WRITE (output_unit, '(T3,A,T47,A)') &
1211 'Type of constraint :', adjustr('Beta spin density constraint')
1212 END IF
1213 CASE DEFAULT
1214 cpabort("Unknown constraint type.")
1215 END SELECT
1216 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1217 'Target value of constraint :', cdft_control%target(igroup)
1218 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1219 'Current value of constraint :', cdft_control%value(igroup)
1220 WRITE (output_unit, '(T3,A,T59,(3X,ES13.3))') &
1221 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup)
1222 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') &
1223 'Strength of constraint :', cdft_control%strength(igroup)
1224 END DO
1225 WRITE (output_unit, '(T3,A)') &
1226 '------------------------------------------------------------------------'
1227 END IF
1228
1229 END SUBROUTINE qs_scf_cdft_constraint_info
1230
1231END 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, tb_tblite)
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.