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