(git:0de0cc2)
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
12  USE atomic_kind_types, ONLY: atomic_kind_type
13  USE cp_blacs_env, ONLY: cp_blacs_env_type
14  USE cp_control_types, ONLY: dft_control_type
18  cp_fm_struct_type
19  USE cp_fm_types, ONLY: cp_fm_init_random,&
20  cp_fm_type
22  cp_logger_type
23  USE cp_output_handling, ONLY: cp_p_file,&
27  USE cp_units, ONLY: cp_unit_from_cp2k
28  USE dbcsr_api, ONLY: dbcsr_p_type,&
29  dbcsr_type
30  USE input_constants, ONLY: &
40  section_vals_type,&
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
47  USE message_passing, ONLY: mp_para_env_type
48  USE particle_types, ONLY: particle_type
49  USE physcon, ONLY: evolt,&
50  kcalmol
51  USE preconditioner_types, ONLY: preconditioner_type
52  USE ps_implicit_types, ONLY: mixed_bc,&
54  neumann_bc,&
56  USE pw_env_types, ONLY: pw_env_type
60  USE qs_cdft_types, ONLY: cdft_control_type
61  USE qs_charges_types, ONLY: qs_charges_type
62  USE qs_energy_types, ONLY: qs_energy_type
63  USE qs_environment_types, ONLY: get_qs_env,&
64  qs_environment_type
65  USE qs_kind_types, ONLY: qs_kind_type
69  calculate_subspace_eigenvalues
70  USE qs_mo_occupation, ONLY: set_mo_occupation
71  USE qs_mo_types, ONLY: allocate_mo_set,&
73  get_mo_set,&
74  init_mo_set,&
75  mo_set_type
77  USE qs_rho_types, ONLY: qs_rho_get,&
78  qs_rho_type
79  USE qs_sccs, ONLY: print_sccs_results
80  USE qs_scf_types, ONLY: ot_method_nr,&
81  qs_scf_env_type,&
83  USE scf_control_types, ONLY: scf_control_type
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 
102 CONTAINS
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)
1010  CASE (outer_scf_optimizer_sd)
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
1069  CASE (shape_function_density)
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"
1081  CASE (becke_cutoff_element)
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)
1153  CASE (cdft_charge_constraint)
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
1169  CASE (cdft_alpha_constraint)
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
1177  CASE (cdft_beta_constraint)
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 
1203 END 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
Definition: cp_blacs_env.F:15
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
Definition: cp_fm_struct.F:14
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
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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
Definition: cp_fm_types.F:452
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.
Definition: kpoint_types.F:15
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
Definition: pw_env_types.F:14
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.
Definition: qs_cdft_types.F:14
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: qs_kind_types.F:23
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
Definition: qs_mo_methods.F:14
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...
Definition: qs_mo_types.F:206
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
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.
Definition: qs_mo_types.F:397
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 ...
Definition: qs_mo_types.F:245
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...
Definition: qs_rho_types.F:18
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...
Definition: qs_rho_types.F:229
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
Definition: qs_scf_types.F:14
integer, parameter, public ot_method_nr
Definition: qs_scf_types.F:51
integer, parameter, public special_diag_method_nr
Definition: qs_scf_types.F:51
parameters that control an scf iteration