1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate Hirshfeld charges and related functions 8!> \par History 9!> 11.2014 created [JGH] 10!> \author JGH 11! ************************************************************************************************** 12MODULE hirshfeld_methods 13 USE atom_kind_orbitals, ONLY: calculate_atomic_density 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind 16 USE cell_types, ONLY: cell_type,& 17 pbc 18 USE cp_control_types, ONLY: dft_control_type 19 USE cp_para_types, ONLY: cp_para_env_type 20 USE cp_units, ONLY: cp_unit_to_cp2k 21 USE cube_utils, ONLY: cube_info_type 22 USE hirshfeld_types, ONLY: get_hirshfeld_info,& 23 hirshfeld_type,& 24 set_hirshfeld_info 25 USE input_constants, ONLY: radius_covalent,& 26 radius_default,& 27 radius_single,& 28 radius_user,& 29 radius_vdw,& 30 shape_function_density,& 31 shape_function_gaussian 32 USE kinds, ONLY: dp,& 33 int_8 34 USE mathconstants, ONLY: pi 35 USE message_passing, ONLY: mp_sum 36 USE particle_types, ONLY: particle_type 37 USE periodic_table, ONLY: get_ptable_info 38 USE pw_env_types, ONLY: pw_env_get,& 39 pw_env_type 40 USE pw_methods, ONLY: pw_integrate_function 41 USE pw_pool_types, ONLY: pw_pool_create_pw,& 42 pw_pool_give_back_pw,& 43 pw_pool_type 44 USE pw_types, ONLY: REALDATA3D,& 45 REALSPACE,& 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_integrate_potential_low, ONLY: integrate_pgf_product_rspace 52 USE qs_kind_types, ONLY: get_qs_kind,& 53 qs_kind_type 54 USE qs_modify_pab_block, ONLY: FUNC_AB 55 USE qs_rho_types, ONLY: qs_rho_get,& 56 qs_rho_type 57 USE realspace_grid_types, ONLY: pw2rs,& 58 realspace_grid_desc_type,& 59 realspace_grid_type,& 60 rs2pw,& 61 rs_grid_release,& 62 rs_grid_retain,& 63 rs_grid_zero,& 64 rs_pw_transfer 65#include "./base/base_uses.f90" 66 67 IMPLICIT NONE 68 PRIVATE 69 70 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hirshfeld_methods' 71 72 PUBLIC :: create_shape_function, comp_hirshfeld_charges, & 73 comp_hirshfeld_i_charges, write_hirshfeld_charges 74 75! ************************************************************************************************** 76 77CONTAINS 78 79! ************************************************************************************************** 80!> \brief ... 81!> \param charges ... 82!> \param hirshfeld_env ... 83!> \param particle_set ... 84!> \param qs_kind_set ... 85!> \param unit_nr ... 86! ************************************************************************************************** 87 SUBROUTINE write_hirshfeld_charges(charges, hirshfeld_env, particle_set, & 88 qs_kind_set, unit_nr) 89 REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges 90 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 91 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 92 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 93 INTEGER, INTENT(IN) :: unit_nr 94 95 CHARACTER(len=*), PARAMETER :: routineN = 'write_hirshfeld_charges', & 96 routineP = moduleN//':'//routineN 97 98 CHARACTER(len=2) :: element_symbol 99 INTEGER :: iatom, ikind, natom, nspin 100 REAL(KIND=dp) :: refc, tc1, zeff 101 102 natom = SIZE(charges, 1) 103 nspin = SIZE(charges, 2) 104 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!' 105 WRITE (UNIT=unit_nr, FMT="(T28,A)") "Hirshfeld Charges" 106 IF (nspin == 1) THEN 107 WRITE (UNIT=unit_nr, FMT="(/,T3,A,A)") & 108 "#Atom Element Kind ", " Ref Charge Population Net charge" 109 ELSE 110 WRITE (UNIT=unit_nr, FMT="(/,T3,A,A)") & 111 "#Atom Element Kind ", " Ref Charge Population Spin moment Net charge" 112 END IF 113 tc1 = 0.0_dp 114 DO iatom = 1, natom 115 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 116 element_symbol=element_symbol, kind_number=ikind) 117 refc = hirshfeld_env%charges(iatom) 118 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 119 IF (nspin == 1) THEN 120 WRITE (UNIT=unit_nr, FMT="(i7,T15,A2,T20,i3,T27,F8.3,T42,F8.3,T72,F8.3)") & 121 iatom, element_symbol, ikind, refc, charges(iatom, 1), zeff - charges(iatom, 1) 122 ELSE 123 WRITE (UNIT=unit_nr, FMT="(i7,T15,A2,T20,i3,T27,F8.3,T36,2F8.3,T61,F8.3,T72,F8.3)") & 124 iatom, element_symbol, ikind, refc, charges(iatom, 1), charges(iatom, 2), & 125 charges(iatom, 1) - charges(iatom, 2), zeff - SUM(charges(iatom, :)) 126 END IF 127 tc1 = tc1 + (zeff - SUM(charges(iatom, :))) 128 END DO 129 WRITE (UNIT=unit_nr, FMT="(/,T3,A,T72,F8.3)") "Total Charge ", tc1 130 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!' 131 132 END SUBROUTINE write_hirshfeld_charges 133 134! ************************************************************************************************** 135!> \brief creates kind specific shape functions for Hirshfeld charges 136!> \param hirshfeld_env the env that holds information about Hirshfeld 137!> \param qs_kind_set the qs_kind_set 138!> \param atomic_kind_set the atomic_kind_set 139!> \param radius optional radius parameter to use for all atomic kinds 140!> \param radii_list optional list of radii to use for different atomic kinds 141! ************************************************************************************************** 142 SUBROUTINE create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, radius, radii_list) 143 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 144 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 145 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 146 REAL(KIND=dp), OPTIONAL :: radius 147 REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii_list 148 149 CHARACTER(len=*), PARAMETER :: routineN = 'create_shape_function', & 150 routineP = moduleN//':'//routineN 151 INTEGER, PARAMETER :: ngto = 8 152 153 CHARACTER(len=2) :: esym 154 INTEGER :: ikind, nkind 155 LOGICAL :: found 156 REAL(KIND=dp) :: al, rco, zeff 157 REAL(KIND=dp), DIMENSION(ngto, 2) :: ppdens 158 TYPE(atomic_kind_type), POINTER :: atomic_kind 159 TYPE(qs_kind_type), POINTER :: qs_kind 160 161 CPASSERT(ASSOCIATED(hirshfeld_env)) 162 163 nkind = SIZE(qs_kind_set) 164 ALLOCATE (hirshfeld_env%kind_shape_fn(nkind)) 165 166 SELECT CASE (hirshfeld_env%shape_function_type) 167 CASE (shape_function_gaussian) 168 DO ikind = 1, nkind 169 hirshfeld_env%kind_shape_fn(ikind)%numexp = 1 170 ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(1)) 171 ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(1)) 172 CALL get_qs_kind(qs_kind_set(ikind), element_symbol=esym) 173 rco = 2.0_dp 174 SELECT CASE (hirshfeld_env%radius_type) 175 CASE (radius_default) 176 CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found) 177 rco = MAX(rco, 1.0_dp) 178 CASE (radius_user) 179 CPASSERT(PRESENT(radii_list)) 180 CPASSERT(ASSOCIATED(radii_list)) 181 CPASSERT(SIZE(radii_list) == nkind) 182 ! Note we assume that radii_list is correctly ordered 183 rco = radii_list(ikind) 184 CASE (radius_vdw) 185 CALL get_ptable_info(symbol=esym, vdw_radius=rco, found=found) 186 IF (.NOT. found) THEN 187 rco = MAX(rco, 1.0_dp) 188 ELSE 189 IF (hirshfeld_env%use_bohr) & 190 rco = cp_unit_to_cp2k(rco, "angstrom") 191 END IF 192 CASE (radius_covalent) 193 CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found) 194 IF (.NOT. found) THEN 195 rco = MAX(rco, 1.0_dp) 196 ELSE 197 IF (hirshfeld_env%use_bohr) & 198 rco = cp_unit_to_cp2k(rco, "angstrom") 199 END IF 200 CASE (radius_single) 201 CPASSERT(PRESENT(radius)) 202 rco = radius 203 END SELECT 204 al = 0.5_dp/rco**2 205 hirshfeld_env%kind_shape_fn(ikind)%zet(1) = al 206 hirshfeld_env%kind_shape_fn(ikind)%coef(1) = (al/pi)**1.5_dp 207 END DO 208 CASE (shape_function_density) 209 ! calculate atomic density 210 DO ikind = 1, nkind 211 atomic_kind => atomic_kind_set(ikind) 212 qs_kind => qs_kind_set(ikind) 213 CALL calculate_atomic_density(ppdens(:, :), atomic_kind, qs_kind, ngto, & 214 confine=.FALSE.) 215 hirshfeld_env%kind_shape_fn(ikind)%numexp = ngto 216 ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(ngto)) 217 ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(ngto)) 218 hirshfeld_env%kind_shape_fn(ikind)%zet(:) = ppdens(:, 1) 219 CALL get_qs_kind(qs_kind, zeff=zeff) 220 hirshfeld_env%kind_shape_fn(ikind)%coef(:) = ppdens(:, 2)/zeff 221 END DO 222 223 CASE DEFAULT 224 CPABORT("Unknown shape function") 225 END SELECT 226 227 END SUBROUTINE create_shape_function 228 229! ************************************************************************************************** 230!> \brief ... 231!> \param qs_env ... 232!> \param hirshfeld_env ... 233!> \param charges ... 234! ************************************************************************************************** 235 SUBROUTINE comp_hirshfeld_charges(qs_env, hirshfeld_env, charges) 236 TYPE(qs_environment_type), POINTER :: qs_env 237 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 238 REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges 239 240 CHARACTER(len=*), PARAMETER :: routineN = 'comp_hirshfeld_charges', & 241 routineP = moduleN//':'//routineN 242 243 INTEGER :: is 244 LOGICAL :: rho_r_valid 245 REAL(KIND=dp) :: tnfun 246 TYPE(pw_env_type), POINTER :: pw_env 247 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 248 TYPE(pw_p_type), POINTER :: rhonorm 249 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 250 TYPE(qs_rho_type), POINTER :: rho 251 252 NULLIFY (rho_r) 253 ! normalization function on grid 254 CALL calculate_hirshfeld_normalization(qs_env, hirshfeld_env) 255 ! check normalization 256 tnfun = pw_integrate_function(hirshfeld_env%fnorm%pw) 257 tnfun = ABS(tnfun - SUM(hirshfeld_env%charges)) 258 ! 259 ALLOCATE (rhonorm) 260 ! 261 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho) 262 CALL qs_rho_get(rho, rho_r=rho_r, rho_r_valid=rho_r_valid) 263 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool) 264 CALL pw_pool_create_pw(auxbas_pw_pool, rhonorm%pw, use_data=REALDATA3D) 265 ! loop over spins 266 DO is = 1, SIZE(rho_r) 267 IF (rho_r_valid) THEN 268 CALL hfun_scale(rhonorm%pw%cr3d, rho_r(is)%pw%cr3d, & 269 hirshfeld_env%fnorm%pw%cr3d) 270 ELSE 271 CPABORT("We need rho in real space") 272 END IF 273 CALL hirshfeld_integration(qs_env, hirshfeld_env, rhonorm, charges(:, is)) 274 charges(:, is) = charges(:, is)*hirshfeld_env%charges(:) 275 END DO 276 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhonorm%pw) 277 ! 278 DEALLOCATE (rhonorm) 279 280 END SUBROUTINE comp_hirshfeld_charges 281! ************************************************************************************************** 282!> \brief Calculate fout = fun1/fun2 283!> \param fout ... 284!> \param fun1 ... 285!> \param fun2 ... 286! ************************************************************************************************** 287 SUBROUTINE hfun_scale(fout, fun1, fun2) 288 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout 289 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2 290 291 CHARACTER(len=*), PARAMETER :: routineN = 'hfun_scale', routineP = moduleN//':'//routineN 292 REAL(KIND=dp), PARAMETER :: small = 1.0e-12_dp 293 294 INTEGER :: i1, i2, i3, n1, n2, n3 295 296 n1 = SIZE(fout, 1) 297 n2 = SIZE(fout, 2) 298 n3 = SIZE(fout, 3) 299 CPASSERT(n1 == SIZE(fun1, 1)) 300 CPASSERT(n2 == SIZE(fun1, 2)) 301 CPASSERT(n3 == SIZE(fun1, 3)) 302 CPASSERT(n1 == SIZE(fun2, 1)) 303 CPASSERT(n2 == SIZE(fun2, 2)) 304 CPASSERT(n3 == SIZE(fun2, 3)) 305 306 DO i3 = 1, n3 307 DO i2 = 1, n2 308 DO i1 = 1, n1 309 IF (fun2(i1, i2, i3) > small) THEN 310 fout(i1, i2, i3) = fun1(i1, i2, i3)/fun2(i1, i2, i3) 311 ELSE 312 fout(i1, i2, i3) = 0.0_dp 313 END IF 314 END DO 315 END DO 316 END DO 317 318 END SUBROUTINE hfun_scale 319 320! ************************************************************************************************** 321!> \brief ... 322!> \param qs_env ... 323!> \param hirshfeld_env ... 324!> \param charges ... 325!> \param ounit ... 326! ************************************************************************************************** 327 SUBROUTINE comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, ounit) 328 TYPE(qs_environment_type), POINTER :: qs_env 329 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 330 REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges 331 INTEGER, INTENT(IN) :: ounit 332 333 CHARACTER(len=*), PARAMETER :: routineN = 'comp_hirshfeld_i_charges', & 334 routineP = moduleN//':'//routineN 335 INTEGER, PARAMETER :: maxloop = 100 336 REAL(KIND=dp), PARAMETER :: maxres = 1.0e-2_dp 337 338 CHARACTER(len=3) :: yesno 339 INTEGER :: iat, iloop, is, natom 340 LOGICAL :: rho_r_valid 341 REAL(KIND=dp) :: res, tnfun 342 TYPE(pw_env_type), POINTER :: pw_env 343 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 344 TYPE(pw_p_type), POINTER :: rhonorm 345 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 346 TYPE(qs_rho_type), POINTER :: rho 347 348 NULLIFY (rho_r) 349 350 natom = SIZE(charges, 1) 351 352 IF (ounit > 0) WRITE (ounit, "(/,T2,A)") "Hirshfeld charge iterations: Residuals ..." 353 ! 354 ALLOCATE (rhonorm) 355 ! 356 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho) 357 CALL qs_rho_get(rho, rho_r=rho_r, rho_r_valid=rho_r_valid) 358 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool) 359 CALL pw_pool_create_pw(auxbas_pw_pool, rhonorm%pw, use_data=REALDATA3D) 360 ! 361 DO iloop = 1, maxloop 362 363 ! normalization function on grid 364 CALL calculate_hirshfeld_normalization(qs_env, hirshfeld_env) 365 ! check normalization 366 tnfun = pw_integrate_function(hirshfeld_env%fnorm%pw) 367 tnfun = ABS(tnfun - SUM(hirshfeld_env%charges)) 368 ! loop over spins 369 DO is = 1, SIZE(rho_r) 370 IF (rho_r_valid) THEN 371 CALL hfun_scale(rhonorm%pw%cr3d, rho_r(is)%pw%cr3d, & 372 hirshfeld_env%fnorm%pw%cr3d) 373 ELSE 374 CPABORT("We need rho in real space") 375 END IF 376 CALL hirshfeld_integration(qs_env, hirshfeld_env, rhonorm, charges(:, is)) 377 charges(:, is) = charges(:, is)*hirshfeld_env%charges(:) 378 END DO 379 ! residual 380 res = 0.0_dp 381 DO iat = 1, natom 382 res = res + (SUM(charges(iat, :)) - hirshfeld_env%charges(iat))**2 383 END DO 384 res = SQRT(res/REAL(natom, KIND=dp)) 385 IF (ounit > 0) THEN 386 yesno = "NO " 387 IF (MOD(iloop, 10) == 0) yesno = "YES" 388 WRITE (ounit, FMT="(F8.3)", ADVANCE=yesno) res 389 END IF 390 ! update 391 DO iat = 1, natom 392 hirshfeld_env%charges(iat) = SUM(charges(iat, :)) 393 END DO 394 IF (res < maxres) EXIT 395 396 END DO 397 ! 398 CALL pw_pool_give_back_pw(auxbas_pw_pool, rhonorm%pw) 399 ! 400 DEALLOCATE (rhonorm) 401 402 END SUBROUTINE comp_hirshfeld_i_charges 403 404! ************************************************************************************************** 405!> \brief ... 406!> \param qs_env ... 407!> \param hirshfeld_env ... 408! ************************************************************************************************** 409 SUBROUTINE calculate_hirshfeld_normalization(qs_env, hirshfeld_env) 410 411 TYPE(qs_environment_type), POINTER :: qs_env 412 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 413 414 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_hirshfeld_normalization', & 415 routineP = moduleN//':'//routineN 416 417 INTEGER :: atom_a, handle, iatom, iex, ikind, & 418 ithread, j, natom, npme, nthread, & 419 numexp 420 INTEGER(KIND=int_8) :: subpatch_pattern 421 INTEGER, DIMENSION(:), POINTER :: atom_list, cores 422 REAL(KIND=dp) :: alpha, coef, eps_rho_rspace 423 REAL(KIND=dp), DIMENSION(3) :: ra 424 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 425 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 426 TYPE(cell_type), POINTER :: cell 427 TYPE(cube_info_type) :: cube_info 428 TYPE(dft_control_type), POINTER :: dft_control 429 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 430 TYPE(pw_env_type), POINTER :: pw_env 431 TYPE(pw_p_type), POINTER :: fnorm 432 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 433 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 434 TYPE(realspace_grid_type), POINTER :: rs_rho 435 436 CALL timeset(routineN, handle) 437 438 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, & 439 dft_control=dft_control, particle_set=particle_set, pw_env=pw_env) 440 CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, auxbas_rs_grid=rs_rho, & 441 auxbas_pw_pool=auxbas_pw_pool) 442 cube_info = pw_env%cube_info(1) 443 ! be careful in parallel nsmax is choosen with multigrid in mind! 444 CALL rs_grid_retain(rs_rho) 445 CALL rs_grid_zero(rs_rho) 446 447 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 448 ALLOCATE (pab(1, 1)) 449 nthread = 1 450 ithread = 0 451 452 DO ikind = 1, SIZE(atomic_kind_set) 453 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp 454 IF (numexp <= 0) CYCLE 455 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list) 456 ALLOCATE (cores(natom)) 457 458 DO iex = 1, numexp 459 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex) 460 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex) 461 npme = 0 462 cores = 0 463 DO iatom = 1, natom 464 atom_a = atom_list(iatom) 465 ra(:) = pbc(particle_set(atom_a)%r, cell) 466 IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN 467 ! replicated realspace grid, split the atoms up between procs 468 IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN 469 npme = npme + 1 470 cores(npme) = iatom 471 ENDIF 472 ELSE 473 npme = npme + 1 474 cores(npme) = iatom 475 ENDIF 476 END DO 477 DO j = 1, npme 478 iatom = cores(j) 479 atom_a = atom_list(iatom) 480 pab(1, 1) = hirshfeld_env%charges(atom_a)*coef 481 ra(:) = pbc(particle_set(atom_a)%r, cell) 482 subpatch_pattern = 0 483 ! la_max==0 so set lmax_global to 0 484 CALL collocate_pgf_product_rspace(0, alpha, 0, 0, 0.0_dp, 0, ra, & 485 (/0.0_dp, 0.0_dp, 0.0_dp/), 0.0_dp, 1.0_dp, pab, 0, 0, rs_rho, & 486 cell, cube_info, eps_rho_rspace, ga_gb_function=FUNC_AB, & 487 ithread=ithread, use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern, & 488 lmax_global=0) 489 END DO 490 END DO 491 492 DEALLOCATE (cores) 493 END DO 494 DEALLOCATE (pab) 495 496 NULLIFY (fnorm) 497 CALL get_hirshfeld_info(hirshfeld_env, fnorm=fnorm) 498 IF (ASSOCIATED(fnorm)) THEN 499 CALL pw_release(fnorm%pw) 500 DEALLOCATE (fnorm) 501 ENDIF 502 ALLOCATE (fnorm) 503 CALL pw_pool_create_pw(auxbas_pw_pool, fnorm%pw, use_data=REALDATA3D) 504 fnorm%pw%in_space = REALSPACE 505 CALL set_hirshfeld_info(hirshfeld_env, fnorm=fnorm) 506 507 CALL rs_pw_transfer(rs_rho, fnorm%pw, rs2pw) 508 CALL rs_grid_release(rs_rho) 509 510 CALL timestop(handle) 511 512 END SUBROUTINE calculate_hirshfeld_normalization 513 514! ************************************************************************************************** 515!> \brief ... 516!> \param qs_env ... 517!> \param hirshfeld_env ... 518!> \param rfun ... 519!> \param fval ... 520!> \param fderiv ... 521! ************************************************************************************************** 522 SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv) 523 524 TYPE(qs_environment_type), POINTER :: qs_env 525 TYPE(hirshfeld_type), POINTER :: hirshfeld_env 526 TYPE(pw_p_type), POINTER :: rfun 527 REAL(KIND=dp), DIMENSION(:), INTENT(inout) :: fval 528 REAL(KIND=dp), DIMENSION(:, :), INTENT(inout), & 529 OPTIONAL :: fderiv 530 531 CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_integration', & 532 routineP = moduleN//':'//routineN 533 534 INTEGER :: atom_a, handle, iatom, iex, ikind, & 535 ithread, j, natom, npme, nthread, & 536 numexp 537 INTEGER, ALLOCATABLE, DIMENSION(:) :: cores 538 INTEGER, DIMENSION(:), POINTER :: atom_list 539 LOGICAL :: do_force 540 REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace 541 REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra 542 REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab 543 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 544 TYPE(cell_type), POINTER :: cell 545 TYPE(cp_para_env_type), POINTER :: para_env 546 TYPE(dft_control_type), POINTER :: dft_control 547 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 548 TYPE(pw_env_type), POINTER :: pw_env 549 TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc 550 TYPE(realspace_grid_type), POINTER :: rs_v 551 552 CALL timeset(routineN, handle) 553 554 do_force = PRESENT(fderiv) 555 fval = 0.0_dp 556 dvol = rfun%pw%pw_grid%dvol 557 558 NULLIFY (pw_env, auxbas_rs_desc) 559 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 560 CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, & 561 auxbas_rs_grid=rs_v) 562 CALL rs_grid_retain(rs_v) 563 CALL rs_pw_transfer(rs_v, rfun%pw, pw2rs) 564 565 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, & 566 dft_control=dft_control, particle_set=particle_set) 567 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace 568 569 nthread = 1 570 ithread = 0 571 ALLOCATE (hab(1, 1), pab(1, 1)) 572 573 DO ikind = 1, SIZE(atomic_kind_set) 574 numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp 575 IF (numexp <= 0) CYCLE 576 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list) 577 ALLOCATE (cores(natom)) 578 579 DO iex = 1, numexp 580 alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex) 581 coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex) 582 npme = 0 583 cores = 0 584 DO iatom = 1, natom 585 atom_a = atom_list(iatom) 586 ra(:) = pbc(particle_set(atom_a)%r, cell) 587 IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN 588 ! replicated realspace grid, split the atoms up between procs 589 IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN 590 npme = npme + 1 591 cores(npme) = iatom 592 ENDIF 593 ELSE 594 npme = npme + 1 595 cores(npme) = iatom 596 ENDIF 597 END DO 598 599 DO j = 1, npme 600 iatom = cores(j) 601 atom_a = atom_list(iatom) 602 ra(:) = pbc(particle_set(atom_a)%r, cell) 603 pab(1, 1) = coef 604 hab(1, 1) = 0.0_dp 605 force_a(:) = 0.0_dp 606 force_b(:) = 0.0_dp 607 ! 608 CALL integrate_pgf_product_rspace(0, alpha, 0, & 609 0, 0.0_dp, 0, ra, (/0.0_dp, 0.0_dp, 0.0_dp/), 0.0_dp, & 610 rs_v, cell, pw_env%cube_info(1), hab, pab=pab, o1=0, o2=0, & 611 eps_gvg_rspace=eps_rho_rspace, calculate_forces=do_force, & 612 force_a=force_a, force_b=force_b, use_virial=.FALSE., & 613 use_subpatch=.TRUE., subpatch_pattern=0_int_8) 614 fval(atom_a) = fval(atom_a) + hab(1, 1)*dvol*coef 615 IF (do_force) THEN 616 fderiv(:, atom_a) = fderiv(:, atom_a) + force_a(:)*dvol 617 END IF 618 END DO 619 620 END DO 621 DEALLOCATE (cores) 622 623 END DO 624 625 CALL rs_grid_release(rs_v) 626 DEALLOCATE (hab, pab) 627 628 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 629 CALL mp_sum(fval, para_env%group) 630 631 CALL timestop(handle) 632 633 END SUBROUTINE hirshfeld_integration 634 635END MODULE hirshfeld_methods 636