1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief routines that build the Kohn-Sham matrix (i.e calculate the coulomb 8!> and xc parts 9!> \par History 10!> 05.2002 moved from qs_scf (see there the history) [fawzi] 11!> JGH [30.08.02] multi-grid arrays independent from density and potential 12!> 10.2002 introduced pools, uses updated rho as input, 13!> removed most temporary variables, renamed may vars, 14!> began conversion to LSD [fawzi] 15!> 10.2004 moved calculate_w_matrix here [Joost VandeVondele] 16!> introduced energy derivative wrt MOs [Joost VandeVondele] 17!> \author Fawzi Mohamed 18! ************************************************************************************************** 19 20MODULE qs_ks_utils 21 USE admm_types, ONLY: admm_type 22 USE atomic_kind_types, ONLY: atomic_kind_type 23 USE cell_types, ONLY: cell_type 24 USE cp_control_types, ONLY: dft_control_type 25 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 26 copy_fm_to_dbcsr,& 27 cp_dbcsr_plus_fm_fm_t,& 28 cp_dbcsr_sm_fm_multiply,& 29 dbcsr_allocate_matrix_set,& 30 dbcsr_deallocate_matrix_set 31 USE cp_ddapc, ONLY: cp_ddapc_apply_CD 32 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 33 cp_fm_struct_release,& 34 cp_fm_struct_type 35 USE cp_fm_types, ONLY: cp_fm_create,& 36 cp_fm_get_info,& 37 cp_fm_p_type,& 38 cp_fm_release,& 39 cp_fm_set_all,& 40 cp_fm_to_fm,& 41 cp_fm_type 42 USE cp_log_handling, ONLY: cp_get_default_logger,& 43 cp_logger_type,& 44 cp_to_string 45 USE cp_output_handling, ONLY: cp_p_file,& 46 cp_print_key_finished_output,& 47 cp_print_key_should_output,& 48 cp_print_key_unit_nr 49 USE cp_para_types, ONLY: cp_para_env_type 50 USE dbcsr_api, ONLY: & 51 dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_dot, dbcsr_get_info, dbcsr_init_p, & 52 dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_scale_by_vector, & 53 dbcsr_set, dbcsr_type 54 USE input_constants, ONLY: & 55 cdft_alpha_constraint, cdft_beta_constraint, cdft_charge_constraint, & 56 cdft_magnetization_constraint, do_admm_aux_exch_func_none, do_ppl_grid, sic_ad, sic_eo, & 57 sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none 58 USE input_section_types, ONLY: section_vals_get_subs_vals,& 59 section_vals_type,& 60 section_vals_val_get 61 USE kahan_sum, ONLY: accurate_sum 62 USE kinds, ONLY: default_string_length,& 63 dp 64 USE kpoint_types, ONLY: get_kpoint_info,& 65 kpoint_type 66 USE lri_environment_methods, ONLY: v_int_ppl_update 67 USE lri_environment_types, ONLY: lri_density_type,& 68 lri_environment_type,& 69 lri_kind_type 70 USE lri_forces, ONLY: calculate_lri_forces,& 71 calculate_ri_forces 72 USE lri_ks_methods, ONLY: calculate_lri_ks_matrix,& 73 calculate_ri_ks_matrix 74 USE message_passing, ONLY: mp_sum 75 USE ps_implicit_types, ONLY: MIXED_BC,& 76 MIXED_PERIODIC_BC,& 77 NEUMANN_BC,& 78 PERIODIC_BC 79 USE pw_env_types, ONLY: pw_env_get,& 80 pw_env_type 81 USE pw_grid_types, ONLY: PW_MODE_DISTRIBUTED 82 USE pw_methods, ONLY: pw_axpy,& 83 pw_copy,& 84 pw_integrate_function,& 85 pw_scale,& 86 pw_transfer,& 87 pw_zero 88 USE pw_poisson_methods, ONLY: pw_poisson_solve 89 USE pw_poisson_types, ONLY: pw_poisson_implicit,& 90 pw_poisson_type 91 USE pw_pool_types, ONLY: pw_pool_create_pw,& 92 pw_pool_give_back_pw,& 93 pw_pool_type 94 USE pw_types, ONLY: COMPLEXDATA1D,& 95 REALDATA3D,& 96 REALSPACE,& 97 RECIPROCALSPACE,& 98 pw_p_type 99 USE qs_cdft_types, ONLY: cdft_control_type 100 USE qs_charges_types, ONLY: qs_charges_type 101 USE qs_collocate_density, ONLY: calculate_rho_elec 102 USE qs_energy_types, ONLY: qs_energy_type 103 USE qs_environment_types, ONLY: get_qs_env,& 104 qs_environment_type 105 USE qs_force_types, ONLY: qs_force_type 106 USE qs_integrate_potential, ONLY: integrate_v_rspace,& 107 integrate_v_rspace_diagonal,& 108 integrate_v_rspace_one_center 109 USE qs_kind_types, ONLY: get_qs_kind_set,& 110 qs_kind_type 111 USE qs_ks_qmmm_methods, ONLY: qmmm_modify_hartree_pot 112 USE qs_ks_types, ONLY: get_ks_env,& 113 qs_ks_env_type 114 USE qs_mo_types, ONLY: get_mo_set,& 115 mo_set_p_type 116 USE qs_rho_types, ONLY: qs_rho_get,& 117 qs_rho_type 118 USE task_list_types, ONLY: task_list_type 119 USE virial_types, ONLY: virial_type 120 USE xc, ONLY: xc_exc_calc,& 121 xc_vxc_pw_create1 122#include "./base/base_uses.f90" 123 124 IMPLICIT NONE 125 126 PRIVATE 127 128 LOGICAL, PARAMETER :: debug_this_module = .TRUE. 129 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_utils' 130 131 PUBLIC :: low_spin_roks, sic_explicit_orbitals, calc_v_sic_rspace, print_densities, & 132 print_detailed_energy, compute_matrix_vxc, sum_up_and_integrate, & 133 calculate_zmp_potential, get_embed_potential_energy 134 135CONTAINS 136 137! ************************************************************************************************** 138!> \brief do ROKS calculations yielding low spin states 139!> \param energy ... 140!> \param qs_env ... 141!> \param dft_control ... 142!> \param just_energy ... 143!> \param calculate_forces ... 144!> \param auxbas_pw_pool ... 145! ************************************************************************************************** 146 SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, & 147 calculate_forces, auxbas_pw_pool) 148 149 TYPE(qs_energy_type), POINTER :: energy 150 TYPE(qs_environment_type), POINTER :: qs_env 151 TYPE(dft_control_type), POINTER :: dft_control 152 LOGICAL, INTENT(IN) :: just_energy, calculate_forces 153 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 154 155 CHARACTER(*), PARAMETER :: routineN = 'low_spin_roks', routineP = moduleN//':'//routineN 156 157 INTEGER :: handle, ispin, iterm, k, k_alpha, & 158 k_beta, n_rep, Nelectron, Nspin, Nterms 159 INTEGER, DIMENSION(:), POINTER :: ivec 160 INTEGER, DIMENSION(:, :, :), POINTER :: occupations 161 LOGICAL :: compute_virial, in_range, & 162 uniform_occupation 163 REAL(KIND=dp) :: exc, total_rho 164 REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp 165 REAL(KIND=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling 166 TYPE(cell_type), POINTER :: cell 167 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_p, mo_derivs, rho_ao 168 TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, & 169 mo_coeff 170 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 171 TYPE(pw_env_type), POINTER :: pw_env 172 TYPE(pw_p_type) :: work_v_rspace 173 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r, tau, vxc, vxc_tau 174 TYPE(pw_pool_type), POINTER :: xc_pw_pool 175 TYPE(qs_ks_env_type), POINTER :: ks_env 176 TYPE(qs_rho_type), POINTER :: rho 177 TYPE(section_vals_type), POINTER :: input, low_spin_roks_section, xc_section 178 TYPE(virial_type), POINTER :: virial 179 180 IF (.NOT. dft_control%low_spin_roks) RETURN 181 NULLIFY (ks_env, rho_ao) 182 183 CALL timeset(routineN, handle) 184 185 CALL get_qs_env(qs_env, & 186 ks_env=ks_env, & 187 mo_derivs=mo_derivs, & 188 mos=mo_array, & 189 rho=rho, & 190 pw_env=pw_env, & 191 input=input, & 192 cell=cell, & 193 virial=virial) 194 195 CALL qs_rho_get(rho, rho_ao=rho_ao) 196 197 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer) 198 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 199 200 ! some assumptions need to be checked 201 ! we have two spins 202 CPASSERT(SIZE(mo_array, 1) == 2) 203 Nspin = 2 204 ! we want uniform occupations 205 CALL get_mo_set(mo_set=mo_array(1)%mo_set, uniform_occupation=uniform_occupation) 206 CPASSERT(uniform_occupation) 207 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation) 208 CPASSERT(uniform_occupation) 209 210 NULLIFY (dbcsr_deriv) 211 CALL dbcsr_init_p(dbcsr_deriv) 212 CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix) 213 CALL dbcsr_set(dbcsr_deriv, 0.0_dp) 214 215 ! basic info 216 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff_b=mo_coeff) 217 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha) 218 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff_b=mo_coeff) 219 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta) 220 221 ! read the input 222 low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS") 223 224 CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec) 225 Nterms = SIZE(rvec) 226 ALLOCATE (energy_scaling(Nterms)) 227 energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp? 228 229 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep) 230 CPASSERT(n_rep == Nterms) 231 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec) 232 Nelectron = SIZE(ivec) 233 CPASSERT(Nelectron == k_alpha - k_beta) 234 ALLOCATE (occupations(2, Nelectron, Nterms)) 235 occupations = 0 236 DO iterm = 1, Nterms 237 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec) 238 CPASSERT(Nelectron == SIZE(ivec)) 239 in_range = ALL(ivec >= 1) .AND. ALL(ivec <= 2) 240 CPASSERT(in_range) 241 DO k = 1, Nelectron 242 occupations(ivec(k), k, iterm) = 1 243 ENDDO 244 ENDDO 245 246 ! set up general data structures 247 ! density matrices, kohn-sham matrices 248 249 NULLIFY (matrix_p) 250 CALL dbcsr_allocate_matrix_set(matrix_p, Nspin) 251 DO ispin = 1, Nspin 252 ALLOCATE (matrix_p(ispin)%matrix) 253 CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, & 254 name="density matrix low spin roks") 255 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp) 256 ENDDO 257 258 NULLIFY (matrix_h) 259 CALL dbcsr_allocate_matrix_set(matrix_h, Nspin) 260 DO ispin = 1, Nspin 261 ALLOCATE (matrix_h(ispin)%matrix) 262 CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, & 263 name="KS matrix low spin roks") 264 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp) 265 ENDDO 266 267 ! grids in real and g space for rho and vxc 268 ! tau functionals are not supported 269 NULLIFY (tau, vxc_tau, vxc) 270 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool) 271 272 ALLOCATE (rho_r(Nspin)) 273 ALLOCATE (rho_g(Nspin)) 274 DO ispin = 1, Nspin 275 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(ispin)%pw, & 276 use_data=REALDATA3D, & 277 in_space=REALSPACE) 278 CALL pw_pool_create_pw(auxbas_pw_pool, rho_g(ispin)%pw, & 279 use_data=COMPLEXDATA1D, & 280 in_space=RECIPROCALSPACE) 281 ENDDO 282 CALL pw_pool_create_pw(auxbas_pw_pool, work_v_rspace%pw, & 283 use_data=REALDATA3D, & 284 in_space=REALSPACE) 285 286 ! get mo matrices needed to construct the density matrices 287 ! we will base all on the alpha spin matrix, obviously possible in ROKS 288 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff_b=mo_coeff) 289 NULLIFY (fm_scaled, fm_deriv) 290 CALL dbcsr_init_p(fm_scaled) 291 CALL dbcsr_init_p(fm_deriv) 292 CALL dbcsr_copy(fm_scaled, mo_coeff) 293 CALL dbcsr_copy(fm_deriv, mo_coeff) 294 295 ALLOCATE (scaling(k_alpha)) 296 297 ! for each term, add it with the given scaling factor to the energy, and compute the required derivatives 298 DO iterm = 1, Nterms 299 300 DO ispin = 1, Nspin 301 ! compute the proper density matrices with the required occupations 302 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp) 303 scaling = 1.0_dp 304 scaling(k_alpha - Nelectron + 1:k_alpha) = occupations(ispin, :, iterm) 305 CALL dbcsr_copy(fm_scaled, mo_coeff) 306 CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right') 307 CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, & 308 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.TRUE.) 309 ! compute the densities on the grid 310 CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, & 311 rho=rho_r(ispin), rho_gspace=rho_g(ispin), total_rho=total_rho, & 312 ks_env=ks_env) 313 ENDDO 314 315 ! compute the exchange energies / potential if needed 316 IF (just_energy) THEN 317 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, & 318 pw_pool=xc_pw_pool) 319 ELSE 320 CPASSERT(.NOT. compute_virial) 321 CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, & 322 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, & 323 pw_pool=xc_pw_pool, compute_virial=.FALSE., virial_xc=virial_xc_tmp) 324 END IF 325 326 energy%exc = energy%exc + energy_scaling(iterm)*exc 327 328 ! add the corresponding derivatives to the MO derivatives 329 IF (.NOT. just_energy) THEN 330 ! get the potential in matrix form 331 DO ispin = 1, Nspin 332 ! use a work_v_rspace 333 work_v_rspace%pw%cr3d = (energy_scaling(iterm)*vxc(ispin)%pw%pw_grid%dvol)* & 334 vxc(ispin)%pw%cr3d 335 ! zero first ?! 336 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp) 337 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), & 338 qs_env=qs_env, calculate_forces=calculate_forces) 339 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc(ispin)%pw) 340 ENDDO 341 DEALLOCATE (vxc) 342 343 ! add this to the mo_derivs, again based on the alpha mo_coeff 344 DO ispin = 1, Nspin 345 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, & 346 0.0_dp, dbcsr_deriv, last_column=k_alpha) 347 348 scaling = 1.0_dp 349 scaling(k_alpha - Nelectron + 1:k_alpha) = occupations(ispin, :, iterm) 350 CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right') 351 CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp) 352 ENDDO 353 354 ENDIF 355 356 ENDDO 357 358 ! release allocated memory 359 DO ispin = 1, Nspin 360 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_r(ispin)%pw) 361 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_g(ispin)%pw) 362 ENDDO 363 DEALLOCATE (rho_r, rho_g) 364 CALL dbcsr_deallocate_matrix_set(matrix_p) 365 CALL dbcsr_deallocate_matrix_set(matrix_h) 366 367 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_v_rspace%pw) 368 369 CALL dbcsr_release_p(fm_deriv) 370 CALL dbcsr_release_p(fm_scaled) 371 372 DEALLOCATE (occupations) 373 DEALLOCATE (energy_scaling) 374 DEALLOCATE (scaling) 375 376 CALL dbcsr_release_p(dbcsr_deriv) 377 378 CALL timestop(handle) 379 380 END SUBROUTINE low_spin_roks 381! ************************************************************************************************** 382!> \brief do sic calculations on explicit orbitals 383!> \param energy ... 384!> \param qs_env ... 385!> \param dft_control ... 386!> \param poisson_env ... 387!> \param just_energy ... 388!> \param calculate_forces ... 389!> \param auxbas_pw_pool ... 390! ************************************************************************************************** 391 SUBROUTINE sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, & 392 calculate_forces, auxbas_pw_pool) 393 394 TYPE(qs_energy_type), POINTER :: energy 395 TYPE(qs_environment_type), POINTER :: qs_env 396 TYPE(dft_control_type), POINTER :: dft_control 397 TYPE(pw_poisson_type), POINTER :: poisson_env 398 LOGICAL, INTENT(IN) :: just_energy, calculate_forces 399 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 400 401 CHARACTER(*), PARAMETER :: routineN = 'sic_explicit_orbitals', & 402 routineP = moduleN//':'//routineN 403 404 INTEGER :: handle, i, Iorb, k_alpha, k_beta, Norb 405 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: sic_orbital_list 406 LOGICAL :: compute_virial, uniform_occupation 407 REAL(KIND=dp) :: ener, exc, total_rho 408 REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp 409 TYPE(cell_type), POINTER :: cell 410 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: mo_derivs_local 411 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 412 TYPE(cp_fm_type), POINTER :: matrix_hv, matrix_v, mo_coeff 413 TYPE(dbcsr_p_type) :: orb_density_matrix_p, orb_h_p 414 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, rho_ao, tmp_dbcsr 415 TYPE(dbcsr_type), POINTER :: orb_density_matrix, orb_h 416 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 417 TYPE(pw_env_type), POINTER :: pw_env 418 TYPE(pw_p_type) :: orb_rho_g, orb_rho_r, tmp_g, tmp_r, & 419 work_v_gspace, work_v_rspace 420 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r, tau, vxc, vxc_tau 421 TYPE(pw_pool_type), POINTER :: xc_pw_pool 422 TYPE(qs_ks_env_type), POINTER :: ks_env 423 TYPE(qs_rho_type), POINTER :: rho 424 TYPE(section_vals_type), POINTER :: input, xc_section 425 TYPE(virial_type), POINTER :: virial 426 427 IF (dft_control%sic_method_id .NE. sic_eo) RETURN 428 429 CALL timeset(routineN, handle) 430 431 NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao) 432 433 ! generate the lists of orbitals that need sic treatment 434 CALL get_qs_env(qs_env, & 435 ks_env=ks_env, & 436 mo_derivs=mo_derivs, & 437 mos=mo_array, & 438 rho=rho, & 439 pw_env=pw_env, & 440 input=input, & 441 cell=cell, & 442 virial=virial) 443 444 CALL qs_rho_get(rho, rho_ao=rho_ao) 445 446 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer) 447 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 448 449 DO i = 1, SIZE(mo_array) !fm->dbcsr 450 IF (mo_array(i)%mo_set%use_mo_coeff_b) THEN !fm->dbcsr 451 CALL copy_dbcsr_to_fm(mo_array(i)%mo_set%mo_coeff_b, & 452 mo_array(i)%mo_set%mo_coeff) !fm->dbcsr 453 ENDIF !fm->dbcsr 454 ENDDO !fm->dbcsr 455 456 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool) 457 458 ! we have two spins 459 CPASSERT(SIZE(mo_array, 1) == 2) 460 ! we want uniform occupations 461 CALL get_mo_set(mo_set=mo_array(1)%mo_set, uniform_occupation=uniform_occupation) 462 CPASSERT(uniform_occupation) 463 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff=mo_coeff, uniform_occupation=uniform_occupation) 464 CPASSERT(uniform_occupation) 465 466 NULLIFY (tmp_dbcsr) 467 CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1)) 468 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr 469 ! 470 NULLIFY (tmp_dbcsr(i)%matrix) 471 CALL dbcsr_init_p(tmp_dbcsr(i)%matrix) 472 CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix) 473 CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp) 474 ENDDO !fm->dbcsr 475 476 k_alpha = 0; k_beta = 0 477 SELECT CASE (dft_control%sic_list_id) 478 CASE (sic_list_all) 479 480 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff=mo_coeff) 481 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha) 482 483 IF (SIZE(mo_array, 1) > 1) THEN 484 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff=mo_coeff) 485 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta) 486 ENDIF 487 488 Norb = k_alpha + k_beta 489 ALLOCATE (sic_orbital_list(3, Norb)) 490 491 iorb = 0 492 DO i = 1, k_alpha 493 iorb = iorb + 1 494 sic_orbital_list(1, iorb) = 1 495 sic_orbital_list(2, iorb) = i 496 sic_orbital_list(3, iorb) = 1 497 ENDDO 498 DO i = 1, k_beta 499 iorb = iorb + 1 500 sic_orbital_list(1, iorb) = 2 501 sic_orbital_list(2, iorb) = i 502 IF (SIZE(mo_derivs, 1) == 1) THEN 503 sic_orbital_list(3, iorb) = 1 504 ELSE 505 sic_orbital_list(3, iorb) = 2 506 ENDIF 507 ENDDO 508 509 CASE (sic_list_unpaired) 510 ! we have two spins 511 CPASSERT(SIZE(mo_array, 1) == 2) 512 ! we have them restricted 513 CPASSERT(SIZE(mo_derivs, 1) == 1) 514 CPASSERT(dft_control%restricted) 515 516 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff=mo_coeff) 517 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha) 518 519 CALL get_mo_set(mo_set=mo_array(2)%mo_set, mo_coeff=mo_coeff) 520 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta) 521 522 Norb = k_alpha - k_beta 523 ALLOCATE (sic_orbital_list(3, Norb)) 524 525 iorb = 0 526 DO i = k_beta + 1, k_alpha 527 iorb = iorb + 1 528 sic_orbital_list(1, iorb) = 1 529 sic_orbital_list(2, iorb) = i 530 ! we are guaranteed to be restricted 531 sic_orbital_list(3, iorb) = 1 532 ENDDO 533 534 CASE DEFAULT 535 CPABORT("") 536 END SELECT 537 538 ! data needed for each of the orbs 539 CALL pw_pool_create_pw(auxbas_pw_pool, orb_rho_r%pw, & 540 use_data=REALDATA3D, & 541 in_space=REALSPACE) 542 CALL pw_pool_create_pw(auxbas_pw_pool, tmp_r%pw, & 543 use_data=REALDATA3D, & 544 in_space=REALSPACE) 545 CALL pw_pool_create_pw(auxbas_pw_pool, orb_rho_g%pw, & 546 use_data=COMPLEXDATA1D, & 547 in_space=RECIPROCALSPACE) 548 CALL pw_pool_create_pw(auxbas_pw_pool, tmp_g%pw, & 549 use_data=COMPLEXDATA1D, & 550 in_space=RECIPROCALSPACE) 551 CALL pw_pool_create_pw(auxbas_pw_pool, work_v_gspace%pw, & 552 use_data=COMPLEXDATA1D, & 553 in_space=RECIPROCALSPACE) 554 CALL pw_pool_create_pw(auxbas_pw_pool, work_v_rspace%pw, & 555 use_data=REALDATA3D, & 556 in_space=REALSPACE) 557 558 ALLOCATE (orb_density_matrix) 559 CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, & 560 name="orb_density_matrix") 561 CALL dbcsr_set(orb_density_matrix, 0.0_dp) 562 orb_density_matrix_p%matrix => orb_density_matrix 563 564 ALLOCATE (orb_h) 565 CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, & 566 name="orb_density_matrix") 567 CALL dbcsr_set(orb_h, 0.0_dp) 568 orb_h_p%matrix => orb_h 569 570 CALL get_mo_set(mo_set=mo_array(1)%mo_set, mo_coeff=mo_coeff) 571 572 CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=1, & 573 template_fmstruct=mo_coeff%matrix_struct) 574 CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v") 575 CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv") 576 CALL cp_fm_struct_release(fm_struct_tmp) 577 578 ALLOCATE (mo_derivs_local(SIZE(mo_array, 1))) 579 DO I = 1, SIZE(mo_array, 1) 580 CALL get_mo_set(mo_set=mo_array(i)%mo_set, mo_coeff=mo_coeff) 581 CALL cp_fm_create(mo_derivs_local(I)%matrix, mo_coeff%matrix_struct) 582 ENDDO 583 584 ALLOCATE (rho_r(2)) 585 rho_r(1)%pw => orb_rho_r%pw 586 rho_r(2)%pw => tmp_r%pw 587 CALL pw_zero(tmp_r%pw) 588 589 ALLOCATE (rho_g(2)) 590 rho_g(1)%pw => orb_rho_g%pw 591 rho_g(2)%pw => tmp_g%pw 592 CALL pw_zero(tmp_g%pw) 593 594 NULLIFY (vxc) 595 ! ALLOCATE(vxc(2),stat=stat) 596 ! CPPostcondition(stat==0,cp_failure_level,routineP,failure) 597 ! CALL pw_pool_create_pw(xc_pw_pool,vxc(1)%pw,& 598 ! in_space=REALSPACE, use_data=REALDATA3D) 599 ! CALL pw_pool_create_pw(xc_pw_pool,vxc(2)%pw,& 600 ! in_space=REALSPACE, use_data=REALDATA3D) 601 602 ! now apply to SIC correction to each selected orbital 603 DO iorb = 1, Norb 604 ! extract the proper orbital from the mo_coeff 605 CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb))%mo_set, mo_coeff=mo_coeff) 606 CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1) 607 608 ! construct the density matrix and the corresponding density 609 CALL dbcsr_set(orb_density_matrix, 0.0_dp) 610 CALL cp_dbcsr_plus_fm_fm_t(orb_density_matrix, matrix_v=matrix_v, ncol=1, & 611 alpha=1.0_dp) 612 613 CALL calculate_rho_elec(matrix_p=orb_density_matrix, & 614 rho=orb_rho_r, rho_gspace=orb_rho_g, total_rho=total_rho, & 615 ks_env=ks_env) 616 617 ! write(*,*) 'Orbital ',sic_orbital_list(1,iorb),sic_orbital_list(2,iorb) 618 ! write(*,*) 'Total orbital rho= ',total_rho 619 620 ! compute the energy functional for this orbital and its derivative 621 622 CALL pw_poisson_solve(poisson_env, orb_rho_g%pw, ener, work_v_gspace%pw) 623 ! no PBC correction is done here, see "calc_v_sic_rspace" for SIC methods 624 ! with PBC aware corrections 625 energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener 626 IF (.NOT. just_energy) THEN 627 CALL pw_transfer(work_v_gspace%pw, work_v_rspace%pw) 628 CALL pw_scale(work_v_rspace%pw, -dft_control%sic_scaling_a*work_v_rspace%pw%pw_grid%dvol) 629 CALL dbcsr_set(orb_h, 0.0_dp) 630 ENDIF 631 632 IF (just_energy) THEN 633 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, & 634 pw_pool=xc_pw_pool) 635 ELSE 636 CPASSERT(.NOT. compute_virial) 637 CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, & 638 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, & 639 pw_pool=xc_pw_pool, compute_virial=compute_virial, virial_xc=virial_xc_tmp) 640 ! add to the existing work_v_rspace 641 work_v_rspace%pw%cr3d = work_v_rspace%pw%cr3d- & 642 dft_control%sic_scaling_b*vxc(1)%pw%pw_grid%dvol*vxc(1)%pw%cr3d 643 END IF 644 energy%exc = energy%exc - dft_control%sic_scaling_b*exc 645 646 IF (.NOT. just_energy) THEN 647 ! note, orb_h (which is being pointed to with orb_h_p) is zeroed above 648 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, & 649 qs_env=qs_env, calculate_forces=calculate_forces) 650 651 ! add this to the mo_derivs 652 CALL cp_dbcsr_sm_fm_multiply(orb_h, matrix_v, matrix_hv, 1) 653 ! silly trick, copy to an array of the right size and add to mo_derivs 654 CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb))%matrix, 0.0_dp) 655 CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb))%matrix, 1, 1, sic_orbital_list(2, iorb)) 656 CALL copy_fm_to_dbcsr(mo_derivs_local(sic_orbital_list(3, iorb))%matrix, & 657 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix) 658 CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, & 659 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp) 660 ! 661 ! need to deallocate vxc 662 CALL pw_pool_give_back_pw(xc_pw_pool, vxc(1)%pw) 663 CALL pw_pool_give_back_pw(xc_pw_pool, vxc(2)%pw) 664 DEALLOCATE (vxc) 665 666 ENDIF 667 668 ENDDO 669 670 CALL pw_pool_give_back_pw(auxbas_pw_pool, orb_rho_r%pw) 671 CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp_r%pw) 672 CALL pw_pool_give_back_pw(auxbas_pw_pool, orb_rho_g%pw) 673 CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp_g%pw) 674 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_v_gspace%pw) 675 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_v_rspace%pw) 676 677 CALL dbcsr_deallocate_matrix(orb_density_matrix) 678 CALL dbcsr_deallocate_matrix(orb_h) 679 CALL cp_fm_release(matrix_v) 680 CALL cp_fm_release(matrix_hv) 681 DO I = 1, SIZE(mo_derivs_local, 1) 682 CALL cp_fm_release(mo_derivs_local(I)%matrix) 683 ENDDO 684 DEALLOCATE (mo_derivs_local) 685 DEALLOCATE (rho_r) 686 DEALLOCATE (rho_g) 687 688 CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr 689 690 CALL timestop(handle) 691 692 END SUBROUTINE sic_explicit_orbitals 693 694! ************************************************************************************************** 695!> \brief do sic calculations on the spin density 696!> \param v_sic_rspace ... 697!> \param energy ... 698!> \param qs_env ... 699!> \param dft_control ... 700!> \param rho ... 701!> \param poisson_env ... 702!> \param just_energy ... 703!> \param calculate_forces ... 704!> \param auxbas_pw_pool ... 705! ************************************************************************************************** 706 SUBROUTINE calc_v_sic_rspace(v_sic_rspace, energy, & 707 qs_env, dft_control, rho, poisson_env, just_energy, & 708 calculate_forces, auxbas_pw_pool) 709 710 TYPE(pw_p_type), INTENT(INOUT) :: v_sic_rspace 711 TYPE(qs_energy_type), POINTER :: energy 712 TYPE(qs_environment_type), POINTER :: qs_env 713 TYPE(dft_control_type), POINTER :: dft_control 714 TYPE(qs_rho_type), POINTER :: rho 715 TYPE(pw_poisson_type), POINTER :: poisson_env 716 LOGICAL, INTENT(IN) :: just_energy, calculate_forces 717 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 718 719 CHARACTER(*), PARAMETER :: routineN = 'calc_v_sic_rspace', routineP = moduleN//':'//routineN 720 721 INTEGER :: i, nelec, nelec_a, nelec_b, nforce 722 REAL(kind=dp) :: ener, full_scaling, scaling 723 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: store_forces 724 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mo_array 725 TYPE(pw_p_type) :: work_rho, work_v 726 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g 727 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 728 729 NULLIFY (mo_array, rho_g) 730 731 IF (dft_control%sic_method_id == sic_none) RETURN 732 IF (dft_control%sic_method_id == sic_eo) RETURN 733 734 IF (dft_control%qs_control%gapw) & 735 CPABORT("sic and GAPW not yet compatible") 736 737 ! OK, right now we like two spins to do sic, could be relaxed for AD 738 CPASSERT(dft_control%nspins == 2) 739 740 CALL pw_pool_create_pw(auxbas_pw_pool, work_rho%pw, & 741 use_data=COMPLEXDATA1D, & 742 in_space=RECIPROCALSPACE) 743 CALL pw_pool_create_pw(auxbas_pw_pool, work_v%pw, & 744 use_data=COMPLEXDATA1D, & 745 in_space=RECIPROCALSPACE) 746 747 CALL qs_rho_get(rho, rho_g=rho_g) 748 749 ! Hartree sic corrections 750 SELECT CASE (dft_control%sic_method_id) 751 CASE (sic_mauri_us, sic_mauri_spz) 752 CALL pw_copy(rho_g(1)%pw, work_rho%pw) 753 CALL pw_axpy(rho_g(2)%pw, work_rho%pw, alpha=-1._dp) 754 CALL pw_poisson_solve(poisson_env, work_rho%pw, ener, work_v%pw) 755 CASE (sic_ad) 756 ! find out how many elecs we have 757 CALL get_qs_env(qs_env, mos=mo_array) 758 CALL get_mo_set(mo_set=mo_array(1)%mo_set, nelectron=nelec_a) 759 CALL get_mo_set(mo_set=mo_array(2)%mo_set, nelectron=nelec_b) 760 nelec = nelec_a + nelec_b 761 CALL pw_copy(rho_g(1)%pw, work_rho%pw) 762 CALL pw_axpy(rho_g(2)%pw, work_rho%pw) 763 scaling = 1.0_dp/REAL(nelec, KIND=dp) 764 CALL pw_scale(work_rho%pw, scaling) 765 CALL pw_poisson_solve(poisson_env, work_rho%pw, ener, work_v%pw) 766 CASE DEFAULT 767 CPABORT("Unknown sic method id") 768 END SELECT 769 770 ! Correct for DDAP charges (if any) 771 ! storing whatever force might be there from previous decoupling 772 IF (calculate_forces) THEN 773 CALL get_qs_env(qs_env=qs_env, force=force) 774 nforce = 0 775 DO i = 1, SIZE(force) 776 nforce = nforce + SIZE(force(i)%ch_pulay, 2) 777 ENDDO 778 ALLOCATE (store_forces(3, nforce)) 779 nforce = 0 780 DO i = 1, SIZE(force) 781 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :) 782 force(i)%ch_pulay(:, :) = 0.0_dp 783 nforce = nforce + SIZE(force(i)%ch_pulay, 2) 784 ENDDO 785 ENDIF 786 787 CALL cp_ddapc_apply_CD(qs_env, & 788 work_rho, & 789 ener, & 790 v_hartree_gspace=work_v, & 791 calculate_forces=calculate_forces, & 792 Itype_of_density="SPIN") 793 794 SELECT CASE (dft_control%sic_method_id) 795 CASE (sic_mauri_us, sic_mauri_spz) 796 full_scaling = -dft_control%sic_scaling_a 797 CASE (sic_ad) 798 full_scaling = -dft_control%sic_scaling_a*nelec 799 CASE DEFAULT 800 CPABORT("Unknown sic method id") 801 END SELECT 802 energy%hartree = energy%hartree + full_scaling*ener 803 804 ! add scaled forces, restoring the old 805 IF (calculate_forces) THEN 806 nforce = 0 807 DO i = 1, SIZE(force) 808 force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + & 809 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) 810 nforce = nforce + SIZE(force(i)%ch_pulay, 2) 811 ENDDO 812 ENDIF 813 814 IF (.NOT. just_energy) THEN 815 CALL pw_pool_create_pw(auxbas_pw_pool, v_sic_rspace%pw, & 816 use_data=REALDATA3D, in_space=REALSPACE) 817 CALL pw_transfer(work_v%pw, v_sic_rspace%pw) 818 ! also take into account the scaling (in addition to the volume element) 819 CALL pw_scale(v_sic_rspace%pw, & 820 dft_control%sic_scaling_a*v_sic_rspace%pw%pw_grid%dvol) 821 ENDIF 822 823 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_rho%pw) 824 CALL pw_pool_give_back_pw(auxbas_pw_pool, work_v%pw) 825 826 END SUBROUTINE calc_v_sic_rspace 827 828! ************************************************************************************************** 829!> \brief ... 830!> \param qs_env ... 831!> \param rho ... 832! ************************************************************************************************** 833 SUBROUTINE print_densities(qs_env, rho) 834 TYPE(qs_environment_type), POINTER :: qs_env 835 TYPE(qs_rho_type), POINTER :: rho 836 837 INTEGER :: img, ispin, n_electrons, output_unit 838 REAL(dp) :: tot1_h, tot1_s, tot_rho_r, trace, & 839 trace_tmp 840 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_arr 841 TYPE(cell_type), POINTER :: cell 842 TYPE(cp_logger_type), POINTER :: logger 843 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, rho_ao 844 TYPE(dft_control_type), POINTER :: dft_control 845 TYPE(qs_charges_type), POINTER :: qs_charges 846 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 847 TYPE(section_vals_type), POINTER :: input, scf_section 848 849 NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, & 850 dft_control, tot_rho_r_arr, rho_ao) 851 852 logger => cp_get_default_logger() 853 854 CALL get_qs_env(qs_env, & 855 qs_kind_set=qs_kind_set, & 856 cell=cell, qs_charges=qs_charges, & 857 input=input, & 858 matrix_s_kp=matrix_s, & 859 dft_control=dft_control) 860 861 CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons) 862 863 scf_section => section_vals_get_subs_vals(input, "DFT%SCF") 864 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", & 865 extension=".scfLog") 866 867 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao) 868 n_electrons = n_electrons - dft_control%charge 869 tot_rho_r = accurate_sum(tot_rho_r_arr) 870 871 trace = 0 872 IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN 873 DO ispin = 1, dft_control%nspins 874 DO img = 1, dft_control%nimages 875 CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp) 876 trace = trace + trace_tmp 877 END DO 878 END DO 879 ENDIF 880 881 IF (output_unit > 0) THEN 882 WRITE (UNIT=output_unit, FMT="(/,T3,A,T41,F20.10)") "Trace(PS):", trace 883 WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") & 884 "Electronic density on regular grids: ", & 885 tot_rho_r, & 886 tot_rho_r + & 887 REAL(n_electrons, dp), & 888 "Core density on regular grids:", & 889 qs_charges%total_rho_core_rspace, & 890 qs_charges%total_rho_core_rspace - REAL(n_electrons + dft_control%charge, dp) 891 END IF 892 IF (dft_control%qs_control%gapw) THEN 893 tot1_h = qs_charges%total_rho1_hard(1) 894 tot1_s = qs_charges%total_rho1_soft(1) 895 DO ispin = 2, dft_control%nspins 896 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin) 897 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin) 898 END DO 899 IF (output_unit > 0) THEN 900 WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") & 901 "Hard and soft densities (Lebedev):", & 902 tot1_h, tot1_s 903 WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") & 904 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", & 905 tot_rho_r + tot1_h - tot1_s, & 906 "Total charge density (r-space): ", & 907 tot_rho_r + tot1_h - tot1_s & 908 + qs_charges%total_rho_core_rspace, & 909 "Total Rho_soft + Rho0_soft (g-space):", & 910 qs_charges%total_rho_gspace 911 END IF 912 qs_charges%background = tot_rho_r + tot1_h - tot1_s + & 913 qs_charges%total_rho_core_rspace 914 ELSE IF (dft_control%qs_control%gapw_xc) THEN 915 tot1_h = qs_charges%total_rho1_hard(1) 916 tot1_s = qs_charges%total_rho1_soft(1) 917 DO ispin = 2, dft_control%nspins 918 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin) 919 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin) 920 END DO 921 IF (output_unit > 0) THEN 922 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T41,2F20.10))") & 923 "Hard and soft densities (Lebedev):", & 924 tot1_h, tot1_s 925 WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") & 926 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", & 927 accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s 928 END IF 929 qs_charges%background = tot_rho_r + & 930 qs_charges%total_rho_core_rspace 931 ELSE 932 IF (output_unit > 0) THEN 933 WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") & 934 "Total charge density on r-space grids: ", & 935 tot_rho_r + & 936 qs_charges%total_rho_core_rspace, & 937 "Total charge density g-space grids: ", & 938 qs_charges%total_rho_gspace 939 END IF 940 qs_charges%background = tot_rho_r + & 941 qs_charges%total_rho_core_rspace 942 END IF 943 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="()") 944 qs_charges%background = qs_charges%background/cell%deth 945 946 CALL cp_print_key_finished_output(output_unit, logger, scf_section, & 947 "PRINT%TOTAL_DENSITIES") 948 949 END SUBROUTINE print_densities 950 951! ************************************************************************************************** 952!> \brief Print detailed energies 953!> 954!> \param qs_env ... 955!> \param dft_control ... 956!> \param input ... 957!> \param energy ... 958!> \param mulliken_order_p ... 959!> \par History 960!> refactoring 04.03.2011 [MI] 961!> \author 962! ************************************************************************************************** 963 SUBROUTINE print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p) 964 965 TYPE(qs_environment_type), POINTER :: qs_env 966 TYPE(dft_control_type), POINTER :: dft_control 967 TYPE(section_vals_type), POINTER :: input 968 TYPE(qs_energy_type), POINTER :: energy 969 REAL(KIND=dp), INTENT(IN) :: mulliken_order_p 970 971 CHARACTER(LEN=*), PARAMETER :: routineN = 'print_detailed_energy', & 972 routineP = moduleN//':'//routineN 973 974 INTEGER :: bc, n, output_unit, psolver 975 REAL(KIND=dp) :: ddapc_order_p, implicit_ps_ehartree, & 976 s2_order_p 977 TYPE(cp_logger_type), POINTER :: logger 978 TYPE(pw_env_type), POINTER :: pw_env 979 980 logger => cp_get_default_logger() 981 982 NULLIFY (pw_env) 983 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 984 psolver = pw_env%poisson_env%parameters%solver 985 986 output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", & 987 extension=".scfLog") 988 IF (output_unit > 0) THEN 989 IF (dft_control%do_admm) THEN 990 WRITE (UNIT=output_unit, FMT="((T3,A,T60,F20.10))") & 991 "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit 992 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN 993 WRITE (UNIT=output_unit, FMT="((T3,A,T60,F20.10))") & 994 "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit 995 END IF 996 END IF 997 IF (dft_control%do_admm) THEN 998 IF (psolver .EQ. pw_poisson_implicit) THEN 999 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree 1000 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 1001 SELECT CASE (bc) 1002 CASE (MIXED_PERIODIC_BC, MIXED_BC) 1003 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1004 "Core Hamiltonian energy: ", energy%core, & 1005 "Hartree energy: ", implicit_ps_ehartree, & 1006 "Electric enthalpy: ", energy%hartree, & 1007 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit 1008 CASE (PERIODIC_BC, NEUMANN_BC) 1009 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1010 "Core Hamiltonian energy: ", energy%core, & 1011 "Hartree energy: ", energy%hartree, & 1012 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit 1013 END SELECT 1014 ELSE 1015 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1016 "Core Hamiltonian energy: ", energy%core, & 1017 "Hartree energy: ", energy%hartree, & 1018 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit 1019 END IF 1020 ELSE 1021!ZMP to print some variables at each step 1022 IF (dft_control%apply_external_density) THEN 1023 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1024 "DOING ZMP CALCULATION FROM EXTERNAL DENSITY " 1025 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1026 "Core Hamiltonian energy: ", energy%core, & 1027 "Hartree energy: ", energy%hartree 1028 ELSE IF (dft_control%apply_external_vxc) THEN 1029 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1030 "DOING ZMP READING EXTERNAL VXC " 1031 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1032 "Core Hamiltonian energy: ", energy%core, & 1033 "Hartree energy: ", energy%hartree 1034 ELSE 1035 IF (psolver .EQ. pw_poisson_implicit) THEN 1036 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree 1037 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition 1038 SELECT CASE (bc) 1039 CASE (MIXED_PERIODIC_BC, MIXED_BC) 1040 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1041 "Core Hamiltonian energy: ", energy%core, & 1042 "Hartree energy: ", implicit_ps_ehartree, & 1043 "Electric enthalpy: ", energy%hartree, & 1044 "Exchange-correlation energy: ", energy%exc 1045 CASE (PERIODIC_BC, NEUMANN_BC) 1046 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1047 "Core Hamiltonian energy: ", energy%core, & 1048 "Hartree energy: ", energy%hartree, & 1049 "Exchange-correlation energy: ", energy%exc 1050 END SELECT 1051 ELSE 1052 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1053 "Core Hamiltonian energy: ", energy%core, & 1054 "Hartree energy: ", energy%hartree, & 1055 "Exchange-correlation energy: ", energy%exc 1056 END IF 1057 END IF 1058 END IF 1059 1060 IF (dft_control%apply_external_density) THEN 1061 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1062 "Integral of the (density * v_xc): ", energy%exc 1063 END IF 1064 1065 IF (energy%e_hartree /= 0.0_dp) & 1066 WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") & 1067 "Coulomb (electron-electron) energy: ", energy%e_hartree 1068 IF (energy%dispersion /= 0.0_dp) & 1069 WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") & 1070 "Dispersion energy: ", energy%dispersion 1071 IF (energy%gcp /= 0.0_dp) & 1072 WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") & 1073 "gCP energy: ", energy%gcp 1074 IF (dft_control%qs_control%gapw) THEN 1075 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1076 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, & 1077 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c 1078 END IF 1079 IF (dft_control%qs_control%gapw_xc) THEN 1080 WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") & 1081 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit 1082 END IF 1083 IF (dft_control%dft_plus_u) THEN 1084 WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") & 1085 "DFT+U energy:", energy%dft_plus_u 1086 END IF 1087 IF (qs_env%qmmm) THEN 1088 WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") & 1089 "QM/MM Electrostatic energy: ", energy%qmmm_el 1090 IF (qs_env%qmmm_env_qm%image_charge) THEN 1091 WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") & 1092 "QM/MM image charge energy: ", energy%image_charge 1093 ENDIF 1094 END IF 1095 IF (dft_control%qs_control%mulliken_restraint) THEN 1096 WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") & 1097 "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken 1098 ENDIF 1099 IF (dft_control%qs_control%ddapc_restraint) THEN 1100 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control) 1101 ddapc_order_p = & 1102 dft_control%qs_control%ddapc_restraint_control(n)%ddapc_restraint_control%ddapc_order_p 1103 WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") & 1104 "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n) 1105 END DO 1106 ENDIF 1107 IF (dft_control%qs_control%s2_restraint) THEN 1108 s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p 1109 WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") & 1110 "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint 1111 ENDIF 1112 1113 END IF ! output_unit 1114 CALL cp_print_key_finished_output(output_unit, logger, input, & 1115 "DFT%SCF%PRINT%DETAILED_ENERGY") 1116 1117 END SUBROUTINE print_detailed_energy 1118 1119! ************************************************************************************************** 1120!> \brief compute matrix_vxc, defined via the potential created by qs_vxc_create 1121!> ignores things like tau functional, gapw, sic, ... 1122!> so only OK for GGA & GPW right now 1123!> \param qs_env ... 1124!> \param v_rspace ... 1125!> \param matrix_vxc ... 1126!> \par History 1127!> created 23.10.2012 [Joost VandeVondele] 1128!> \author 1129! ************************************************************************************************** 1130 SUBROUTINE compute_matrix_vxc(qs_env, v_rspace, matrix_vxc) 1131 TYPE(qs_environment_type), POINTER :: qs_env 1132 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace 1133 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc 1134 1135 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_matrix_vxc', & 1136 routineP = moduleN//':'//routineN 1137 1138 INTEGER :: handle, ispin 1139 LOGICAL :: gapw 1140 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks 1141 TYPE(dft_control_type), POINTER :: dft_control 1142 1143 CALL timeset(routineN, handle) 1144 1145 ! create the matrix using matrix_ks as a template 1146 IF (ASSOCIATED(matrix_vxc)) THEN 1147 CALL dbcsr_deallocate_matrix_set(matrix_vxc) 1148 ENDIF 1149 CALL get_qs_env(qs_env, matrix_ks=matrix_ks) 1150 ALLOCATE (matrix_vxc(SIZE(matrix_ks))) 1151 DO ispin = 1, SIZE(matrix_ks) 1152 NULLIFY (matrix_vxc(ispin)%matrix) 1153 CALL dbcsr_init_p(matrix_vxc(ispin)%matrix) 1154 CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, & 1155 name="Matrix VXC of spin "//cp_to_string(ispin)) 1156 CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp) 1157 ENDDO 1158 1159 ! and integrate 1160 CALL get_qs_env(qs_env, dft_control=dft_control) 1161 gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc 1162 DO ispin = 1, SIZE(matrix_ks) 1163 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), & 1164 hmat=matrix_vxc(ispin), & 1165 qs_env=qs_env, & 1166 calculate_forces=.FALSE., & 1167 gapw=gapw) 1168 ! scale by the volume element... should really become part of integrate_v_rspace 1169 CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw%pw_grid%dvol) 1170 ENDDO 1171 1172 CALL timestop(handle) 1173 1174 END SUBROUTINE compute_matrix_vxc 1175 1176! ************************************************************************************************** 1177!> \brief Sum up all potentials defined on the grid and integrate 1178!> 1179!> \param qs_env ... 1180!> \param ks_matrix ... 1181!> \param rho ... 1182!> \param my_rho ... 1183!> \param vppl_rspace ... 1184!> \param v_rspace_new ... 1185!> \param v_rspace_new_aux_fit ... 1186!> \param v_tau_rspace ... 1187!> \param v_tau_rspace_aux_fit ... 1188!> \param v_efield_rspace ... 1189!> \param v_sic_rspace ... 1190!> \param v_spin_ddapc_rest_r ... 1191!> \param v_sccs_rspace ... 1192!> \param v_rspace_embed ... 1193!> \param cdft_control ... 1194!> \param calculate_forces ... 1195!> \par History 1196!> - refactoring 04.03.2011 [MI] 1197!> - SCCS implementation (16.10.2013,MK) 1198!> \author 1199! ************************************************************************************************** 1200 SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, & 1201 vppl_rspace, v_rspace_new, & 1202 v_rspace_new_aux_fit, v_tau_rspace, & 1203 v_tau_rspace_aux_fit, v_efield_rspace, & 1204 v_sic_rspace, v_spin_ddapc_rest_r, & 1205 v_sccs_rspace, v_rspace_embed, cdft_control, & 1206 calculate_forces) 1207 1208 TYPE(qs_environment_type), POINTER :: qs_env 1209 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 1210 TYPE(qs_rho_type), POINTER :: rho 1211 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: my_rho 1212 TYPE(pw_p_type), POINTER :: vppl_rspace 1213 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, v_rspace_new_aux_fit, & 1214 v_tau_rspace, v_tau_rspace_aux_fit 1215 TYPE(pw_p_type) :: v_efield_rspace, v_sic_rspace, & 1216 v_spin_ddapc_rest_r, v_sccs_rspace 1217 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_embed 1218 TYPE(cdft_control_type), POINTER :: cdft_control 1219 LOGICAL, INTENT(in) :: calculate_forces 1220 1221 CHARACTER(LEN=*), PARAMETER :: routineN = 'sum_up_and_integrate', & 1222 routineP = moduleN//':'//routineN 1223 1224 CHARACTER(LEN=default_string_length) :: basis_type 1225 INTEGER :: handle, igroup, ikind, ispin, nkind, & 1226 nspins 1227 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 1228 LOGICAL :: do_ppl, gapw, gapw_xc, lrigpw, rigpw 1229 REAL(dp) :: dvol, sign 1230 TYPE(admm_type), POINTER :: admm_env 1231 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1232 TYPE(cp_para_env_type), POINTER :: para_env 1233 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_ks_aux_fit, & 1234 matrix_ks_aux_fit_dft, rho_ao, & 1235 rho_ao_aux, rho_ao_nokp, smat 1236 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp 1237 TYPE(dft_control_type), POINTER :: dft_control 1238 TYPE(kpoint_type), POINTER :: kpoints 1239 TYPE(lri_density_type), POINTER :: lri_density 1240 TYPE(lri_environment_type), POINTER :: lri_env 1241 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int 1242 TYPE(pw_env_type), POINTER :: pw_env 1243 TYPE(pw_p_type) :: v_rspace 1244 TYPE(pw_p_type), POINTER :: vee 1245 TYPE(pw_poisson_type), POINTER :: poisson_env 1246 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1247 TYPE(qs_ks_env_type), POINTER :: ks_env 1248 TYPE(qs_rho_type), POINTER :: rho_aux_fit 1249 TYPE(task_list_type), POINTER :: task_list 1250 1251 CALL timeset(routineN, handle) 1252 1253 NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, & 1254 v_rspace%pw, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, & 1255 ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, & 1256 rho_ao_nokp, ks_env, admm_env, task_list) 1257 1258 CALL get_qs_env(qs_env, & 1259 dft_control=dft_control, & 1260 pw_env=pw_env, & 1261 matrix_ks_aux_fit=matrix_ks_aux_fit, & 1262 matrix_ks_aux_fit_dft=matrix_ks_aux_fit_dft, & 1263 v_hartree_rspace=v_rspace%pw, & 1264 rho_aux_fit=rho_aux_fit, & 1265 vee=vee) 1266 1267 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 1268 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux) 1269 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool) 1270 gapw = dft_control%qs_control%gapw 1271 gapw_xc = dft_control%qs_control%gapw_xc 1272 do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid 1273 1274 rigpw = dft_control%qs_control%rigpw 1275 lrigpw = dft_control%qs_control%lrigpw 1276 IF (lrigpw .OR. rigpw) THEN 1277 CALL get_qs_env(qs_env, & 1278 lri_env=lri_env, & 1279 lri_density=lri_density, & 1280 atomic_kind_set=atomic_kind_set) 1281 ENDIF 1282 1283 nspins = dft_control%nspins 1284 1285 ! sum up potentials and integrate 1286 IF (ASSOCIATED(v_rspace_new)) THEN 1287 DO ispin = 1, nspins 1288 IF (gapw_xc) THEN 1289 ! SIC not implemented (or at least not tested) 1290 CPASSERT(dft_control%sic_method_id == sic_none) 1291 !Only the xc potential, because it has to be integrated with the soft basis 1292 v_rspace_new(ispin)%pw%cr3d = & 1293 v_rspace_new(ispin)%pw%pw_grid%dvol* & 1294 v_rspace_new(ispin)%pw%cr3d 1295 1296 ! add the xc part due to v_rspace soft 1297 rho_ao => rho_ao_kp(ispin, :) 1298 ksmat => ks_matrix(ispin, :) 1299 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1300 pmat_kp=rho_ao, hmat_kp=ksmat, & 1301 qs_env=qs_env, & 1302 calculate_forces=calculate_forces, & 1303 gapw=gapw_xc) 1304 1305 ! Now the Hartree potential to be integrated with the full basis 1306 v_rspace_new(ispin)%pw%cr3d = v_rspace%pw%cr3d 1307 ELSE 1308 ! Add v_hartree + v_xc = v_rspace_new 1309 v_rspace_new(ispin)%pw%cr3d = & 1310 v_rspace_new(ispin)%pw%pw_grid%dvol* & 1311 v_rspace_new(ispin)%pw%cr3d+v_rspace%pw%cr3d 1312 END IF ! gapw_xc 1313 IF (dft_control%qs_control%ddapc_explicit_potential) THEN 1314 IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN 1315 IF (ispin == 1) THEN 1316 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d & 1317 +v_spin_ddapc_rest_r%pw%cr3d 1318 ELSE 1319 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d & 1320 -v_spin_ddapc_rest_r%pw%cr3d 1321 ENDIF 1322 ELSE 1323 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d & 1324 +v_spin_ddapc_rest_r%pw%cr3d 1325 END IF 1326 END IF 1327 ! CDFT constraint contribution 1328 IF (dft_control%qs_control%cdft) THEN 1329 DO igroup = 1, SIZE(cdft_control%group) 1330 SELECT CASE (cdft_control%group(igroup)%constraint_type) 1331 CASE (cdft_charge_constraint) 1332 sign = 1.0_dp 1333 CASE (cdft_magnetization_constraint) 1334 IF (ispin == 1) THEN 1335 sign = 1.0_dp 1336 ELSE 1337 sign = -1.0_dp 1338 END IF 1339 CASE (cdft_alpha_constraint) 1340 sign = 1.0_dp 1341 IF (ispin == 2) CYCLE 1342 CASE (cdft_beta_constraint) 1343 sign = 1.0_dp 1344 IF (ispin == 1) CYCLE 1345 CASE DEFAULT 1346 CPABORT("Unknown constraint type.") 1347 END SELECT 1348 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d & 1349 +sign*cdft_control%group(igroup)%weight%pw%cr3d* & 1350 cdft_control%strength(igroup) 1351 END DO 1352 END IF 1353 ! The efield contribution 1354 IF (dft_control%apply_efield_field) THEN 1355 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1356 v_efield_rspace%pw%cr3d 1357 END IF 1358 ! functional derivative of the Hartree energy wrt the density in the presence of dielectric 1359 ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function 1360 ! of the charge density 1361 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN 1362 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol 1363 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1364 poisson_env%implicit_env%v_eps%cr3d*dvol 1365 END IF 1366 ! Add SCCS contribution 1367 IF (dft_control%do_sccs) THEN 1368 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1369 v_sccs_rspace%pw%cr3d 1370 END IF 1371 ! External electrostatic potential 1372 IF (dft_control%apply_external_potential) THEN 1373 CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), & 1374 v_qmmm=vee, scale=-1.0_dp) 1375 END IF 1376 IF (do_ppl) THEN 1377 CPASSERT(.NOT. gapw) 1378 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1379 vppl_rspace%pw%cr3d*vppl_rspace%pw%pw_grid%dvol 1380 END IF 1381 ! the electrostatic sic contribution 1382 SELECT CASE (dft_control%sic_method_id) 1383 CASE (sic_none) 1384 ! 1385 CASE (sic_mauri_us, sic_mauri_spz) 1386 IF (ispin == 1) THEN 1387 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d- & 1388 v_sic_rspace%pw%cr3d 1389 ELSE 1390 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1391 v_sic_rspace%pw%cr3d 1392 ENDIF 1393 CASE (sic_ad) 1394 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d-v_sic_rspace%pw%cr3d 1395 CASE (sic_eo) 1396 ! NOTHING TO BE DONE 1397 END SELECT 1398 ! DFT embedding 1399 IF (dft_control%apply_embed_pot) THEN 1400 v_rspace_new(ispin)%pw%cr3d = v_rspace_new(ispin)%pw%cr3d+ & 1401 v_rspace_embed(ispin)%pw%pw_grid%dvol* & 1402 v_rspace_embed(ispin)%pw%cr3d 1403 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_embed(ispin)%pw) 1404 ENDIF 1405 IF (lrigpw) THEN 1406 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 1407 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env) 1408 DO ikind = 1, nkind 1409 lri_v_int(ikind)%v_int = 0.0_dp 1410 END DO 1411 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, & 1412 lri_v_int, calculate_forces, "LRI_AUX") 1413 DO ikind = 1, nkind 1414 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 1415 END DO 1416 IF (lri_env%exact_1c_terms) THEN 1417 rho_ao => my_rho(ispin, :) 1418 ksmat => ks_matrix(ispin, :) 1419 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, & 1420 rho_ao(1)%matrix, qs_env, & 1421 calculate_forces, "ORB") 1422 END IF 1423 IF (lri_env%ppl_ri) THEN 1424 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces) 1425 END IF 1426 ELSEIF (rigpw) THEN 1427 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 1428 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env) 1429 DO ikind = 1, nkind 1430 lri_v_int(ikind)%v_int = 0.0_dp 1431 END DO 1432 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, & 1433 lri_v_int, calculate_forces, "RI_HXC") 1434 DO ikind = 1, nkind 1435 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 1436 END DO 1437 ELSE 1438 rho_ao => my_rho(ispin, :) 1439 ksmat => ks_matrix(ispin, :) 1440 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 1441 pmat_kp=rho_ao, hmat_kp=ksmat, & 1442 qs_env=qs_env, & 1443 calculate_forces=calculate_forces, & 1444 gapw=gapw) 1445 ENDIF 1446 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(ispin)%pw) 1447 END DO ! ispin 1448 1449 SELECT CASE (dft_control%sic_method_id) 1450 CASE (sic_none) 1451 CASE (sic_mauri_us, sic_mauri_spz, sic_ad) 1452 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_sic_rspace%pw) 1453 END SELECT 1454 DEALLOCATE (v_rspace_new) 1455 1456 ELSE 1457 ! not implemented (or at least not tested) 1458 CPASSERT(dft_control%sic_method_id == sic_none) 1459 CPASSERT(.NOT. dft_control%qs_control%ddapc_restraint_is_spin) 1460 DO ispin = 1, nspins 1461 ! the efield contribution 1462 IF (dft_control%apply_efield_field) THEN 1463 v_rspace%pw%cr3d = v_rspace%pw%cr3d+v_efield_rspace%pw%cr3d 1464 END IF 1465 ! extra contribution attributed to the dependency of the dielectric constant to the charge density 1466 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN 1467 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol 1468 v_rspace%pw%cr3d = v_rspace%pw%cr3d+poisson_env%implicit_env%v_eps%cr3d*dvol 1469 END IF 1470 ! Add SCCS contribution 1471 IF (dft_control%do_sccs) THEN 1472 v_rspace%pw%cr3d = v_rspace%pw%cr3d+v_sccs_rspace%pw%cr3d 1473 END IF 1474 ! DFT embedding 1475 IF (dft_control%apply_embed_pot) THEN 1476 v_rspace%pw%cr3d = v_rspace%pw%cr3d+ & 1477 v_rspace_embed(ispin)%pw%pw_grid%dvol* & 1478 v_rspace_embed(ispin)%pw%cr3d 1479 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_embed(ispin)%pw) 1480 ENDIF 1481 IF (lrigpw) THEN 1482 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 1483 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env) 1484 DO ikind = 1, nkind 1485 lri_v_int(ikind)%v_int = 0.0_dp 1486 END DO 1487 CALL integrate_v_rspace_one_center(v_rspace, qs_env, & 1488 lri_v_int, calculate_forces, "LRI_AUX") 1489 DO ikind = 1, nkind 1490 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 1491 END DO 1492 IF (lri_env%exact_1c_terms) THEN 1493 rho_ao => my_rho(ispin, :) 1494 ksmat => ks_matrix(ispin, :) 1495 CALL integrate_v_rspace_diagonal(v_rspace, ksmat(1)%matrix, & 1496 rho_ao(1)%matrix, qs_env, & 1497 calculate_forces, "ORB") 1498 END IF 1499 IF (lri_env%ppl_ri) THEN 1500 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces) 1501 END IF 1502 ELSEIF (rigpw) THEN 1503 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 1504 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env) 1505 DO ikind = 1, nkind 1506 lri_v_int(ikind)%v_int = 0.0_dp 1507 END DO 1508 CALL integrate_v_rspace_one_center(v_rspace, qs_env, & 1509 lri_v_int, calculate_forces, "RI_HXC") 1510 DO ikind = 1, nkind 1511 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 1512 END DO 1513 ELSE 1514 rho_ao => my_rho(ispin, :) 1515 ksmat => ks_matrix(ispin, :) 1516 CALL integrate_v_rspace(v_rspace=v_rspace, & 1517 pmat_kp=rho_ao, & 1518 hmat_kp=ksmat, & 1519 qs_env=qs_env, & 1520 calculate_forces=calculate_forces, & 1521 gapw=gapw) 1522 ENDIF 1523 END DO 1524 END IF ! ASSOCIATED(v_rspace_new) 1525 1526 ! **** LRIGPW: KS matrix from integrated potential 1527 IF (lrigpw) THEN 1528 CALL get_qs_env(qs_env, ks_env=ks_env) 1529 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 1530 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 1531 DO ispin = 1, nspins 1532 ksmat => ks_matrix(ispin, :) 1533 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, & 1534 cell_to_index=cell_to_index) 1535 ENDDO 1536 IF (calculate_forces) THEN 1537 CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set) 1538 ENDIF 1539 ELSEIF (rigpw) THEN 1540 CALL get_qs_env(qs_env, matrix_s=smat) 1541 DO ispin = 1, nspins 1542 CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, & 1543 smat(1)%matrix, atomic_kind_set, ispin) 1544 ENDDO 1545 IF (calculate_forces) THEN 1546 rho_ao_nokp => rho_ao_kp(:, 1) 1547 CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set) 1548 ENDIF 1549 ENDIF 1550 1551 IF (ASSOCIATED(v_tau_rspace)) THEN 1552 IF (lrigpw .OR. rigpw) THEN 1553 CPABORT("LRIGPW/RIGPW not implemented for meta-GGAs") 1554 ENDIF 1555 DO ispin = 1, nspins 1556 v_tau_rspace(ispin)%pw%cr3d = & 1557 v_tau_rspace(ispin)%pw%pw_grid%dvol* & 1558 v_tau_rspace(ispin)%pw%cr3d 1559 1560 rho_ao => rho_ao_kp(ispin, :) 1561 ksmat => ks_matrix(ispin, :) 1562 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), & 1563 pmat_kp=rho_ao, hmat_kp=ksmat, & 1564 qs_env=qs_env, & 1565 calculate_forces=calculate_forces, compute_tau=.TRUE., & 1566 gapw=gapw) 1567 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw) 1568 END DO 1569 DEALLOCATE (v_tau_rspace) 1570 ENDIF 1571 1572 ! Add contributions from ADMM if requested 1573 IF (dft_control%do_admm) THEN 1574 CALL get_qs_env(qs_env, admm_env=admm_env) 1575 IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN 1576 DO ispin = 1, nspins 1577 ! Calculate the xc potential 1578 v_rspace_new_aux_fit(ispin)%pw%cr3d = & 1579 v_rspace_new_aux_fit(ispin)%pw%pw_grid%dvol* & 1580 v_rspace_new_aux_fit(ispin)%pw%cr3d 1581 1582 ! set matrix_ks_aux_fit_dft = matrix_ks_aux_fit(k_HF) 1583 CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin)%matrix, matrix_ks_aux_fit(ispin)%matrix, & 1584 name="DFT exch. part of matrix_ks_aux_fit") 1585 1586 ! Add potential to ks_matrix aux_fit, skip integration if no DFT correction 1587 1588 IF (admm_env%aux_exch_func .NE. do_admm_aux_exch_func_none) THEN 1589 1590 !GPW by default. IF GAPW, then take relevant task list and basis 1591 CALL get_qs_env(qs_env, task_list_aux_fit=task_list) 1592 basis_type = "AUX_FIT" 1593 IF (admm_env%do_gapw) THEN 1594 task_list => admm_env%admm_gapw_env%task_list 1595 basis_type = "AUX_FIT_SOFT" 1596 END IF 1597 1598 CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), & 1599 pmat=rho_ao_aux(ispin), & 1600 hmat=matrix_ks_aux_fit(ispin), & 1601 qs_env=qs_env, & 1602 calculate_forces=calculate_forces, & 1603 force_adm=.TRUE., & 1604 ispin=ispin, & 1605 gapw=.FALSE., & !even if actual GAPW calculation, want to use AUX_FIT_SOFT 1606 basis_type=basis_type, & 1607 task_list_external=task_list) 1608 END IF 1609 1610 ! matrix_ks_aux_fit_dft(x_DFT)=matrix_ks_aux_fit_dft(old,k_HF)-matrix_ks_aux_fit(k_HF-x_DFT) 1611 CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin)%matrix, & 1612 matrix_ks_aux_fit(ispin)%matrix, 1.0_dp, -1.0_dp) 1613 1614 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new_aux_fit(ispin)%pw) 1615 END DO 1616 DEALLOCATE (v_rspace_new_aux_fit) 1617 END IF 1618 ! Clean up v_tau_rspace_aux_fit, which is actually not needed 1619 IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN 1620 DO ispin = 1, nspins 1621 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace_aux_fit(ispin)%pw) 1622 END DO 1623 DEALLOCATE (v_tau_rspace_aux_fit) 1624 END IF 1625 END IF 1626 1627 IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed) 1628 1629 CALL timestop(handle) 1630 1631 END SUBROUTINE sum_up_and_integrate 1632 1633!************************************************************************** 1634!> \brief Calculate the ZMP potential and energy as in Zhao, Morrison Parr 1635!> PRA 50i, 2138 (1994) 1636!> V_c^\lambda defined as int_rho-rho_0/r-r' or rho-rho_0 times a Lagrange 1637!> multiplier, plus Fermi-Amaldi potential that should give the V_xc in the 1638!> limit \lambda --> \infty 1639!> 1640!> \param qs_env ... 1641!> \param v_rspace_new ... 1642!> \param rho ... 1643!> \param exc ... 1644!> \author D. Varsano [daniele.varsano@nano.cnr.it] 1645! ************************************************************************************************** 1646 SUBROUTINE calculate_zmp_potential(qs_env, v_rspace_new, rho, exc) 1647 1648 TYPE(qs_environment_type), POINTER :: qs_env 1649 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new 1650 TYPE(qs_rho_type), POINTER :: rho 1651 REAL(KIND=dp) :: exc 1652 1653 CHARACTER(*), PARAMETER :: routineN = 'calculate_zmp_potential', & 1654 routineP = moduleN//':'//routineN 1655 1656 INTEGER :: handle, my_val, nelectron, nspins, stat 1657 INTEGER, DIMENSION(2) :: nelectron_spin 1658 LOGICAL :: do_zmp_read, fermi_amaldi 1659 REAL(KIND=dp) :: factor, lambda, total_rho 1660 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_ext_r 1661 TYPE(dft_control_type), POINTER :: dft_control 1662 TYPE(pw_env_type), POINTER :: pw_env 1663 TYPE(pw_p_type) :: rho_eff_gspace, v_xc_gspace, v_xc_rspace 1664 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_ext_g, rho_g, rho_r 1665 TYPE(pw_poisson_type), POINTER :: poisson_env 1666 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1667 TYPE(qs_ks_env_type), POINTER :: ks_env 1668 TYPE(section_vals_type), POINTER :: ext_den_section, input 1669 1670!, v_h_gspace, & 1671 1672 CALL timeset(routineN, handle) 1673 NULLIFY (auxbas_pw_pool) 1674 NULLIFY (pw_env) 1675 NULLIFY (poisson_env) 1676 NULLIFY (v_rspace_new) 1677 NULLIFY (dft_control) 1678 NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g) 1679 CALL get_qs_env(qs_env=qs_env, & 1680 pw_env=pw_env, & 1681 ks_env=ks_env, & 1682 rho=rho, & 1683 input=input, & 1684 nelectron_spin=nelectron_spin, & 1685 dft_control=dft_control) 1686 CALL pw_env_get(pw_env=pw_env, & 1687 auxbas_pw_pool=auxbas_pw_pool, & 1688 poisson_env=poisson_env) 1689 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g) 1690 nspins = 1 1691 ALLOCATE (v_rspace_new(nspins), stat=stat) 1692 CPASSERT(stat == 0) 1693 CALL pw_pool_create_pw(pool=auxbas_pw_pool, pw=v_rspace_new(1)%pw, & 1694 use_data=REALDATA3D, in_space=REALSPACE) 1695 CALL pw_pool_create_pw(pool=auxbas_pw_pool, pw=v_xc_rspace%pw, & 1696 use_data=REALDATA3D, in_space=REALSPACE) 1697 1698 CALL pw_zero(v_rspace_new(1)%pw) 1699 do_zmp_read = dft_control%apply_external_vxc 1700 IF (do_zmp_read) THEN 1701 CALL pw_copy(qs_env%external_vxc%pw, v_rspace_new(1)%pw) 1702 exc = 0.0_dp 1703 exc = accurate_sum(v_rspace_new(1)%pw%cr3d*rho_r(1)%pw%cr3d)* & 1704 v_rspace_new(1)%pw%pw_grid%dvol 1705 ELSE 1706 CALL pw_pool_create_pw(pool=auxbas_pw_pool, & 1707 pw=rho_eff_gspace%pw, & 1708 use_data=COMPLEXDATA1D, & 1709 in_space=RECIPROCALSPACE) 1710 CALL pw_pool_create_pw(pool=auxbas_pw_pool, & 1711 pw=v_xc_gspace%pw, & 1712 use_data=COMPLEXDATA1D, & 1713 in_space=RECIPROCALSPACE) 1714 CALL pw_zero(rho_eff_gspace%pw) 1715 CALL pw_zero(v_xc_gspace%pw) 1716 CALL pw_zero(v_xc_rspace%pw) 1717 factor = pw_integrate_function(rho_g(1)%pw) 1718 CALL qs_rho_get(qs_env%rho_external, & 1719 rho_g=rho_ext_g, & 1720 tot_rho_r=tot_rho_ext_r) 1721 factor = tot_rho_ext_r(1)/factor 1722 1723 CALL pw_axpy(rho_g(1)%pw, rho_eff_gspace%pw, alpha=factor) 1724 CALL pw_axpy(rho_ext_g(1)%pw, rho_eff_gspace%pw, alpha=-1.0_dp) 1725 total_rho = pw_integrate_function(rho_eff_gspace%pw, isign=1) 1726 ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY") 1727 CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda) 1728 CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val) 1729 CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi) 1730 1731 CALL pw_scale(rho_eff_gspace%pw, a=lambda) 1732 nelectron = nelectron_spin(1) 1733 factor = -1.0_dp/nelectron 1734 CALL pw_axpy(rho_g(1)%pw, rho_eff_gspace%pw, alpha=factor) 1735 1736 CALL pw_poisson_solve(poisson_env, rho_eff_gspace%pw, vhartree=v_xc_gspace%pw) 1737 CALL pw_transfer(v_xc_gspace%pw, v_rspace_new(1)%pw) 1738 CALL pw_copy(v_rspace_new(1)%pw, v_xc_rspace%pw) 1739 1740 exc = 0.0_dp 1741 exc = accurate_sum(v_rspace_new(1)%pw%cr3d*rho_r(1)%pw%cr3d)* & 1742 v_rspace_new(1)%pw%pw_grid%dvol 1743 IF (v_rspace_new(1)%pw%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN 1744 CALL mp_sum(exc, v_rspace_new(1)%pw%pw_grid%para%group) 1745 END IF 1746 1747! IF (v_rspace_new(1)%pw%pw_grid%para%my_pos==0) & 1748! WRITE(*,FMT="(T3,A,T61,F20.10)") "ZMP| Integral of (electron density)*(v_xc): " , exc 1749!Note that this is not the xc energy but \int(\rho*v_xc) 1750!Vxc---> v_rspace_new 1751!Exc---> energy%exc 1752 CALL pw_pool_give_back_pw(auxbas_pw_pool, & 1753 rho_eff_gspace%pw) 1754 CALL pw_pool_give_back_pw(auxbas_pw_pool, & 1755 v_xc_gspace%pw) 1756 ENDIF 1757 1758 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_xc_rspace%pw) 1759 1760 CALL timestop(handle) 1761 1762 END SUBROUTINE calculate_zmp_potential 1763 1764! ************************************************************************************************** 1765!> \brief ... 1766!> \param qs_env ... 1767!> \param rho ... 1768!> \param v_rspace_embed ... 1769!> \param dft_control ... 1770!> \param embed_corr ... 1771!> \param just_energy ... 1772! ************************************************************************************************** 1773 SUBROUTINE get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, & 1774 just_energy) 1775 TYPE(qs_environment_type), POINTER :: qs_env 1776 TYPE(qs_rho_type), POINTER :: rho 1777 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_embed 1778 TYPE(dft_control_type), POINTER :: dft_control 1779 REAL(KIND=dp) :: embed_corr 1780 LOGICAL :: just_energy 1781 1782 CHARACTER(*), PARAMETER :: routineN = 'get_embed_potential_energy', & 1783 routineP = moduleN//':'//routineN 1784 1785 INTEGER :: handle, ispin, ns1, ns2 1786 REAL(KIND=dp) :: embed_corr_local 1787 TYPE(pw_env_type), POINTER :: pw_env 1788 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 1789 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 1790 1791 CALL timeset(routineN, handle) 1792 1793 NULLIFY (auxbas_pw_pool) 1794 NULLIFY (pw_env) 1795 NULLIFY (rho_r) 1796 CALL get_qs_env(qs_env=qs_env, & 1797 pw_env=pw_env, & 1798 rho=rho) 1799 CALL pw_env_get(pw_env=pw_env, & 1800 auxbas_pw_pool=auxbas_pw_pool) 1801 CALL qs_rho_get(rho, rho_r=rho_r) 1802 ALLOCATE (v_rspace_embed(dft_control%nspins)) 1803 1804 embed_corr = 0.0_dp 1805 1806 DO ispin = 1, dft_control%nspins 1807 CALL pw_pool_create_pw(pool=auxbas_pw_pool, pw=v_rspace_embed(ispin)%pw, & 1808 use_data=REALDATA3D, in_space=REALSPACE) 1809 CALL pw_zero(v_rspace_embed(ispin)%pw) 1810 ns1 = SIZE(qs_env%embed_pot%pw%cr3d) 1811 ns2 = SIZE(v_rspace_embed(ispin)%pw%cr3d) 1812 IF (ns1 .NE. ns2) CPABORT("Subsystem grids must be identical") 1813 1814 v_rspace_embed(ispin)%pw%cr3d(:, :, :) = qs_env%embed_pot%pw%cr3d 1815 embed_corr_local = 0.0_dp 1816 1817 ! Spin embedding potential in open-shell case 1818 IF (dft_control%nspins .EQ. 2) THEN 1819 ns1 = SIZE(qs_env%spin_embed_pot%pw%cr3d) 1820 IF (ns1 .NE. ns2) CPABORT("Subsystem grids must be identical") 1821 IF (ispin .EQ. 1) v_rspace_embed(ispin)%pw%cr3d(:, :, :) = & 1822 v_rspace_embed(ispin)%pw%cr3d(:, :, :) + qs_env%spin_embed_pot%pw%cr3d 1823 IF (ispin .EQ. 2) v_rspace_embed(ispin)%pw%cr3d(:, :, :) = & 1824 v_rspace_embed(ispin)%pw%cr3d(:, :, :) - qs_env%spin_embed_pot%pw%cr3d 1825 ENDIF 1826 ! Integrate the density*potential 1827 embed_corr_local = accurate_sum(v_rspace_embed(ispin)%pw%cr3d*rho_r(ispin)%pw%cr3d)* & 1828 v_rspace_embed(ispin)%pw%pw_grid%dvol 1829 1830 IF (v_rspace_embed(ispin)%pw%pw_grid%para%mode == PW_MODE_DISTRIBUTED) THEN 1831 CALL mp_sum(embed_corr_local, v_rspace_embed(ispin)%pw%pw_grid%para%group) 1832 END IF 1833 1834 embed_corr = embed_corr + embed_corr_local 1835 1836 ENDDO 1837 1838 ! If only energy requiested we delete the potential 1839 IF (just_energy) THEN 1840 DO ispin = 1, dft_control%nspins 1841 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_embed(ispin)%pw) 1842 ENDDO 1843 DEALLOCATE (v_rspace_embed) 1844 ENDIF 1845 1846 CALL timestop(handle) 1847 1848 END SUBROUTINE get_embed_potential_energy 1849 1850END MODULE qs_ks_utils 1851