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