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