1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6MODULE hartree_local_methods 7 USE atomic_kind_types, ONLY: atomic_kind_type,& 8 get_atomic_kind 9 USE atprop_types, ONLY: atprop_array_init,& 10 atprop_type 11 USE basis_set_types, ONLY: get_gto_basis_set,& 12 gto_basis_set_type 13 USE cell_types, ONLY: cell_type 14 USE cp_control_types, ONLY: dft_control_type 15 USE cp_para_types, ONLY: cp_para_env_type 16 USE hartree_local_types, ONLY: allocate_ecoul_1center,& 17 ecoul_1center_type,& 18 hartree_local_type,& 19 set_ecoul_1c 20 USE input_constants, ONLY: tddfpt_singlet 21 USE kinds, ONLY: dp 22 USE mathconstants, ONLY: fourpi,& 23 pi 24 USE message_passing, ONLY: mp_sum 25 USE orbital_pointers, ONLY: indso,& 26 nsoset 27 USE pw_env_types, ONLY: pw_env_get,& 28 pw_env_type 29 USE pw_poisson_types, ONLY: pw_poisson_periodic,& 30 pw_poisson_type 31 USE qs_charges_types, ONLY: qs_charges_type 32 USE qs_environment_types, ONLY: get_qs_env,& 33 qs_environment_type 34 USE qs_grid_atom, ONLY: grid_atom_type 35 USE qs_harmonics_atom, ONLY: get_none0_cg_list,& 36 harmonics_atom_type 37 USE qs_kind_types, ONLY: get_qs_kind,& 38 get_qs_kind_set,& 39 qs_kind_type 40 USE qs_local_rho_types, ONLY: rhoz_type 41 USE qs_p_env_types, ONLY: qs_p_env_type 42 USE qs_rho0_types, ONLY: get_rho0_mpole,& 43 rho0_atom_type,& 44 rho0_mpole_type 45 USE qs_rho_atom_types, ONLY: get_rho_atom,& 46 rho_atom_coeff,& 47 rho_atom_type 48 USE util, ONLY: get_limit 49#include "./base/base_uses.f90" 50 51 IMPLICIT NONE 52 53 PRIVATE 54 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods' 56 57 ! Public Subroutine 58 59 PUBLIC :: init_coulomb_local, calculate_Vh_1center, Vh_1c_gg_integrals 60 61CONTAINS 62 63! ************************************************************************************************** 64!> \brief ... 65!> \param hartree_local ... 66!> \param natom ... 67! ************************************************************************************************** 68 SUBROUTINE init_coulomb_local(hartree_local, natom) 69 70 TYPE(hartree_local_type), POINTER :: hartree_local 71 INTEGER, INTENT(IN) :: natom 72 73 CHARACTER(len=*), PARAMETER :: routineN = 'init_coulomb_local', & 74 routineP = moduleN//':'//routineN 75 76 INTEGER :: handle 77 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c 78 79 CALL timeset(routineN, handle) 80 81 NULLIFY (ecoul_1c) 82 ! Allocate and Initialize 1-center Potentials and Integrals 83 CALL allocate_ecoul_1center(ecoul_1c, natom) 84 hartree_local%ecoul_1c => ecoul_1c 85 86 CALL timestop(handle) 87 88 END SUBROUTINE init_coulomb_local 89 90! ************************************************************************************************** 91!> \brief Calculates Hartree potential for hard and soft densities (including 92!> nuclear charge and compensation charges) using numerical integration 93!> \param vrad_h ... 94!> \param vrad_s ... 95!> \param rrad_h ... 96!> \param rrad_s ... 97!> \param rrad_0 ... 98!> \param rrad_z ... 99!> \param grid_atom ... 100!> \par History 101!> 05.2012 JGH refactoring 102!> \author ?? 103! ************************************************************************************************** 104 SUBROUTINE calculate_Vh_1center(vrad_h, vrad_s, rrad_h, rrad_s, rrad_0, rrad_z, grid_atom) 105 106 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s 107 TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: rrad_h, rrad_s 108 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0 109 REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z 110 TYPE(grid_atom_type), POINTER :: grid_atom 111 112 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_Vh_1center', & 113 routineP = moduleN//':'//routineN 114 115 INTEGER :: handle, ir, iso, ispin, l_ang, & 116 max_s_harm, nchannels, nr, nspins 117 REAL(dp) :: I1_down, I1_up, I2_down, I2_up, prefactor 118 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rho_1, rho_2 119 REAL(dp), DIMENSION(:), POINTER :: wr 120 REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l 121 122 CALL timeset(routineN, handle) 123 124 nr = grid_atom%nr 125 max_s_harm = SIZE(vrad_h, 2) 126 nspins = SIZE(rrad_h, 1) 127 nchannels = SIZE(rrad_0, 2) 128 129 r2l => grid_atom%rad2l 130 oor2l => grid_atom%oorad2l 131 wr => grid_atom%wr 132 133 ALLOCATE (rho_1(nr, max_s_harm), rho_2(nr, max_s_harm)) 134 rho_1 = 0.0_dp 135 rho_2 = 0.0_dp 136 137 ! Case lm = 0 138 rho_1(:, 1) = rrad_z(:) 139 rho_2(:, 1) = rrad_0(:, 1) 140 141 DO iso = 2, nchannels 142 rho_2(:, iso) = rrad_0(:, iso) 143 END DO 144 145 DO iso = 1, max_s_harm 146 DO ispin = 1, nspins 147 rho_1(:, iso) = rho_1(:, iso) + rrad_h(ispin)%r_coef(:, iso) 148 rho_2(:, iso) = rho_2(:, iso) + rrad_s(ispin)%r_coef(:, iso) 149 END DO 150 151 l_ang = indso(1, iso) 152 prefactor = fourpi/(2._dp*l_ang + 1._dp) 153 154 rho_1(:, iso) = rho_1(:, iso)*wr(:) 155 rho_2(:, iso) = rho_2(:, iso)*wr(:) 156 157 I1_up = 0.0_dp 158 I1_down = 0.0_dp 159 I2_up = 0.0_dp 160 I2_down = 0.0_dp 161 162 I1_up = r2l(nr, l_ang)*rho_1(nr, iso) 163 I2_up = r2l(nr, l_ang)*rho_2(nr, iso) 164 165 DO ir = nr - 1, 1, -1 166 I1_down = I1_down + oor2l(ir, l_ang + 1)*rho_1(ir, iso) 167 I2_down = I2_down + oor2l(ir, l_ang + 1)*rho_2(ir, iso) 168 END DO 169 170 vrad_h(nr, iso) = vrad_h(nr, iso) + prefactor* & 171 (oor2l(nr, l_ang + 1)*I1_up + r2l(nr, l_ang)*I1_down) 172 vrad_s(nr, iso) = vrad_s(nr, iso) + prefactor* & 173 (oor2l(nr, l_ang + 1)*I2_up + r2l(nr, l_ang)*I2_down) 174 175 DO ir = nr - 1, 1, -1 176 I1_up = I1_up + r2l(ir, l_ang)*rho_1(ir, iso) 177 I1_down = I1_down - oor2l(ir, l_ang + 1)*rho_1(ir, iso) 178 I2_up = I2_up + r2l(ir, l_ang)*rho_2(ir, iso) 179 I2_down = I2_down - oor2l(ir, l_ang + 1)*rho_2(ir, iso) 180 181 vrad_h(ir, iso) = vrad_h(ir, iso) + prefactor* & 182 (oor2l(ir, l_ang + 1)*I1_up + r2l(ir, l_ang)*I1_down) 183 vrad_s(ir, iso) = vrad_s(ir, iso) + prefactor* & 184 (oor2l(ir, l_ang + 1)*I2_up + r2l(ir, l_ang)*I2_down) 185 186 END DO 187 188 END DO 189 190 DEALLOCATE (rho_1, rho_2) 191 192 CALL timestop(handle) 193 194 END SUBROUTINE calculate_Vh_1center 195 196! ************************************************************************************************** 197!> \brief Calculates one center GAPW Hartree energies and matrix elements 198!> Hartree potentials are input 199!> Takes possible background charge into account 200!> Special case for densities without core charge 201!> \param qs_env ... 202!> \param energy_hartree_1c ... 203!> \param tddft ... 204!> \param do_triplet ... 205!> \param p_env ... 206!> \par History 207!> 05.2012 JGH refactoring 208!> \author ?? 209! ************************************************************************************************** 210 SUBROUTINE Vh_1c_gg_integrals(qs_env, energy_hartree_1c, tddft, do_triplet, p_env) 211 212 TYPE(qs_environment_type), POINTER :: qs_env 213 REAL(kind=dp), INTENT(out) :: energy_hartree_1c 214 LOGICAL, INTENT(IN), OPTIONAL :: tddft, do_triplet 215 TYPE(qs_p_env_type), OPTIONAL, POINTER :: p_env 216 217 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_gg_integrals', & 218 routineP = moduleN//':'//routineN 219 220 INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, lmax0, & 221 lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, natom, & 222 nchan_0, nkind, nr, nset, nsotot, nspins, num_pe 223 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list 224 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list 225 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf 226 LOGICAL :: atenergy, core_charge, my_periodic, & 227 my_tddft, paw_atom 228 REAL(dp) :: back_ch, factor 229 REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr 230 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_hh, aVh1b_ss, g0_h_w 231 REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z 232 REAL(dp), DIMENSION(:, :), POINTER :: g0_h, gsph, rrad_0, Vh1_h, Vh1_s, & 233 vrrad_0, zet 234 REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg 235 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 236 TYPE(atprop_type), POINTER :: atprop 237 TYPE(cell_type), POINTER :: cell 238 TYPE(cp_para_env_type), POINTER :: para_env 239 TYPE(dft_control_type), POINTER :: dft_control 240 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c 241 TYPE(grid_atom_type), POINTER :: grid_atom 242 TYPE(gto_basis_set_type), POINTER :: orb_basis 243 TYPE(harmonics_atom_type), POINTER :: harmonics 244 TYPE(pw_env_type), POINTER :: pw_env 245 TYPE(pw_poisson_type), POINTER :: poisson_env 246 TYPE(qs_charges_type), POINTER :: qs_charges 247 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 248 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set 249 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 250 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 251 TYPE(rho_atom_type), POINTER :: rho_atom 252 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set 253 254 CALL timeset(routineN, handle) 255 256 NULLIFY (cell, dft_control, para_env, poisson_env, pw_env, qs_charges) 257 NULLIFY (atomic_kind_set, qs_kind_set, rho_atom_set, rho0_atom_set) 258 NULLIFY (rho0_mpole, rhoz_set, ecoul_1c) 259 NULLIFY (atom_list, grid_atom, harmonics) 260 NULLIFY (orb_basis, lmin, lmax, npgf, zet) 261 NULLIFY (gsph, atprop) 262 263 my_tddft = .FALSE. 264 IF (PRESENT(tddft)) my_tddft = tddft 265 core_charge = .NOT. my_tddft 266 267 CALL get_qs_env(qs_env=qs_env, & 268 cell=cell, dft_control=dft_control, & 269 para_env=para_env, & 270 atomic_kind_set=atomic_kind_set, & 271 qs_kind_set=qs_kind_set, atprop=atprop, & 272 pw_env=pw_env, qs_charges=qs_charges) 273 274 CALL pw_env_get(pw_env, poisson_env=poisson_env) 275 my_periodic = (poisson_env%method == pw_poisson_periodic) 276 277 back_ch = qs_charges%background*cell%deth 278 279 IF (my_tddft) THEN 280 rho_atom_set => p_env%local_rho_set%rho_atom_set 281 rho0_atom_set => p_env%local_rho_set%rho0_atom_set 282 rho0_mpole => p_env%local_rho_set%rho0_mpole 283 ecoul_1c => p_env%hartree_local%ecoul_1c 284 ELSE 285 CALL get_qs_env(qs_env=qs_env, & 286 rho_atom_set=rho_atom_set, & 287 rho0_atom_set=rho0_atom_set, & 288 rho0_mpole=rho0_mpole, & 289 rhoz_set=rhoz_set, & 290 ecoul_1c=ecoul_1c) 291 END IF 292 293 nkind = SIZE(atomic_kind_set, 1) 294 nspins = dft_control%nspins 295 296 IF (my_tddft) THEN 297 IF (PRESENT(do_triplet)) THEN 298 IF (nspins == 1 .AND. do_triplet) RETURN 299 ELSE 300 IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN 301 ENDIF 302 END IF 303 304 CALL get_qs_kind_set(qs_kind_set, maxg_iso_not0=max_iso) 305 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax_0) 306 307 atenergy = .FALSE. 308 IF (ASSOCIATED(atprop)) THEN 309 atenergy = atprop%energy 310 IF (atenergy) THEN 311 CALL get_qs_env(qs_env=qs_env, natom=natom) 312 CALL atprop_array_init(atprop%ate1c, natom) 313 END IF 314 END IF 315 316 ! Put to 0 the local hartree energy contribution from 1 center integrals 317 energy_hartree_1c = 0.0_dp 318 319 ! Here starts the loop over all the atoms 320 DO ikind = 1, nkind 321 322 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat) 323 CALL get_qs_kind(qs_kind_set(ikind), & 324 basis_set=orb_basis, & 325 grid_atom=grid_atom, & 326 harmonics=harmonics, ngrid_rad=nr, & 327 max_iso_not0=max_iso_not0, paw_atom=paw_atom) 328 329 IF (paw_atom) THEN 330!=========== PAW =============== 331 CALL get_gto_basis_set(gto_basis_set=orb_basis, lmax=lmax, lmin=lmin, & 332 maxso=maxso, npgf=npgf, maxl=maxl, & 333 nset=nset, zet=zet) 334 335 max_s_harm = harmonics%max_s_harm 336 llmax = harmonics%llmax 337 338 nsotot = maxso*nset 339 ALLOCATE (gsph(nr, nsotot)) 340 ALLOCATE (gexp(nr)) 341 ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0)) 342 343 NULLIFY (Vh1_h, Vh1_s) 344 ALLOCATE (Vh1_h(nr, max_iso_not0)) 345 ALLOCATE (Vh1_s(nr, max_iso_not0)) 346 347 ALLOCATE (aVh1b_hh(nsotot, nsotot)) 348 ALLOCATE (aVh1b_ss(nsotot, nsotot)) 349 ALLOCATE (aVh1b_00(nsotot, nsotot)) 350 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm)) 351 352 NULLIFY (Qlm_gg, g0_h) 353 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, & 354 l0_ikind=lmax0, & 355 Qlm_gg=Qlm_gg, g0_h=g0_h) 356 357 nchan_0 = nsoset(lmax0) 358 359 IF (nchan_0 > max_iso_not0) CPABORT("channels for rho0 > # max of spherical harmonics") 360 361 NULLIFY (rrad_z, my_CG) 362 my_CG => harmonics%my_CG 363 364 ! set to zero temporary arrays 365 sqrtwr = 0.0_dp 366 g0_h_w = 0.0_dp 367 gexp = 0.0_dp 368 gsph = 0.0_dp 369 370 sqrtwr(1:nr) = SQRT(grid_atom%wr(1:nr)) 371 DO l_ang = 0, lmax0 372 g0_h_w(1:nr, l_ang) = g0_h(1:nr, l_ang)*grid_atom%wr(1:nr) 373 END DO 374 375 m1 = 0 376 DO iset1 = 1, nset 377 n1 = nsoset(lmax(iset1)) 378 DO ipgf1 = 1, npgf(iset1) 379 gexp(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr) 380 DO is1 = nsoset(lmin(iset1) - 1) + 1, nsoset(lmax(iset1)) 381 iso = is1 + (ipgf1 - 1)*n1 + m1 382 l_ang = indso(1, is1) 383 gsph(1:nr, iso) = grid_atom%rad2l(1:nr, l_ang)*gexp(1:nr) 384 END DO ! is1 385 END DO ! ipgf1 386 m1 = m1 + maxso 387 END DO ! iset1 388 389 ! Distribute the atoms of this kind 390 num_pe = para_env%num_pe 391 mepos = para_env%mepos 392 bo = get_limit(nat, num_pe, mepos) 393 394 DO iat = bo(1), bo(2) !1,nat 395 iatom = atom_list(iat) 396 rho_atom => rho_atom_set(iatom) 397 398 NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0) 399 IF (core_charge) THEN 400 rrad_z => rhoz_set(ikind)%r_coef 401 vrrad_z => rhoz_set(ikind)%vr_coef 402 END IF 403 rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef 404 vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef 405 IF (my_periodic .AND. back_ch .GT. 1.E-3_dp) THEN 406 factor = -2.0_dp*pi/3.0_dp*SQRT(fourpi)*qs_charges%background 407 ELSE 408 factor = 0._dp 409 END IF 410 411 CALL Vh_1c_atom_potential(rho_atom, vrrad_0, & 412 grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, & 413 nchan_0, nspins, max_iso_not0, factor) 414 415 CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, & 416 grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, & 417 nchan_0, nspins, max_iso_not0, atenergy, atprop%ate1c) 418 419 CALL Vh_1c_atom_integrals(rho_atom, & 420 aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, & 421 max_s_harm, llmax, cg_list, cg_n_list, & 422 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, g0_h_w, my_CG, Qlm_gg) 423 424 END DO ! iat 425 426 DEALLOCATE (aVh1b_hh) 427 DEALLOCATE (aVh1b_ss) 428 DEALLOCATE (aVh1b_00) 429 DEALLOCATE (Vh1_h, Vh1_s) 430 DEALLOCATE (cg_list, cg_n_list) 431 DEALLOCATE (gsph) 432 DEALLOCATE (gexp) 433 DEALLOCATE (sqrtwr, g0_h_w) 434 435 ELSE 436!=========== NO PAW =============== 437! This term is taken care of using the core density as in GPW 438 CYCLE 439 END IF ! paw 440 END DO ! ikind 441 442 CALL mp_sum(energy_hartree_1c, para_env%group) 443 444 CALL timestop(handle) 445 446 END SUBROUTINE Vh_1c_gg_integrals 447 448! ************************************************************************************************** 449 450! ************************************************************************************************** 451!> \brief ... 452!> \param rho_atom ... 453!> \param vrrad_0 ... 454!> \param grid_atom ... 455!> \param core_charge ... 456!> \param vrrad_z ... 457!> \param Vh1_h ... 458!> \param Vh1_s ... 459!> \param nchan_0 ... 460!> \param nspins ... 461!> \param max_iso_not0 ... 462!> \param bfactor ... 463! ************************************************************************************************** 464 SUBROUTINE Vh_1c_atom_potential(rho_atom, vrrad_0, & 465 grid_atom, core_charge, vrrad_z, Vh1_h, Vh1_s, & 466 nchan_0, nspins, max_iso_not0, bfactor) 467 468 TYPE(rho_atom_type), POINTER :: rho_atom 469 REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0 470 TYPE(grid_atom_type), POINTER :: grid_atom 471 LOGICAL, INTENT(IN) :: core_charge 472 REAL(dp), DIMENSION(:), POINTER :: vrrad_z 473 REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s 474 INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0 475 REAL(dp), INTENT(IN) :: bfactor 476 477 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_atom_potential', & 478 routineP = moduleN//':'//routineN 479 480 INTEGER :: ir, iso, ispin, nr 481 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: vr_h, vr_s 482 483 nr = grid_atom%nr 484 485 NULLIFY (vr_h, vr_s) 486 CALL get_rho_atom(rho_atom=rho_atom, vrho_rad_h=vr_h, vrho_rad_s=vr_s) 487 488 Vh1_h = 0.0_dp 489 Vh1_s = 0.0_dp 490 491 IF (core_charge) Vh1_h(:, 1) = vrrad_z(:) 492 493 DO iso = 1, nchan_0 494 Vh1_s(:, iso) = vrrad_0(:, iso) 495 END DO 496 497 DO ispin = 1, nspins 498 DO iso = 1, max_iso_not0 499 Vh1_h(:, iso) = Vh1_h(:, iso) + vr_h(ispin)%r_coef(:, iso) 500 Vh1_s(:, iso) = Vh1_s(:, iso) + vr_s(ispin)%r_coef(:, iso) 501 END DO 502 END DO 503 504 IF (bfactor /= 0._dp) THEN 505 DO ir = 1, nr 506 Vh1_h(ir, 1) = Vh1_h(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir) 507 Vh1_s(ir, 1) = Vh1_s(ir, 1) + bfactor*grid_atom%rad2(ir)*grid_atom%wr(ir) 508 END DO 509 END IF 510 511 END SUBROUTINE Vh_1c_atom_potential 512 513! ************************************************************************************************** 514 515! ************************************************************************************************** 516!> \brief ... 517!> \param energy_hartree_1c ... 518!> \param ecoul_1c ... 519!> \param rho_atom ... 520!> \param rrad_0 ... 521!> \param grid_atom ... 522!> \param iatom ... 523!> \param core_charge ... 524!> \param rrad_z ... 525!> \param Vh1_h ... 526!> \param Vh1_s ... 527!> \param nchan_0 ... 528!> \param nspins ... 529!> \param max_iso_not0 ... 530!> \param atenergy ... 531!> \param ate1c ... 532! ************************************************************************************************** 533 SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, & 534 grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, & 535 nchan_0, nspins, max_iso_not0, atenergy, ate1c) 536 537 REAL(dp), INTENT(INOUT) :: energy_hartree_1c 538 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c 539 TYPE(rho_atom_type), POINTER :: rho_atom 540 REAL(dp), DIMENSION(:, :), POINTER :: rrad_0 541 TYPE(grid_atom_type), POINTER :: grid_atom 542 INTEGER, INTENT(IN) :: iatom 543 LOGICAL, INTENT(IN) :: core_charge 544 REAL(dp), DIMENSION(:), POINTER :: rrad_z 545 REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s 546 INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0 547 LOGICAL, INTENT(IN) :: atenergy 548 REAL(dp), DIMENSION(:), POINTER :: ate1c 549 550 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_atom_energy', & 551 routineP = moduleN//':'//routineN 552 553 INTEGER :: iso, ispin, nr 554 REAL(dp) :: ecoul_1_0, ecoul_1_h, ecoul_1_s, & 555 ecoul_1_z 556 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s 557 558 nr = grid_atom%nr 559 560 NULLIFY (r_h, r_s) 561 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s) 562 563 ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz 564 ecoul_1_z = 0.0_dp 565 IF (core_charge) THEN 566 ecoul_1_z = 0.5_dp*SUM(Vh1_h(:, 1)*rrad_z(:)*grid_atom%wr(:)) 567 END IF 568 569 ! Calculate the contributions to Ecoul coming from Vh1_s*rho0 570 ecoul_1_0 = 0.0_dp 571 DO iso = 1, nchan_0 572 ecoul_1_0 = ecoul_1_0 + 0.5_dp*SUM(Vh1_s(:, iso)*rrad_0(:, iso)*grid_atom%wr(:)) 573 END DO 574 575 ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s 576 ecoul_1_s = 0.0_dp 577 ecoul_1_h = 0.0_dp 578 DO ispin = 1, nspins 579 DO iso = 1, max_iso_not0 580 ecoul_1_s = ecoul_1_s + 0.5_dp*SUM(Vh1_s(:, iso)*r_s(ispin)%r_coef(:, iso)*grid_atom%wr(:)) 581 ecoul_1_h = ecoul_1_h + 0.5_dp*SUM(Vh1_h(:, iso)*r_h(ispin)%r_coef(:, iso)*grid_atom%wr(:)) 582 END DO 583 END DO 584 585 CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1_z, ecoul_1_0=ecoul_1_0) 586 CALL set_ecoul_1c(ecoul_1c=ecoul_1c, iatom=iatom, ecoul_1_h=ecoul_1_h, ecoul_1_s=ecoul_1_s) 587 588 energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0 589 energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s 590 591 ! atomic energy contribution 592 IF (atenergy) THEN 593 ate1c(iatom) = ate1c(iatom) + ecoul_1_z - ecoul_1_0 594 END IF 595 596 END SUBROUTINE Vh_1c_atom_energy 597 598!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 599 600! ************************************************************************************************** 601!> \brief ... 602!> \param rho_atom ... 603!> \param aVh1b_hh ... 604!> \param aVh1b_ss ... 605!> \param aVh1b_00 ... 606!> \param Vh1_h ... 607!> \param Vh1_s ... 608!> \param max_iso_not0 ... 609!> \param max_s_harm ... 610!> \param llmax ... 611!> \param cg_list ... 612!> \param cg_n_list ... 613!> \param nset ... 614!> \param npgf ... 615!> \param lmin ... 616!> \param lmax ... 617!> \param nsotot ... 618!> \param maxso ... 619!> \param nspins ... 620!> \param nchan_0 ... 621!> \param gsph ... 622!> \param g0_h_w ... 623!> \param my_CG ... 624!> \param Qlm_gg ... 625! ************************************************************************************************** 626 SUBROUTINE Vh_1c_atom_integrals(rho_atom, & 627 aVh1b_hh, aVh1b_ss, aVh1b_00, Vh1_h, Vh1_s, max_iso_not0, & 628 max_s_harm, llmax, cg_list, cg_n_list, & 629 nset, npgf, lmin, lmax, nsotot, maxso, nspins, nchan_0, gsph, g0_h_w, my_CG, Qlm_gg) 630 631 TYPE(rho_atom_type), POINTER :: rho_atom 632 REAL(dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00 633 REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s 634 INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, llmax 635 INTEGER, DIMENSION(:, :, :) :: cg_list 636 INTEGER, DIMENSION(:) :: cg_n_list 637 INTEGER, INTENT(IN) :: nset 638 INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax 639 INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0 640 REAL(dp), DIMENSION(:, :), POINTER :: gsph 641 REAL(dp), DIMENSION(:, 0:) :: g0_h_w 642 REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg 643 644 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_atom_integrals', & 645 routineP = moduleN//':'//routineN 646 647 INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, & 648 iset2, iso, iso1, iso2, ispin, l_ang, & 649 m1, m2, max_iso_not0_local, n1, n2, nr 650 REAL(dp) :: gVg_0, gVg_h, gVg_s 651 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s 652 653 NULLIFY (int_local_h, int_local_s) 654 CALL get_rho_atom(rho_atom=rho_atom, & 655 ga_Vlocal_gb_h=int_local_h, & 656 ga_Vlocal_gb_s=int_local_s) 657 658 ! Calculate the integrals of the potential with 2 primitives 659 aVh1b_hh = 0.0_dp 660 aVh1b_ss = 0.0_dp 661 aVh1b_00 = 0.0_dp 662 663 nr = SIZE(gsph, 1) 664 665 m1 = 0 666 DO iset1 = 1, nset 667 m2 = 0 668 DO iset2 = 1, nset 669 CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), & 670 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local) 671 672 n1 = nsoset(lmax(iset1)) 673 DO ipgf1 = 1, npgf(iset1) 674 n2 = nsoset(lmax(iset2)) 675 DO ipgf2 = 1, npgf(iset2) 676 ! with contributions to V1_s*rho0 677 DO iso = 1, nchan_0 678 l_ang = indso(1, iso) 679 gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang)) 680 DO icg = 1, cg_n_list(iso) 681 is1 = cg_list(1, icg, iso) 682 is2 = cg_list(2, icg, iso) 683 684 iso1 = is1 + n1*(ipgf1 - 1) + m1 685 iso2 = is2 + n2*(ipgf2 - 1) + m2 686 gVg_h = 0.0_dp 687 gVg_s = 0.0_dp 688 689 DO ir = 1, nr 690 gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso) 691 gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso) 692 END DO ! ir 693 694 aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso) 695 aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso) 696 aVh1b_00(iso1, iso2) = aVh1b_00(iso1, iso2) + gVg_0*Qlm_gg(iso1, iso2, iso) 697 698 END DO !icg 699 END DO ! iso 700 ! without contributions to V1_s*rho0 701 DO iso = nchan_0 + 1, max_iso_not0 702 DO icg = 1, cg_n_list(iso) 703 is1 = cg_list(1, icg, iso) 704 is2 = cg_list(2, icg, iso) 705 706 iso1 = is1 + n1*(ipgf1 - 1) + m1 707 iso2 = is2 + n2*(ipgf2 - 1) + m2 708 gVg_h = 0.0_dp 709 gVg_s = 0.0_dp 710 711 DO ir = 1, nr 712 gVg_h = gVg_h + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_h(ir, iso) 713 gVg_s = gVg_s + gsph(ir, iso1)*gsph(ir, iso2)*Vh1_s(ir, iso) 714 END DO ! ir 715 716 aVh1b_hh(iso1, iso2) = aVh1b_hh(iso1, iso2) + gVg_h*my_CG(is1, is2, iso) 717 aVh1b_ss(iso1, iso2) = aVh1b_ss(iso1, iso2) + gVg_s*my_CG(is1, is2, iso) 718 719 END DO !icg 720 END DO ! iso 721 END DO ! ipgf2 722 END DO ! ipgf1 723 m2 = m2 + maxso 724 END DO ! iset2 725 m1 = m1 + maxso 726 END DO !iset1 727 728 DO ispin = 1, nspins 729 CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_hh, 1, int_local_h(ispin)%r_coef, 1) 730 CALL daxpy(nsotot*nsotot, 1.0_dp, aVh1b_ss, 1, int_local_s(ispin)%r_coef, 1) 731 CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_h(ispin)%r_coef, 1) 732 CALL daxpy(nsotot*nsotot, -1.0_dp, aVh1b_00, 1, int_local_s(ispin)%r_coef, 1) 733 END DO ! ispin 734 735 END SUBROUTINE Vh_1c_atom_integrals 736 737!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 738 739END MODULE hartree_local_methods 740 741