1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of the core Hamiltonian integral matrix <a|H|b> over 8!> Cartesian Gaussian-type functions. 9!> 10!> <a|H|b> = <a|T|b> + <a|V|b> 11!> 12!> Kinetic energy: 13!> 14!> <a|T|b> = <a|-nabla**2/2|b> 15!> \_______________/ 16!> | 17!> kinetic 18!> 19!> Nuclear potential energy: 20!> 21!> a) Allelectron calculation: 22!> 23!> erfc(r) 24!> <a|V|b> = -Z*<a|---------|b> 25!> r 26!> 27!> 1 - erf(r) 28!> = -Z*<a|------------|b> 29!> r 30!> 31!> 1 erf(r) 32!> = -Z*(<a|---|b> - <a|--------|b>) 33!> r r 34!> 35!> 1 36!> = -Z*(<a|---|b> - N*<ab||c>) 37!> r 38!> 39!> -Z 40!> = <a|---|b> + Z*N*<ab||c> 41!> r 42!> \_______/ \_____/ 43!> | | 44!> nuclear coulomb 45!> 46!> b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH): 47!> 48!> <a|V|b> = <a|(V(local) + V(non-local))|b> 49!> 50!> = <a|(V(local)|b> + <a|V(non-local))|b> 51!> 52!> <a|V(local)|b> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r + 53!> (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 + 54!> C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b> 55!> 56!> <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b> 57!> \par Literature 58!> S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996) 59!> C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998) 60!> M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000) 61!> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986) 62!> \par History 63!> - Joost VandeVondele (April 2003) : added LSD forces 64!> - Non-redundant calculation of the non-local part of the GTH PP 65!> (22.05.2003,MK) 66!> - New parallelization scheme (27.06.2003,MK) 67!> - OpenMP version (07.12.2003,JGH) 68!> - Binary search loop for VPPNL operators (09.01.2004,JGH,MK) 69!> - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH) 70!> - General refactoring (01.10.2010,JGH) 71!> - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH) 72!> - k-point functionality (07.2015,JGH) 73!> \author Matthias Krack (14.09.2000,21.03.02) 74! ************************************************************************************************** 75MODULE qs_core_hamiltonian 76 USE atomic_kind_types, ONLY: atomic_kind_type,& 77 get_atomic_kind_set 78 USE core_ae, ONLY: build_core_ae 79 USE core_ppl, ONLY: build_core_ppl 80 USE core_ppnl, ONLY: build_core_ppnl 81 USE cp_blacs_env, ONLY: cp_blacs_env_type 82 USE cp_control_types, ONLY: dft_control_type 83 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 84 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& 85 dbcsr_deallocate_matrix_set 86 USE cp_dbcsr_output, ONLY: cp_dbcsr_write_matrix_dist,& 87 cp_dbcsr_write_sparse_matrix 88 USE cp_log_handling, ONLY: cp_get_default_logger,& 89 cp_logger_type 90 USE cp_output_handling, ONLY: cp_p_file,& 91 cp_print_key_finished_output,& 92 cp_print_key_should_output,& 93 cp_print_key_unit_nr 94 USE cp_para_types, ONLY: cp_para_env_type 95 USE dbcsr_api, ONLY: & 96 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_iterator_blocks_left, & 97 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 98 dbcsr_p_type, dbcsr_type 99 USE input_constants, ONLY: do_admm_purify_none,& 100 do_ppl_analytic,& 101 kg_tnadd_atomic,& 102 rel_none,& 103 rel_trans_atom 104 USE input_section_types, ONLY: section_vals_val_get 105 USE kg_environment_types, ONLY: kg_environment_type 106 USE kg_tnadd_mat, ONLY: build_tnadd_mat 107 USE kinds, ONLY: dp 108 USE kpoint_types, ONLY: get_kpoint_info,& 109 kpoint_type 110 USE lri_environment_types, ONLY: lri_environment_type 111 USE particle_types, ONLY: particle_type 112 USE qs_condnum, ONLY: overlap_condnum 113 USE qs_environment_types, ONLY: get_qs_env,& 114 qs_environment_type,& 115 set_qs_env 116 USE qs_force_types, ONLY: qs_force_type 117 USE qs_kind_types, ONLY: get_qs_kind,& 118 qs_kind_type 119 USE qs_kinetic, ONLY: build_kinetic_matrix 120 USE qs_ks_types, ONLY: get_ks_env,& 121 qs_ks_env_type,& 122 set_ks_env 123 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 124 USE qs_oce_methods, ONLY: build_oce_matrices 125 USE qs_oce_types, ONLY: allocate_oce_set,& 126 create_oce_set,& 127 oce_matrix_type 128 USE qs_overlap, ONLY: build_overlap_matrix 129 USE qs_rho_types, ONLY: qs_rho_get,& 130 qs_rho_type 131 USE virial_types, ONLY: virial_type 132#include "./base/base_uses.f90" 133 134 IMPLICIT NONE 135 136 PRIVATE 137 138 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian' 139 140 PUBLIC :: build_core_hamiltonian_matrix 141 PUBLIC :: dump_info_core_hamiltonian 142 143CONTAINS 144 145! ************************************************************************************************** 146!> \brief Cosntruction of the QS Core Hamiltonian Matrix 147!> \param qs_env ... 148!> \param calculate_forces ... 149!> \author Creation (11.03.2002,MK) 150!> Non-redundant calculation of the non-local part of the GTH PP (22.05.2003,MK) 151!> New parallelization scheme (27.06.2003,MK) 152! ************************************************************************************************** 153 SUBROUTINE build_core_hamiltonian_matrix(qs_env, calculate_forces) 154 155 TYPE(qs_environment_type), POINTER :: qs_env 156 LOGICAL, INTENT(IN) :: calculate_forces 157 158 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix', & 159 routineP = moduleN//':'//routineN 160 161 INTEGER :: handle, ic, img, iw, nder, nders, & 162 nimages, nkind 163 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 164 LOGICAL :: norml1, norml2, ofdft, use_arnoldi, & 165 use_virial 166 REAL(KIND=dp) :: eps_filter, eps_fit 167 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 168 TYPE(cp_blacs_env_type), POINTER :: blacs_env 169 TYPE(cp_logger_type), POINTER :: logger 170 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 171 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_aux_fit_vs_orb 172 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s, matrix_t, & 173 matrix_w 174 TYPE(dft_control_type), POINTER :: dft_control 175 TYPE(kg_environment_type), POINTER :: kg_env 176 TYPE(kpoint_type), POINTER :: kpoints 177 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 178 POINTER :: sab_aux_fit, sab_aux_fit_vs_orb, & 179 sab_orb, sap_oce 180 TYPE(oce_matrix_type), POINTER :: oce 181 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 182 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 183 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 184 TYPE(qs_ks_env_type), POINTER :: ks_env 185 TYPE(qs_rho_type), POINTER :: rho 186 TYPE(virial_type), POINTER :: virial 187 188 IF (calculate_forces) THEN 189 CALL timeset(routineN//"_forces", handle) 190 ELSE 191 CALL timeset(routineN, handle) 192 ENDIF 193 194 NULLIFY (logger) 195 logger => cp_get_default_logger() 196 197 NULLIFY (dft_control) 198 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 199 200 ! is this a orbital-free method calculation 201 ofdft = dft_control%qs_control%ofgpw 202 203 nimages = dft_control%nimages 204 IF (ofdft) THEN 205 CPASSERT(nimages == 1) 206 END IF 207 208 nders = 0 209 IF (calculate_forces) THEN 210 nder = 1 211 ELSE 212 IF (cp_print_key_should_output(logger%iter_info, qs_env%input, & 213 "DFT%PRINT%AO_MATRICES/DERIVATIVES") /= 0) THEN 214 nder = 1 215 ELSE 216 nder = 0 217 END IF 218 END IF 219 220 IF ((cp_print_key_should_output(logger%iter_info, qs_env%input, & 221 "DFT%PRINT%AO_MATRICES/OVERLAP") /= 0 .AND. & 222 BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, & 223 "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file))) THEN 224 nders = 1 225 END IF 226 227 ! the delta pulse in the periodic case needs the momentum operator, 228 ! which is equivalent to the derivative of the overlap matrix 229 IF (ASSOCIATED(dft_control%rtp_control)) THEN 230 IF (dft_control%rtp_control%apply_delta_pulse .AND. & 231 dft_control%rtp_control%periodic) THEN 232 nders = 1 233 ENDIF 234 ENDIF 235 236 IF (dft_control%tddfpt2_control%enabled) THEN 237 nders = 1 238 IF (dft_control%do_admm) THEN 239 IF (dft_control%admm_control%purification_method /= do_admm_purify_none) & 240 CALL cp_abort(__LOCATION__, & 241 "Only purification method NONE is possible with TDDFT at the moment") 242 END IF 243 END IF 244 245 ! filter for new matrices 246 eps_filter = dft_control%qs_control%eps_filter_matrix 247 ! 248 NULLIFY (ks_env) 249 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) 250 NULLIFY (matrix_s, matrix_t) 251 CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s) 252 NULLIFY (sab_orb) 253 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb) 254 NULLIFY (rho, force, matrix_p, matrix_w) 255 IF (calculate_forces) THEN 256 CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w) 257 CALL get_qs_env(qs_env=qs_env, rho=rho) 258 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 259 ! *** If LSD, then combine alpha density and beta density to 260 ! *** total density: alpha <- alpha + beta and 261 ! *** spin density: beta <- alpha - beta 262 ! (since all things can be computed based on the sum of these matrices anyway) 263 ! (matrix_p is restored at the end of the run, matrix_w is left in its modified state 264 ! (as it should not be needed afterwards) 265 IF (SIZE(matrix_p, 1) == 2) THEN 266 DO img = 1, nimages 267 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 268 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 269 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, & 270 alpha_scalar=-2.0_dp, beta_scalar=1.0_dp) 271 CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, & 272 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 273 END DO 274 END IF 275 ! S matrix 276 CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, & 277 matrix_name="OVERLAP MATRIX", & 278 basis_type_a="ORB", & 279 basis_type_b="ORB", & 280 sab_nl=sab_orb, calculate_forces=.TRUE., & 281 matrixkp_p=matrix_w) 282 ! T matrix 283 IF (.NOT. ofdft) & 284 CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, & 285 matrix_name="KINETIC ENERGY MATRIX", & 286 basis_type="ORB", & 287 sab_nl=sab_orb, calculate_forces=.TRUE., & 288 matrixkp_p=matrix_p, & 289 eps_filter=eps_filter) 290 IF (calculate_forces) THEN 291 ! *** If LSD, then recover alpha density and beta density *** 292 ! *** from the total density (1) and the spin density (2) *** 293 ! *** The W matrix is neglected, since it will be destroyed *** 294 ! *** in the calling force routine after leaving this routine *** 295 IF (SIZE(matrix_p, 1) == 2) THEN 296 DO img = 1, nimages 297 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, & 298 alpha_scalar=0.5_dp, beta_scalar=0.5_dp) 299 CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, & 300 alpha_scalar=-1.0_dp, beta_scalar=1.0_dp) 301 END DO 302 END IF 303 END IF 304 ELSE 305 ! S matrix 306 CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, & 307 matrix_name="OVERLAP MATRIX", & 308 basis_type_a="ORB", & 309 basis_type_b="ORB", & 310 sab_nl=sab_orb) 311 ! T matrix 312 IF (.NOT. ofdft) & 313 CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, & 314 matrix_name="KINETIC ENERGY MATRIX", & 315 basis_type="ORB", & 316 sab_nl=sab_orb, & 317 eps_filter=eps_filter) 318 END IF 319 320 ! ADMM 321 IF (.NOT. calculate_forces) THEN 322 IF (dft_control%do_admm) THEN 323 NULLIFY (matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, & 324 sab_aux_fit, sab_aux_fit_vs_orb) 325 CALL get_qs_env(qs_env=qs_env, matrix_s_aux_fit=matrix_s_aux_fit, & 326 sab_aux_fit=sab_aux_fit) 327 CALL build_overlap_matrix(ks_env, matrix_s=matrix_s_aux_fit, & 328 matrix_name="AUX_FIT_OVERLAP", & 329 basis_type_a="AUX_FIT", & 330 basis_type_b="AUX_FIT", & 331 sab_nl=sab_aux_fit) 332 CALL set_ks_env(ks_env, matrix_s_aux_fit=matrix_s_aux_fit) 333 CALL get_qs_env(qs_env=qs_env, matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, & 334 sab_aux_fit_vs_orb=sab_aux_fit_vs_orb) 335 CALL build_overlap_matrix(ks_env, matrix_s=matrix_s_aux_fit_vs_orb, & 336 matrix_name="MIXED_OVERLAP", & 337 basis_type_a="AUX_FIT", & 338 basis_type_b="ORB", & 339 sab_nl=sab_aux_fit_vs_orb) 340 CALL set_ks_env(ks_env, matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb) 341 END IF 342 END IF 343 344 ! initialize H matrix 345 NULLIFY (matrix_h) 346 CALL get_qs_env(qs_env=qs_env, matrix_h_kp=matrix_h) 347 CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages) 348 DO img = 1, nimages 349 ALLOCATE (matrix_h(1, img)%matrix) 350 CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix) 351 CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb) 352 IF (.NOT. ofdft) THEN 353 CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, & 354 keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX") 355 END IF 356 END DO 357 358 NULLIFY (qs_kind_set, atomic_kind_set, particle_set) 359 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & 360 particle_set=particle_set) 361 362 IF (.NOT. ofdft) THEN 363 ! relativistic atomic correction to kinetic energy 364 IF (qs_env%rel_control%rel_method /= rel_none) THEN 365 IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN 366 IF (nimages == 1) THEN 367 ic = 1 368 ELSE 369 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 370 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 371 ic = cell_to_index(0, 0, 0) 372 END IF 373 CALL build_atomic_relmat(matrix_h(1, 1)%matrix, & 374 atomic_kind_set, qs_kind_set, particle_set) 375 ELSE 376 CPABORT("Relativistic corrections of this type are currently not implemented") 377 END IF 378 END IF 379 END IF 380 381 ! *** core and pseudopotentials 382 CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder) 383 384 ! *** GAPW one-center-expansion (oce) matrices 385 NULLIFY (sap_oce) 386 CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce) 387 NULLIFY (oce) 388 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN 389 CALL get_qs_env(qs_env=qs_env, oce=oce) 390 CALL create_oce_set(oce) 391 nkind = SIZE(atomic_kind_set) 392 CALL allocate_oce_set(oce, nkind) 393 eps_fit = dft_control%qs_control%gapw_control%eps_fit 394 IF (ASSOCIATED(sap_oce)) & 395 CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, & 396 sap_oce, eps_fit) 397 END IF 398 399 ! *** KG atomic potentials for nonadditive kinetic energy 400 IF (dft_control%qs_control%do_kg) THEN 401 IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN 402 CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, virial=virial, dbcsr_dist=dbcsr_dist) 403 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 404 CALL build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, & 405 qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist) 406 END IF 407 END IF 408 409 ! *** Put the core Hamiltonian matrix in the QS environment *** 410 CALL set_qs_env(qs_env, oce=oce) 411 CALL set_ks_env(ks_env, matrix_s_kp=matrix_s, kinetic_kp=matrix_t, matrix_h_kp=matrix_h) 412 413 ! *** Print matrices if requested 414 CALL dump_info_core_hamiltonian(qs_env, calculate_forces) 415 416 ! *** Overlap condition number 417 IF (.NOT. calculate_forces) THEN 418 IF (cp_print_key_should_output(logger%iter_info, qs_env%input, & 419 "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN 420 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", & 421 extension=".Log") 422 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1) 423 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2) 424 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi) 425 CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env) 426 CALL overlap_condnum(matrix_s, iw, norml1, norml2, use_arnoldi, blacs_env) 427 END IF 428 END IF 429 430 CALL timestop(handle) 431 432 END SUBROUTINE build_core_hamiltonian_matrix 433 434! ************************************************************************************************** 435!> \brief ... 436!> \param qs_env ... 437!> \param matrix_h ... 438!> \param matrix_p ... 439!> \param calculate_forces ... 440!> \param nder ... 441! ************************************************************************************************** 442 SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder) 443 444 TYPE(qs_environment_type), POINTER :: qs_env 445 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p 446 LOGICAL, INTENT(IN) :: calculate_forces 447 INTEGER, INTENT(IN) :: nder 448 449 CHARACTER(LEN=*), PARAMETER :: routineN = 'core_matrices', routineP = moduleN//':'//routineN 450 451 INTEGER :: nimages 452 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 453 LOGICAL :: all_present, ppl_present, ppnl_present, & 454 use_virial 455 REAL(KIND=dp) :: eps_ppnl 456 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 457 TYPE(dft_control_type), POINTER :: dft_control 458 TYPE(kpoint_type), POINTER :: kpoints 459 TYPE(lri_environment_type), POINTER :: lri_env 460 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 461 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl 462 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 463 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 464 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 465 TYPE(qs_ks_env_type), POINTER :: ks_env 466 TYPE(virial_type), POINTER :: virial 467 468 NULLIFY (dft_control) 469 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control) 470 nimages = dft_control%nimages 471 ! prepare for k-points 472 NULLIFY (cell_to_index) 473 IF (nimages > 1) THEN 474 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 475 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 476 dft_control%qs_control%do_ppl_method = do_ppl_analytic 477 END IF 478 ! force analytic ppl calculation for GAPW methods 479 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN 480 dft_control%qs_control%do_ppl_method = do_ppl_analytic 481 ENDIF 482 483 ! force 484 NULLIFY (force) 485 IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force) 486 ! virial 487 CALL get_qs_env(qs_env=qs_env, virial=virial) 488 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 489 490 NULLIFY (qs_kind_set, atomic_kind_set, particle_set) 491 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, & 492 particle_set=particle_set) 493 494 NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl) 495 CALL get_qs_env(qs_env=qs_env, & 496 sab_orb=sab_orb, & 497 sac_ae=sac_ae, & 498 sac_ppl=sac_ppl, & 499 sap_ppnl=sap_ppnl) 500 501 ! *** compute the nuclear attraction contribution to the core hamiltonian *** 502 all_present = ASSOCIATED(sac_ae) 503 IF (all_present) THEN 504 CALL build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & 505 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index) 506 END IF 507 508 ! *** compute the ppl contribution to the core hamiltonian *** 509 ppl_present = ASSOCIATED(sac_ppl) 510 IF (ppl_present) THEN 511 IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN 512 IF (dft_control%qs_control%lrigpw) THEN 513 CALL get_qs_env(qs_env, lri_env=lri_env) 514 IF (lri_env%ppl_ri) THEN 515 IF (lri_env%exact_1c_terms) THEN 516 CPABORT("not implemented") 517 END IF 518 ELSE 519 CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & 520 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & 521 nimages, cell_to_index, "ORB") 522 END IF 523 ELSE 524 CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & 525 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & 526 nimages, cell_to_index, "ORB") 527 END IF 528 END IF 529 END IF 530 531 ! *** compute the ppnl contribution to the core hamiltonian *** 532 eps_ppnl = dft_control%qs_control%eps_ppnl 533 ppnl_present = ASSOCIATED(sap_ppnl) 534 IF (ppnl_present) THEN 535 CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, & 536 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, & 537 nimages, cell_to_index, "ORB") 538 END IF 539 540 END SUBROUTINE core_matrices 541 542! ************************************************************************************************** 543!> \brief Adds atomic blocks of relativistic correction for the kinetic energy 544!> \param matrix_h ... 545!> \param atomic_kind_set ... 546!> \param qs_kind_set ... 547!> \param particle_set ... 548! ************************************************************************************************** 549 SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set, particle_set) 550 TYPE(dbcsr_type), POINTER :: matrix_h 551 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 552 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 553 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 554 555 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_atomic_relmat', & 556 routineP = moduleN//':'//routineN 557 558 INTEGER :: blk, iatom, ikind, jatom, natom 559 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of 560 REAL(KIND=dp), DIMENSION(:, :), POINTER :: hblock, reltmat 561 TYPE(dbcsr_iterator_type) :: iter 562 563 natom = SIZE(particle_set) 564 ALLOCATE (kind_of(natom)) 565 566 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of) 567 568 CALL dbcsr_iterator_start(iter, matrix_h) 569 DO WHILE (dbcsr_iterator_blocks_left(iter)) 570 CALL dbcsr_iterator_next_block(iter, iatom, jatom, hblock, blk) 571 IF (iatom == jatom) THEN 572 ikind = kind_of(iatom) 573 CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat) 574 IF (ASSOCIATED(reltmat)) hblock = hblock + reltmat 575 END IF 576 END DO 577 CALL dbcsr_iterator_stop(iter) 578 579 DEALLOCATE (kind_of) 580 581 END SUBROUTINE build_atomic_relmat 582 583! ************************************************************************************************** 584!> \brief Possibly prints matrices after the construction of the Core 585!> Hamiltonian Matrix 586!> \param qs_env ... 587!> \param calculate_forces ... 588! ************************************************************************************************** 589 SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces) 590 TYPE(qs_environment_type), POINTER :: qs_env 591 LOGICAL, INTENT(IN) :: calculate_forces 592 593 CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian', & 594 routineP = moduleN//':'//routineN 595 596 INTEGER :: after, handle, i, ic, iw, output_unit 597 LOGICAL :: omit_headers 598 TYPE(cp_logger_type), POINTER :: logger 599 TYPE(cp_para_env_type), POINTER :: para_env 600 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_v 601 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_s, matrixkp_t 602 603 CALL timeset(routineN, handle) 604 605 NULLIFY (logger, matrix_v, matrix_s, para_env) 606 logger => cp_get_default_logger() 607 CALL get_qs_env(qs_env, para_env=para_env) 608 609 ! Print the distribution of the overlap matrix blocks 610 ! this duplicates causes duplicate printing at the force calc 611 IF (.NOT. calculate_forces) THEN 612 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 613 qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN 614 output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", & 615 extension=".distribution") 616 CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s) 617 CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env) 618 CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION") 619 END IF 620 END IF 621 622 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers) 623 ! Print the overlap integral matrix, if requested 624 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 625 qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN 626 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", & 627 extension=".Log") 628 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 629 after = MIN(MAX(after, 1), 16) 630 CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s) 631 IF (ASSOCIATED(matrixkp_s)) THEN 632 DO ic = 1, SIZE(matrixkp_s, 2) 633 CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, & 634 output_unit=iw, omit_headers=omit_headers) 635 END DO 636 IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, & 637 "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file) & 638 .AND. ASSOCIATED(matrix_s)) THEN 639 DO ic = 1, SIZE(matrixkp_s, 2) 640 DO i = 2, SIZE(matrix_s) 641 CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, & 642 output_unit=iw, omit_headers=omit_headers) 643 END DO 644 END DO 645 END IF 646 END IF 647 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 648 "DFT%PRINT%AO_MATRICES/OVERLAP") 649 END IF 650 651 ! Print the kinetic energy integral matrix, if requested 652 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 653 qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN 654 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", & 655 extension=".Log") 656 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 657 after = MIN(MAX(after, 1), 16) 658 CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t) 659 IF (ASSOCIATED(matrixkp_t)) THEN 660 DO ic = 1, SIZE(matrixkp_t, 2) 661 CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, & 662 output_unit=iw, omit_headers=omit_headers) 663 END DO 664 END IF 665 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 666 "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY") 667 END IF 668 669 ! Print the potential energy matrix, if requested 670 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 671 qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN 672 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", & 673 extension=".Log") 674 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 675 after = MIN(MAX(after, 1), 16) 676 CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t) 677 IF (ASSOCIATED(matrixkp_h)) THEN 678 IF (SIZE(matrixkp_h, 2) == 1) THEN 679 CALL dbcsr_allocate_matrix_set(matrix_v, 1) 680 ALLOCATE (matrix_v(1)%matrix) 681 CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX") 682 CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, & 683 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) 684 CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, & 685 para_env, output_unit=iw, omit_headers=omit_headers) 686 CALL dbcsr_deallocate_matrix_set(matrix_v) 687 ELSE 688 CPWARN("Printing of potential energy matrix not implemented for k-points") 689 END IF 690 END IF 691 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 692 "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY") 693 END IF 694 695 ! Print the core Hamiltonian matrix, if requested 696 IF (BTEST(cp_print_key_should_output(logger%iter_info, & 697 qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN 698 iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", & 699 extension=".Log") 700 CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after) 701 after = MIN(MAX(after, 1), 16) 702 CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h) 703 IF (ASSOCIATED(matrixkp_h)) THEN 704 DO ic = 1, SIZE(matrixkp_h, 2) 705 CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, & 706 output_unit=iw, omit_headers=omit_headers) 707 END DO 708 END IF 709 CALL cp_print_key_finished_output(iw, logger, qs_env%input, & 710 "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN") 711 END IF 712 713 CALL timestop(handle) 714 715 END SUBROUTINE dump_info_core_hamiltonian 716 717END MODULE qs_core_hamiltonian 718