1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6MODULE qs_scf_output 7 USE atomic_kind_types, ONLY: atomic_kind_type 8 USE cp_control_types, ONLY: dft_control_type 9 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_sparse_matrix 10 USE cp_log_handling, ONLY: cp_get_default_logger,& 11 cp_logger_type 12 USE cp_output_handling, ONLY: cp_p_file,& 13 cp_print_key_finished_output,& 14 cp_print_key_should_output,& 15 cp_print_key_unit_nr 16 USE cp_para_types, ONLY: cp_para_env_type 17 USE cp_units, ONLY: cp_unit_from_cp2k 18 USE dbcsr_api, ONLY: dbcsr_p_type 19 USE input_constants, ONLY: & 20 becke_cutoff_element, becke_cutoff_global, cdft_alpha_constraint, cdft_beta_constraint, & 21 cdft_charge_constraint, cdft_magnetization_constraint, outer_scf_becke_constraint, & 22 outer_scf_hirshfeld_constraint, outer_scf_optimizer_bisect, outer_scf_optimizer_broyden, & 23 outer_scf_optimizer_diis, outer_scf_optimizer_newton, outer_scf_optimizer_newton_ls, & 24 outer_scf_optimizer_sd, outer_scf_optimizer_secant, radius_covalent, radius_default, & 25 radius_single, radius_user, radius_vdw, shape_function_density, shape_function_gaussian 26 USE input_section_types, ONLY: section_vals_get_subs_vals,& 27 section_vals_type,& 28 section_vals_val_get 29 USE kahan_sum, ONLY: accurate_sum 30 USE kinds, ONLY: dp 31 USE machine, ONLY: m_flush 32 USE particle_types, ONLY: particle_type 33 USE physcon, ONLY: evolt,& 34 kcalmol 35 USE ps_implicit_types, ONLY: MIXED_BC,& 36 MIXED_PERIODIC_BC,& 37 NEUMANN_BC,& 38 PERIODIC_BC 39 USE pw_env_types, ONLY: pw_env_type 40 USE pw_poisson_types, ONLY: pw_poisson_implicit 41 USE qmmm_image_charge, ONLY: print_image_coefficients 42 USE qs_cdft_opt_types, ONLY: cdft_opt_type_write 43 USE qs_cdft_types, ONLY: cdft_control_type 44 USE qs_charges_types, ONLY: qs_charges_type 45 USE qs_energy_types, ONLY: qs_energy_type 46 USE qs_environment_types, ONLY: get_qs_env,& 47 qs_environment_type 48 USE qs_kind_types, ONLY: qs_kind_type 49 USE qs_mo_io, ONLY: write_mo_set 50 USE qs_mo_methods, ONLY: calculate_magnitude,& 51 calculate_orthonormality 52 USE qs_mo_types, ONLY: get_mo_set,& 53 mo_set_p_type 54 USE qs_rho_types, ONLY: qs_rho_get,& 55 qs_rho_type 56 USE qs_scf_types, ONLY: qs_scf_env_type,& 57 special_diag_method_nr 58 USE scf_control_types, ONLY: scf_control_type 59#include "./base/base_uses.f90" 60 61 IMPLICIT NONE 62 63 PRIVATE 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_output' 66 67 PUBLIC :: qs_scf_loop_info, & 68 qs_scf_print_summary, & 69 qs_scf_loop_print, & 70 qs_scf_outer_loop_info, & 71 qs_scf_initial_info, & 72 qs_scf_write_mos, & 73 qs_scf_cdft_info, & 74 qs_scf_cdft_initial_info, & 75 qs_scf_cdft_constraint_info 76 77CONTAINS 78 79! ************************************************************************************************** 80!> \brief writes a summary of information after scf 81!> \param output_unit ... 82!> \param qs_env ... 83! ************************************************************************************************** 84 SUBROUTINE qs_scf_print_summary(output_unit, qs_env) 85 INTEGER, INTENT(IN) :: output_unit 86 TYPE(qs_environment_type), POINTER :: qs_env 87 88 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_print_summary', & 89 routineP = moduleN//':'//routineN 90 91 INTEGER :: nelectron_total 92 LOGICAL :: gapw, gapw_xc, qmmm 93 TYPE(dft_control_type), POINTER :: dft_control 94 TYPE(qs_charges_type), POINTER :: qs_charges 95 TYPE(qs_energy_type), POINTER :: energy 96 TYPE(qs_rho_type), POINTER :: rho 97 TYPE(qs_scf_env_type), POINTER :: scf_env 98 99 NULLIFY (rho, energy, dft_control, scf_env, qs_charges) 100 CALL get_qs_env(qs_env=qs_env, rho=rho, energy=energy, dft_control=dft_control, & 101 scf_env=scf_env, qs_charges=qs_charges) 102 103 gapw = dft_control%qs_control%gapw 104 gapw_xc = dft_control%qs_control%gapw_xc 105 qmmm = qs_env%qmmm 106 nelectron_total = scf_env%nelectron 107 108 CALL qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, & 109 dft_control, qmmm, qs_env, gapw, gapw_xc) 110 111 END SUBROUTINE qs_scf_print_summary 112 113! ************************************************************************************************** 114!> \brief writes basic information at the beginning of an scf run 115!> \param output_unit ... 116!> \param mos ... 117!> \param dft_control ... 118! ************************************************************************************************** 119 SUBROUTINE qs_scf_initial_info(output_unit, mos, dft_control) 120 INTEGER :: output_unit 121 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 122 TYPE(dft_control_type), POINTER :: dft_control 123 124 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_initial_info', & 125 routineP = moduleN//':'//routineN 126 127 INTEGER :: homo, ispin, nao, nelectron_spin, nmo 128 129! print occupation numbers 130 131 IF (output_unit > 0) THEN 132 DO ispin = 1, dft_control%nspins 133 CALL get_mo_set(mo_set=mos(ispin)%mo_set, & 134 homo=homo, & 135 nelectron=nelectron_spin, & 136 nao=nao, & 137 nmo=nmo) 138 IF (dft_control%nspins > 1) THEN 139 WRITE (UNIT=output_unit, FMT="(/,T2,A,I2)") "Spin", ispin 140 END IF 141 WRITE (UNIT=output_unit, FMT="(/,(T2,A,T71,I10))") & 142 "Number of electrons:", nelectron_spin, & 143 "Number of occupied orbitals:", homo, & 144 "Number of molecular orbitals:", nmo 145 END DO 146 WRITE (UNIT=output_unit, FMT="(/,T2,A,T71,I10)") & 147 "Number of orbital functions:", nao 148 END IF 149 150 END SUBROUTINE qs_scf_initial_info 151 152! ************************************************************************************************** 153!> \brief writes out the mos in the scf loop if needed 154!> \param mos ... 155!> \param atomic_kind_set ... 156!> \param qs_kind_set ... 157!> \param particle_set ... 158!> \param dft_section ... 159! ************************************************************************************************** 160 SUBROUTINE qs_scf_write_mos(mos, atomic_kind_set, qs_kind_set, particle_set, dft_section) 161 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 162 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 163 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 164 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 165 TYPE(section_vals_type), POINTER :: dft_section 166 167 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_write_mos', & 168 routineP = moduleN//':'//routineN 169 170 IF (SIZE(mos) > 1) THEN 171 CALL write_mo_set(mos(1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 172 dft_section, spin="ALPHA", last=.FALSE.) 173 CALL write_mo_set(mos(2)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 174 dft_section, spin="BETA", last=.FALSE.) 175 ELSE 176 CALL write_mo_set(mos(1)%mo_set, atomic_kind_set, qs_kind_set, particle_set, 4, & 177 dft_section, last=.FALSE.) 178 END IF 179 180 END SUBROUTINE qs_scf_write_mos 181! ************************************************************************************************** 182!> \brief writes basic information obtained in a scf outer loop step 183!> \param output_unit ... 184!> \param scf_control ... 185!> \param scf_env ... 186!> \param energy ... 187!> \param total_steps ... 188!> \param should_stop ... 189!> \param outer_loop_converged ... 190! ************************************************************************************************** 191 192 SUBROUTINE qs_scf_outer_loop_info(output_unit, scf_control, scf_env, & 193 energy, total_steps, should_stop, outer_loop_converged) 194 INTEGER :: output_unit 195 TYPE(scf_control_type), POINTER :: scf_control 196 TYPE(qs_scf_env_type), POINTER :: scf_env 197 TYPE(qs_energy_type), POINTER :: energy 198 INTEGER :: total_steps 199 LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged 200 201 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_outer_loop_info', & 202 routineP = moduleN//':'//routineN 203 204 REAL(KIND=dp) :: outer_loop_eps 205 206 outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2)) 207 IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') & 208 "outer SCF iter = ", scf_env%outer_scf%iter_count, & 209 " RMS gradient = ", outer_loop_eps, " energy =", energy%total 210 211 IF (outer_loop_converged) THEN 212 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') & 213 "outer SCF loop converged in", scf_env%outer_scf%iter_count, & 214 " iterations or ", total_steps, " steps" 215 ELSE IF (scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf & 216 .OR. should_stop) THEN 217 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') & 218 "outer SCF loop FAILED to converge after ", & 219 scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps" 220 END IF 221 222 END SUBROUTINE qs_scf_outer_loop_info 223 224! ************************************************************************************************** 225!> \brief writes basic information obtained in a scf step 226!> \param scf_env ... 227!> \param output_unit ... 228!> \param just_energy ... 229!> \param t1 ... 230!> \param t2 ... 231!> \param energy ... 232! ************************************************************************************************** 233 234 SUBROUTINE qs_scf_loop_info(scf_env, output_unit, just_energy, t1, t2, energy) 235 236 TYPE(qs_scf_env_type), POINTER :: scf_env 237 INTEGER :: output_unit 238 LOGICAL :: just_energy 239 REAL(KIND=dp) :: t1, t2 240 TYPE(qs_energy_type), POINTER :: energy 241 242 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_info', & 243 routineP = moduleN//':'//routineN 244 245 IF ((output_unit > 0) .AND. scf_env%print_iter_line) THEN 246 IF (just_energy) THEN 247 WRITE (UNIT=output_unit, & 248 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,16X,F20.10)") & 249 scf_env%iter_count, TRIM(scf_env%iter_method), scf_env%iter_param, & 250 t2 - t1, energy%total 251 ELSE 252 IF ((ABS(scf_env%iter_delta) < 1.0E-8_dp) .OR. & 253 (ABS(scf_env%iter_delta) >= 1.0E5_dp)) THEN 254 WRITE (UNIT=output_unit, & 255 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,ES14.4,1X,F20.10,1X,ES9.2)") & 256 scf_env%iter_count, TRIM(scf_env%iter_method), scf_env%iter_param, & 257 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old 258 ELSE 259 WRITE (UNIT=output_unit, & 260 FMT="(T2,I5,1X,A,T20,E8.2,1X,F6.1,1X,F14.8,1X,F20.10,1X,ES9.2)") & 261 scf_env%iter_count, TRIM(scf_env%iter_method), scf_env%iter_param, & 262 t2 - t1, scf_env%iter_delta, energy%total, energy%total - energy%tot_old 263 END IF 264 END IF 265 END IF 266 267 END SUBROUTINE qs_scf_loop_info 268 269! ************************************************************************************************** 270!> \brief writes rather detailed summary of densities and energies 271!> after the SCF 272!> \param output_unit ... 273!> \param rho ... 274!> \param qs_charges ... 275!> \param energy ... 276!> \param nelectron_total ... 277!> \param dft_control ... 278!> \param qmmm ... 279!> \param qs_env ... 280!> \param gapw ... 281!> \param gapw_xc ... 282!> \par History 283!> 03.2006 created [Joost VandeVondele] 284! ************************************************************************************************** 285 SUBROUTINE qs_scf_print_scf_summary(output_unit, rho, qs_charges, energy, nelectron_total, & 286 dft_control, qmmm, qs_env, gapw, gapw_xc) 287 INTEGER, INTENT(IN) :: output_unit 288 TYPE(qs_rho_type), POINTER :: rho 289 TYPE(qs_charges_type), POINTER :: qs_charges 290 TYPE(qs_energy_type), POINTER :: energy 291 INTEGER, INTENT(IN) :: nelectron_total 292 TYPE(dft_control_type), POINTER :: dft_control 293 LOGICAL, INTENT(IN) :: qmmm 294 TYPE(qs_environment_type), POINTER :: qs_env 295 LOGICAL, INTENT(IN) :: gapw, gapw_xc 296 297 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_print_scf_summary', & 298 routineP = moduleN//':'//routineN 299 300 INTEGER :: bc, handle, ispin, psolver 301 REAL(kind=dp) :: exc_energy, implicit_ps_ehartree, & 302 tot1_h, tot1_s 303 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r 304 TYPE(pw_env_type), POINTER :: pw_env 305 306 NULLIFY (tot_rho_r, pw_env) 307 CALL timeset(routineN, handle) 308 309 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 310 psolver = pw_env%poisson_env%parameters%solver 311 312 IF (output_unit > 0) THEN 313 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r) 314 IF (.NOT. (dft_control%qs_control%semi_empirical .OR. & 315 dft_control%qs_control%xtb .OR. & 316 dft_control%qs_control%dftb)) THEN 317 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T41,2F20.10))") & 318 "Electronic density on regular grids: ", & 319 accurate_sum(tot_rho_r), & 320 accurate_sum(tot_rho_r) + nelectron_total, & 321 "Core density on regular grids:", & 322 qs_charges%total_rho_core_rspace, & 323 qs_charges%total_rho_core_rspace - REAL(nelectron_total + dft_control%charge, dp) 324 IF (gapw) THEN 325 tot1_h = qs_charges%total_rho1_hard(1) 326 tot1_s = qs_charges%total_rho1_soft(1) 327 DO ispin = 2, dft_control%nspins 328 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin) 329 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin) 330 END DO 331 WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") & 332 "Hard and soft densities (Lebedev):", & 333 tot1_h, tot1_s 334 WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") & 335 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", & 336 accurate_sum(tot_rho_r) + tot1_h - tot1_s, & 337 "Total charge density (r-space): ", & 338 accurate_sum(tot_rho_r) + tot1_h - tot1_s & 339 + qs_charges%total_rho_core_rspace, & 340 "Total Rho_soft + Rho0_soft (g-space):", & 341 qs_charges%total_rho_gspace 342 ELSE 343 WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") & 344 "Total charge density on r-space grids: ", & 345 accurate_sum(tot_rho_r) + & 346 qs_charges%total_rho_core_rspace, & 347 "Total charge density g-space grids: ", & 348 qs_charges%total_rho_gspace 349 END IF 350 END IF 351 IF (dft_control%qs_control%semi_empirical) THEN 352 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 353 "Core-core repulsion energy [eV]: ", energy%core_overlap*evolt, & 354 "Core Hamiltonian energy [eV]: ", energy%core*evolt, & 355 "Two-electron integral energy [eV]: ", energy%hartree*evolt, & 356 "Electronic energy [eV]: ", & 357 (energy%core + 0.5_dp*energy%hartree)*evolt 358 IF (energy%dispersion /= 0.0_dp) & 359 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 360 "Dispersion energy [eV]: ", energy%dispersion*evolt 361 ELSEIF (dft_control%qs_control%dftb) THEN 362 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 363 "Core Hamiltonian energy: ", energy%core, & 364 "Repulsive potential energy: ", energy%repulsive, & 365 "Electronic energy: ", energy%hartree, & 366 "Dispersion energy: ", energy%dispersion 367 IF (energy%dftb3 /= 0.0_dp) & 368 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 369 "DFTB3 3rd order energy: ", energy%dftb3 370 IF (energy%efield /= 0.0_dp) & 371 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 372 "Electric field interaction energy: ", energy%efield 373 ELSEIF (dft_control%qs_control%xtb) THEN 374 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 375 "Core Hamiltonian energy: ", energy%core, & 376 "Repulsive potential energy: ", energy%repulsive, & 377 "Electronic energy: ", energy%hartree, & 378 "DFTB3 3rd order energy: ", energy%dftb3, & 379 "Dispersion energy: ", energy%dispersion 380 IF (energy%efield /= 0.0_dp) & 381 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 382 "Electric field interaction energy: ", energy%efield 383 ELSE 384 IF (dft_control%do_admm) THEN 385 exc_energy = energy%exc + energy%exc_aux_fit 386 ELSE 387 exc_energy = energy%exc 388 END IF 389 390 IF (psolver .EQ. pw_poisson_implicit) THEN 391 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree 392 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 393 SELECT CASE (bc) 394 CASE (MIXED_PERIODIC_BC, MIXED_BC) 395 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 396 "Overlap energy of the core charge distribution:", energy%core_overlap, & 397 "Self energy of the core charge distribution: ", energy%core_self, & 398 "Core Hamiltonian energy: ", energy%core, & 399 "Hartree energy: ", implicit_ps_ehartree, & 400 "Electric enthalpy: ", energy%hartree, & 401 "Exchange-correlation energy: ", exc_energy 402 CASE (PERIODIC_BC, NEUMANN_BC) 403 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 404 "Overlap energy of the core charge distribution:", energy%core_overlap, & 405 "Self energy of the core charge distribution: ", energy%core_self, & 406 "Core Hamiltonian energy: ", energy%core, & 407 "Hartree energy: ", energy%hartree, & 408 "Exchange-correlation energy: ", exc_energy 409 END SELECT 410 ELSE 411 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 412 "Overlap energy of the core charge distribution:", energy%core_overlap, & 413 "Self energy of the core charge distribution: ", energy%core_self, & 414 "Core Hamiltonian energy: ", energy%core, & 415 "Hartree energy: ", energy%hartree, & 416 "Exchange-correlation energy: ", exc_energy 417 END IF 418 IF (energy%e_hartree /= 0.0_dp) & 419 WRITE (UNIT=output_unit, FMT="(T3,A,/,T3,A,T56,F25.14)") & 420 "Coulomb Electron-Electron Interaction Energy ", & 421 "- Already included in the total Hartree term ", energy%e_hartree 422 IF (energy%ex /= 0.0_dp) & 423 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 424 "Hartree-Fock Exchange energy: ", energy%ex 425 IF (energy%dispersion /= 0.0_dp) & 426 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 427 "Dispersion energy: ", energy%dispersion 428 IF (energy%gcp /= 0.0_dp) & 429 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 430 "gCP energy: ", energy%gcp 431 IF (gapw) THEN 432 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 433 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1, & 434 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c 435 END IF 436 IF (gapw_xc) THEN 437 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 438 "GAPW_XC| Exc from hard and soft atomic rho1: ", energy%exc1 439 END IF 440 END IF 441 IF (dft_control%smear) THEN 442 WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & 443 "Electronic entropic energy:", energy%kTS 444 WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") & 445 "Fermi energy:", energy%efermi 446 END IF 447 IF (dft_control%dft_plus_u) THEN 448 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 449 "DFT+U energy:", energy%dft_plus_u 450 END IF 451 IF (dft_control%do_sccs) THEN 452 WRITE (UNIT=output_unit, FMT="(/,T3,A,T56,F25.14)") & 453 "SCCS| Hartree energy of solute and solvent [Hartree]", energy%sccs_hartree 454 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14,/,T3,A,T61,F20.3)") & 455 "SCCS| Polarisation energy [Hartree]", energy%sccs_pol, & 456 "SCCS| [kcal/mol]", & 457 cp_unit_from_cp2k(energy%sccs_pol, "kcalmol") 458 END IF 459 IF (qmmm) THEN 460 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 461 "QM/MM Electrostatic energy: ", energy%qmmm_el 462 IF (qs_env%qmmm_env_qm%image_charge) THEN 463 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 464 "QM/MM image charge energy: ", energy%image_charge 465 ENDIF 466 END IF 467 IF (dft_control%qs_control%mulliken_restraint) THEN 468 WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") & 469 "Mulliken restraint energy: ", energy%mulliken 470 END IF 471 IF (dft_control%qs_control%semi_empirical) THEN 472 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 473 "Total energy [eV]: ", energy%total*evolt 474 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 475 "Atomic reference energy [eV]: ", energy%core_self*evolt, & 476 "Heat of formation [kcal/mol]: ", & 477 (energy%total + energy%core_self)*kcalmol 478 ELSE 479 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") & 480 "Total energy: ", energy%total 481 END IF 482 IF (qmmm) THEN 483 IF (qs_env%qmmm_env_qm%image_charge) THEN 484 CALL print_image_coefficients(qs_env%image_coeff, qs_env) 485 ENDIF 486 ENDIF 487 CALL m_flush(output_unit) 488 END IF 489 490 CALL timestop(handle) 491 492 END SUBROUTINE qs_scf_print_scf_summary 493 494! ************************************************************************************************** 495!> \brief collects the 'heavy duty' printing tasks out of the SCF loop 496!> \param qs_env ... 497!> \param scf_env ... 498!> \param para_env ... 499!> \par History 500!> 03.2006 created [Joost VandeVondele] 501! ************************************************************************************************** 502 SUBROUTINE qs_scf_loop_print(qs_env, scf_env, para_env) 503 TYPE(qs_environment_type), POINTER :: qs_env 504 TYPE(qs_scf_env_type), POINTER :: scf_env 505 TYPE(cp_para_env_type), POINTER :: para_env 506 507 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_loop_print', & 508 routineP = moduleN//':'//routineN 509 510 INTEGER :: after, handle, ic, ispin, iw 511 LOGICAL :: do_kpoints, omit_headers 512 REAL(KIND=dp) :: mo_mag_max, mo_mag_min, orthonormality 513 TYPE(cp_logger_type), POINTER :: logger 514 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s 515 TYPE(dft_control_type), POINTER :: dft_control 516 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 517 TYPE(qs_rho_type), POINTER :: rho 518 TYPE(section_vals_type), POINTER :: dft_section, input, scf_section 519 520 logger => cp_get_default_logger() 521 CALL timeset(routineN, handle) 522 523 CALL get_qs_env(qs_env=qs_env, input=input, dft_control=dft_control, & 524 do_kpoints=do_kpoints) 525 526 dft_section => section_vals_get_subs_vals(input, "DFT") 527 scf_section => section_vals_get_subs_vals(dft_section, "SCF") 528 529 CALL section_vals_val_get(input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 530 DO ispin = 1, dft_control%nspins 531 532 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 533 dft_section, "PRINT%AO_MATRICES/DENSITY"), cp_p_file)) THEN 534 CALL get_qs_env(qs_env, rho=rho) 535 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 536 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/DENSITY", & 537 extension=".Log") 538 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after) 539 after = MIN(MAX(after, 1), 16) 540 DO ic = 1, SIZE(matrix_p, 2) 541 CALL cp_dbcsr_write_sparse_matrix(matrix_p(ispin, ic)%matrix, 4, after, qs_env, para_env, & 542 output_unit=iw, omit_headers=omit_headers) 543 END DO 544 CALL cp_print_key_finished_output(iw, logger, dft_section, & 545 "PRINT%AO_MATRICES/DENSITY") 546 END IF 547 548 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 549 dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX"), cp_p_file)) THEN 550 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX", & 551 extension=".Log") 552 CALL section_vals_val_get(dft_section, "PRINT%AO_MATRICES%NDIGITS", i_val=after) 553 after = MIN(MAX(after, 1), 16) 554 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks) 555 DO ic = 1, SIZE(matrix_ks, 2) 556 IF (dft_control%qs_control%semi_empirical) THEN 557 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, & 558 scale=evolt, output_unit=iw, omit_headers=omit_headers) 559 ELSE 560 CALL cp_dbcsr_write_sparse_matrix(matrix_ks(ispin, ic)%matrix, 4, after, qs_env, para_env, & 561 output_unit=iw, omit_headers=omit_headers) 562 END IF 563 END DO 564 CALL cp_print_key_finished_output(iw, logger, dft_section, & 565 "PRINT%AO_MATRICES/KOHN_SHAM_MATRIX") 566 END IF 567 568 ENDDO 569 570 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 571 scf_section, "PRINT%MO_ORTHONORMALITY"), cp_p_file)) THEN 572 IF (do_kpoints) THEN 573 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", & 574 extension=".scfLog") 575 IF (iw > 0) THEN 576 WRITE (iw, '(T8,A)') & 577 " K-points: Maximum deviation from MO S-orthonormality not determined" 578 ENDIF 579 CALL cp_print_key_finished_output(iw, logger, scf_section, & 580 "PRINT%MO_ORTHONORMALITY") 581 ELSE 582 CALL get_qs_env(qs_env, mos=mos) 583 IF (scf_env%method == special_diag_method_nr) THEN 584 CALL calculate_orthonormality(orthonormality, mos) 585 ELSE 586 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s) 587 CALL calculate_orthonormality(orthonormality, mos, matrix_s(1, 1)%matrix) 588 END IF 589 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_ORTHONORMALITY", & 590 extension=".scfLog") 591 IF (iw > 0) THEN 592 WRITE (iw, '(T8,A,T61,E20.4)') & 593 " Maximum deviation from MO S-orthonormality", orthonormality 594 ENDIF 595 CALL cp_print_key_finished_output(iw, logger, scf_section, & 596 "PRINT%MO_ORTHONORMALITY") 597 END IF 598 ENDIF 599 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 600 scf_section, "PRINT%MO_MAGNITUDE"), cp_p_file)) THEN 601 IF (do_kpoints) THEN 602 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", & 603 extension=".scfLog") 604 IF (iw > 0) THEN 605 WRITE (iw, '(T8,A)') & 606 " K-points: Minimum/Maximum MO magnitude not determined" 607 ENDIF 608 CALL cp_print_key_finished_output(iw, logger, scf_section, & 609 "PRINT%MO_MAGNITUDE") 610 ELSE 611 CALL get_qs_env(qs_env, mos=mos) 612 CALL calculate_magnitude(mos, mo_mag_min, mo_mag_max) 613 iw = cp_print_key_unit_nr(logger, scf_section, "PRINT%MO_MAGNITUDE", & 614 extension=".scfLog") 615 IF (iw > 0) THEN 616 WRITE (iw, '(T8,A,T41,2E20.4)') & 617 " Minimum/Maximum MO magnitude ", mo_mag_min, mo_mag_max 618 ENDIF 619 CALL cp_print_key_finished_output(iw, logger, scf_section, & 620 "PRINT%MO_MAGNITUDE") 621 END IF 622 ENDIF 623 624 CALL timestop(handle) 625 626 END SUBROUTINE qs_scf_loop_print 627 628! ************************************************************************************************** 629!> \brief writes CDFT constraint information and optionally CDFT scf loop info 630!> \param output_unit where to write the information 631!> \param scf_control settings of the SCF loop 632!> \param scf_env the env which holds convergence data 633!> \param cdft_control the env which holds information about the constraint 634!> \param energy the total energy 635!> \param total_steps the total number of performed SCF iterations 636!> \param should_stop if the calculation should stop 637!> \param outer_loop_converged logical which determines if the CDFT SCF loop converged 638!> \param cdft_loop logical which determines a CDFT SCF loop is active 639!> \par History 640!> 12.2015 created [Nico Holmberg] 641! ************************************************************************************************** 642 SUBROUTINE qs_scf_cdft_info(output_unit, scf_control, scf_env, cdft_control, & 643 energy, total_steps, should_stop, outer_loop_converged, & 644 cdft_loop) 645 INTEGER :: output_unit 646 TYPE(scf_control_type), POINTER :: scf_control 647 TYPE(qs_scf_env_type), POINTER :: scf_env 648 TYPE(cdft_control_type), POINTER :: cdft_control 649 TYPE(qs_energy_type), POINTER :: energy 650 INTEGER :: total_steps 651 LOGICAL, INTENT(IN) :: should_stop, outer_loop_converged, & 652 cdft_loop 653 654 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_cdft_info', & 655 routineP = moduleN//':'//routineN 656 657 REAL(KIND=dp) :: outer_loop_eps 658 659 IF (cdft_loop) THEN 660 outer_loop_eps = SQRT(MAXVAL(scf_env%outer_scf%gradient(:, scf_env%outer_scf%iter_count)**2)) 661 IF (output_unit > 0) WRITE (output_unit, '(/,T3,A,I4,A,E10.2,A,F22.10)') & 662 "CDFT SCF iter = ", scf_env%outer_scf%iter_count, & 663 " RMS gradient = ", outer_loop_eps, " energy =", energy%total 664 IF (outer_loop_converged) THEN 665 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') & 666 "CDFT SCF loop converged in", scf_env%outer_scf%iter_count, & 667 " iterations or ", total_steps, " steps" 668 END IF 669 IF ((scf_env%outer_scf%iter_count > scf_control%outer_scf%max_scf .OR. should_stop) & 670 .AND. .NOT. outer_loop_converged) THEN 671 IF (output_unit > 0) WRITE (output_unit, '(T3,A,I4,A,I4,A,/)') & 672 "CDFT SCF loop FAILED to converge after ", & 673 scf_env%outer_scf%iter_count, " iterations or ", total_steps, " steps" 674 END IF 675 END IF 676 CALL qs_scf_cdft_constraint_info(output_unit, cdft_control) 677 678 END SUBROUTINE qs_scf_cdft_info 679 680! ************************************************************************************************** 681!> \brief writes information about the CDFT env 682!> \param output_unit where to write the information 683!> \param cdft_control the CDFT env that stores information about the constraint calculation 684!> \par History 685!> 12.2015 created [Nico Holmberg] 686! ************************************************************************************************** 687 SUBROUTINE qs_scf_cdft_initial_info(output_unit, cdft_control) 688 INTEGER :: output_unit 689 TYPE(cdft_control_type), POINTER :: cdft_control 690 691 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_cdft_initial_info', & 692 routineP = moduleN//':'//routineN 693 694 IF (output_unit > 0) THEN 695 WRITE (output_unit, '(/,A)') & 696 " ---------------------------------- CDFT --------------------------------------" 697 WRITE (output_unit, '(A)') & 698 " Optimizing a density constraint in an external SCF loop " 699 WRITE (output_unit, '(A)') " " 700 SELECT CASE (cdft_control%type) 701 CASE (outer_scf_hirshfeld_constraint) 702 WRITE (output_unit, '(A)') " Type of constraint: Hirshfeld" 703 CASE (outer_scf_becke_constraint) 704 WRITE (output_unit, '(A)') " Type of constraint: Becke" 705 END SELECT 706 WRITE (output_unit, '(A,I8)') " Number of constraints: ", SIZE(cdft_control%group) 707 WRITE (output_unit, '(A,L8)') " Using fragment densities:", cdft_control%fragment_density 708 WRITE (output_unit, '(A)') " " 709 IF (cdft_control%atomic_charges) WRITE (output_unit, '(A,/)') " Calculating atomic CDFT charges" 710 SELECT CASE (cdft_control%constraint_control%optimizer) 711 CASE (outer_scf_optimizer_sd) 712 WRITE (output_unit, '(A)') & 713 " Minimizer : SD : steepest descent" 714 CASE (outer_scf_optimizer_diis) 715 WRITE (output_unit, '(A)') & 716 " Minimizer : DIIS : direct inversion" 717 WRITE (output_unit, '(A)') & 718 " in the iterative subspace" 719 WRITE (output_unit, '(A,I3,A)') & 720 " using ", & 721 cdft_control%constraint_control%diis_buffer_length, " DIIS vectors" 722 CASE (outer_scf_optimizer_bisect) 723 WRITE (output_unit, '(A)') & 724 " Minimizer : BISECT : gradient bisection" 725 WRITE (output_unit, '(A,I3)') & 726 " using a trust count of", & 727 cdft_control%constraint_control%bisect_trust_count 728 CASE (outer_scf_optimizer_broyden, outer_scf_optimizer_newton, & 729 outer_scf_optimizer_newton_ls) 730 CALL cdft_opt_type_write(cdft_control%constraint_control%cdft_opt_control, & 731 cdft_control%constraint_control%optimizer, output_unit) 732 CASE (outer_scf_optimizer_secant) 733 WRITE (output_unit, '(A)') " Minimizer : Secant" 734 CASE DEFAULT 735 CPABORT("") 736 END SELECT 737 WRITE (output_unit, '(/,A,L7)') & 738 " Reusing OT preconditioner: ", cdft_control%reuse_precond 739 IF (cdft_control%reuse_precond) THEN 740 WRITE (output_unit, '(A,I3,A,I3,A)') & 741 " using old preconditioner for upto ", & 742 cdft_control%max_reuse, " subsequent CDFT SCF" 743 WRITE (output_unit, '(A,I3,A,I3,A)') & 744 " iterations if the relevant loop converged in less than ", & 745 cdft_control%precond_freq, " steps" 746 END IF 747 SELECT CASE (cdft_control%type) 748 CASE (outer_scf_hirshfeld_constraint) 749 WRITE (output_unit, '(/,A)') " Hirshfeld constraint settings" 750 WRITE (output_unit, '(A)') " " 751 SELECT CASE (cdft_control%hirshfeld_control%shape_function) 752 CASE (shape_function_gaussian) 753 WRITE (output_unit, '(A, A8)') & 754 " Shape function type: ", "Gaussian" 755 WRITE (output_unit, '(A)', ADVANCE='NO') & 756 " Type of Gaussian: " 757 SELECT CASE (cdft_control%hirshfeld_control%gaussian_shape) 758 CASE (radius_default) 759 WRITE (output_unit, '(A13)') "Default" 760 CASE (radius_covalent) 761 WRITE (output_unit, '(A13)') "Covalent" 762 CASE (radius_single) 763 WRITE (output_unit, '(A13)') "Fixed radius" 764 CASE (radius_vdw) 765 WRITE (output_unit, '(A13)') "Van der Waals" 766 CASE (radius_user) 767 WRITE (output_unit, '(A13)') "User-defined" 768 769 END SELECT 770 CASE (shape_function_density) 771 WRITE (output_unit, '(A, A8)') & 772 " Shape function type: ", "Density" 773 END SELECT 774 CASE (outer_scf_becke_constraint) 775 WRITE (output_unit, '(/, A)') " Becke constraint settings" 776 WRITE (output_unit, '(A)') " " 777 SELECT CASE (cdft_control%becke_control%cutoff_type) 778 CASE (becke_cutoff_global) 779 WRITE (output_unit, '(A,F8.3,A)') & 780 " Cutoff for partitioning :", cp_unit_from_cp2k(cdft_control%becke_control%rglobal, & 781 "angstrom"), " angstrom" 782 CASE (becke_cutoff_element) 783 WRITE (output_unit, '(A)') & 784 " Using element specific cutoffs for partitioning" 785 END SELECT 786 WRITE (output_unit, '(A,L7)') & 787 " Skipping distant gpoints: ", cdft_control%becke_control%should_skip 788 WRITE (output_unit, '(A,L7)') & 789 " Precompute gradients : ", cdft_control%becke_control%in_memory 790 WRITE (output_unit, '(A)') " " 791 IF (cdft_control%becke_control%adjust) & 792 WRITE (output_unit, '(A)') & 793 " Using atomic radii to generate a heteronuclear charge partitioning" 794 WRITE (output_unit, '(A)') " " 795 IF (.NOT. cdft_control%becke_control%cavity_confine) THEN 796 WRITE (output_unit, '(A)') & 797 " No confinement is active" 798 ELSE 799 WRITE (output_unit, '(A)') " Confinement using a Gaussian shaped cavity is active" 800 SELECT CASE (cdft_control%becke_control%cavity_shape) 801 CASE (radius_single) 802 WRITE (output_unit, '(A,F8.4, A)') & 803 " Type of Gaussian : Fixed radius: ", & 804 cp_unit_from_cp2k(cdft_control%becke_control%rcavity, "angstrom"), " angstrom" 805 CASE (radius_covalent) 806 WRITE (output_unit, '(A)') & 807 " Type of Gaussian : Covalent radius " 808 CASE (radius_vdw) 809 WRITE (output_unit, '(A)') & 810 " Type of Gaussian : vdW radius " 811 CASE (radius_user) 812 WRITE (output_unit, '(A)') & 813 " Type of Gaussian : User radius " 814 END SELECT 815 WRITE (output_unit, '(A,ES12.4)') & 816 " Cavity threshold : ", cdft_control%becke_control%eps_cavity 817 END IF 818 END SELECT 819 WRITE (output_unit, '(/,A)') & 820 " ---------------------------------- CDFT --------------------------------------" 821 END IF 822 823 END SUBROUTINE qs_scf_cdft_initial_info 824 825! ************************************************************************************************** 826!> \brief writes CDFT constraint information 827!> \param output_unit where to write the information 828!> \param cdft_control the env which holds information about the constraint 829!> \par History 830!> 08.2018 separated from qs_scf_cdft_info to make code callable elsewhere [Nico Holmberg] 831! ************************************************************************************************** 832 SUBROUTINE qs_scf_cdft_constraint_info(output_unit, cdft_control) 833 INTEGER :: output_unit 834 TYPE(cdft_control_type), POINTER :: cdft_control 835 836 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_scf_cdft_constraint_info', & 837 routineP = moduleN//':'//routineN 838 839 INTEGER :: igroup 840 841 IF (output_unit > 0) THEN 842 SELECT CASE (cdft_control%type) 843 CASE (outer_scf_hirshfeld_constraint) 844 WRITE (output_unit, '(/,T3,A,T60)') & 845 '------------------- Hirshfeld constraint information -------------------' 846 CASE (outer_scf_becke_constraint) 847 WRITE (output_unit, '(/,T3,A,T60)') & 848 '--------------------- Becke constraint information ---------------------' 849 CASE DEFAULT 850 CPABORT("Unknown CDFT constraint.") 851 END SELECT 852 DO igroup = 1, SIZE(cdft_control%target) 853 IF (igroup > 1) WRITE (output_unit, '(T3,A)') ' ' 854 WRITE (output_unit, '(T3,A,T54,(3X,I18))') & 855 'Atomic group :', igroup 856 SELECT CASE (cdft_control%group(igroup)%constraint_type) 857 CASE (cdft_charge_constraint) 858 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN 859 WRITE (output_unit, '(T3,A,T42,A)') & 860 'Type of constraint :', ADJUSTR('Charge density constraint (frag.)') 861 ELSE 862 WRITE (output_unit, '(T3,A,T50,A)') & 863 'Type of constraint :', ADJUSTR('Charge density constraint') 864 END IF 865 CASE (cdft_magnetization_constraint) 866 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN 867 WRITE (output_unit, '(T3,A,T35,A)') & 868 'Type of constraint :', ADJUSTR('Magnetization density constraint (frag.)') 869 ELSE 870 WRITE (output_unit, '(T3,A,T43,A)') & 871 'Type of constraint :', ADJUSTR('Magnetization density constraint') 872 END IF 873 CASE (cdft_alpha_constraint) 874 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN 875 WRITE (output_unit, '(T3,A,T38,A)') & 876 'Type of constraint :', ADJUSTR('Alpha spin density constraint (frag.)') 877 ELSE 878 WRITE (output_unit, '(T3,A,T46,A)') & 879 'Type of constraint :', ADJUSTR('Alpha spin density constraint') 880 END IF 881 CASE (cdft_beta_constraint) 882 IF (cdft_control%group(igroup)%is_fragment_constraint) THEN 883 WRITE (output_unit, '(T3,A,T39,A)') & 884 'Type of constraint :', ADJUSTR('Beta spin density constraint (frag.)') 885 ELSE 886 WRITE (output_unit, '(T3,A,T47,A)') & 887 'Type of constraint :', ADJUSTR('Beta spin density constraint') 888 END IF 889 CASE DEFAULT 890 CPABORT("Unknown constraint type.") 891 END SELECT 892 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') & 893 'Target value of constraint :', cdft_control%target(igroup) 894 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') & 895 'Current value of constraint :', cdft_control%value(igroup) 896 WRITE (output_unit, '(T3,A,T59,(3X,ES13.3))') & 897 'Deviation from target :', cdft_control%value(igroup) - cdft_control%target(igroup) 898 WRITE (output_unit, '(T3,A,T54,(3X,F18.12))') & 899 'Strength of constraint :', cdft_control%strength(igroup) 900 END DO 901 WRITE (output_unit, '(T3,A)') & 902 '------------------------------------------------------------------------' 903 END IF 904 905 END SUBROUTINE qs_scf_cdft_constraint_info 906 907END MODULE qs_scf_output 908