1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief 8!> 9!> 10!> \par History 11!> refactoring 03-2011 [MI] 12!> \author MI 13! ************************************************************************************************** 14MODULE qs_vxc 15 16 USE cell_types, ONLY: cell_type 17 USE cp_control_types, ONLY: dft_control_type 18 USE cp_para_types, ONLY: cp_para_env_type 19 USE input_constants, ONLY: sic_ad,& 20 sic_eo,& 21 sic_mauri_spz,& 22 sic_mauri_us,& 23 sic_none,& 24 xc_none,& 25 xc_vdw_fun_nonloc 26 USE input_section_types, ONLY: section_vals_type,& 27 section_vals_val_get 28 USE kinds, ONLY: dp 29 USE pw_env_types, ONLY: pw_env_get,& 30 pw_env_type 31 USE pw_grids, ONLY: pw_grid_compare 32 USE pw_methods, ONLY: pw_axpy,& 33 pw_copy,& 34 pw_multiply,& 35 pw_scale,& 36 pw_transfer,& 37 pw_zero 38 USE pw_pool_types, ONLY: pw_pool_create_pw,& 39 pw_pool_give_back_pw,& 40 pw_pool_type 41 USE pw_types, ONLY: COMPLEXDATA1D,& 42 REALDATA3D,& 43 REALSPACE,& 44 RECIPROCALSPACE,& 45 pw_p_type,& 46 pw_release,& 47 pw_type 48 USE qs_dispersion_nonloc, ONLY: calculate_dispersion_nonloc 49 USE qs_dispersion_types, ONLY: qs_dispersion_type 50 USE qs_ks_types, ONLY: get_ks_env,& 51 qs_ks_env_type 52 USE qs_rho_types, ONLY: qs_rho_get,& 53 qs_rho_type 54 USE virial_types, ONLY: virial_type 55 USE xc, ONLY: xc_exc_calc,& 56 xc_vxc_pw_create1 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 61 PRIVATE 62 63 ! *** Public subroutines *** 64 PUBLIC :: qs_vxc_create, qs_xc_density 65 66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc' 67 68CONTAINS 69 70! ************************************************************************************************** 71!> \brief calculates and allocates the xc potential, already reducing it to 72!> the dependence on rho and the one on tau 73!> \param ks_env to get all the needed things 74!> \param rho_struct density for which v_xc is calculated 75!> \param xc_section ... 76!> \param vxc_rho will contain the v_xc part that depend on rho 77!> (if one of the choosen xc functionals has it it is allocated and you 78!> are responsible for it) 79!> \param vxc_tau will contain the kinetic (tau) part of v_xc 80!> (if one of the choosen xc functionals has it it is allocated and you 81!> are responsible for it) 82!> \param exc ... 83!> \param just_energy if true calculates just the energy, and does not 84!> allocate v_*_rspace 85!> \param edisp ... 86!> \param dispersion_env ... 87!> \param adiabatic_rescale_factor ... 88!> \param pw_env_external external plane wave environment 89!> \par History 90!> - 05.2002 modified to use the mp_allgather function each pe 91!> computes only part of the grid and this is broadcasted to all 92!> instead of summed. 93!> This scales significantly better (e.g. factor 3 on 12 cpus 94!> 32 H2O) [Joost VdV] 95!> - moved to qs_ks_methods [fawzi] 96!> - sic alterations [Joost VandeVondele] 97!> \author Fawzi Mohamed 98! ************************************************************************************************** 99 SUBROUTINE qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, & 100 just_energy, edisp, dispersion_env, adiabatic_rescale_factor, & 101 pw_env_external) 102 103 TYPE(qs_ks_env_type), POINTER :: ks_env 104 TYPE(qs_rho_type), POINTER :: rho_struct 105 TYPE(section_vals_type), POINTER :: xc_section 106 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau 107 REAL(KIND=dp), INTENT(out) :: exc 108 LOGICAL, INTENT(in), OPTIONAL :: just_energy 109 REAL(KIND=dp), INTENT(out), OPTIONAL :: edisp 110 TYPE(qs_dispersion_type), OPTIONAL, POINTER :: dispersion_env 111 REAL(KIND=dp), INTENT(in), OPTIONAL :: adiabatic_rescale_factor 112 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external 113 114 CHARACTER(len=*), PARAMETER :: routineN = 'qs_vxc_create', routineP = moduleN//':'//routineN 115 116 INTEGER :: handle, ispin, myfun, nelec_spin(2), vdw 117 LOGICAL :: compute_virial, do_adiabatic_rescaling, my_just_energy, rho_g_valid, & 118 sic_scaling_b_zero, tau_r_valid, uf_grid, vdW_nl 119 REAL(KIND=dp) :: exc_m, factor, & 120 my_adiabatic_rescale_factor, & 121 my_scaling, nelec_s_inv 122 REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp 123 TYPE(cell_type), POINTER :: cell 124 TYPE(cp_para_env_type), POINTER :: para_env 125 TYPE(dft_control_type), POINTER :: dft_control 126 TYPE(pw_env_type), POINTER :: pw_env 127 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, rho_g, & 128 rho_m_gspace, rho_m_rspace, rho_r, & 129 rho_struct_g, rho_struct_r, tau, & 130 tau_struct_r 131 TYPE(pw_p_type), POINTER :: rho_nlcc, rho_nlcc_g 132 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, xc_pw_pool 133 TYPE(pw_type), POINTER :: tmp_g, tmp_g2, tmp_pw 134 TYPE(virial_type), POINTER :: virial 135 136 CALL timeset(routineN, handle) 137 138 CPASSERT(.NOT. ASSOCIATED(vxc_rho)) 139 CPASSERT(.NOT. ASSOCIATED(vxc_tau)) 140 NULLIFY (dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, & 141 tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, & 142 rho_m_gspace, rho_nlcc, rho_nlcc_g, rho_struct_r, rho_struct_g, tau_struct_r) 143 144 exc = 0.0_dp 145 my_just_energy = .FALSE. 146 IF (PRESENT(just_energy)) my_just_energy = just_energy 147 my_adiabatic_rescale_factor = 1.0_dp 148 do_adiabatic_rescaling = .FALSE. 149 IF (PRESENT(adiabatic_rescale_factor)) THEN 150 my_adiabatic_rescale_factor = adiabatic_rescale_factor 151 do_adiabatic_rescaling = .TRUE. 152 END IF 153 154 CALL get_ks_env(ks_env, & 155 dft_control=dft_control, & 156 pw_env=pw_env, & 157 cell=cell, & 158 virial=virial, & 159 rho_nlcc=rho_nlcc, & 160 rho_nlcc_g=rho_nlcc_g) 161 162 CALL qs_rho_get(rho_struct, & 163 tau_r_valid=tau_r_valid, & 164 rho_g_valid=rho_g_valid, & 165 rho_r=rho_struct_r, & 166 rho_g=rho_struct_g, & 167 tau_r=tau_struct_r) 168 169 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer) 170 171 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", & 172 i_val=myfun) 173 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", & 174 i_val=vdw) 175 176 vdW_nl = (vdw == xc_vdw_fun_nonloc) 177 ! this combination has not been investigated 178 CPASSERT(.NOT. (do_adiabatic_rescaling .AND. vdW_nl)) 179 ! are the necessary inputs available 180 IF (.NOT. (PRESENT(dispersion_env) .AND. PRESENT(edisp))) THEN 181 vdW_nl = .FALSE. 182 END IF 183 IF (PRESENT(edisp)) edisp = 0.0_dp 184 185 IF (myfun /= xc_none .OR. vdW_nl) THEN 186 187 ! test if the real space density is available 188 CPASSERT(ASSOCIATED(rho_struct)) 189 IF (dft_control%nspins /= 1 .AND. dft_control%nspins /= 2) & 190 CPABORT("nspins must be 1 or 2") 191 ! there are some options related to SIC here. 192 ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD) 193 ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta)) 194 ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0) 195 196 ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term 197 my_scaling = 1.0_dp 198 SELECT CASE (dft_control%sic_method_id) 199 CASE (sic_none) 200 ! all fine 201 CASE (sic_mauri_spz, sic_ad) 202 ! no idea yet what to do here in that case 203 CPASSERT(.NOT. tau_r_valid) 204 CASE (sic_mauri_us) 205 my_scaling = 1.0_dp - dft_control%sic_scaling_b 206 ! no idea yet what to do here in that case 207 CPASSERT(.NOT. tau_r_valid) 208 CASE (sic_eo) 209 ! NOTHING TO BE DONE 210 CASE DEFAULT 211 ! this case has not yet been treated here 212 CPABORT("NYI") 213 END SELECT 214 215 IF (dft_control%sic_scaling_b .EQ. 0.0_dp) THEN 216 sic_scaling_b_zero = .TRUE. 217 ELSE 218 sic_scaling_b_zero = .FALSE. 219 ENDIF 220 221 IF (PRESENT(pw_env_external)) & 222 pw_env => pw_env_external 223 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool) 224 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid) 225 226 ALLOCATE (rho_r(dft_control%nspins)) 227 IF (.NOT. uf_grid) THEN 228 DO ispin = 1, dft_control%nspins 229 rho_r(ispin)%pw => rho_struct_r(ispin)%pw 230 END DO 231 232 IF (tau_r_valid) THEN 233 ALLOCATE (tau(dft_control%nspins)) 234 DO ispin = 1, dft_control%nspins 235 tau(ispin)%pw => tau_struct_r(ispin)%pw 236 END DO 237 END IF 238 239 ! for gradient corrected functional the density in g space might 240 ! be useful so if we have it, we pass it in 241 IF (rho_g_valid) THEN 242 ALLOCATE (rho_g(dft_control%nspins)) 243 DO ispin = 1, dft_control%nspins 244 rho_g(ispin)%pw => rho_struct_g(ispin)%pw 245 END DO 246 END IF 247 ELSE 248 CPASSERT(rho_g_valid) 249 ALLOCATE (rho_g(dft_control%nspins)) 250 DO ispin = 1, dft_control%nspins 251 CALL pw_pool_create_pw(xc_pw_pool, rho_g(ispin)%pw, & 252 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D) 253 CALL pw_transfer(rho_struct_g(ispin)%pw, rho_g(ispin)%pw) 254 END DO 255 DO ispin = 1, dft_control%nspins 256 CALL pw_pool_create_pw(xc_pw_pool, rho_r(ispin)%pw, & 257 in_space=REALSPACE, use_data=REALDATA3D) 258 CALL pw_transfer(rho_g(ispin)%pw, rho_r(ispin)%pw) 259 END DO 260 IF (tau_r_valid) THEN 261 ! tau with finer grids is not implemented (at least not correctly), which this asserts 262 CPABORT("tau with finer grids") 263! ALLOCATE(tau(dft_control%nspins),stat=stat) 264! CPPostcondition(stat==0,cp_failure_level,routineP,failure) 265! DO ispin=1,dft_control%nspins 266! CALL pw_pool_create_pw(xc_pw_pool,tau(ispin)%pw,& 267! in_space=REALSPACE, use_data=REALDATA3D) 268! 269! CALL pw_pool_create_pw(xc_pw_pool,tmp_g,& 270! in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D) 271! CALL pw_pool_create_pw(auxbas_pw_pool,tmp_g2,& 272! in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D) 273! CALL pw_transfer(tau(ispin)%pw,tmp_g) 274! CALL pw_transfer(tmp_g,tmp_g2) 275! CALL pw_transfer(tmp_g2,tmp_pw) 276! CALL pw_pool_give_back_pw(auxbas_pw_pool,tmp_g2) 277! CALL pw_pool_give_back_pw(xc_pw_pool,tmp_g) 278! END DO 279 END IF 280 END IF 281 282 ! add the nlcc densities 283 IF (ASSOCIATED(rho_nlcc)) THEN 284 factor = 1.0_dp 285 DO ispin = 1, dft_control%nspins 286 CALL pw_axpy(rho_nlcc%pw, rho_r(ispin)%pw, factor) 287 CALL pw_axpy(rho_nlcc_g%pw, rho_g(ispin)%pw, factor) 288 ENDDO 289 ENDIF 290 291 ! 292 ! here the rho_r, rho_g, tau is what it should be 293 ! we get back the right my_vxc_rho and my_vxc_tau as required 294 ! 295 IF (my_just_energy) THEN 296 exc = xc_exc_calc(rho_r=rho_r, tau=tau, & 297 rho_g=rho_g, xc_section=xc_section, & 298 pw_pool=xc_pw_pool) 299 300 ELSE 301 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, & 302 rho_g=rho_g, tau=tau, exc=exc, & 303 xc_section=xc_section, & 304 pw_pool=xc_pw_pool, & 305 compute_virial=compute_virial, & 306 virial_xc=virial%pv_xc) 307 END IF 308 309 ! remove the nlcc densities (keep stuff in original state) 310 IF (ASSOCIATED(rho_nlcc)) THEN 311 factor = -1.0_dp 312 DO ispin = 1, dft_control%nspins 313 CALL pw_axpy(rho_nlcc%pw, rho_r(ispin)%pw, factor) 314 CALL pw_axpy(rho_nlcc_g%pw, rho_g(ispin)%pw, factor) 315 ENDDO 316 ENDIF 317 318 ! calclulate non-local vdW functional 319 ! only if this XC_SECTION has it 320 ! if yes, we use the dispersion_env from ks_env 321 ! this is dangerous, as it assumes a special connection xc_section -> qs_env 322 IF (vdW_nl) THEN 323 CALL get_ks_env(ks_env=ks_env, para_env=para_env) 324 ! no SIC functionals allowed 325 CPASSERT(dft_control%sic_method_id == sic_none) 326 ! 327 CALL pw_env_get(pw_env, vdw_pw_pool=vdw_pw_pool) 328 IF (my_just_energy) THEN 329 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, & 330 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env) 331 ELSE 332 CALL calculate_dispersion_nonloc(my_vxc_rho, rho_r, rho_g, edisp, dispersion_env, & 333 my_just_energy, vdw_pw_pool, xc_pw_pool, para_env, virial=virial) 334 END IF 335 END IF 336 337 !! Apply rescaling to the potential if requested 338 IF (.NOT. my_just_energy) THEN 339 IF (do_adiabatic_rescaling) THEN 340 IF (ASSOCIATED(my_vxc_rho)) THEN 341 DO ispin = 1, SIZE(my_vxc_rho) 342 my_vxc_rho(ispin)%pw%cr3d = my_vxc_rho(ispin)%pw%cr3d*my_adiabatic_rescale_factor 343 END DO 344 END IF 345 END IF 346 END IF 347 348 IF (my_scaling .NE. 1.0_dp) THEN 349 exc = exc*my_scaling 350 IF (ASSOCIATED(my_vxc_rho)) THEN 351 DO ispin = 1, SIZE(my_vxc_rho) 352 my_vxc_rho(ispin)%pw%cr3d = my_vxc_rho(ispin)%pw%cr3d*my_scaling 353 ENDDO 354 ENDIF 355 IF (ASSOCIATED(my_vxc_tau)) THEN 356 DO ispin = 1, SIZE(my_vxc_tau) 357 my_vxc_tau(ispin)%pw%cr3d = my_vxc_tau(ispin)%pw%cr3d*my_scaling 358 ENDDO 359 ENDIF 360 ENDIF 361 362 ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer 363 ! pw -> coeff 364 IF (ASSOCIATED(my_vxc_rho)) THEN 365 ALLOCATE (vxc_rho(dft_control%nspins)) 366 DO ispin = 1, dft_control%nspins 367 vxc_rho(ispin)%pw => my_vxc_rho(ispin)%pw 368 END DO 369 DEALLOCATE (my_vxc_rho) 370 END IF 371 IF (ASSOCIATED(my_vxc_tau)) THEN 372 ALLOCATE (vxc_tau(dft_control%nspins)) 373 DO ispin = 1, dft_control%nspins 374 vxc_tau(ispin)%pw => my_vxc_tau(ispin)%pw 375 END DO 376 DEALLOCATE (my_vxc_tau) 377 END IF 378 379 ! compute again the xc but now for Exc(m,o) and the opposite sign 380 IF (dft_control%sic_method_id .EQ. sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN 381 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2)) 382 CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(1)%pw, & 383 use_data=COMPLEXDATA1D, & 384 in_space=RECIPROCALSPACE) 385 CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(1)%pw, & 386 use_data=REALDATA3D, & 387 in_space=REALSPACE) 388 CALL pw_copy(rho_struct_r(1)%pw, rho_m_rspace(1)%pw) 389 CALL pw_axpy(rho_struct_r(2)%pw, rho_m_rspace(1)%pw, alpha=-1._dp) 390 CALL pw_copy(rho_struct_g(1)%pw, rho_m_gspace(1)%pw) 391 CALL pw_axpy(rho_struct_g(2)%pw, rho_m_gspace(1)%pw, alpha=-1._dp) 392 ! bit sad, these will be just zero... 393 CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(2)%pw, & 394 use_data=COMPLEXDATA1D, & 395 in_space=RECIPROCALSPACE) 396 CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(2)%pw, & 397 use_data=REALDATA3D, & 398 in_space=REALSPACE) 399 CALL pw_zero(rho_m_rspace(2)%pw) 400 CALL pw_zero(rho_m_gspace(2)%pw) 401 402 rho_g(1)%pw => rho_m_gspace(1)%pw 403 rho_g(2)%pw => rho_m_gspace(2)%pw 404 rho_r(1)%pw => rho_m_rspace(1)%pw 405 rho_r(2)%pw => rho_m_rspace(2)%pw 406 407 IF (my_just_energy) THEN 408 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, & 409 rho_g=rho_g, xc_section=xc_section, & 410 pw_pool=xc_pw_pool) 411 ELSE 412 ! virial untested 413 CPASSERT(.NOT. compute_virial) 414 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, & 415 rho_g=rho_g, tau=tau, exc=exc_m, & 416 xc_section=xc_section, & 417 pw_pool=xc_pw_pool, & 418 compute_virial=.FALSE., & 419 virial_xc=virial_xc_tmp) 420 END IF 421 422 exc = exc - dft_control%sic_scaling_b*exc_m 423 424 ! and take care of the potential only vxc_rho is taken into account 425 IF (.NOT. my_just_energy) THEN 426 vxc_rho(1)%pw%cr3d = vxc_rho(1)%pw%cr3d-dft_control%sic_scaling_b* & 427 my_vxc_rho(1)%pw%cr3d 428 vxc_rho(2)%pw%cr3d = vxc_rho(2)%pw%cr3d+dft_control%sic_scaling_b* & 429 my_vxc_rho(1)%pw%cr3d ! 1=m 430 CALL pw_release(my_vxc_rho(1)%pw) 431 CALL pw_release(my_vxc_rho(2)%pw) 432 DEALLOCATE (my_vxc_rho) 433 ENDIF 434 435 DO ispin = 1, 2 436 CALL pw_pool_give_back_pw(xc_pw_pool, rho_m_rspace(ispin)%pw) 437 CALL pw_pool_give_back_pw(xc_pw_pool, rho_m_gspace(ispin)%pw) 438 ENDDO 439 DEALLOCATE (rho_m_rspace) 440 DEALLOCATE (rho_m_gspace) 441 442 ENDIF 443 444 ! now we have - sum_s N_s * Exc(rho_s/N_s,0) 445 IF (dft_control%sic_method_id .EQ. sic_ad .AND. .NOT. sic_scaling_b_zero) THEN 446 447 ! find out how many elecs we have 448 CALL get_ks_env(ks_env, nelectron_spin=nelec_spin) 449 450 ALLOCATE (rho_m_rspace(2), rho_m_gspace(2)) 451 DO ispin = 1, 2 452 CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(ispin)%pw, & 453 use_data=COMPLEXDATA1D, & 454 in_space=RECIPROCALSPACE) 455 CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(ispin)%pw, & 456 use_data=REALDATA3D, & 457 in_space=REALSPACE) 458 ENDDO 459 460 rho_g(1)%pw => rho_m_gspace(1)%pw 461 rho_g(2)%pw => rho_m_gspace(2)%pw 462 rho_r(1)%pw => rho_m_rspace(1)%pw 463 rho_r(2)%pw => rho_m_rspace(2)%pw 464 465 DO ispin = 1, 2 466 IF (nelec_spin(ispin) .GT. 0.0_dp) THEN 467 nelec_s_inv = 1.0_dp/nelec_spin(ispin) 468 ELSE 469 ! does it matter if there are no electrons with this spin (H) ? 470 nelec_s_inv = 0.0_dp 471 ENDIF 472 CALL pw_copy(rho_struct_r(ispin)%pw, rho_m_rspace(1)%pw) 473 CALL pw_copy(rho_struct_g(ispin)%pw, rho_m_gspace(1)%pw) 474 CALL pw_scale(rho_m_rspace(1)%pw, nelec_s_inv) 475 CALL pw_scale(rho_m_gspace(1)%pw, nelec_s_inv) 476 CALL pw_zero(rho_m_rspace(2)%pw) 477 CALL pw_zero(rho_m_gspace(2)%pw) 478 479 IF (my_just_energy) THEN 480 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, & 481 rho_g=rho_g, xc_section=xc_section, & 482 pw_pool=xc_pw_pool) 483 ELSE 484 ! virial untested 485 CPASSERT(.NOT. compute_virial) 486 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, & 487 rho_g=rho_g, tau=tau, exc=exc_m, & 488 xc_section=xc_section, & 489 pw_pool=xc_pw_pool, & 490 compute_virial=.FALSE., & 491 virial_xc=virial_xc_tmp) 492 END IF 493 494 exc = exc - dft_control%sic_scaling_b*nelec_spin(ispin)*exc_m 495 496 ! and take care of the potential only vxc_rho is taken into account 497 IF (.NOT. my_just_energy) THEN 498 vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d-dft_control%sic_scaling_b* & 499 my_vxc_rho(1)%pw%cr3d 500 CALL pw_release(my_vxc_rho(1)%pw) 501 CALL pw_release(my_vxc_rho(2)%pw) 502 DEALLOCATE (my_vxc_rho) 503 ENDIF 504 ENDDO 505 506 DO ispin = 1, 2 507 CALL pw_pool_give_back_pw(xc_pw_pool, rho_m_rspace(ispin)%pw) 508 CALL pw_pool_give_back_pw(xc_pw_pool, rho_m_gspace(ispin)%pw) 509 ENDDO 510 DEALLOCATE (rho_m_rspace) 511 DEALLOCATE (rho_m_gspace) 512 513 ENDIF 514 515 ! compute again the xc but now for Exc(n_down,n_down) 516 IF (dft_control%sic_method_id .EQ. sic_mauri_us .AND. .NOT. sic_scaling_b_zero) THEN 517 rho_r(1)%pw => rho_struct_r(2)%pw 518 rho_r(2)%pw => rho_struct_r(2)%pw 519 IF (rho_g_valid) THEN 520 rho_g(1)%pw => rho_struct_g(2)%pw 521 rho_g(2)%pw => rho_struct_g(2)%pw 522 ENDIF 523 524 IF (my_just_energy) THEN 525 exc_m = xc_exc_calc(rho_r=rho_r, tau=tau, & 526 rho_g=rho_g, xc_section=xc_section, & 527 pw_pool=xc_pw_pool) 528 ELSE 529 ! virial untested 530 CPASSERT(.NOT. compute_virial) 531 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, & 532 rho_g=rho_g, tau=tau, exc=exc_m, & 533 xc_section=xc_section, & 534 pw_pool=xc_pw_pool, & 535 compute_virial=.FALSE., & 536 virial_xc=virial_xc_tmp) 537 END IF 538 539 exc = exc + dft_control%sic_scaling_b*exc_m 540 541 ! and take care of the potential 542 IF (.NOT. my_just_energy) THEN 543 ! both go to minority spin 544 vxc_rho(2)%pw%cr3d = vxc_rho(2)%pw%cr3d+ & 545 2.0_dp*dft_control%sic_scaling_b*my_vxc_rho(1)%pw%cr3d 546 CALL pw_release(my_vxc_rho(1)%pw) 547 CALL pw_release(my_vxc_rho(2)%pw) 548 DEALLOCATE (my_vxc_rho) 549 ENDIF 550 551 ENDIF 552 553 ! 554 ! cleanups 555 ! 556 IF (uf_grid) THEN 557 DO ispin = 1, SIZE(rho_r) 558 CALL pw_pool_give_back_pw(xc_pw_pool, rho_r(ispin)%pw) 559 END DO 560 IF (ASSOCIATED(vxc_rho)) THEN 561 DO ispin = 1, SIZE(vxc_rho) 562 CALL pw_pool_create_pw(auxbas_pw_pool, tmp_pw, & 563 in_space=REALSPACE, use_data=REALDATA3D) 564 565 CALL pw_pool_create_pw(xc_pw_pool, tmp_g, & 566 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D) 567 CALL pw_pool_create_pw(auxbas_pw_pool, tmp_g2, & 568 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D) 569 CALL pw_transfer(vxc_rho(ispin)%pw, tmp_g) 570 CALL pw_transfer(tmp_g, tmp_g2) 571 CALL pw_transfer(tmp_g2, tmp_pw) 572 CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp_g2) 573 CALL pw_pool_give_back_pw(xc_pw_pool, tmp_g) 574 !FM CALL pw_zero(tmp_pw) 575 !FM CALL pw_restrict_s3(vxc_rho(ispin)%pw,tmp_pw,& 576 !FM auxbas_pw_pool,param_section=interp_section) 577 CALL pw_pool_give_back_pw(xc_pw_pool, vxc_rho(ispin)%pw) 578 vxc_rho(ispin)%pw => tmp_pw 579 NULLIFY (tmp_pw) 580 END DO 581 END IF 582 IF (ASSOCIATED(vxc_tau)) THEN 583 DO ispin = 1, SIZE(vxc_tau) 584 CALL pw_pool_create_pw(auxbas_pw_pool, tmp_pw, & 585 in_space=REALSPACE, use_data=REALDATA3D) 586 587 CALL pw_pool_create_pw(xc_pw_pool, tmp_g, & 588 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D) 589 CALL pw_pool_create_pw(auxbas_pw_pool, tmp_g2, & 590 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D) 591 CALL pw_transfer(vxc_tau(ispin)%pw, tmp_g) 592 CALL pw_transfer(tmp_g, tmp_g2) 593 CALL pw_transfer(tmp_g2, tmp_pw) 594 CALL pw_pool_give_back_pw(auxbas_pw_pool, tmp_g2) 595 CALL pw_pool_give_back_pw(xc_pw_pool, tmp_g) 596 !FM CALL pw_zero(tmp_pw) 597 !FM CALL pw_restrict_s3(vxc_rho(ispin)%pw,tmp_pw,& 598 !FM auxbas_pw_pool,param_section=interp_section) 599 CALL pw_pool_give_back_pw(xc_pw_pool, vxc_tau(ispin)%pw) 600 vxc_tau(ispin)%pw => tmp_pw 601 NULLIFY (tmp_pw) 602 END DO 603 END IF 604 605 END IF 606 DEALLOCATE (rho_r) 607 IF (ASSOCIATED(rho_g)) THEN 608 IF (uf_grid) THEN 609 DO ispin = 1, SIZE(rho_g) 610 CALL pw_pool_give_back_pw(xc_pw_pool, rho_g(ispin)%pw) 611 END DO 612 END IF 613 DEALLOCATE (rho_g) 614 END IF 615 IF (ASSOCIATED(tau)) THEN 616 IF (uf_grid) THEN 617 DO ispin = 1, SIZE(tau) 618 CALL pw_pool_give_back_pw(xc_pw_pool, tau(ispin)%pw) 619 END DO 620 END IF 621 DEALLOCATE (tau) 622 END IF 623 624 END IF 625 626 CALL timestop(handle) 627 628 END SUBROUTINE qs_vxc_create 629 630! ************************************************************************************************** 631!> \brief calculates the XC density: E_xc(r) - V_xc(r)*rho(r) 632!> \param ks_env to get all the needed things 633!> \param rho_struct density 634!> \param xc_section ... 635!> \param xc_den will contain the xc energy density 636!> \author JGH 637! ************************************************************************************************** 638 SUBROUTINE qs_xc_density(ks_env, rho_struct, xc_section, xc_den) 639 640 TYPE(qs_ks_env_type), POINTER :: ks_env 641 TYPE(qs_rho_type), POINTER :: rho_struct 642 TYPE(section_vals_type), POINTER :: xc_section 643 TYPE(pw_p_type), INTENT(INOUT) :: xc_den 644 645 CHARACTER(len=*), PARAMETER :: routineN = 'qs_xc_density', routineP = moduleN//':'//routineN 646 647 INTEGER :: handle, ispin, myfun, nspins, vdw 648 LOGICAL :: rho_g_valid, tau_r_valid, uf_grid, vdW_nl 649 REAL(KIND=dp) :: exc, factor, ovol 650 REAL(KIND=dp), DIMENSION(3, 3) :: vdum 651 TYPE(cell_type), POINTER :: cell 652 TYPE(dft_control_type), POINTER :: dft_control 653 TYPE(pw_env_type), POINTER :: pw_env 654 TYPE(pw_p_type) :: exc_r 655 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g, rho_r, tau_r, vxc_rho, vxc_tau 656 TYPE(pw_p_type), POINTER :: rho_nlcc, rho_nlcc_g 657 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool 658 659 CALL timeset(routineN, handle) 660 661 CALL pw_zero(xc_den%pw) 662 663 CALL get_ks_env(ks_env, & 664 dft_control=dft_control, & 665 pw_env=pw_env, & 666 cell=cell, & 667 rho_nlcc=rho_nlcc, & 668 rho_nlcc_g=rho_nlcc_g) 669 670 CALL qs_rho_get(rho_struct, & 671 tau_r_valid=tau_r_valid, & 672 rho_g_valid=rho_g_valid, & 673 rho_r=rho_r, & 674 rho_g=rho_g, & 675 tau_r=tau_r) 676 677 CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun) 678 CALL section_vals_val_get(xc_section, "VDW_POTENTIAL%POTENTIAL_TYPE", i_val=vdw) 679 vdW_nl = (vdw == xc_vdw_fun_nonloc) 680 IF (tau_r_valid) THEN 681 CALL cp_warn(__LOCATION__, "Tau contribution will not be correctly handled") 682 END IF 683 IF (vdW_nl) THEN 684 CALL cp_warn(__LOCATION__, "vdW functional contribution will be ignored") 685 END IF 686 687 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool) 688 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid) 689 IF (uf_grid) THEN 690 CALL cp_warn(__LOCATION__, "Fine grid option not possible with local energy") 691 CPABORT("Fine Grid in Local Energy") 692 END IF 693 694 IF (myfun /= xc_none) THEN 695 696 CPASSERT(ASSOCIATED(rho_struct)) 697 CPASSERT(dft_control%sic_method_id == sic_none) 698 699 nspins = dft_control%nspins 700 ! add the nlcc densities 701 IF (ASSOCIATED(rho_nlcc)) THEN 702 factor = 1.0_dp 703 DO ispin = 1, nspins 704 CALL pw_axpy(rho_nlcc%pw, rho_r(ispin)%pw, factor) 705 CALL pw_axpy(rho_nlcc_g%pw, rho_g(ispin)%pw, factor) 706 ENDDO 707 ENDIF 708 NULLIFY (vxc_rho, vxc_tau, exc_r%pw) 709 CALL xc_vxc_pw_create1(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, & 710 rho_g=rho_g, tau=tau_r, exc=exc, & 711 xc_section=xc_section, & 712 pw_pool=xc_pw_pool, & 713 compute_virial=.FALSE., & 714 virial_xc=vdum, & 715 exc_r=exc_r) 716 ! remove the nlcc densities (keep stuff in original state) 717 IF (ASSOCIATED(rho_nlcc)) THEN 718 factor = -1.0_dp 719 DO ispin = 1, dft_control%nspins 720 CALL pw_axpy(rho_nlcc%pw, rho_r(ispin)%pw, factor) 721 CALL pw_axpy(rho_nlcc_g%pw, rho_g(ispin)%pw, factor) 722 ENDDO 723 ENDIF 724 ! 725 ovol = 1.0_dp/xc_den%pw%pw_grid%dvol 726 CALL pw_copy(exc_r%pw, xc_den%pw) 727 DO ispin = 1, nspins 728!deb CALL pw_multiply(xc_den%pw, vxc_rho(ispin)%pw, rho_r(ispin)%pw, alpha=-ovol) 729 CALL pw_multiply(xc_den%pw, vxc_rho(ispin)%pw, rho_r(ispin)%pw, alpha=-1.0_dp) 730 END DO 731 ! remove arrays 732 DO ispin = 1, nspins 733 CALL pw_release(vxc_rho(ispin)%pw) 734 END DO 735 DEALLOCATE (vxc_rho) 736 CALL pw_release(exc_r%pw) 737 ! 738 END IF 739 740 CALL timestop(handle) 741 742 END SUBROUTINE qs_xc_density 743 744END MODULE qs_vxc 745