1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6MODULE qs_rho0_types 7 8 USE cp_units, ONLY: cp_unit_from_cp2k 9 USE kinds, ONLY: dp 10 USE mathconstants, ONLY: fourpi,& 11 pi,& 12 rootpi 13 USE memory_utilities, ONLY: reallocate 14 USE pw_types, ONLY: pw_p_type,& 15 pw_release 16 USE qs_grid_atom, ONLY: grid_atom_type 17 USE qs_rho_atom_types, ONLY: rho_atom_coeff 18 USE whittaker, ONLY: whittaker_c0a,& 19 whittaker_ci 20#include "./base/base_uses.f90" 21 22 IMPLICIT NONE 23 24 PRIVATE 25 26! *** Global parameters (only in this module) 27 28 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_types' 29 30! *** Define multipole type *** 31 32! ************************************************************************************************** 33 TYPE mpole_rho_atom 34 REAL(dp), DIMENSION(:), POINTER :: Qlm_h, & 35 Qlm_s, & 36 Qlm_tot, & 37 Qlm_car 38 REAL(dp) :: Qlm_z 39 REAL(dp), DIMENSION(2) :: Q0 40 END TYPE mpole_rho_atom 41 42! ************************************************************************************************** 43 TYPE mpole_gau_overlap 44 REAL(dp), DIMENSION(:, :, :), POINTER :: Qlm_gg 45 REAL(dp), DIMENSION(:, :), POINTER :: g0_h, vg0_h 46 REAL(dp) :: rpgf0_h, rpgf0_s 47 END TYPE mpole_gau_overlap 48 49! ************************************************************************************************** 50 TYPE rho0_mpole_type 51 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho 52 TYPE(mpole_gau_overlap), DIMENSION(:), & 53 POINTER :: mp_gau 54 REAL(dp) :: zet0_h, & 55 total_rho0_h 56 REAL(dp) :: max_rpgf0_s 57 REAL(dp), DIMENSION(:), POINTER :: norm_g0l_h 58 INTEGER, DIMENSION(:), POINTER :: lmax0_kind 59 INTEGER :: lmax_0, igrid_zet0_s 60 TYPE(pw_p_type), POINTER :: rho0_s_rs, & 61 rho0_s_gs 62 END TYPE rho0_mpole_type 63 64! ************************************************************************************************** 65 TYPE rho0_atom_type 66 TYPE(rho_atom_coeff), POINTER :: rho0_rad_h, & 67 vrho0_rad_h 68 END TYPE rho0_atom_type 69 70! Public Types 71 72 PUBLIC :: mpole_rho_atom, mpole_gau_overlap, & 73 rho0_atom_type, rho0_mpole_type 74 75! Public Subroutine 76 77 PUBLIC :: allocate_multipoles, allocate_rho0_mpole, & 78 allocate_rho0_atom, allocate_rho0_atom_rad, & 79 deallocate_rho0_atom, deallocate_rho0_mpole, & 80 calculate_g0, get_rho0_mpole, initialize_mpole_rho, & 81 write_rho0_info 82 83CONTAINS 84 85! ************************************************************************************************** 86!> \brief ... 87!> \param mp_rho ... 88!> \param natom ... 89!> \param mp_gau ... 90!> \param nkind ... 91! ************************************************************************************************** 92 SUBROUTINE allocate_multipoles(mp_rho, natom, mp_gau, nkind) 93 94 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho 95 INTEGER, INTENT(IN) :: natom 96 TYPE(mpole_gau_overlap), DIMENSION(:), POINTER :: mp_gau 97 INTEGER, INTENT(IN) :: nkind 98 99 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_multipoles', & 100 routineP = moduleN//':'//routineN 101 102 INTEGER :: iat, ikind 103 104 IF (ASSOCIATED(mp_rho)) THEN 105 CALL deallocate_mpole_rho(mp_rho) 106 END IF 107 108 ALLOCATE (mp_rho(natom)) 109 110 DO iat = 1, natom 111 NULLIFY (mp_rho(iat)%Qlm_h) 112 NULLIFY (mp_rho(iat)%Qlm_s) 113 NULLIFY (mp_rho(iat)%Qlm_tot) 114 NULLIFY (mp_rho(iat)%Qlm_car) 115 END DO 116 117 IF (ASSOCIATED(mp_gau)) THEN 118 CALL deallocate_mpole_gau(mp_gau) 119 END IF 120 121 ALLOCATE (mp_gau(nkind)) 122 123 DO ikind = 1, nkind 124 NULLIFY (mp_gau(ikind)%Qlm_gg) 125 NULLIFY (mp_gau(ikind)%g0_h) 126 NULLIFY (mp_gau(ikind)%vg0_h) 127 mp_gau(ikind)%rpgf0_h = 0.0_dp 128 mp_gau(ikind)%rpgf0_s = 0.0_dp 129 END DO 130 131 END SUBROUTINE allocate_multipoles 132 133! ************************************************************************************************** 134!> \brief ... 135!> \param rho0_set ... 136!> \param natom ... 137! ************************************************************************************************** 138 SUBROUTINE allocate_rho0_atom(rho0_set, natom) 139 140 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_set 141 INTEGER, INTENT(IN) :: natom 142 143 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho0_atom', & 144 routineP = moduleN//':'//routineN 145 146 INTEGER :: iat 147 148 IF (ASSOCIATED(rho0_set)) THEN 149 CALL deallocate_rho0_atom(rho0_set) 150 END IF 151 152 ALLOCATE (rho0_set(natom)) 153 154 DO iat = 1, natom 155 NULLIFY (rho0_set(iat)%rho0_rad_h) 156 NULLIFY (rho0_set(iat)%vrho0_rad_h) 157 ENDDO 158 159 END SUBROUTINE allocate_rho0_atom 160 161! ************************************************************************************************** 162!> \brief ... 163!> \param rho0_atom ... 164!> \param nr ... 165!> \param nchannels ... 166! ************************************************************************************************** 167 SUBROUTINE allocate_rho0_atom_rad(rho0_atom, nr, nchannels) 168 169 TYPE(rho0_atom_type), INTENT(OUT) :: rho0_atom 170 INTEGER, INTENT(IN) :: nr, nchannels 171 172 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho0_atom_rad', & 173 routineP = moduleN//':'//routineN 174 175 ALLOCATE (rho0_atom%rho0_rad_h) 176 177 NULLIFY (rho0_atom%rho0_rad_h%r_coef) 178 ALLOCATE (rho0_atom%rho0_rad_h%r_coef(1:nr, 1:nchannels)) 179 rho0_atom%rho0_rad_h%r_coef = 0.0_dp 180 181 ALLOCATE (rho0_atom%vrho0_rad_h) 182 183 NULLIFY (rho0_atom%vrho0_rad_h%r_coef) 184 ALLOCATE (rho0_atom%vrho0_rad_h%r_coef(1:nr, 1:nchannels)) 185 rho0_atom%vrho0_rad_h%r_coef = 0.0_dp 186 187 END SUBROUTINE allocate_rho0_atom_rad 188 189! ************************************************************************************************** 190!> \brief ... 191!> \param rho0 ... 192! ************************************************************************************************** 193 SUBROUTINE allocate_rho0_mpole(rho0) 194 195 TYPE(rho0_mpole_type), POINTER :: rho0 196 197 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_rho0_mpole', & 198 routineP = moduleN//':'//routineN 199 200 IF (ASSOCIATED(rho0)) THEN 201 CALL deallocate_rho0_mpole(rho0) 202 END IF 203 204 ALLOCATE (rho0) 205 206 NULLIFY (rho0%mp_rho) 207 NULLIFY (rho0%mp_gau) 208 NULLIFY (rho0%norm_g0l_h) 209 NULLIFY (rho0%lmax0_kind) 210 NULLIFY (rho0%rho0_s_rs) 211 NULLIFY (rho0%rho0_s_gs) 212 213 END SUBROUTINE allocate_rho0_mpole 214 215! ************************************************************************************************** 216!> \brief ... 217!> \param rho0_mpole ... 218!> \param grid_atom ... 219!> \param ik ... 220! ************************************************************************************************** 221 SUBROUTINE calculate_g0(rho0_mpole, grid_atom, ik) 222 223 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 224 TYPE(grid_atom_type), POINTER :: grid_atom 225 INTEGER, INTENT(IN) :: ik 226 227 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_g0', routineP = moduleN//':'//routineN 228 229 INTEGER :: ir, l, lmax, nr 230 REAL(dp) :: c1, prefactor, root_z_h, z_h 231 REAL(dp), ALLOCATABLE, DIMENSION(:) :: erf_z_h, gexp, gh_tmp, int1, int2 232 233 nr = grid_atom%nr 234 lmax = rho0_mpole%lmax0_kind(ik) 235 z_h = rho0_mpole%zet0_h 236 root_z_h = SQRT(z_h) 237 238! Allocate g0 239 CALL reallocate(rho0_mpole%mp_gau(ik)%g0_h, 1, nr, 0, lmax) 240 CALL reallocate(rho0_mpole%mp_gau(ik)%vg0_h, 1, nr, 0, lmax) 241 242 ALLOCATE (gexp(nr), gh_tmp(nr), erf_z_h(nr), int1(nr), int2(nr)) 243 244 gh_tmp(1:nr) = EXP(-z_h*grid_atom%rad2(1:nr)) 245 246 DO ir = 1, nr 247 erf_z_h(ir) = erf(grid_atom%rad(ir)*root_z_h) 248 END DO 249 250 DO ir = 1, nr 251 IF (gh_tmp(ir) < 1.0E-30_dp) gh_tmp(ir) = 0.0_dp 252 END DO 253 254 gexp(1:nr) = gh_tmp(1:nr) 255 rho0_mpole%mp_gau(ik)%g0_h(1:nr, 0) = gh_tmp(1:nr)* & 256 rho0_mpole%norm_g0l_h(0) 257 CALL whittaker_c0a(int1, grid_atom%rad, gh_tmp, erf_z_h, z_h, 0, 0, nr) 258 CALL whittaker_ci(int2, grid_atom%rad, gh_tmp, z_h, 0, nr) 259 260 prefactor = fourpi*rho0_mpole%norm_g0l_h(0) 261 262 c1 = SQRT(pi*pi*pi/(z_h*z_h*z_h))*rho0_mpole%norm_g0l_h(0) 263 264 DO ir = 1, nr 265 rho0_mpole%mp_gau(ik)%vg0_h(ir, 0) = c1*erf_z_h(ir)*grid_atom%oorad2l(ir, 1) 266 END DO 267 268 DO l = 1, lmax 269 gh_tmp(1:nr) = gh_tmp(1:nr)*grid_atom%rad(1:nr) 270 rho0_mpole%mp_gau(ik)%g0_h(1:nr, l) = gh_tmp(1:nr)* & 271 rho0_mpole%norm_g0l_h(l) 272 273 prefactor = fourpi/(2.0_dp*l + 1.0_dp)*rho0_mpole%norm_g0l_h(l) 274 CALL whittaker_c0a(int1, grid_atom%rad, gexp, erf_z_h, z_h, l, l, nr) 275 DO ir = 1, nr 276 rho0_mpole%mp_gau(ik)%vg0_h(ir, l) = prefactor*(int1(ir) + & 277 int2(ir)*grid_atom%rad2l(ir, l)) 278 END DO 279 280 END DO ! l 281 282 DEALLOCATE (gexp, erf_z_h, gh_tmp, int1, int2) 283 END SUBROUTINE calculate_g0 284 285! ************************************************************************************************** 286!> \brief ... 287!> \param mp_gau ... 288! ************************************************************************************************** 289 SUBROUTINE deallocate_mpole_gau(mp_gau) 290 291 TYPE(mpole_gau_overlap), DIMENSION(:), POINTER :: mp_gau 292 293 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_mpole_gau', & 294 routineP = moduleN//':'//routineN 295 296 INTEGER :: ikind, nkind 297 298 nkind = SIZE(mp_gau) 299 300 DO ikind = 1, nkind 301 302 IF (ASSOCIATED(mp_gau(ikind)%Qlm_gg)) THEN 303 DEALLOCATE (mp_gau(ikind)%Qlm_gg) 304 END IF 305 306 DEALLOCATE (mp_gau(ikind)%g0_h) 307 308 DEALLOCATE (mp_gau(ikind)%vg0_h) 309 END DO 310 311 DEALLOCATE (mp_gau) 312 313 END SUBROUTINE deallocate_mpole_gau 314 315! ************************************************************************************************** 316!> \brief ... 317!> \param mp_rho ... 318! ************************************************************************************************** 319 SUBROUTINE deallocate_mpole_rho(mp_rho) 320 321 TYPE(mpole_rho_atom), DIMENSION(:), POINTER :: mp_rho 322 323 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_mpole_rho', & 324 routineP = moduleN//':'//routineN 325 326 INTEGER :: iat, natom 327 328 natom = SIZE(mp_rho) 329 330 DO iat = 1, natom 331 DEALLOCATE (mp_rho(iat)%Qlm_h) 332 DEALLOCATE (mp_rho(iat)%Qlm_s) 333 DEALLOCATE (mp_rho(iat)%Qlm_tot) 334 DEALLOCATE (mp_rho(iat)%Qlm_car) 335 END DO 336 337 DEALLOCATE (mp_rho) 338 339 END SUBROUTINE deallocate_mpole_rho 340 341! ************************************************************************************************** 342!> \brief ... 343!> \param rho0_atom_set ... 344! ************************************************************************************************** 345 SUBROUTINE deallocate_rho0_atom(rho0_atom_set) 346 347 TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set 348 349 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_rho0_atom', & 350 routineP = moduleN//':'//routineN 351 352 INTEGER :: iat, natom 353 354 IF (ASSOCIATED(rho0_atom_set)) THEN 355 356 natom = SIZE(rho0_atom_set) 357 358 DO iat = 1, natom 359 IF (ASSOCIATED(rho0_atom_set(iat)%rho0_rad_h)) THEN 360 DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h%r_coef) 361 DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h) 362 ENDIF 363 IF (ASSOCIATED(rho0_atom_set(iat)%vrho0_rad_h)) THEN 364 DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h%r_coef) 365 DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h) 366 ENDIF 367 ENDDO 368 369 DEALLOCATE (rho0_atom_set) 370 ELSE 371 CALL cp_abort(__LOCATION__, & 372 "The pointer rho0_atom_set is not associated and "// & 373 "cannot be deallocated") 374 END IF 375 376 END SUBROUTINE deallocate_rho0_atom 377! ************************************************************************************************** 378!> \brief ... 379!> \param rho0 ... 380! ************************************************************************************************** 381 SUBROUTINE deallocate_rho0_mpole(rho0) 382 383 TYPE(rho0_mpole_type), POINTER :: rho0 384 385 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_rho0_mpole', & 386 routineP = moduleN//':'//routineN 387 388 IF (ASSOCIATED(rho0)) THEN 389 390 IF (ASSOCIATED(rho0%mp_gau)) CALL deallocate_mpole_gau(rho0%mp_gau) 391 392 IF (ASSOCIATED(rho0%mp_rho)) CALL deallocate_mpole_rho(rho0%mp_rho) 393 394 IF (ASSOCIATED(rho0%lmax0_kind)) THEN 395 DEALLOCATE (rho0%lmax0_kind) 396 END IF 397 398 IF (ASSOCIATED(rho0%norm_g0l_h)) THEN 399 DEALLOCATE (rho0%norm_g0l_h) 400 END IF 401 402 IF (ASSOCIATED(rho0%rho0_s_rs)) THEN 403 CALL pw_release(rho0%rho0_s_rs%pw) 404 DEALLOCATE (rho0%rho0_s_rs) 405 ENDIF 406 407 IF (ASSOCIATED(rho0%rho0_s_gs)) THEN 408 CALL pw_release(rho0%rho0_s_gs%pw) 409 DEALLOCATE (rho0%rho0_s_gs) 410 411 ENDIF 412 DEALLOCATE (rho0) 413 ELSE 414 CALL cp_abort(__LOCATION__, & 415 "The pointer rho0 is not associated and "// & 416 "cannot be deallocated") 417 END IF 418 419 END SUBROUTINE deallocate_rho0_mpole 420 421! ************************************************************************************************** 422!> \brief ... 423!> \param rho0_mpole ... 424!> \param g0_h ... 425!> \param vg0_h ... 426!> \param iat ... 427!> \param ikind ... 428!> \param lmax_0 ... 429!> \param l0_ikind ... 430!> \param mp_gau_ikind ... 431!> \param mp_rho ... 432!> \param norm_g0l_h ... 433!> \param Qlm_gg ... 434!> \param Qlm_car ... 435!> \param Qlm_tot ... 436!> \param zet0_h ... 437!> \param igrid_zet0_s ... 438!> \param rpgf0_h ... 439!> \param rpgf0_s ... 440!> \param max_rpgf0_s ... 441!> \param rho0_s_rs ... 442!> \param rho0_s_gs ... 443! ************************************************************************************************** 444 SUBROUTINE get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, & 445 mp_gau_ikind, mp_rho, norm_g0l_h, & 446 Qlm_gg, Qlm_car, Qlm_tot, & 447 zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, & 448 max_rpgf0_s, rho0_s_rs, rho0_s_gs) 449 450 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 451 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: g0_h, vg0_h 452 INTEGER, INTENT(IN), OPTIONAL :: iat, ikind 453 INTEGER, INTENT(OUT), OPTIONAL :: lmax_0, l0_ikind 454 TYPE(mpole_gau_overlap), OPTIONAL, POINTER :: mp_gau_ikind 455 TYPE(mpole_rho_atom), DIMENSION(:), OPTIONAL, & 456 POINTER :: mp_rho 457 REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: norm_g0l_h 458 REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg 459 REAL(dp), DIMENSION(:), OPTIONAL, POINTER :: Qlm_car, Qlm_tot 460 REAL(dp), INTENT(OUT), OPTIONAL :: zet0_h 461 INTEGER, INTENT(OUT), OPTIONAL :: igrid_zet0_s 462 REAL(dp), INTENT(OUT), OPTIONAL :: rpgf0_h, rpgf0_s, max_rpgf0_s 463 TYPE(pw_p_type), OPTIONAL, POINTER :: rho0_s_rs, rho0_s_gs 464 465 CHARACTER(len=*), PARAMETER :: routineN = 'get_rho0_mpole', routineP = moduleN//':'//routineN 466 467 IF (ASSOCIATED(rho0_mpole)) THEN 468 469 IF (PRESENT(lmax_0)) lmax_0 = rho0_mpole%lmax_0 470 IF (PRESENT(mp_rho)) mp_rho => rho0_mpole%mp_rho 471 IF (PRESENT(norm_g0l_h)) norm_g0l_h => rho0_mpole%norm_g0l_h 472 IF (PRESENT(zet0_h)) zet0_h = rho0_mpole%zet0_h 473 IF (PRESENT(igrid_zet0_s)) igrid_zet0_s = rho0_mpole%igrid_zet0_s 474 IF (PRESENT(max_rpgf0_s)) max_rpgf0_s = rho0_mpole%max_rpgf0_s 475 IF (PRESENT(rho0_s_rs)) rho0_s_rs => rho0_mpole%rho0_s_rs 476 IF (PRESENT(rho0_s_gs)) rho0_s_gs => rho0_mpole%rho0_s_gs 477 478 IF (PRESENT(ikind)) THEN 479 IF (PRESENT(l0_ikind)) l0_ikind = rho0_mpole%lmax0_kind(ikind) 480 IF (PRESENT(mp_gau_ikind)) mp_gau_ikind => rho0_mpole%mp_gau(ikind) 481 IF (PRESENT(g0_h)) g0_h => rho0_mpole%mp_gau(ikind)%g0_h 482 IF (PRESENT(vg0_h)) vg0_h => rho0_mpole%mp_gau(ikind)%vg0_h 483 IF (PRESENT(Qlm_gg)) Qlm_gg => rho0_mpole%mp_gau(ikind)%Qlm_gg 484 IF (PRESENT(rpgf0_h)) rpgf0_h = rho0_mpole%mp_gau(ikind)%rpgf0_h 485 IF (PRESENT(rpgf0_s)) rpgf0_s = rho0_mpole%mp_gau(ikind)%rpgf0_s 486 END IF 487 IF (PRESENT(iat)) THEN 488 IF (PRESENT(Qlm_car)) Qlm_car => rho0_mpole%mp_rho(iat)%Qlm_car 489 IF (PRESENT(Qlm_tot)) Qlm_tot => rho0_mpole%mp_rho(iat)%Qlm_tot 490 END IF 491 492 ELSE 493 CPABORT("The pointer rho0_mpole is not associated") 494 END IF 495 496 END SUBROUTINE get_rho0_mpole 497 498! ************************************************************************************************** 499!> \brief ... 500!> \param mp_rho ... 501!> \param nchan_s ... 502!> \param nchan_c ... 503!> \param zeff ... 504!> \param tddft ... 505! ************************************************************************************************** 506 SUBROUTINE initialize_mpole_rho(mp_rho, nchan_s, nchan_c, zeff, tddft) 507 508 TYPE(mpole_rho_atom) :: mp_rho 509 INTEGER, INTENT(IN) :: nchan_s, nchan_c 510 REAL(KIND=dp), INTENT(IN) :: zeff 511 LOGICAL, OPTIONAL :: tddft 512 513 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_mpole_rho', & 514 routineP = moduleN//':'//routineN 515 516 LOGICAL :: my_tddft 517 518 my_tddft = .FALSE. 519 IF (PRESENT(tddft)) my_tddft = tddft 520 521 CALL reallocate(mp_rho%Qlm_h, 1, nchan_s) 522 CALL reallocate(mp_rho%Qlm_s, 1, nchan_s) 523 CALL reallocate(mp_rho%Qlm_tot, 1, nchan_s) 524 CALL reallocate(mp_rho%Qlm_car, 1, nchan_c) 525 526 mp_rho%Qlm_h = 0.0_dp 527 mp_rho%Qlm_s = 0.0_dp 528 mp_rho%Qlm_tot = 0.0_dp 529 mp_rho%Qlm_car = 0.0_dp 530 IF (.NOT. my_tddft) THEN 531 mp_rho%Qlm_z = -2.0_dp*rootpi*Zeff 532 ELSE 533 mp_rho%Qlm_z = 0.0_dp 534 END IF 535 mp_rho%Q0 = 0.0_dp 536 537 END SUBROUTINE initialize_mpole_rho 538 539! ************************************************************************************************** 540!> \brief ... 541!> \param rho0_mpole ... 542!> \param unit_str ... 543!> \param output_unit ... 544! ************************************************************************************************** 545 SUBROUTINE write_rho0_info(rho0_mpole, unit_str, output_unit) 546 547 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 548 CHARACTER(LEN=*), INTENT(IN) :: unit_str 549 INTEGER, INTENT(in) :: output_unit 550 551 CHARACTER(len=*), PARAMETER :: routineN = 'write_rho0_info', & 552 routineP = moduleN//':'//routineN 553 554 INTEGER :: ikind, l, nkind 555 REAL(dp) :: conv 556 557 IF (ASSOCIATED(rho0_mpole)) THEN 558 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 559 560 WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & 561 "*** Compensation density charges data set ***" 562 WRITE (UNIT=output_unit, FMT="(T2,A,T35,f16.10)") & 563 "- Rho0 exponent :", rho0_mpole%zet0_h 564 WRITE (UNIT=output_unit, FMT="(T2,A,T35,I5)") & 565 "- Global max l :", rho0_mpole%lmax_0 566 567 WRITE (UNIT=output_unit, FMT="(T2,A)") & 568 "- Normalization constants for g0" 569 DO l = 0, rho0_mpole%lmax_0 570 WRITE (UNIT=output_unit, FMT="(T20,A,T31,I2,T38,A,f15.5)") & 571 "ang. mom.= ", l, " hard= ", rho0_mpole%norm_g0l_h(l) 572 END DO 573 574 nkind = SIZE(rho0_mpole%lmax0_kind, 1) 575 DO ikind = 1, nkind 576 WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,I2)") & 577 "- rho0 max L and radii in "//TRIM(unit_str)// & 578 " for the atom kind :", ikind 579 580 WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,I5)") & 581 "=> l max :", rho0_mpole%lmax0_kind(ikind) 582 583 WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,f20.10)") & 584 "=> max radius of g0: ", & 585 rho0_mpole%mp_gau(ikind)%rpgf0_h*conv 586 END DO ! ikind 587 588 ELSE 589 WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") & 590 ' WARNING: I cannot print rho0, it is not associated' 591 END IF 592 593 END SUBROUTINE write_rho0_info 594END MODULE qs_rho0_types 595