1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6MODULE qs_rho0_ggrid 7 8 USE atomic_kind_types, ONLY: atomic_kind_type,& 9 get_atomic_kind 10 USE basis_set_types, ONLY: get_gto_basis_set,& 11 gto_basis_set_type 12 USE cell_types, ONLY: cell_type,& 13 pbc 14 USE cp_control_types, ONLY: dft_control_type 15 USE cp_para_types, ONLY: cp_para_env_type 16 USE cube_utils, ONLY: cube_info_type 17 USE input_constants, ONLY: tddfpt_singlet 18 USE kinds, ONLY: dp,& 19 int_8 20 USE mathconstants, ONLY: dfac,& 21 fourpi 22 USE memory_utilities, ONLY: reallocate 23 USE message_passing, ONLY: mp_sum 24 USE orbital_pointers, ONLY: indco,& 25 nco,& 26 ncoset,& 27 nso,& 28 nsoset 29 USE orbital_transformation_matrices, ONLY: orbtramat 30 USE particle_types, ONLY: particle_type 31 USE pw_env_types, ONLY: pw_env_get,& 32 pw_env_type 33 USE pw_methods, ONLY: pw_axpy,& 34 pw_copy,& 35 pw_integrate_function,& 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_p_type,& 41 pw_pool_type 42 USE pw_types, ONLY: COMPLEXDATA1D,& 43 REALDATA3D,& 44 REALSPACE,& 45 RECIPROCALSPACE,& 46 pw_p_type,& 47 pw_release 48 USE qs_collocate_density, ONLY: collocate_pgf_product_rspace 49 USE qs_environment_types, ONLY: get_qs_env,& 50 qs_environment_type 51 USE qs_force_types, ONLY: qs_force_type 52 USE qs_harmonics_atom, ONLY: get_none0_cg_list,& 53 harmonics_atom_type 54 USE qs_integrate_potential, ONLY: integrate_pgf_product_rspace 55 USE qs_kind_types, ONLY: get_qs_kind,& 56 qs_kind_type 57 USE qs_p_env_types, ONLY: qs_p_env_type 58 USE qs_rho0_types, ONLY: get_rho0_mpole,& 59 rho0_mpole_type 60 USE qs_rho_atom_types, ONLY: get_rho_atom,& 61 rho_atom_coeff,& 62 rho_atom_type 63 USE realspace_grid_types, ONLY: & 64 pw2rs, realspace_grid_desc_p_type, realspace_grid_desc_type, realspace_grid_p_type, & 65 realspace_grid_type, rs2pw, rs_grid_create, rs_grid_release, rs_grid_retain, rs_grid_zero, & 66 rs_pw_transfer 67 USE util, ONLY: get_limit 68 USE virial_types, ONLY: virial_type 69#include "./base/base_uses.f90" 70 71 IMPLICIT NONE 72 73 PRIVATE 74 75! *** Global parameters (only in this module) 76 77 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_ggrid' 78 79! *** Public subroutines *** 80 81 PUBLIC :: put_rho0_on_grid, rho0_s_grid_create, integrate_vhg0_rspace 82 83CONTAINS 84 85! ************************************************************************************************** 86!> \brief ... 87!> \param qs_env ... 88!> \param rho0 ... 89!> \param tot_rs_int ... 90! ************************************************************************************************** 91 SUBROUTINE put_rho0_on_grid(qs_env, rho0, tot_rs_int) 92 93 TYPE(qs_environment_type), POINTER :: qs_env 94 TYPE(rho0_mpole_type), POINTER :: rho0 95 REAL(dp), INTENT(OUT) :: tot_rs_int 96 97 CHARACTER(LEN=*), PARAMETER :: routineN = 'put_rho0_on_grid', & 98 routineP = moduleN//':'//routineN 99 100 INTEGER :: auxbas_grid, handle, iat, iatom, igrid, & 101 ikind, ithread, j, l0_ikind, lmax0, & 102 lmax_global, nat, nch_ik, nch_max, npme 103 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 104 LOGICAL :: paw_atom 105 REAL(KIND=dp) :: eps_rho_rspace, rpgf0, zet0 106 REAL(KIND=dp), DIMENSION(3) :: ra 107 REAL(KIND=dp), DIMENSION(:), POINTER :: Qlm_c 108 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 109 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 110 TYPE(cell_type), POINTER :: cell 111 TYPE(cp_para_env_type), POINTER :: para_env 112 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info 113 TYPE(dft_control_type), POINTER :: dft_control 114 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 115 TYPE(pw_env_type), POINTER :: pw_env 116 TYPE(pw_p_type) :: rho0_r_tmp 117 TYPE(pw_p_type), POINTER :: coeff_gspace, coeff_rspace, rho0_s_gs, & 118 rho0_s_rs 119 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 120 TYPE(pw_pool_type), POINTER :: pw_pool 121 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 122 TYPE(realspace_grid_desc_p_type), DIMENSION(:), & 123 POINTER :: descs 124 TYPE(realspace_grid_desc_type), POINTER :: desc 125 TYPE(realspace_grid_p_type), DIMENSION(:), POINTER :: grids 126 TYPE(realspace_grid_type), POINTER :: rs_grid 127 128 CALL timeset(routineN, handle) 129 130 NULLIFY (atomic_kind_set, qs_kind_set, cores, pab, Qlm_c) 131 132 NULLIFY (dft_control, pw_env, particle_set, para_env, cell, rho0_s_gs, rho0_s_rs) 133 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 134 particle_set=particle_set, & 135 atomic_kind_set=atomic_kind_set, & 136 qs_kind_set=qs_kind_set, & 137 para_env=para_env, & 138 pw_env=pw_env, cell=cell) 139 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 140 141 NULLIFY (descs, pw_pools) 142 CALL pw_env_get(pw_env=pw_env, rs_descs=descs, rs_grids=grids, pw_pools=pw_pools) 143 cube_info => pw_env%cube_info 144 auxbas_grid = pw_env%auxbas_grid 145 146 NULLIFY (rho0_s_gs, rho0_s_rs) 147 CALL get_rho0_mpole(rho0_mpole=rho0, lmax_0=lmax0, & 148 zet0_h=zet0, igrid_zet0_s=igrid, & 149 rho0_s_gs=rho0_s_gs, & 150 rho0_s_rs=rho0_s_rs) 151 152 ! *** set up the rs grid at level igrid 153 NULLIFY (rs_grid, desc, pw_pool, coeff_rspace, coeff_gspace) 154 desc => descs(igrid)%rs_desc 155 rs_grid => grids(igrid)%rs_grid 156 pw_pool => pw_pools(igrid)%pool 157 158 CPASSERT(ASSOCIATED(desc)) 159 CPASSERT(ASSOCIATED(pw_pool)) 160 161 IF (igrid /= auxbas_grid) THEN 162 ALLOCATE (coeff_rspace) 163 ALLOCATE (coeff_gspace) 164 165 CALL pw_pool_create_pw(pw_pool, coeff_rspace%pw, use_data=REALDATA3D, & 166 in_space=REALSPACE) 167 CALL pw_pool_create_pw(pw_pool, coeff_gspace%pw, & 168 use_data=COMPLEXDATA1D, & 169 in_space=RECIPROCALSPACE) 170 END IF 171 CALL rs_grid_retain(rs_grid) 172 CALL rs_grid_zero(rs_grid) 173 174 nch_max = ncoset(lmax0) 175 176 ALLOCATE (pab(nch_max, 1)) 177 178 !Compute global lmax for use in collocate_pgf_product_rspace 179 lmax_global = 0 180 DO ikind = 1, SIZE(atomic_kind_set) 181 CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind) 182 lmax_global = MAX(lmax_global, l0_ikind) 183 END DO 184 185 DO ikind = 1, SIZE(atomic_kind_set) 186 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat) 187 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom) 188 189 IF (.NOT. paw_atom .AND. dft_control%qs_control%gapw_control%nopaw_as_gpw) CYCLE 190 191 CALL get_rho0_mpole(rho0_mpole=rho0, ikind=ikind, l0_ikind=l0_ikind, & 192 rpgf0_s=rpgf0) 193 194 nch_ik = ncoset(l0_ikind) 195 pab = 0.0_dp 196 197 CALL reallocate(cores, 1, nat) 198 npme = 0 199 cores = 0 200 201 DO iat = 1, nat 202 iatom = atom_list(iat) 203 ra(:) = pbc(particle_set(iatom)%r, cell) 204 IF (rs_grid%desc%parallel .AND. .NOT. rs_grid%desc%distributed) THEN 205 ! replicated realspace grid, split the atoms up between procs 206 IF (MODULO(nat, rs_grid%desc%group_size) == rs_grid%desc%my_pos) THEN 207 npme = npme + 1 208 cores(npme) = iat 209 ENDIF 210 ELSE 211 npme = npme + 1 212 cores(npme) = iat 213 ENDIF 214 215 END DO 216 217 ithread = 0 218 DO j = 1, npme 219 220 iat = cores(j) 221 iatom = atom_list(iat) 222 223 CALL get_rho0_mpole(rho0_mpole=rho0, iat=iatom, Qlm_car=Qlm_c) 224 225 pab(1:nch_ik, 1) = Qlm_c(1:nch_ik) 226 227 ra(:) = pbc(particle_set(iatom)%r, cell) 228 229 CALL collocate_pgf_product_rspace( & 230 l0_ikind, zet0, 0, 0, 0.0_dp, 0, & 231 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 0.0_dp, 1.0_dp, pab, 0, 0, & 232 rs_grid, cell, cube_info(igrid), eps_rho_rspace, ga_gb_function=401, & 233 ithread=ithread, collocate_rho0=.TRUE., rpgf0_s=rpgf0, & 234 use_subpatch=.TRUE., subpatch_pattern=0_int_8, lmax_global=lmax_global) 235 236 END DO ! j 237 238 END DO ! ikind 239 240 IF (ASSOCIATED(cores)) THEN 241 DEALLOCATE (cores) 242 END IF 243 244 DEALLOCATE (pab) 245 246 IF (igrid /= auxbas_grid) THEN 247 CALL rs_pw_transfer(rs_grid, coeff_rspace%pw, rs2pw) 248 CALL rs_grid_release(rs_grid) 249 CALL pw_zero(rho0_s_gs%pw) 250 CALL pw_transfer(coeff_rspace%pw, coeff_gspace%pw) 251 CALL pw_axpy(coeff_gspace%pw, rho0_s_gs%pw) 252 253 tot_rs_int = pw_integrate_function(coeff_rspace%pw, isign=-1) 254 255 CALL pw_pool_give_back_pw(pw_pool, coeff_rspace%pw) 256 CALL pw_pool_give_back_pw(pw_pool, coeff_gspace%pw) 257 258 DEALLOCATE (coeff_rspace, coeff_gspace) 259 260 ELSE 261 262 CALL pw_pool_create_pw(pw_pool, rho0_r_tmp%pw, & 263 use_data=REALDATA3D, in_space=REALSPACE) 264 265 CALL rs_pw_transfer(rs_grid, rho0_r_tmp%pw, rs2pw) 266 CALL rs_grid_release(rs_grid) 267 268 tot_rs_int = pw_integrate_function(rho0_r_tmp%pw, isign=-1) 269 270 CALL pw_transfer(rho0_r_tmp%pw, rho0_s_rs%pw) 271 CALL pw_pool_give_back_pw(pw_pool, rho0_r_tmp%pw) 272 273 CALL pw_zero(rho0_s_gs%pw) 274 CALL pw_transfer(rho0_s_rs%pw, rho0_s_gs%pw) 275 END IF 276 CALL timestop(handle) 277 278 END SUBROUTINE put_rho0_on_grid 279 280! ************************************************************************************************** 281!> \brief ... 282!> \param qs_env ... 283!> \param rho0_mpole ... 284!> \param tddft ... 285! ************************************************************************************************** 286 SUBROUTINE rho0_s_grid_create(qs_env, rho0_mpole, tddft) 287 288 TYPE(qs_environment_type), POINTER :: qs_env 289 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 290 LOGICAL, INTENT(IN), OPTIONAL :: tddft 291 292 CHARACTER(len=*), PARAMETER :: routineN = 'rho0_s_grid_create', & 293 routineP = moduleN//':'//routineN 294 295 LOGICAL :: my_tddft 296 TYPE(pw_env_type), POINTER :: new_pw_env 297 TYPE(pw_p_type), POINTER :: rho_gs, rho_rs 298 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 299 300 my_tddft = .FALSE. 301 IF (PRESENT(tddft)) my_tddft = tddft 302 303 NULLIFY (new_pw_env) 304 CALL get_qs_env(qs_env, pw_env=new_pw_env) 305 CPASSERT(ASSOCIATED(new_pw_env)) 306 307 NULLIFY (auxbas_pw_pool) 308 CALL pw_env_get(new_pw_env, auxbas_pw_pool=auxbas_pw_pool) 309 CPASSERT(ASSOCIATED(auxbas_pw_pool)) 310 311 ! reallocate rho0 on the global grid in real and reciprocal space 312 NULLIFY (rho_rs, rho_gs) 313 CPASSERT(ASSOCIATED(rho0_mpole)) 314 rho_rs => rho0_mpole%rho0_s_rs 315 rho_gs => rho0_mpole%rho0_s_gs 316 317 ! rho0 density in real space 318 IF (ASSOCIATED(rho_rs)) THEN 319 CALL pw_release(rho_rs%pw) 320 DEALLOCATE (rho_rs) 321 END IF 322 ALLOCATE (rho_rs) 323 CALL pw_pool_create_pw(auxbas_pw_pool, rho_rs%pw, & 324 use_data=REALDATA3D, in_space=REALSPACE) 325 rho0_mpole%rho0_s_rs => rho_rs 326 327 ! rho0 density in reciprocal space 328 IF (ASSOCIATED(rho_gs)) THEN 329 CALL pw_release(rho_gs%pw) 330 DEALLOCATE (rho_gs) 331 END IF 332 ALLOCATE (rho_gs) 333 CALL pw_pool_create_pw(auxbas_pw_pool, rho_gs%pw, & 334 use_data=COMPLEXDATA1D) 335 rho_gs%pw%in_space = RECIPROCALSPACE 336 rho0_mpole%rho0_s_gs => rho_gs 337 338 IF (my_tddft) THEN 339 rho0_mpole%igrid_zet0_s = qs_env%local_rho_set%rho0_mpole%igrid_zet0_s 340 END IF 341 342 END SUBROUTINE rho0_s_grid_create 343 344! ************************************************************************************************** 345!> \brief ... 346!> \param qs_env ... 347!> \param v_rspace ... 348!> \param calculate_forces ... 349!> \param tddft ... 350!> \param do_triplet ... 351!> \param p_env ... 352! ************************************************************************************************** 353 SUBROUTINE integrate_vhg0_rspace(qs_env, v_rspace, calculate_forces, & 354 tddft, do_triplet, p_env) 355 356 TYPE(qs_environment_type), POINTER :: qs_env 357 TYPE(pw_p_type) :: v_rspace 358 LOGICAL, INTENT(IN) :: calculate_forces 359 LOGICAL, INTENT(IN), OPTIONAL :: tddft, do_triplet 360 TYPE(qs_p_env_type), OPTIONAL, POINTER :: p_env 361 362 CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace', & 363 routineP = moduleN//':'//routineN 364 365 INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, & 366 ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, lmax0, & 367 lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1, n2, & 368 nat, nch_ik, nch_max, ncurr, nset, nsotot, nspins, num_pe 369 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list 370 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list 371 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf 372 LOGICAL :: grid_distributed, my_tddft, paw_atom, & 373 use_virial 374 REAL(dp) :: c4pi, eps_rho_rspace, force_tmp(3), & 375 ra(3), rpgf0, scale, zet0 376 REAL(dp), DIMENSION(:), POINTER :: hab_sph, norm_l, Qlm 377 REAL(dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, pab 378 REAL(dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hadb, hdab, Qlm_gg 379 REAL(dp), DIMENSION(:, :, :, :), POINTER :: a_hdab 380 REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b 381 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 382 TYPE(cell_type), POINTER :: cell 383 TYPE(cp_para_env_type), POINTER :: para_env 384 TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info 385 TYPE(dft_control_type), POINTER :: dft_control 386 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 387 TYPE(harmonics_atom_type), POINTER :: harmonics 388 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 389 TYPE(pw_env_type), POINTER :: pw_env 390 TYPE(pw_p_type), POINTER :: coeff_gaux, coeff_gspace, coeff_raux, & 391 coeff_rspace 392 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 393 TYPE(pw_pool_type), POINTER :: pw_aux, pw_pool 394 TYPE(qs_force_type), DIMENSION(:), POINTER :: force 395 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 396 TYPE(realspace_grid_desc_p_type), DIMENSION(:), & 397 POINTER :: rs_descs 398 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 399 TYPE(realspace_grid_type), POINTER :: rs_v 400 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 401 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s 402 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 403 TYPE(rho_atom_type), POINTER :: rho_atom 404 TYPE(virial_type), POINTER :: virial 405 406 c4pi = fourpi 407 408 CALL timeset(routineN, handle) 409 410 NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, particle_set) 411 NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set) 412 413 my_tddft = .FALSE. 414 IF (PRESENT(tddft)) my_tddft = tddft 415 416 CALL get_qs_env(qs_env=qs_env, & 417 atomic_kind_set=atomic_kind_set, & 418 qs_kind_set=qs_kind_set, & 419 cell=cell, & 420 dft_control=dft_control, & 421 para_env=para_env, & 422 force=force, pw_env=pw_env, & 423 rho0_mpole=rho0_mpole, & 424 rho_atom_set=rho_atom_set, & 425 particle_set=particle_set, & 426 virial=virial) 427 428 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 429 430 nspins = dft_control%nspins 431 432 IF (my_tddft) THEN 433 IF (PRESENT(do_triplet)) THEN 434 IF (nspins == 1 .AND. do_triplet) RETURN 435 ELSE 436 IF (nspins == 1 .AND. dft_control%tddfpt_control%res_etype /= tddfpt_singlet) RETURN 437 ENDIF 438 END IF 439 440 IF (my_tddft) THEN 441 rho0_mpole => p_env%local_rho_set%rho0_mpole 442 rho_atom_set => p_env%local_rho_set%rho_atom_set 443 END IF 444 445 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, lmax_0=lmax0, & 446 zet0_h=zet0, igrid_zet0_s=igrid, & 447 norm_g0l_h=norm_l) 448 449 ! *** set up of the potential on the multigrids 450 NULLIFY (rs_descs, pw_pools) 451 CPASSERT(ASSOCIATED(pw_env)) 452 CALL pw_env_get(pw_env, rs_descs=rs_descs, pw_pools=pw_pools) 453 454 ! *** assign from pw_env 455 auxbas_grid = pw_env%auxbas_grid 456 cube_info => pw_env%cube_info 457 458 ! *** Get the potential on the right grid 459 NULLIFY (rs_v, rs_desc, pw_pool, pw_aux, coeff_rspace, coeff_gspace) 460 rs_desc => rs_descs(igrid)%rs_desc 461 pw_pool => pw_pools(igrid)%pool 462 463 ALLOCATE (coeff_gspace) 464 ALLOCATE (coeff_rspace) 465 CALL pw_pool_create_pw(pw_pool, coeff_gspace%pw, & 466 use_data=COMPLEXDATA1D, & 467 in_space=RECIPROCALSPACE) 468 469 CALL pw_pool_create_pw(pw_pool, coeff_rspace%pw, use_data=REALDATA3D, & 470 in_space=REALSPACE) 471 472 IF (igrid /= auxbas_grid) THEN 473 pw_aux => pw_pools(auxbas_grid)%pool 474 ALLOCATE (coeff_gaux) 475 CALL pw_pool_create_pw(pw_aux, coeff_gaux%pw, & 476 use_data=COMPLEXDATA1D, & 477 in_space=RECIPROCALSPACE) 478 CALL pw_transfer(v_rspace%pw, coeff_gaux%pw) 479 CALL pw_copy(coeff_gaux%pw, coeff_gspace%pw) 480 CALL pw_transfer(coeff_gspace%pw, coeff_rspace%pw) 481 CALL pw_pool_give_back_pw(pw_aux, coeff_gaux%pw) 482 DEALLOCATE (coeff_gaux) 483 ALLOCATE (coeff_raux) 484 CALL pw_pool_create_pw(pw_aux, coeff_raux%pw, use_data=REALDATA3D, & 485 in_space=REALSPACE) 486 scale = coeff_rspace%pw%pw_grid%dvol/coeff_raux%pw%pw_grid%dvol 487 coeff_rspace%pw%cr3d = scale*coeff_rspace%pw%cr3d 488 CALL pw_pool_give_back_pw(pw_aux, coeff_raux%pw) 489 DEALLOCATE (coeff_raux) 490 ELSE 491 492 IF (coeff_gspace%pw%pw_grid%spherical) THEN 493 CALL pw_transfer(v_rspace%pw, coeff_gspace%pw) 494 CALL pw_transfer(coeff_gspace%pw, coeff_rspace%pw) 495 ELSE 496 CALL pw_copy(v_rspace%pw, coeff_rspace%pw) 497 END IF 498 END IF 499 CALL pw_pool_give_back_pw(pw_pool, coeff_gspace%pw) 500 DEALLOCATE (coeff_gspace) 501 502 ! *** set up the rs grid at level igrid 503 CALL rs_grid_create(rs_v, rs_desc) 504 CALL rs_grid_zero(rs_v) 505 CALL rs_pw_transfer(rs_v, coeff_rspace%pw, pw2rs) 506 507 CALL pw_pool_give_back_pw(pw_pool, coeff_rspace%pw) 508 DEALLOCATE (coeff_rspace) 509 510 ! *** Now the potential is on the right grid => integration 511 512 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 513 514 ! *** Allocate work storage *** 515 516 NULLIFY (hab, hab_sph, hdab, hdab_sph, hadb, pab, a_hdab, a_hdab_sph) 517 nch_max = ncoset(lmax0) 518 CALL reallocate(hab, 1, nch_max, 1, 1) 519 CALL reallocate(hab_sph, 1, nch_max) 520 CALL reallocate(hdab, 1, 3, 1, nch_max, 1, 1) 521 CALL reallocate(hadb, 1, 3, 1, nch_max, 1, 1) 522 CALL reallocate(hdab_sph, 1, 3, 1, nch_max) 523 CALL reallocate(a_hdab, 1, 3, 1, 3, 1, nch_max, 1, 1) 524 CALL reallocate(a_hdab_sph, 1, 3, 1, 3, 1, nch_max) 525 CALL reallocate(pab, 1, nch_max, 1, 1) 526 527 ncurr = -1 528 529 grid_distributed = rs_v%desc%distributed 530 531 DO ikind = 1, SIZE(atomic_kind_set, 1) 532 NULLIFY (orb_basis_set, atom_list, harmonics) 533 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat) 534 CALL get_qs_kind(qs_kind_set(ikind), & 535 basis_set=orb_basis_set, & 536 paw_atom=paw_atom, & 537 harmonics=harmonics) 538 539 IF (.NOT. paw_atom) CYCLE 540 541 NULLIFY (Qlm_gg, lmax, npgf) 542 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, & 543 l0_ikind=l0_ikind, Qlm_gg=Qlm_gg, & 544 rpgf0_s=rpgf0) 545 546 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 547 lmax=lmax, lmin=lmin, & 548 maxso=maxso, maxl=maxl, & 549 nset=nset, npgf=npgf) 550 551 nsotot = maxso*nset 552 ALLOCATE (intloc(nsotot, nsotot)) 553 554 ! Initialize the local KS integrals 555 556 nch_ik = ncoset(l0_ikind) 557 pab = 1.0_dp 558 max_s_harm = harmonics%max_s_harm 559 llmax = harmonics%llmax 560 561 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm)) 562 563 num_pe = para_env%num_pe 564 mepos = para_env%mepos 565 DO j = 0, num_pe - 1 566 bo = get_limit(nat, num_pe, j) 567 IF (.NOT. grid_distributed .AND. j /= mepos) CYCLE 568 569 DO iat = bo(1), bo(2) 570 iatom = atom_list(iat) 571 ra(:) = pbc(particle_set(iatom)%r, cell) 572 573 NULLIFY (Qlm) 574 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, iat=iatom, Qlm_tot=Qlm) 575 576 hab = 0.0_dp 577 hdab = 0.0_dp 578 intloc = 0._dp 579 IF (use_virial) THEN 580 my_virial_a = 0.0_dp 581 my_virial_b = 0.0_dp 582 a_hdab = 0.0_dp 583 END IF 584 585 CALL integrate_pgf_product_rspace( & 586 l0_ikind, zet0, 0, 0, 0.0_dp, 0, & 587 ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 0.0_dp, rs_v, cell, & 588 cube_info(igrid), hab, pab, o1=0, o2=0, & 589 eps_gvg_rspace=eps_rho_rspace, & 590 calculate_forces=calculate_forces, & 591 hdab=hdab, hadb=hadb, & 592 collocate_rho0=.TRUE., rpgf0_s=rpgf0, & 593 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b, & 594 a_hdab=a_hdab, use_subpatch=.TRUE., subpatch_pattern=0_int_8) 595 596 ! Convert from cartesian to spherical 597 DO lshell = 0, l0_ikind 598 DO is = 1, nso(lshell) 599 iso = is + nsoset(lshell - 1) 600 hab_sph(iso) = 0.0_dp 601 hdab_sph(1:3, iso) = 0.0_dp 602 a_hdab_sph(1:3, 1:3, iso) = 0.0_dp 603 DO ic = 1, nco(lshell) 604 ico = ic + ncoset(lshell - 1) 605 lx = indco(1, ico) 606 ly = indco(2, ico) 607 lz = indco(3, ico) 608 609 hab_sph(iso) = hab_sph(iso) + & 610 orbtramat(lshell)%c2s(is, ic)*hab(ico, 1)* & 611 norm_l(lshell)/ & 612 SQRT(c4pi*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)/ & 613 dfac(2*lshell + 1)) 614 615 IF (calculate_forces) THEN 616 hdab_sph(1:3, iso) = hdab_sph(1:3, iso) + & 617 orbtramat(lshell)%c2s(is, ic)*hdab(1:3, ico, 1)* & 618 norm_l(lshell)/ & 619 SQRT(c4pi*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)/ & 620 dfac(2*lshell + 1)) 621 END IF 622 IF (use_virial) THEN 623 DO ii = 1, 3 624 DO i = 1, 3 625 a_hdab_sph(i, ii, iso) = a_hdab_sph(i, ii, iso) + & 626 orbtramat(lshell)%c2s(is, ic)*a_hdab(i, ii, ico, 1)* & 627 norm_l(lshell)/ & 628 SQRT(c4pi*dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)/ & 629 dfac(2*lshell + 1)) 630 END DO 631 END DO 632 END IF 633 634 END DO ! ic 635 END DO ! is 636 END DO ! lshell 637 638 m1 = 0 639 DO iset1 = 1, nset 640 641 m2 = 0 642 DO iset2 = 1, nset 643 CALL get_none0_cg_list(harmonics%my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), & 644 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local) 645 n1 = nsoset(lmax(iset1)) 646 DO ipgf1 = 1, npgf(iset1) 647 n2 = nsoset(lmax(iset2)) 648 DO ipgf2 = 1, npgf(iset2) 649 650 DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local) 651 DO icg = 1, cg_n_list(iso) 652 iso1 = cg_list(1, icg, iso) 653 iso2 = cg_list(2, icg, iso) 654 655 ig1 = iso1 + n1*(ipgf1 - 1) + m1 656 ig2 = iso2 + n2*(ipgf2 - 1) + m2 657 658 intloc(ig1, ig2) = intloc(ig1, ig2) + Qlm_gg(ig1, ig2, iso)*hab_sph(iso) 659 660 END DO ! icg 661 END DO ! iso 662 663 END DO ! ipgf2 664 END DO ! ipgf1 665 m2 = m2 + maxso 666 END DO ! iset2 667 m1 = m1 + maxso 668 END DO ! iset1 669 670 IF (grid_distributed) THEN 671 ! sum result over all processors 672 CALL mp_sum(intloc, para_env%group) 673 END IF 674 675 IF (j == mepos) THEN 676 rho_atom => rho_atom_set(iatom) 677 CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_local_h, ga_Vlocal_gb_s=int_local_s) 678 DO ispin = 1, nspins 679 int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc 680 int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc 681 END DO 682 END IF 683 684 IF (calculate_forces) THEN 685 force_tmp(1:3) = 0.0_dp 686 DO iso = 1, nsoset(l0_ikind) 687 force_tmp(1) = force_tmp(1) + Qlm(iso)*hdab_sph(1, iso) 688 force_tmp(2) = force_tmp(2) + Qlm(iso)*hdab_sph(2, iso) 689 force_tmp(3) = force_tmp(3) + Qlm(iso)*hdab_sph(3, iso) 690 END DO 691 force(ikind)%g0s_Vh_elec(1:3, iat) = force(ikind)%g0s_Vh_elec(1:3, iat) + force_tmp(1:3) 692 END IF 693 IF (use_virial) THEN 694 my_virial_a = 0.0_dp 695 DO iso = 1, nsoset(l0_ikind) 696 DO ii = 1, 3 697 DO i = 1, 3 698 virial%pv_virial(i, ii) = virial%pv_virial(i, ii) + Qlm(iso)*a_hdab_sph(i, ii, iso) 699 END DO 700 END DO 701 END DO 702 END IF 703 704 END DO 705 END DO 706 707 DEALLOCATE (intloc) 708 DEALLOCATE (cg_list, cg_n_list) 709 710 END DO ! ikind 711 712 CALL rs_grid_release(rs_v) 713 714 DEALLOCATE (hab, hdab, hadb, hab_sph, hdab_sph, pab, a_hdab, a_hdab_sph) 715 716 CALL timestop(handle) 717 718 END SUBROUTINE integrate_vhg0_rspace 719 720END MODULE qs_rho0_ggrid 721