1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE qs_rho0_methods 8 9 USE ao_util, ONLY: exp_radius,& 10 gaussint_sph,& 11 trace_r_AxB 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind,& 14 get_atomic_kind_set 15 USE basis_set_types, ONLY: get_gto_basis_set,& 16 gto_basis_set_type 17 USE cp_control_types, ONLY: dft_control_type,& 18 gapw_control_type 19 USE cp_log_handling, ONLY: cp_get_default_logger,& 20 cp_logger_type 21 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 22 cp_print_key_unit_nr 23 USE input_section_types, ONLY: section_vals_get_subs_vals,& 24 section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: default_string_length,& 27 dp 28 USE mathconstants, ONLY: dfac,& 29 fourpi 30 USE memory_utilities, ONLY: reallocate 31 USE orbital_pointers, ONLY: indco,& 32 indso,& 33 nco,& 34 ncoset,& 35 nso,& 36 nsoset 37 USE orbital_transformation_matrices, ONLY: orbtramat 38 USE qs_environment_types, ONLY: get_qs_env,& 39 qs_environment_type,& 40 set_qs_env 41 USE qs_grid_atom, ONLY: grid_atom_type 42 USE qs_harmonics_atom, ONLY: get_none0_cg_list,& 43 harmonics_atom_type 44 USE qs_kind_types, ONLY: get_qs_kind,& 45 qs_kind_type,& 46 set_qs_kind 47 USE qs_local_rho_types, ONLY: allocate_rhoz,& 48 calculate_rhoz,& 49 local_rho_type,& 50 rhoz_type 51 USE qs_oce_methods, ONLY: prj_scatter 52 USE qs_rho0_ggrid, ONLY: rho0_s_grid_create 53 USE qs_rho0_types, ONLY: & 54 allocate_multipoles, allocate_rho0_atom, allocate_rho0_atom_rad, allocate_rho0_mpole, & 55 calculate_g0, get_rho0_mpole, initialize_mpole_rho, mpole_gau_overlap, mpole_rho_atom, & 56 rho0_atom_type, rho0_mpole_type, write_rho0_info 57 USE qs_rho_atom_types, ONLY: get_rho_atom,& 58 rho_atom_coeff,& 59 rho_atom_type 60#include "./base/base_uses.f90" 61 62 IMPLICIT NONE 63 64 PRIVATE 65 66! *** Global parameters (only in this module) 67 68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_methods' 69 70! *** Public subroutines *** 71 72 PUBLIC :: calculate_rho0_atom, init_rho0 73 74CONTAINS 75 76! ************************************************************************************************** 77!> \brief ... 78!> \param mp_gau ... 79!> \param orb_basis ... 80!> \param harmonics ... 81!> \param nchannels ... 82!> \param nsotot ... 83! ************************************************************************************************** 84 SUBROUTINE calculate_mpole_gau(mp_gau, orb_basis, harmonics, nchannels, nsotot) 85 86 TYPE(mpole_gau_overlap) :: mp_gau 87 TYPE(gto_basis_set_type), POINTER :: orb_basis 88 TYPE(harmonics_atom_type), POINTER :: harmonics 89 INTEGER, INTENT(IN) :: nchannels, nsotot 90 91 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_mpole_gau', & 92 routineP = moduleN//':'//routineN 93 94 INTEGER :: handle, icg, ig1, ig2, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso2, l, l1, l2, & 95 llmax, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, nset 96 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list 97 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list 98 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf 99 REAL(dp) :: zet1, zet2 100 REAL(dp), DIMENSION(:, :), POINTER :: zet 101 REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG 102 103 CALL timeset(routineN, handle) 104 105 NULLIFY (lmax, lmin, npgf, my_CG, zet) 106 107 CALL reallocate(mp_gau%Qlm_gg, 1, nsotot, 1, nsotot, 1, nchannels) 108 109 CALL get_gto_basis_set(gto_basis_set=orb_basis, & 110 lmax=lmax, lmin=lmin, maxso=maxso, & 111 npgf=npgf, nset=nset, zet=zet, maxl=maxl) 112 113 max_s_harm = harmonics%max_s_harm 114 llmax = harmonics%llmax 115 116 ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm)) 117 118 my_CG => harmonics%my_CG 119 120 m1 = 0 121 DO iset1 = 1, nset 122 m2 = 0 123 DO iset2 = 1, nset 124 125 CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), & 126 max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local) 127 128 n1 = nsoset(lmax(iset1)) 129 DO ipgf1 = 1, npgf(iset1) 130 zet1 = zet(ipgf1, iset1) 131 132 n2 = nsoset(lmax(iset2)) 133 DO ipgf2 = 1, npgf(iset2) 134 zet2 = zet(ipgf2, iset2) 135 136 DO iso = 1, MIN(nchannels, max_iso_not0_local) 137 l = indso(1, iso) 138 DO icg = 1, cg_n_list(iso) 139 iso1 = cg_list(1, icg, iso) 140 iso2 = cg_list(2, icg, iso) 141 142 l1 = indso(1, iso1) 143 l2 = indso(1, iso2) 144 ig1 = iso1 + n1*(ipgf1 - 1) + m1 145 ig2 = iso2 + n2*(ipgf2 - 1) + m2 146 147 mp_gau%Qlm_gg(ig1, ig2, iso) = fourpi/(2._dp*l + 1._dp)* & 148 my_CG(iso1, iso2, iso)*gaussint_sph(zet1 + zet2, l + l1 + l2) 149 END DO ! icg 150 END DO ! iso 151 152 END DO ! ipgf2 153 END DO ! ipgf1 154 m2 = m2 + maxso 155 END DO ! iset2 156 m1 = m1 + maxso 157 END DO ! iset1 158 159 DEALLOCATE (cg_list, cg_n_list) 160 161 CALL timestop(handle) 162 END SUBROUTINE calculate_mpole_gau 163 164! ************************************************************************************************** 165!> \brief ... 166!> \param gapw_control ... 167!> \param rho_atom_set ... 168!> \param rho0_atom_set ... 169!> \param rho0_mp ... 170!> \param a_list ... 171!> \param g_atom ... 172!> \param paw_atom ... 173!> \param natom ... 174!> \param ikind ... 175!> \param qs_kind ... 176!> \param harmonics ... 177!> \param rho0_h_tot ... 178! ************************************************************************************************** 179 SUBROUTINE calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, & 180 rho0_mp, a_list, g_atom, & 181 paw_atom, natom, ikind, qs_kind, harmonics, & 182 rho0_h_tot) 183 184 TYPE(gapw_control_type), POINTER :: gapw_control 185 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set 186 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set 187 TYPE(rho0_mpole_type), POINTER :: rho0_mp 188 INTEGER, DIMENSION(:), INTENT(IN) :: a_list 189 TYPE(grid_atom_type), INTENT(IN) :: g_atom 190 LOGICAL, INTENT(IN) :: paw_atom 191 INTEGER, INTENT(IN) :: natom, ikind 192 TYPE(qs_kind_type), INTENT(IN) :: qs_kind 193 TYPE(harmonics_atom_type), POINTER :: harmonics 194 REAL(dp), INTENT(INOUT) :: rho0_h_tot 195 196 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rho0_atom', & 197 routineP = moduleN//':'//routineN 198 199 INTEGER :: handle, iat, iatom, ic, ico, ir, is, & 200 iso, ispin, l, lmax0, lshell, lx, ly, & 201 lz, nr, nsotot, nspins 202 REAL(dp) :: sum1 203 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc_ah, cpc_as 204 REAL(dp), DIMENSION(:), POINTER :: norm_g0l_h 205 REAL(dp), DIMENSION(:, :), POINTER :: g0_h, vg0_h 206 TYPE(mpole_gau_overlap), POINTER :: mpole_gau 207 TYPE(mpole_rho_atom), POINTER :: mpole_rho 208 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s 209 TYPE(rho_atom_type), POINTER :: rho_atom 210 211 CALL timeset(routineN, handle) 212 213 NULLIFY (mpole_gau) 214 NULLIFY (mpole_rho) 215 NULLIFY (g0_h, vg0_h) 216 NULLIFY (norm_g0l_h) 217 218 CALL get_rho0_mpole(rho0_mpole=rho0_mp, ikind=ikind, & 219 l0_ikind=lmax0, mp_gau_ikind=mpole_gau, & 220 g0_h=g0_h, & 221 vg0_h=vg0_h, & 222 norm_g0l_h=norm_g0l_h) 223 224 nr = g_atom%nr 225 226! Set density coefficient to zero before the calculation 227 DO iat = 1, natom 228 iatom = a_list(iat) 229 rho0_atom_set(iatom)%rho0_rad_h%r_coef = 0.0_dp 230 rho0_mp%mp_rho(iatom)%Qlm_tot = 0.0_dp 231 rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z 232 rho0_mp%mp_rho(iatom)%Q0 = 0.0_dp 233 rho0_mp%mp_rho(iatom)%Qlm_car = 0.0_dp 234 ENDDO 235 236 IF (.NOT. (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw)) THEN 237 DO iat = 1, natom 238 iatom = a_list(iat) 239 mpole_rho => rho0_mp%mp_rho(iatom) 240 rho_atom => rho_atom_set(iatom) 241 242 IF (paw_atom) THEN 243 NULLIFY (cpc_h, cpc_s) 244 CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s) 245 nspins = SIZE(cpc_h) 246 nsotot = SIZE(mpole_gau%Qlm_gg, 1) 247 ALLOCATE (cpc_ah(nsotot, nsotot, nspins)) 248 cpc_ah = 0._dp 249 ALLOCATE (cpc_as(nsotot, nsotot, nspins)) 250 cpc_as = 0._dp 251 DO ispin = 1, nspins 252 CALL prj_scatter(cpc_h(ispin)%r_coef, cpc_ah(:, :, ispin), qs_kind) 253 CALL prj_scatter(cpc_s(ispin)%r_coef, cpc_as(:, :, ispin), qs_kind) 254 END DO 255 END IF 256 257 ! Total charge (hard-soft) at atom 258 IF (paw_atom) THEN 259 DO ispin = 1, nspins 260 mpole_rho%Q0(ispin) = (trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, & 261 cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) & 262 - trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, & 263 cpc_as(:, :, ispin), nsotot, nsotot, nsotot))/SQRT(fourpi) 264 END DO 265 END IF 266 ! Multipoles of local charge distribution 267 DO iso = 1, nsoset(lmax0) 268 l = indso(1, iso) 269 IF (paw_atom) THEN 270 mpole_rho%Qlm_h(iso) = 0.0_dp 271 mpole_rho%Qlm_s(iso) = 0.0_dp 272 273 DO ispin = 1, nspins 274 mpole_rho%Qlm_h(iso) = mpole_rho%Qlm_h(iso) + & 275 trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, & 276 cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) 277 mpole_rho%Qlm_s(iso) = mpole_rho%Qlm_s(iso) + & 278 trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, & 279 cpc_as(:, :, ispin), nsotot, nsotot, nsotot) 280 END DO ! ispin 281 282 mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) + & 283 mpole_rho%Qlm_h(iso) - mpole_rho%Qlm_s(iso) 284 END IF 285 286 rho0_atom_set(iatom)%rho0_rad_h%r_coef(1:nr, iso) = & 287 g0_h(1:nr, l)*mpole_rho%Qlm_tot(iso) 288 rho0_atom_set(iatom)%vrho0_rad_h%r_coef(1:nr, iso) = & 289 vg0_h(1:nr, l)*mpole_rho%Qlm_tot(iso) 290 291 sum1 = 0.0_dp 292 DO ir = 1, nr 293 sum1 = sum1 + g_atom%wr(ir)* & 294 rho0_atom_set(iatom)%rho0_rad_h%r_coef(ir, iso) 295 ENDDO 296 rho0_h_tot = rho0_h_tot + sum1*harmonics%slm_int(iso) 297 END DO ! iso 298 IF (paw_atom) THEN 299 DEALLOCATE (cpc_ah, cpc_as) 300 END IF 301 END DO ! iat 302 END IF 303 304! Transform the coefficinets from spherical to cartesian 305 IF (.NOT. paw_atom .AND. gapw_control%nopaw_as_gpw) THEN 306 DO iat = 1, natom 307 iatom = a_list(iat) 308 mpole_rho => rho0_mp%mp_rho(iatom) 309 310 DO lshell = 0, lmax0 311 DO ic = 1, nco(lshell) 312 ico = ic + ncoset(lshell - 1) 313 mpole_rho%Qlm_car(ico) = 0.0_dp 314 END DO 315 END DO 316 END DO 317 ELSE 318 DO iat = 1, natom 319 iatom = a_list(iat) 320 mpole_rho => rho0_mp%mp_rho(iatom) 321 322 DO lshell = 0, lmax0 323 DO ic = 1, nco(lshell) 324 ico = ic + ncoset(lshell - 1) 325 mpole_rho%Qlm_car(ico) = 0.0_dp 326 lx = indco(1, ico) 327 ly = indco(2, ico) 328 lz = indco(3, ico) 329 330 DO is = 1, nso(lshell) 331 iso = is + nsoset(lshell - 1) 332 333 mpole_rho%Qlm_car(ico) = mpole_rho%Qlm_car(ico) + & 334 orbtramat(lshell)%c2s(is, ic)*mpole_rho%Qlm_tot(iso)* & 335 norm_g0l_h(lshell) & 336 /SQRT(dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)*fourpi/dfac(2*lshell + 1)) 337 338 END DO 339 END DO 340 END DO ! lshell 341 END DO ! iat 342 END IF 343!MI get rid of full gapw 344 345 CALL timestop(handle) 346 347 END SUBROUTINE calculate_rho0_atom 348 349! ************************************************************************************************** 350!> \brief ... 351!> \param qs_env ... 352!> \param gapw_control ... 353!> \param tddft ... 354!> \param tddft_local_rho_set ... 355! ************************************************************************************************** 356 SUBROUTINE init_rho0(qs_env, gapw_control, & 357 tddft, tddft_local_rho_set) 358 359 TYPE(qs_environment_type), POINTER :: qs_env 360 TYPE(gapw_control_type), POINTER :: gapw_control 361 LOGICAL, INTENT(IN), OPTIONAL :: tddft 362 TYPE(local_rho_type), OPTIONAL, POINTER :: tddft_local_rho_set 363 364 CHARACTER(len=*), PARAMETER :: routineN = 'init_rho0', routineP = moduleN//':'//routineN 365 366 CHARACTER(LEN=default_string_length) :: unit_str 367 INTEGER :: handle, iat, iatom, ikind, l, l_rho1_max, laddg, lmaxg, maxl, maxnset, maxso, & 368 nat, natom, nchan_c, nchan_s, nkind, nr, nset, nsotot, output_unit 369 INTEGER, DIMENSION(:), POINTER :: atom_list 370 LOGICAL :: my_tddft, paw_atom 371 REAL(dp) :: alpha_core, eps_Vrho0, max_rpgf0_s, & 372 radius, rc_min, rc_orb, & 373 total_rho_core_rspace 374 REAL(KIND=dp) :: zeff 375 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 376 TYPE(cp_logger_type), POINTER :: logger 377 TYPE(dft_control_type), POINTER :: dft_control 378 TYPE(grid_atom_type), POINTER :: grid_atom 379 TYPE(gto_basis_set_type), POINTER :: orb_basis 380 TYPE(harmonics_atom_type), POINTER :: harmonics 381 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 382 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set 383 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 384 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set 385 TYPE(section_vals_type), POINTER :: dft_section 386 387 CALL timeset(routineN, handle) 388 389 NULLIFY (logger) 390 logger => cp_get_default_logger() 391 392 NULLIFY (qs_kind_set) 393 NULLIFY (atomic_kind_set) 394 NULLIFY (dft_control) 395 NULLIFY (harmonics) 396 NULLIFY (orb_basis) 397 NULLIFY (rho0_mpole) 398 NULLIFY (rho0_atom_set) 399 NULLIFY (rhoz_set) 400 401 my_tddft = .FALSE. 402 IF (PRESENT(tddft)) my_tddft = tddft 403 404 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, & 405 atomic_kind_set=atomic_kind_set, & 406 dft_control=dft_control) 407 408 nkind = SIZE(atomic_kind_set) 409 eps_Vrho0 = gapw_control%eps_Vrho0 410 411! Initialize rhoz total to zero 412! in gapw rhoz is calculated on local the lebedev grids 413 total_rho_core_rspace = 0.0_dp 414 415 CALL get_atomic_kind_set(atomic_kind_set, natom=natom) 416 417! Initialize the multipole and the compensation charge type 418 CALL allocate_rho0_mpole(rho0_mpole) 419 CALL allocate_rho0_atom(rho0_atom_set, natom) 420 421! Allocate the multipole set 422 CALL allocate_multipoles(rho0_mpole%mp_rho, natom, rho0_mpole%mp_gau, nkind) 423 424! Allocate the core density on the radial grid for each kind: rhoz_set 425 CALL allocate_rhoz(rhoz_set, nkind) 426 427! For each kind, determine the max l for the compensation charge density 428 lmaxg = gapw_control%lmax_rho0 429 laddg = gapw_control%ladd_rho0 430 431 CALL reallocate(rho0_mpole%lmax0_kind, 1, nkind) 432 433 rho0_mpole%lmax_0 = 0 434 rc_min = 100.0_dp 435 maxnset = 0 436 DO ikind = 1, nkind 437 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat) 438 CALL get_qs_kind(qs_kind_set(ikind), & 439 basis_set=orb_basis, & 440 ngrid_rad=nr, & 441 grid_atom=grid_atom, & 442 harmonics=harmonics, & 443 paw_atom=paw_atom, & 444 hard0_radius=rc_orb, & 445 zeff=zeff, & 446 alpha_core_charge=alpha_core) 447 448 CALL get_gto_basis_set(gto_basis_set=orb_basis, & 449 maxl=maxl, & 450 maxso=maxso, nset=nset) 451 452 maxnset = MAX(maxnset, nset) 453 454 l_rho1_max = indso(1, harmonics%max_iso_not0) 455 IF (paw_atom) THEN 456 rho0_mpole%lmax0_kind(ikind) = MIN(2*maxl, l_rho1_max, maxl + laddg, lmaxg) 457 ELSE 458 rho0_mpole%lmax0_kind(ikind) = 0 459 END IF 460 461 CALL set_qs_kind(qs_kind_set(ikind), lmax_rho0=rho0_mpole%lmax0_kind(ikind)) 462 463 IF (gapw_control%lrho1_eq_lrho0) harmonics%max_iso_not0 = & 464 nsoset(rho0_mpole%lmax0_kind(ikind)) 465 466 rho0_mpole%lmax_0 = MAX(rho0_mpole%lmax_0, rho0_mpole%lmax0_kind(ikind)) 467 rc_min = MIN(rc_min, rc_orb) 468 469 nchan_s = nsoset(rho0_mpole%lmax0_kind(ikind)) 470 nchan_c = ncoset(rho0_mpole%lmax0_kind(ikind)) 471 nsotot = maxso*nset 472 473 DO iat = 1, nat 474 iatom = atom_list(iat) 475! Allocate the multipole for rho1_h rho1_s and rho_z 476 CALL initialize_mpole_rho(rho0_mpole%mp_rho(iatom), nchan_s, nchan_c, zeff, my_tddft) 477! Allocate the radial part of rho0_h and rho0_s 478! This is calculated on the radial grid centered at the atomic position 479 CALL allocate_rho0_atom_rad(rho0_atom_set(iatom), nr, nchan_s) 480 END DO 481! 482 IF (paw_atom) THEN 483! Calculate multipoles given by the product of 2 primitives Qlm_gg 484 CALL calculate_mpole_gau(rho0_mpole%mp_gau(ikind), & 485 orb_basis, harmonics, nchan_s, nsotot) 486 END IF 487 488! Calculate the core density rhoz 489! exp(-alpha_c**2 r**2)Z(alpha_c**2/pi)**(3/2) 490! on the logarithmic radial grid 491! WARNING: alpha_core_charge = alpha_c**2 492 CALL calculate_rhoz(rhoz_set(ikind), grid_atom, alpha_core, zeff, & 493 nat, total_rho_core_rspace, harmonics) 494 END DO ! ikind 495 total_rho_core_rspace = -total_rho_core_rspace 496 497 IF (gapw_control%alpha0_hard_from_input) THEN 498! The Exponent for the compensation charge rho0_hard is read from input 499 rho0_mpole%zet0_h = gapw_control%alpha0_hard 500 ELSE 501! Calculate the exponent for the compensation charge rho0_hard 502 rho0_mpole%zet0_h = 0.1_dp 503 DO 504 radius = exp_radius(rho0_mpole%lmax_0, rho0_mpole%zet0_h, eps_Vrho0, 1.0_dp) 505 IF (radius <= rc_min) EXIT 506 rho0_mpole%zet0_h = rho0_mpole%zet0_h + 0.1_dp 507 END DO 508 509 END IF 510 511! Allocate and calculate the normalization factors for g0_lm_h and g0_lm_s 512 CALL reallocate(rho0_mpole%norm_g0l_h, 0, rho0_mpole%lmax_0) 513 DO l = 0, rho0_mpole%lmax_0 514 rho0_mpole%norm_g0l_h(l) = (2._dp*l + 1._dp)/ & 515 (fourpi*gaussint_sph(rho0_mpole%zet0_h, 2*l)) 516 END DO 517 518! Allocate and Initialize the g0 gaussians used to build the compensation density 519! and calculate the interaction radii 520 max_rpgf0_s = 0.0_dp 521 DO ikind = 1, nkind 522 CALL get_qs_kind(qs_kind_set(ikind), grid_atom=grid_atom) 523 CALL calculate_g0(rho0_mpole, grid_atom, ikind) 524 CALL interaction_radii_g0(rho0_mpole, ikind, eps_Vrho0, max_rpgf0_s) 525 END DO 526 rho0_mpole%max_rpgf0_s = max_rpgf0_s 527 528 IF (.NOT. my_tddft) THEN 529 CALL set_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole, & 530 rho0_atom_set=rho0_atom_set, & 531 rhoz_set=rhoz_set, & 532 rhoz_tot=total_rho_core_rspace) 533 ELSE 534 tddft_local_rho_set%rho0_mpole => rho0_mpole 535 tddft_local_rho_set%rho0_atom_set => rho0_atom_set 536 tddft_local_rho_set%rhoz_set => rhoz_set 537 tddft_local_rho_set%rhoz_tot = total_rho_core_rspace 538 CALL rho0_s_grid_create(qs_env, rho0_mpole, .TRUE.) 539 END IF 540 541 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") 542 output_unit = cp_print_key_unit_nr(logger, dft_section, "PRINT%GAPW%RHO0_INFORMATION", & 543 extension=".Log") 544 CALL section_vals_val_get(dft_section, "PRINT%GAPW%RHO0_INFORMATION%UNIT", c_val=unit_str) 545 IF (output_unit > 0) THEN 546 CALL write_rho0_info(rho0_mpole, unit_str, output_unit) 547 END IF 548 CALL cp_print_key_finished_output(output_unit, logger, dft_section, & 549 "PRINT%GAPW%RHO0_INFORMATION") 550 551 CALL timestop(handle) 552 553 END SUBROUTINE init_rho0 554 555! ************************************************************************************************** 556!> \brief ... 557!> \param rho0_mpole ... 558!> \param ik ... 559!> \param eps_Vrho0 ... 560!> \param max_rpgf0_s ... 561! ************************************************************************************************** 562 SUBROUTINE interaction_radii_g0(rho0_mpole, ik, eps_Vrho0, max_rpgf0_s) 563 564 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 565 INTEGER, INTENT(IN) :: ik 566 REAL(dp), INTENT(IN) :: eps_Vrho0 567 REAL(dp), INTENT(INOUT) :: max_rpgf0_s 568 569 CHARACTER(len=*), PARAMETER :: routineN = 'interaction_radii_g0', & 570 routineP = moduleN//':'//routineN 571 572 INTEGER :: l, lmax 573 REAL(dp) :: r_h, z0_h 574 REAL(dp), DIMENSION(:), POINTER :: ng0_h 575 576 CALL get_rho0_mpole(rho0_mpole, ikind=ik, l0_ikind=lmax, & 577 zet0_h=z0_h, norm_g0l_h=ng0_h) 578 r_h = 0.0_dp 579 DO l = 0, lmax 580 r_h = MAX(r_h, exp_radius(l, z0_h, eps_Vrho0, ng0_h(l))) 581 END DO 582 583 rho0_mpole%mp_gau(ik)%rpgf0_h = r_h 584 rho0_mpole%mp_gau(ik)%rpgf0_s = r_h 585 max_rpgf0_s = MAX(max_rpgf0_s, r_h) 586 587 END SUBROUTINE interaction_radii_g0 588 589! ************************************************************************************************** 590 591END MODULE qs_rho0_methods 592 593