1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utilities for hfx and admm methods 8!> 9!> 10!> \par History 11!> refactoring 03-2011 [MI] 12!> \author MI 13! ************************************************************************************************** 14MODULE hfx_admm_utils 15 USE admm_dm_types, ONLY: admm_dm_create,& 16 admm_dm_type 17 USE admm_methods, ONLY: scale_dm 18 USE admm_types, ONLY: admm_env_create,& 19 admm_type 20 USE atomic_kind_types, ONLY: atomic_kind_type 21 USE cell_types, ONLY: cell_type 22 USE cp_control_types, ONLY: dft_control_type 23 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm 24 USE cp_fm_types, ONLY: cp_fm_get_info,& 25 cp_fm_type 26 USE cp_log_handling, ONLY: cp_get_default_logger,& 27 cp_logger_type 28 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 29 cp_print_key_unit_nr 30 USE cp_para_types, ONLY: cp_para_env_type 31 USE dbcsr_api, ONLY: dbcsr_add,& 32 dbcsr_p_type,& 33 dbcsr_set,& 34 dbcsr_type 35 USE hfx_derivatives, ONLY: derivatives_four_center 36 USE hfx_energy_potential, ONLY: integrate_four_center 37 USE hfx_types, ONLY: hfx_type 38 USE input_constants, ONLY: & 39 do_admm_aux_exch_func_bee, do_admm_aux_exch_func_default, do_admm_aux_exch_func_none, & 40 do_admm_aux_exch_func_opt, do_admm_aux_exch_func_pbex, do_potential_coulomb, & 41 do_potential_short, do_potential_truncated, xc_funct_no_shortcut 42 USE input_section_types, ONLY: section_vals_duplicate,& 43 section_vals_get,& 44 section_vals_get_subs_vals,& 45 section_vals_get_subs_vals2,& 46 section_vals_remove_values,& 47 section_vals_type,& 48 section_vals_val_get,& 49 section_vals_val_set 50 USE kinds, ONLY: dp 51 USE particle_types, ONLY: particle_type 52 USE pw_env_types, ONLY: pw_env_get,& 53 pw_env_type 54 USE pw_methods, ONLY: pw_transfer 55 USE pw_poisson_methods, ONLY: pw_poisson_solve 56 USE pw_poisson_types, ONLY: pw_poisson_type 57 USE pw_pool_types, ONLY: pw_pool_create_pw,& 58 pw_pool_give_back_pw,& 59 pw_pool_type 60 USE pw_types, ONLY: COMPLEXDATA1D,& 61 REALDATA3D,& 62 REALSPACE,& 63 RECIPROCALSPACE,& 64 pw_create,& 65 pw_p_type,& 66 pw_release 67 USE qs_collocate_density, ONLY: calculate_wavefunction 68 USE qs_energy_types, ONLY: qs_energy_type 69 USE qs_environment_types, ONLY: get_qs_env,& 70 qs_environment_type,& 71 set_qs_env 72 USE qs_force_types, ONLY: qs_force_type 73 USE qs_kind_types, ONLY: qs_kind_type 74 USE qs_ks_types, ONLY: qs_ks_env_type,& 75 set_ks_env 76 USE qs_mo_types, ONLY: get_mo_set,& 77 mo_set_p_type 78 USE qs_rho_types, ONLY: qs_rho_get,& 79 qs_rho_type 80 USE rt_propagation_types, ONLY: rt_prop_type 81 USE virial_types, ONLY: virial_type 82 USE xc_adiabatic_utils, ONLY: rescale_xc_potential 83#include "./base/base_uses.f90" 84 85 IMPLICIT NONE 86 87 PRIVATE 88 89 ! *** Public subroutines *** 90 PUBLIC :: hfx_ks_matrix, hfx_admm_init, create_admm_xc_section, tddft_hfx_matrix 91 92 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_admm_utils' 93 94CONTAINS 95 96! ************************************************************************************************** 97!> \brief ... 98!> \param qs_env ... 99!> \param do_mp2_hfx ... 100! ************************************************************************************************** 101 SUBROUTINE hfx_admm_init(qs_env, do_mp2_hfx) 102 103 TYPE(qs_environment_type), POINTER :: qs_env 104 LOGICAL, INTENT(in), OPTIONAL :: do_mp2_hfx 105 106 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_admm_init', routineP = moduleN//':'//routineN 107 108 INTEGER :: handle, n_rep_hf, natoms, nspins 109 LOGICAL :: do_hfx, my_do_mp2_hfx, s_mstruct_changed 110 TYPE(admm_dm_type), POINTER :: admm_dm 111 TYPE(admm_type), POINTER :: admm_env 112 TYPE(cp_para_env_type), POINTER :: para_env 113 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s_aux_fit, matrix_s_aux_fit_vs_orb 114 TYPE(dft_control_type), POINTER :: dft_control 115 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos, mos_aux_fit 116 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 117 TYPE(qs_ks_env_type), POINTER :: ks_env 118 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit 119 TYPE(section_vals_type), POINTER :: hfx_sections, input, xc_section 120 121 my_do_mp2_hfx = .FALSE. 122 IF (PRESENT(do_mp2_hfx)) my_do_mp2_hfx = do_mp2_hfx 123 124 CALL timeset(routineN, handle) 125 NULLIFY (admm_env, admm_dm, hfx_sections, mos, mos_aux_fit, para_env, & 126 particle_set, xc_section, matrix_s_aux_fit, & 127 matrix_s_aux_fit_vs_orb, rho, rho_aux_fit, ks_env, dft_control, & 128 input) 129 130 CALL get_qs_env(qs_env, & 131 mos_aux_fit=mos_aux_fit, & 132 mos=mos, & 133 admm_env=admm_env, & 134 admm_dm=admm_dm, & 135 matrix_s_aux_fit=matrix_s_aux_fit, & 136 matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, & 137 para_env=para_env, & 138 s_mstruct_changed=s_mstruct_changed, & 139 rho=rho, & 140 rho_aux_fit=rho_aux_fit, & 141 ks_env=ks_env, & 142 dft_control=dft_control, & 143 input=input) 144 145 nspins = dft_control%nspins 146 147 IF (my_do_mp2_hfx) THEN 148 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF") 149 ELSE 150 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") 151 END IF 152 153 CALL section_vals_get(hfx_sections, explicit=do_hfx) 154 155 !! ** ADMM can only be used with HFX 156 IF (.NOT. do_hfx) & 157 CPABORT("Wavefunction fitting requested without Hartree-Fock.") 158 159 ! ** Method runs with GAPW only if no DFT exchange correction 160 IF (dft_control%qs_control%gapw .AND. & 161 dft_control%admm_control%aux_exch_func .NE. do_admm_aux_exch_func_none) THEN 162 CPABORT("ADMM only works with GAPW without DFT exchange correction") 163 END IF 164 165 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) 166 IF (n_rep_hf > 1) & 167 CPABORT("ADMM can handle only one HF section.") 168 169 IF (.NOT. ASSOCIATED(admm_env)) THEN 170 ! setup admm environment 171 CALL get_qs_env(qs_env, input=input, particle_set=particle_set) 172 natoms = SIZE(particle_set, 1) 173 CALL admm_env_create(admm_env, dft_control%admm_control, mos, mos_aux_fit, & 174 para_env, natoms) 175 CALL set_qs_env(qs_env, admm_env=admm_env) 176 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 177 CALL create_admm_xc_section(qs_env=qs_env, xc_section=xc_section, & 178 admm_env=admm_env) 179 ENDIF 180 181 IF (dft_control%do_admm_dm .AND. .NOT. ASSOCIATED(admm_dm)) THEN 182 CALL admm_dm_create(admm_dm, dft_control%admm_control, nspins=nspins, natoms=natoms) 183 CALL set_ks_env(ks_env, admm_dm=admm_dm) 184 ENDIF 185 186 CALL timestop(handle) 187 188 END SUBROUTINE hfx_admm_init 189 190! ************************************************************************************************** 191!> \brief Add the hfx contributions to the Hamiltonian 192!> 193!> \param qs_env ... 194!> \param matrix_ks ... 195!> \param rho ... 196!> \param energy ... 197!> \param calculate_forces ... 198!> \param just_energy ... 199!> \param v_rspace_new ... 200!> \param v_tau_rspace ... 201!> \par History 202!> refactoring 03-2011 [MI] 203! ************************************************************************************************** 204 205 SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, & 206 just_energy, v_rspace_new, v_tau_rspace) 207 208 TYPE(qs_environment_type), POINTER :: qs_env 209 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks 210 TYPE(qs_rho_type), POINTER :: rho 211 TYPE(qs_energy_type), POINTER :: energy 212 LOGICAL, INTENT(in) :: calculate_forces, just_energy 213 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace 214 215 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ks_matrix', routineP = moduleN//':'//routineN 216 217 INTEGER :: handle, ikind, img, irep, ispin, mspin, & 218 n_rep_hf, nimages, ns, nspins 219 LOGICAL :: distribute_fock_matrix, & 220 do_adiabatic_rescaling, & 221 hfx_treat_lsd_in_core, & 222 s_mstruct_changed, use_virial 223 REAL(dp) :: eh1, ehfx, ehfxrt 224 REAL(dp), ALLOCATABLE, DIMENSION(:) :: hf_energy 225 TYPE(cp_para_env_type), POINTER :: para_env 226 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_1d, matrix_ks_aux_fit, & 227 matrix_ks_aux_fit_hfx, & 228 matrix_ks_aux_fit_im, matrix_ks_im, & 229 rho_ao_1d 230 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks_orb, rho_ao_orb 231 TYPE(dft_control_type), POINTER :: dft_control 232 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data 233 TYPE(pw_env_type), POINTER :: pw_env 234 TYPE(pw_poisson_type), POINTER :: poisson_env 235 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 236 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 237 TYPE(qs_rho_type), POINTER :: rho_orb 238 TYPE(rt_prop_type), POINTER :: rtp 239 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, & 240 hfx_sections, input 241 TYPE(virial_type), POINTER :: virial 242 243 CALL timeset(routineN, handle) 244 245 NULLIFY (auxbas_pw_pool, dft_control, force, hfx_sections, input, & 246 para_env, poisson_env, pw_env, virial, & 247 matrix_ks_im, matrix_ks_orb, rho_ao_orb, & 248 matrix_h, matrix_ks_aux_fit, matrix_ks_aux_fit_im, matrix_ks_aux_fit_hfx) 249 250 CALL get_qs_env(qs_env=qs_env, & 251 dft_control=dft_control, & 252 input=input, & 253 matrix_h_kp=matrix_h, & 254 para_env=para_env, & 255 pw_env=pw_env, & 256 virial=virial, & 257 matrix_ks_im=matrix_ks_im, & 258 matrix_ks_aux_fit=matrix_ks_aux_fit, & 259 matrix_ks_aux_fit_im=matrix_ks_aux_fit_im, & 260 matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx, & 261 s_mstruct_changed=s_mstruct_changed, & 262 x_data=x_data) 263 264 nspins = dft_control%nspins 265 nimages = dft_control%nimages 266 267 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 268 269 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") 270 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) 271 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, & 272 i_rep_section=1) 273 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING") 274 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) 275 276 ! *** Initialize the auxiliary ks matrix to zero if required 277 IF (dft_control%do_admm) THEN 278 DO ispin = 1, nspins 279 CALL dbcsr_set(matrix_ks_aux_fit(ispin)%matrix, 0.0_dp) 280 CALL dbcsr_set(matrix_ks_aux_fit_hfx(ispin)%matrix, 0.0_dp) 281 END DO 282 END IF 283 DO ispin = 1, nspins 284 DO img = 1, nimages 285 CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp) 286 END DO 287 END DO 288 289 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) 290 291 IF (calculate_forces) THEN 292 !! initalize force array to zero 293 CALL get_qs_env(qs_env=qs_env, force=force) 294 DO ikind = 1, SIZE(force) 295 force(ikind)%fock_4c(:, :) = 0.0_dp 296 END DO 297 END IF 298 ALLOCATE (hf_energy(n_rep_hf)) 299 300 DO irep = 1, n_rep_hf 301 302 IF (do_adiabatic_rescaling .AND. hfx_treat_lsd_in_core) & 303 CPABORT("HFX_TREAT_LSD_IN_CORE not implemented for adiabatically rescaled hybrids") 304 ! everything is calulated with adiabatic rescaling but the potential is not added in a first step 305 distribute_fock_matrix = .NOT. do_adiabatic_rescaling 306 307 mspin = 1 308 IF (hfx_treat_lsd_in_core) mspin = nspins 309 310 ! fetch the correct matrices for normal HFX or ADMM 311 IF (dft_control%do_admm) THEN 312 CALL get_qs_env(qs_env=qs_env, matrix_ks_aux_fit=matrix_ks_1d, & 313 rho_aux_fit=rho_orb) 314 ns = SIZE(matrix_ks_1d) 315 matrix_ks_orb(1:ns, 1:1) => matrix_ks_1d(1:ns) 316 ELSE 317 CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_orb, rho=rho_orb) 318 END IF 319 CALL qs_rho_get(rho_struct=rho_orb, rho_ao_kp=rho_ao_orb) 320 ! Finally the real hfx calulation 321 ehfx = 0.0_dp 322 DO ispin = 1, mspin 323 CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, & 324 para_env, s_mstruct_changed, irep, distribute_fock_matrix, & 325 ispin=ispin) 326 ehfx = ehfx + eh1 327 END DO 328 329 IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN 330 !Scale auxiliary density matrix for ADMMP (see Merlot2014) with gsi(ispin) to scale force 331 IF (dft_control%do_admm) THEN 332 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.FALSE.) 333 END IF 334 CALL derivatives_four_center(qs_env, rho_ao_orb, hfx_sections, & 335 para_env, irep, use_virial) 336 !Scale auxiliary density matrix for ADMMP back with 1/gsi(ispin) 337 IF (dft_control%do_admm) THEN 338 CALL scale_dm(qs_env, rho_ao_orb, scale_back=.TRUE.) 339 END IF 340 END IF 341 342 !! If required, the calculation of the forces will be done later with adiabatic rescaling 343 IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx 344 345 ! special case RTP/EMD we have a full complex density and HFX has a contrinution from the imaginary part 346 ehfxrt = 0.0_dp 347 IF (qs_env%run_rtp) THEN 348 349 CALL get_qs_env(qs_env=qs_env, rtp=rtp) 350 DO ispin = 1, nspins 351 CALL dbcsr_set(matrix_ks_im(ispin)%matrix, 0.0_dp) 352 END DO 353 IF (dft_control%do_admm) THEN 354 ! matrix_ks_orb => matrix_ks_aux_fit_im 355 ns = SIZE(matrix_ks_aux_fit_im) 356 matrix_ks_orb(1:ns, 1:1) => matrix_ks_aux_fit_im(1:ns) 357 DO ispin = 1, nspins 358 CALL dbcsr_set(matrix_ks_aux_fit_im(ispin)%matrix, 0.0_dp) 359 END DO 360 ELSE 361 ! matrix_ks_orb => matrix_ks_im 362 ns = SIZE(matrix_ks_im) 363 matrix_ks_orb(1:ns, 1:1) => matrix_ks_im(1:ns) 364 END IF 365 366 CALL qs_rho_get(rho_orb, rho_ao_im=rho_ao_1d) 367 ns = SIZE(rho_ao_1d) 368 rho_ao_orb(1:ns, 1:1) => rho_ao_1d(1:ns) 369 370 ehfxrt = 0.0_dp 371 DO ispin = 1, mspin 372 CALL integrate_four_center(qs_env, x_data, matrix_ks_orb, eh1, rho_ao_orb, hfx_sections, & 373 para_env, .FALSE., irep, distribute_fock_matrix, & 374 ispin=ispin) 375 ehfxrt = ehfxrt + eh1 376 END DO 377 378 IF (calculate_forces .AND. .NOT. do_adiabatic_rescaling) THEN 379 CALL derivatives_four_center(qs_env, rho_ao_orb, hfx_sections, & 380 para_env, irep, use_virial) 381 END IF 382 383 !! If required, the calculation of the forces will be done later with adiabatic rescaling 384 IF (do_adiabatic_rescaling) hf_energy(irep) = ehfx + ehfxrt 385 END IF 386 387 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 388 poisson_env=poisson_env) 389 CALL pw_hfx(qs_env, energy, hfx_sections, poisson_env, auxbas_pw_pool, irep) 390 391 END DO 392 393 ! *** Set the total HFX energy 394 energy%ex = ehfx + ehfxrt 395 396 ! *** Add Core-Hamiltonian-Matrix *** 397 DO ispin = 1, nspins 398 DO img = 1, nimages 399 CALL dbcsr_add(matrix_ks(ispin, img)%matrix, matrix_h(1, img)%matrix, & 400 1.0_dp, 1.0_dp) 401 END DO 402 END DO 403 IF (use_virial .AND. calculate_forces) THEN 404 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c 405 virial%pv_calculate = .FALSE. 406 ENDIF 407 408 !! If we perform adiabatic rescaling we are now able to rescale the xc-potential 409 IF (do_adiabatic_rescaling) THEN 410 CALL rescale_xc_potential(qs_env, matrix_ks, rho, energy, v_rspace_new, v_tau_rspace, & 411 hf_energy, just_energy, calculate_forces, use_virial) 412 END IF ! do_adiabatic_rescaling 413 414 CALL timestop(handle) 415 416 END SUBROUTINE hfx_ks_matrix 417 418! ************************************************************************************************** 419!> \brief computes the Hartree-Fock energy brute force in a pw basis 420!> \param qs_env ... 421!> \param energy ... 422!> \param hfx_section ... 423!> \param poisson_env ... 424!> \param auxbas_pw_pool ... 425!> \param irep ... 426!> \par History 427!> 12.2007 created [Joost VandeVondele] 428!> \note 429!> only computes the HFX energy, no derivatives as yet 430! ************************************************************************************************** 431 SUBROUTINE pw_hfx(qs_env, energy, hfx_section, poisson_env, auxbas_pw_pool, irep) 432 TYPE(qs_environment_type), POINTER :: qs_env 433 TYPE(qs_energy_type), POINTER :: energy 434 TYPE(section_vals_type), POINTER :: hfx_section 435 TYPE(pw_poisson_type), POINTER :: poisson_env 436 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 437 INTEGER :: irep 438 439 CHARACTER(*), PARAMETER :: routineN = 'pw_hfx', routineP = moduleN//':'//routineN 440 441 INTEGER :: blocksize, handle, iloc, iorb, & 442 iorb_block, ispin, iw, jloc, jorb, & 443 jorb_block, norb 444 LOGICAL :: do_pw_hfx 445 REAL(KIND=dp) :: exchange_energy, fraction, pair_energy, & 446 scaling 447 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 448 TYPE(cell_type), POINTER :: cell 449 TYPE(cp_fm_type), POINTER :: mo_coeff 450 TYPE(cp_logger_type), POINTER :: logger 451 TYPE(dbcsr_type), POINTER :: mo_coeff_b 452 TYPE(dft_control_type), POINTER :: dft_control 453 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 454 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 455 TYPE(pw_env_type), POINTER :: pw_env 456 TYPE(pw_p_type) :: pot_g, rho_g, rho_r 457 TYPE(pw_p_type), ALLOCATABLE, DIMENSION(:) :: rho_i, rho_j 458 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 459 460 CALL timeset(routineN, handle) 461 logger => cp_get_default_logger() 462 463 CALL section_vals_val_get(hfx_section, "PW_HFX", l_val=do_pw_hfx, i_rep_section=irep) 464 465 IF (do_pw_hfx) THEN 466 CALL section_vals_val_get(hfx_section, "FRACTION", r_val=fraction) 467 CALL section_vals_val_get(hfx_section, "PW_HFX_BLOCKSIZE", i_val=blocksize) 468 469 CALL get_qs_env(qs_env, mos=mo_array, pw_env=pw_env, dft_control=dft_control, & 470 cell=cell, particle_set=particle_set, & 471 atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set) 472 473 ! limit the blocksize by the number of orbitals 474 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff=mo_coeff) 475 CALL cp_fm_get_info(mo_coeff, ncol_global=norb) 476 blocksize = MAX(1, MIN(blocksize, norb)) 477 478 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r%pw, & 479 use_data=REALDATA3D, & 480 in_space=REALSPACE) 481 CALL pw_pool_create_pw(auxbas_pw_pool, rho_g%pw, & 482 use_data=COMPLEXDATA1D, & 483 in_space=RECIPROCALSPACE) 484 CALL pw_pool_create_pw(auxbas_pw_pool, pot_g%pw, & 485 use_data=COMPLEXDATA1D, & 486 in_space=RECIPROCALSPACE) 487 488 ALLOCATE (rho_i(blocksize)) 489 ALLOCATE (rho_j(blocksize)) 490 491 DO iorb_block = 1, blocksize 492 NULLIFY (rho_i(iorb_block)%pw) 493 CALL pw_create(rho_i(iorb_block)%pw, rho_r%pw%pw_grid, & 494 use_data=REALDATA3D, & 495 in_space=REALSPACE) 496 NULLIFY (rho_j(iorb_block)%pw) 497 CALL pw_create(rho_j(iorb_block)%pw, rho_r%pw%pw_grid, & 498 use_data=REALDATA3D, & 499 in_space=REALSPACE) 500 ENDDO 501 502 exchange_energy = 0.0_dp 503 504 DO ispin = 1, SIZE(mo_array) 505 CALL get_mo_set(mo_set=mo_array(ispin)%mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b) 506 507 IF (mo_array(ispin)%mo_set%use_mo_coeff_b) THEN !fm->dbcsr 508 CALL copy_dbcsr_to_fm(mo_coeff_b, mo_coeff) !fm->dbcsr 509 ENDIF !fm->dbcsr 510 511 CALL cp_fm_get_info(mo_coeff, ncol_global=norb) 512 513 DO iorb_block = 1, norb, blocksize 514 515 DO iorb = iorb_block, MIN(iorb_block + blocksize - 1, norb) 516 517 iloc = iorb - iorb_block + 1 518 CALL calculate_wavefunction(mo_coeff, iorb, rho_i(iloc), rho_g, & 519 atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, & 520 pw_env) 521 522 ENDDO 523 524 DO jorb_block = iorb_block, norb, blocksize 525 526 DO jorb = jorb_block, MIN(jorb_block + blocksize - 1, norb) 527 528 jloc = jorb - jorb_block + 1 529 CALL calculate_wavefunction(mo_coeff, jorb, rho_j(jloc), rho_g, & 530 atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, & 531 pw_env) 532 533 ENDDO 534 535 DO iorb = iorb_block, MIN(iorb_block + blocksize - 1, norb) 536 iloc = iorb - iorb_block + 1 537 DO jorb = jorb_block, MIN(jorb_block + blocksize - 1, norb) 538 jloc = jorb - jorb_block + 1 539 IF (jorb < iorb) CYCLE 540 541 ! compute the pair density 542 rho_r%pw%cr3d = rho_i(iloc)%pw%cr3d*rho_j(jloc)%pw%cr3d 543 544 ! go the g-space and compute hartree energy 545 CALL pw_transfer(rho_r%pw, rho_g%pw) 546 CALL pw_poisson_solve(poisson_env, rho_g%pw, pair_energy, pot_g%pw) 547 548 ! sum up to the full energy 549 scaling = fraction 550 IF (SIZE(mo_array) == 1) scaling = scaling*2.0_dp 551 IF (iorb /= jorb) scaling = scaling*2.0_dp 552 553 exchange_energy = exchange_energy - scaling*pair_energy 554 555 ENDDO 556 ENDDO 557 558 ENDDO 559 ENDDO 560 ENDDO 561 562 DO iorb_block = 1, blocksize 563 CALL pw_release(rho_i(iorb_block)%pw) 564 CALL pw_release(rho_j(iorb_block)%pw) 565 ENDDO 566 567 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r%pw) 568 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_g%pw) 569 CALL pw_pool_give_back_pw(auxbas_pw_pool, pot_g%pw) 570 571 iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", & 572 extension=".scfLog") 573 574 IF (iw > 0) THEN 575 WRITE (UNIT=iw, FMT="((T3,A,T61,F20.10))") & 576 "HF_PW_HFX| PW exchange energy:", exchange_energy 577 WRITE (UNIT=iw, FMT="((T3,A,T61,F20.10),/)") & 578 "HF_PW_HFX| Gaussian exchange energy:", energy%ex 579 ENDIF 580 581 CALL cp_print_key_finished_output(iw, logger, hfx_section, & 582 "HF_INFO") 583 584 ENDIF 585 586 CALL timestop(handle) 587 588 END SUBROUTINE pw_hfx 589 590! ************************************************************************************************** 591!> \brief This routine modifies the xc section depending on the potential type 592!> used for the HF exchange and the resulting correction term. Currently 593!> three types of corrections are implemented: 594!> 595!> coulomb: Ex,hf = Ex,hf' + (PBEx-PBEx') 596!> shortrange: Ex,hf = Ex,hf' + (XWPBEX-XWPBEX') 597!> truncated: Ex,hf = Ex,hf' + ( (XWPBEX0-PBE_HOLE_TC_LR) -(XWPBEX0-PBE_HOLE_TC_LR)' ) 598!> 599!> with ' denoting the auxiliary basis set and 600!> 601!> PBEx: PBE exchange functional 602!> XWPBEX: PBE exchange hole for short-range potential (erfc(omega*r)/r) 603!> XWPBEX0: PBE exchange hole for standard coulomb potential 604!> PBE_HOLE_TC_LR: PBE exchange hole for longrange truncated coulomb potential 605!> 606!> 607!> \param qs_env the qs environment 608!> \param xc_section the original xc_section 609!> \param admm_env the ADMM environment 610!> \par History 611!> 12.2009 created [Manuel Guidon] 612!> \author Manuel Guidon 613! ************************************************************************************************** 614 SUBROUTINE create_admm_xc_section(qs_env, xc_section, admm_env) 615 TYPE(qs_environment_type), POINTER :: qs_env 616 TYPE(section_vals_type), POINTER :: xc_section 617 TYPE(admm_type), POINTER :: admm_env 618 619 CHARACTER(LEN=*), PARAMETER :: routineN = 'create_admm_xc_section', & 620 routineP = moduleN//':'//routineN 621 622 CHARACTER(LEN=20) :: name_x_func 623 INTEGER :: hfx_potential_type, ifun, nfun 624 LOGICAL :: funct_found 625 REAL(dp) :: cutoff_radius, hfx_fraction, omega, & 626 scale_x 627 TYPE(section_vals_type), POINTER :: xc_fun, xc_fun_section 628 629 NULLIFY (admm_env%xc_section_aux, admm_env%xc_section_primary) 630 631 CALL get_qs_env(qs_env) 632 633 !! ** Duplicate existing xc-section 634 CALL section_vals_duplicate(xc_section, admm_env%xc_section_aux) 635 CALL section_vals_duplicate(xc_section, admm_env%xc_section_primary) 636 !** Now modify the auxiliary basis 637 !** First remove all functionals 638 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL") 639 640 !* Overwrite possible shortcut 641 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", & 642 i_val=xc_funct_no_shortcut) 643 644 !** Get number of Functionals in the list 645 ifun = 0 646 nfun = 0 647 DO 648 ifun = ifun + 1 649 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 650 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 651 nfun = nfun + 1 652 END DO 653 654 ifun = 0 655 DO ifun = 1, nfun 656 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=1) 657 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 658 CALL section_vals_remove_values(xc_fun) 659 END DO 660 661 hfx_potential_type = qs_env%x_data(1, 1)%potential_parameter%potential_type 662 hfx_fraction = qs_env%x_data(1, 1)%general_parameter%fraction 663 664 !in case of no admm exchange corr., no auxiliary exchange functional needed 665 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN 666 hfx_fraction = 0.0_dp 667 END IF 668 669 ! default PBE Functional 670 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_default .OR. & 671 admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN 672 673 !! ** Add functionals evaluated with auxiliary basis 674 SELECT CASE (hfx_potential_type) 675 CASE (do_potential_coulomb) 676 CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", & 677 l_val=.TRUE.) 678 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", & 679 r_val=-hfx_fraction) 680 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", & 681 r_val=0.0_dp) 682 CASE (do_potential_short) 683 omega = qs_env%x_data(1, 1)%potential_parameter%omega 684 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", & 685 l_val=.TRUE.) 686 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", & 687 r_val=-hfx_fraction) 688 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", & 689 r_val=0.0_dp) 690 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", & 691 r_val=omega) 692 CASE (do_potential_truncated) 693 cutoff_radius = qs_env%x_data(1, 1)%potential_parameter%cutoff_radius 694 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", & 695 l_val=.TRUE.) 696 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", & 697 r_val=hfx_fraction) 698 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", & 699 r_val=cutoff_radius) 700 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", & 701 l_val=.TRUE.) 702 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", & 703 r_val=0.0_dp) 704 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", & 705 r_val=-hfx_fraction) 706 CASE DEFAULT 707 CPABORT("") 708 END SELECT 709 710 !** Now modify the functionals for the primary basis 711 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL") 712 !* Overwrite possible shortcut 713 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", & 714 i_val=xc_funct_no_shortcut) 715 716 SELECT CASE (hfx_potential_type) 717 CASE (do_potential_coulomb) 718 ifun = 0 719 funct_found = .FALSE. 720 DO 721 ifun = ifun + 1 722 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 723 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 724 IF (xc_fun%section%name == "PBE") THEN 725 funct_found = .TRUE. 726 END IF 727 END DO 728 IF (.NOT. funct_found) THEN 729 CALL section_vals_val_set(xc_fun_section, "PBE%_SECTION_PARAMETERS_", & 730 l_val=.TRUE.) 731 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", & 732 r_val=hfx_fraction) 733 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_C", & 734 r_val=0.0_dp) 735 ELSE 736 CALL section_vals_val_get(xc_fun_section, "PBE%SCALE_X", & 737 r_val=scale_x) 738 scale_x = scale_x + hfx_fraction 739 CALL section_vals_val_set(xc_fun_section, "PBE%SCALE_X", & 740 r_val=scale_x) 741 END IF 742 CASE (do_potential_short) 743 omega = qs_env%x_data(1, 1)%potential_parameter%omega 744 ifun = 0 745 funct_found = .FALSE. 746 DO 747 ifun = ifun + 1 748 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 749 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 750 IF (xc_fun%section%name == "XWPBE") THEN 751 funct_found = .TRUE. 752 END IF 753 END DO 754 IF (.NOT. funct_found) THEN 755 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", & 756 l_val=.TRUE.) 757 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", & 758 r_val=hfx_fraction) 759 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", & 760 r_val=0.0_dp) 761 CALL section_vals_val_set(xc_fun_section, "XWPBE%OMEGA", & 762 r_val=omega) 763 ELSE 764 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X", & 765 r_val=scale_x) 766 scale_x = scale_x + hfx_fraction 767 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", & 768 r_val=scale_x) 769 END IF 770 CASE (do_potential_truncated) 771 cutoff_radius = qs_env%x_data(1, 1)%potential_parameter%cutoff_radius 772 ifun = 0 773 funct_found = .FALSE. 774 DO 775 ifun = ifun + 1 776 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 777 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 778 IF (xc_fun%section%name == "PBE_HOLE_T_C_LR") THEN 779 funct_found = .TRUE. 780 END IF 781 END DO 782 IF (.NOT. funct_found) THEN 783 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%_SECTION_PARAMETERS_", & 784 l_val=.TRUE.) 785 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", & 786 r_val=-hfx_fraction) 787 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", & 788 r_val=cutoff_radius) 789 790 ELSE 791 CALL section_vals_val_get(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", & 792 r_val=scale_x) 793 scale_x = scale_x - hfx_fraction 794 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%SCALE_X", & 795 r_val=scale_x) 796 CALL section_vals_val_set(xc_fun_section, "PBE_HOLE_T_C_LR%CUTOFF_RADIUS", & 797 r_val=cutoff_radius) 798 END IF 799 ifun = 0 800 funct_found = .FALSE. 801 DO 802 ifun = ifun + 1 803 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 804 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 805 IF (xc_fun%section%name == "XWPBE") THEN 806 funct_found = .TRUE. 807 END IF 808 END DO 809 IF (.NOT. funct_found) THEN 810 CALL section_vals_val_set(xc_fun_section, "XWPBE%_SECTION_PARAMETERS_", & 811 l_val=.TRUE.) 812 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", & 813 r_val=hfx_fraction) 814 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X", & 815 r_val=0.0_dp) 816 817 ELSE 818 CALL section_vals_val_get(xc_fun_section, "XWPBE%SCALE_X0", & 819 r_val=scale_x) 820 scale_x = scale_x + hfx_fraction 821 CALL section_vals_val_set(xc_fun_section, "XWPBE%SCALE_X0", & 822 r_val=scale_x) 823 END IF 824 825 END SELECT 826 827 ! PBEX (always bare form), OPTX and Becke88 functional 828 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex .OR. & 829 admm_env%aux_exch_func == do_admm_aux_exch_func_opt .OR. & 830 admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN 831 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN 832 name_x_func = 'PBE' 833 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN 834 name_x_func = 'OPTX' 835 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_bee) THEN 836 name_x_func = 'BECKE88' 837 END IF 838 !primary basis 839 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", & 840 l_val=.TRUE.) 841 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", & 842 r_val=-hfx_fraction) 843 844 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN 845 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_C", r_val=0.0_dp) 846 END IF 847 848 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN 849 IF (admm_env%aux_exch_func_param) THEN 850 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A1", & 851 r_val=admm_env%aux_x_param(1)) 852 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A2", & 853 r_val=admm_env%aux_x_param(2)) 854 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%GAMMA", & 855 r_val=admm_env%aux_x_param(3)) 856 END IF 857 END IF 858 859 !** Now modify the functionals for the primary basis 860 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL") 861 !* Overwrite possible L") 862 !* Overwrite possible shortcut 863 CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", & 864 i_val=xc_funct_no_shortcut) 865 866 ifun = 0 867 funct_found = .FALSE. 868 DO 869 ifun = ifun + 1 870 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 871 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 872 IF (xc_fun%section%name == TRIM(name_x_func)) THEN 873 funct_found = .TRUE. 874 END IF 875 END DO 876 IF (.NOT. funct_found) THEN 877 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%_SECTION_PARAMETERS_", & 878 l_val=.TRUE.) 879 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", & 880 r_val=hfx_fraction) 881 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_pbex) THEN 882 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_C", & 883 r_val=0.0_dp) 884 ELSE IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN 885 IF (admm_env%aux_exch_func_param) THEN 886 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A1", & 887 r_val=admm_env%aux_x_param(1)) 888 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%A2", & 889 r_val=admm_env%aux_x_param(2)) 890 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%GAMMA", & 891 r_val=admm_env%aux_x_param(3)) 892 END IF 893 END IF 894 895 ELSE 896 CALL section_vals_val_get(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", & 897 r_val=scale_x) 898 scale_x = scale_x + hfx_fraction 899 CALL section_vals_val_set(xc_fun_section, TRIM(name_x_func)//"%SCALE_X", & 900 r_val=scale_x) 901 IF (admm_env%aux_exch_func == do_admm_aux_exch_func_opt) THEN 902 CPASSERT(.NOT. admm_env%aux_exch_func_param) 903 END IF 904 END IF 905 906 END IF 907 908 IF (1 == 0) THEN 909 WRITE (*, *) "primary" 910 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_primary, "XC_FUNCTIONAL") 911 ifun = 0 912 funct_found = .FALSE. 913 DO 914 ifun = ifun + 1 915 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 916 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 917 918 scale_x = -1000.0_dp 919 IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN 920 CALL section_vals_val_get(xc_fun, "SCALE_X", & 921 r_val=scale_x) 922 END IF 923 IF (xc_fun%section%name == "XWPBE") THEN 924 CALL section_vals_val_get(xc_fun, "SCALE_X0", & 925 r_val=hfx_fraction) 926 927 WRITE (*, *) xc_fun%section%name, scale_x, hfx_fraction 928 ELSE 929 WRITE (*, *) xc_fun%section%name, scale_x 930 END IF 931 END DO 932 933 WRITE (*, *) "auxiliary" 934 xc_fun_section => section_vals_get_subs_vals(admm_env%xc_section_aux, "XC_FUNCTIONAL") 935 ifun = 0 936 funct_found = .FALSE. 937 DO 938 ifun = ifun + 1 939 xc_fun => section_vals_get_subs_vals2(xc_fun_section, i_section=ifun) 940 IF (.NOT. ASSOCIATED(xc_fun)) EXIT 941 scale_x = -1000.0_dp 942 IF (xc_fun%section%name /= "LYP" .AND. xc_fun%section%name /= "VWN") THEN 943 CALL section_vals_val_get(xc_fun, "SCALE_X", & 944 r_val=scale_x) 945 END IF 946 IF (xc_fun%section%name == "XWPBE") THEN 947 CALL section_vals_val_get(xc_fun, "SCALE_X0", & 948 r_val=hfx_fraction) 949 950 WRITE (*, *) xc_fun%section%name, scale_x, hfx_fraction 951 ELSE 952 WRITE (*, *) xc_fun%section%name, scale_x 953 END IF 954 END DO 955 END IF 956 957 END SUBROUTINE create_admm_xc_section 958 959! ************************************************************************************************** 960!> \brief Add the hfx contributions to the Hamiltonian 961!> 962!> \param matrix_ks Kohn-Sham matrix (updated on exit) 963!> \param rho_ao electron density expressed in terms of atomic orbitals 964!> \param qs_env Quickstep environment 965!> \note 966!> Simplified version of subroutine hfx_ks_matrix() 967! ************************************************************************************************** 968 SUBROUTINE tddft_hfx_matrix(matrix_ks, rho_ao, qs_env) 969 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, rho_ao 970 TYPE(qs_environment_type), POINTER :: qs_env 971 972 CHARACTER(LEN=*), PARAMETER :: routineN = 'tddft_hfx_matrix', & 973 routineP = moduleN//':'//routineN 974 975 INTEGER :: handle, irep, ispin, mspin, n_rep_hf, & 976 nspins 977 LOGICAL :: distribute_fock_matrix, do_hfx, & 978 hfx_treat_lsd_in_core, & 979 s_mstruct_changed 980 REAL(KIND=dp) :: eh1, ehfx 981 TYPE(cp_para_env_type), POINTER :: para_env 982 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, rho_ao_kp 983 TYPE(dft_control_type), POINTER :: dft_control 984 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data 985 TYPE(qs_energy_type), POINTER :: energy 986 TYPE(section_vals_type), POINTER :: hfx_sections, input 987 988 CALL timeset(routineN, handle) 989 990 NULLIFY (dft_control, hfx_sections, input, para_env, matrix_ks_kp, rho_ao_kp) 991 992 CALL get_qs_env(qs_env=qs_env, & 993 dft_control=dft_control, & 994 energy=energy, & 995 input=input, & 996 para_env=para_env, & 997 s_mstruct_changed=s_mstruct_changed, & 998 x_data=x_data) 999 1000 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") 1001 CALL section_vals_get(hfx_sections, explicit=do_hfx) 1002 1003 IF (do_hfx) THEN 1004 CPASSERT(dft_control%nimages == 1) 1005 nspins = dft_control%nspins 1006 1007 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) 1008 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, & 1009 i_rep_section=1) 1010 1011 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) 1012 distribute_fock_matrix = .TRUE. 1013 1014 mspin = 1 1015 IF (hfx_treat_lsd_in_core) mspin = nspins 1016 1017 matrix_ks_kp(1:nspins, 1:1) => matrix_ks(1:nspins) 1018 rho_ao_kp(1:nspins, 1:1) => rho_ao(1:nspins) 1019 1020 DO irep = 1, n_rep_hf 1021 ! the real hfx calulation 1022 ehfx = 0.0_dp 1023 DO ispin = 1, mspin 1024 CALL integrate_four_center(qs_env, x_data, matrix_ks_kp, eh1, rho_ao_kp, hfx_sections, para_env, & 1025 s_mstruct_changed, irep, distribute_fock_matrix, ispin=ispin) 1026 ehfx = ehfx + eh1 1027 END DO 1028 END DO 1029 energy%ex = ehfx 1030 END IF 1031 1032 CALL timestop(handle) 1033 END SUBROUTINE tddft_hfx_matrix 1034 1035END MODULE hfx_admm_utils 1036