1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5! ************************************************************************************************** 6MODULE qs_grid_atom 7 8 USE input_constants, ONLY: do_gapw_gcs,& 9 do_gapw_gct,& 10 do_gapw_log 11 USE kinds, ONLY: dp 12 USE lebedev, ONLY: get_number_of_lebedev_grid,& 13 lebedev_grid 14 USE mathconstants, ONLY: pi 15 USE memory_utilities, ONLY: reallocate 16#include "./base/base_uses.f90" 17 18 IMPLICIT NONE 19 20 PRIVATE 21 22 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_grid_atom' 23 24 TYPE grid_batch_type 25 INTEGER :: np 26 REAL(KIND=dp), DIMENSION(3) :: rcenter 27 REAL(KIND=dp) :: rad 28 REAL(dp), DIMENSION(:, :), ALLOCATABLE :: rco 29 REAL(dp), DIMENSION(:), ALLOCATABLE :: weight 30 END TYPE grid_batch_type 31 32 TYPE atom_integration_grid_type 33 INTEGER :: nr, na 34 INTEGER :: np, ntot 35 INTEGER :: lebedev_grid 36 REAL(dp), DIMENSION(:), ALLOCATABLE :: rr 37 REAL(dp), DIMENSION(:), ALLOCATABLE :: wr, wa 38 INTEGER :: nbatch 39 TYPE(grid_batch_type), DIMENSION(:), ALLOCATABLE :: batch 40 END TYPE atom_integration_grid_type 41 42 TYPE grid_atom_type 43 INTEGER :: quadrature 44 INTEGER :: nr, ng_sphere 45 REAL(dp), DIMENSION(:), POINTER :: rad, rad2, & 46 wr, wa, & 47 azi, cos_azi, sin_azi, & 48 pol, cos_pol, sin_pol, usin_azi 49 REAL(dp), DIMENSION(:, :), & 50 POINTER :: rad2l, oorad2l, weight 51 END TYPE grid_atom_type 52 53 PUBLIC :: allocate_grid_atom, create_grid_atom, deallocate_grid_atom 54 PUBLIC :: grid_atom_type 55 PUBLIC :: initialize_atomic_grid 56 PUBLIC :: atom_integration_grid_type, deallocate_atom_int_grid 57 58! ************************************************************************************************** 59 60CONTAINS 61 62! ************************************************************************************************** 63!> \brief Initialize components of the grid_atom_type structure 64!> \param grid_atom ... 65!> \date 03.11.2000 66!> \author MK 67!> \author Matthias Krack (MK) 68!> \version 1.0 69! ************************************************************************************************** 70 SUBROUTINE allocate_grid_atom(grid_atom) 71 72 TYPE(grid_atom_type), POINTER :: grid_atom 73 74 CHARACTER(len=*), PARAMETER :: routineN = 'allocate_grid_atom', & 75 routineP = moduleN//':'//routineN 76 77 IF (ASSOCIATED(grid_atom)) CALL deallocate_grid_atom(grid_atom) 78 79 ALLOCATE (grid_atom) 80 81 NULLIFY (grid_atom%rad) 82 NULLIFY (grid_atom%rad2) 83 NULLIFY (grid_atom%wr) 84 NULLIFY (grid_atom%wa) 85 NULLIFY (grid_atom%weight) 86 NULLIFY (grid_atom%azi) 87 NULLIFY (grid_atom%cos_azi) 88 NULLIFY (grid_atom%sin_azi) 89 NULLIFY (grid_atom%pol) 90 NULLIFY (grid_atom%cos_pol) 91 NULLIFY (grid_atom%sin_pol) 92 NULLIFY (grid_atom%usin_azi) 93 NULLIFY (grid_atom%rad2l) 94 NULLIFY (grid_atom%oorad2l) 95 96 END SUBROUTINE allocate_grid_atom 97 98! ************************************************************************************************** 99!> \brief Deallocate a Gaussian-type orbital (GTO) basis set data set. 100!> \param grid_atom ... 101!> \date 03.11.2000 102!> \author MK 103!> \version 1.0 104! ************************************************************************************************** 105 SUBROUTINE deallocate_grid_atom(grid_atom) 106 TYPE(grid_atom_type), POINTER :: grid_atom 107 108 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_grid_atom', & 109 routineP = moduleN//':'//routineN 110 111 IF (ASSOCIATED(grid_atom)) THEN 112 113 IF (ASSOCIATED(grid_atom%rad)) THEN 114 DEALLOCATE (grid_atom%rad) 115 END IF 116 117 IF (ASSOCIATED(grid_atom%rad2)) THEN 118 DEALLOCATE (grid_atom%rad2) 119 END IF 120 121 IF (ASSOCIATED(grid_atom%wr)) THEN 122 DEALLOCATE (grid_atom%wr) 123 END IF 124 125 IF (ASSOCIATED(grid_atom%wa)) THEN 126 DEALLOCATE (grid_atom%wa) 127 END IF 128 129 IF (ASSOCIATED(grid_atom%weight)) THEN 130 DEALLOCATE (grid_atom%weight) 131 END IF 132 133 IF (ASSOCIATED(grid_atom%azi)) THEN 134 DEALLOCATE (grid_atom%azi) 135 END IF 136 137 IF (ASSOCIATED(grid_atom%cos_azi)) THEN 138 DEALLOCATE (grid_atom%cos_azi) 139 END IF 140 141 IF (ASSOCIATED(grid_atom%sin_azi)) THEN 142 DEALLOCATE (grid_atom%sin_azi) 143 END IF 144 145 IF (ASSOCIATED(grid_atom%pol)) THEN 146 DEALLOCATE (grid_atom%pol) 147 END IF 148 149 IF (ASSOCIATED(grid_atom%cos_pol)) THEN 150 DEALLOCATE (grid_atom%cos_pol) 151 END IF 152 153 IF (ASSOCIATED(grid_atom%sin_pol)) THEN 154 DEALLOCATE (grid_atom%sin_pol) 155 END IF 156 157 IF (ASSOCIATED(grid_atom%usin_azi)) THEN 158 DEALLOCATE (grid_atom%usin_azi) 159 END IF 160 161 IF (ASSOCIATED(grid_atom%rad2l)) THEN 162 DEALLOCATE (grid_atom%rad2l) 163 END IF 164 165 IF (ASSOCIATED(grid_atom%oorad2l)) THEN 166 DEALLOCATE (grid_atom%oorad2l) 167 END IF 168 169 DEALLOCATE (grid_atom) 170 ELSE 171 CALL cp_abort(__LOCATION__, & 172 "The pointer grid_atom is not associated and "// & 173 "cannot be deallocated") 174 END IF 175 END SUBROUTINE deallocate_grid_atom 176 177! ************************************************************************************************** 178!> \brief ... 179!> \param grid_atom ... 180!> \param nr ... 181!> \param na ... 182!> \param llmax ... 183!> \param ll ... 184!> \param quadrature ... 185! ************************************************************************************************** 186 SUBROUTINE create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature) 187 188 TYPE(grid_atom_type), POINTER :: grid_atom 189 INTEGER, INTENT(IN) :: nr, na, llmax, ll, quadrature 190 191 CHARACTER(len=*), PARAMETER :: routineN = 'create_grid_atom', & 192 routineP = moduleN//':'//routineN 193 194 INTEGER :: handle, ia, ir, l 195 REAL(dp) :: cosia, pol 196 REAL(dp), DIMENSION(:), POINTER :: rad, rad2, wr 197 198 CALL timeset(routineN, handle) 199 200 NULLIFY (rad, rad2, wr) 201 202 IF (ASSOCIATED(grid_atom)) THEN 203 204 ! Allocate the radial grid arrays 205 CALL reallocate(grid_atom%rad, 1, nr) 206 CALL reallocate(grid_atom%rad2, 1, nr) 207 CALL reallocate(grid_atom%wr, 1, nr) 208 CALL reallocate(grid_atom%wa, 1, na) 209 CALL reallocate(grid_atom%weight, 1, na, 1, nr) 210 CALL reallocate(grid_atom%azi, 1, na) 211 CALL reallocate(grid_atom%cos_azi, 1, na) 212 CALL reallocate(grid_atom%sin_azi, 1, na) 213 CALL reallocate(grid_atom%pol, 1, na) 214 CALL reallocate(grid_atom%cos_pol, 1, na) 215 CALL reallocate(grid_atom%sin_pol, 1, na) 216 CALL reallocate(grid_atom%usin_azi, 1, na) 217 CALL reallocate(grid_atom%rad2l, 1, nr, 0, llmax + 1) 218 CALL reallocate(grid_atom%oorad2l, 1, nr, 0, llmax + 1) 219 220 ! Calculate the radial grid for this kind 221 rad => grid_atom%rad 222 rad2 => grid_atom%rad2 223 wr => grid_atom%wr 224 225 grid_atom%quadrature = quadrature 226 CALL radial_grid(nr, rad, rad2, wr, quadrature) 227 228 grid_atom%rad2l(:, 0) = 1._dp 229 grid_atom%oorad2l(:, 0) = 1._dp 230 DO l = 1, llmax + 1 231 grid_atom%rad2l(:, l) = grid_atom%rad2l(:, l - 1)*rad(:) 232 grid_atom%oorad2l(:, l) = grid_atom%oorad2l(:, l - 1)/rad(:) 233 ENDDO 234 235 IF (ll > 0) THEN 236 grid_atom%wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:na) 237 DO ir = 1, nr 238 DO ia = 1, na 239 grid_atom%weight(ia, ir) = grid_atom%wr(ir)*grid_atom%wa(ia) 240 END DO 241 END DO 242 243 DO ia = 1, na 244 ! polar angle: pol = acos(r(3)) 245 cosia = lebedev_grid(ll)%r(3, ia) 246 grid_atom%cos_pol(ia) = cosia 247 ! azimuthal angle: pol = atan(r(2)/r(1)) 248 IF (ABS(lebedev_grid(ll)%r(2, ia)) < EPSILON(1.0_dp) .AND. & 249 ABS(lebedev_grid(ll)%r(1, ia)) < EPSILON(1.0_dp)) THEN 250 grid_atom%azi(ia) = 0.0_dp 251 ELSE 252 grid_atom%azi(ia) = ATAN2(lebedev_grid(ll)%r(2, ia), lebedev_grid(ll)%r(1, ia)) 253 END IF 254 grid_atom%cos_azi(ia) = COS(grid_atom%azi(ia)) 255 pol = ACOS(cosia) 256 IF (grid_atom%sin_pol(ia) < 0.0_dp) pol = -pol 257 grid_atom%pol(ia) = pol 258 grid_atom%sin_pol(ia) = SIN(grid_atom%pol(ia)) 259 260 grid_atom%sin_azi(ia) = SIN(grid_atom%azi(ia)) 261 IF (ABS(grid_atom%sin_azi(ia)) > EPSILON(1.0_dp)) THEN 262 grid_atom%usin_azi(ia) = 1.0_dp/grid_atom%sin_azi(ia) 263 ELSE 264 grid_atom%usin_azi(ia) = 1.0_dp 265 END IF 266 267 ENDDO 268 269 END IF 270 271 ELSE 272 CPABORT("The pointer grid_atom is not associated") 273 END IF 274 275 CALL timestop(handle) 276 277 END SUBROUTINE create_grid_atom 278 279! ************************************************************************************************** 280!> \brief Initialize atomic grid 281!> \param int_grid ... 282!> \param nr ... 283!> \param na ... 284!> \param rmax ... 285!> \param quadrature ... 286!> \param iunit ... 287!> \date 02.2018 288!> \author JGH 289!> \version 1.0 290! ************************************************************************************************** 291 SUBROUTINE initialize_atomic_grid(int_grid, nr, na, rmax, quadrature, iunit) 292 TYPE(atom_integration_grid_type), POINTER :: int_grid 293 INTEGER, INTENT(IN) :: nr, na 294 REAL(KIND=dp), INTENT(IN) :: rmax 295 INTEGER, INTENT(IN), OPTIONAL :: quadrature, iunit 296 297 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_atomic_grid', & 298 routineP = moduleN//':'//routineN 299 300 INTEGER :: ia, ig, ir, ix, iy, iz, la, ll, my_quad, & 301 n1, n2, n3, nbatch, ng, no, np, ntot, & 302 nu, nx 303 INTEGER, ALLOCATABLE, DIMENSION(:) :: icell 304 REAL(KIND=dp) :: ag, dd, dmax, r1, r2, r3 305 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rad, rad2, wa, wc, wr 306 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rang, rco 307 REAL(KIND=dp), DIMENSION(10) :: dco 308 REAL(KIND=dp), DIMENSION(3) :: rm 309 TYPE(atom_integration_grid_type), POINTER :: igr 310 311 ALLOCATE (igr) 312 313 ! type of quadrature grid 314 IF (PRESENT(quadrature)) THEN 315 my_quad = quadrature 316 ELSE 317 my_quad = do_gapw_log 318 END IF 319 320 ! radial grid 321 CPASSERT(nr > 1) 322 ALLOCATE (rad(nr), rad2(nr), wr(nr)) 323 CALL radial_grid(nr, rad, rad2, wr, my_quad) 324 ! 325 igr%nr = nr 326 ALLOCATE (igr%rr(nr)) 327 ALLOCATE (igr%wr(nr)) 328 ! store grid points always in ascending order 329 IF (rad(1) > rad(nr)) THEN 330 DO ir = nr, 1, -1 331 igr%rr(nr - ir + 1) = rad(ir) 332 igr%wr(nr - ir + 1) = wr(ir) 333 END DO 334 ELSE 335 igr%rr(1:nr) = rad(1:nr) 336 igr%wr(1:nr) = wr(1:nr) 337 END IF 338 ! only include grid points smaller than rmax 339 np = 0 340 DO ir = 1, nr 341 IF (igr%rr(ir) < rmax) THEN 342 np = np + 1 343 rad(np) = igr%rr(ir) 344 wr(np) = igr%wr(ir) 345 END IF 346 END DO 347 igr%np = np 348 ! 349 ! angular grid 350 CPASSERT(na > 1) 351 ll = get_number_of_lebedev_grid(n=na) 352 np = lebedev_grid(ll)%n 353 la = lebedev_grid(ll)%l 354 ALLOCATE (rang(3, np), wa(np)) 355 wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:np) 356 rang(1:3, 1:np) = lebedev_grid(ll)%r(1:3, 1:np) 357 igr%lebedev_grid = ll 358 ALLOCATE (igr%wa(np)) 359 igr%na = np 360 igr%wa(1:np) = wa(1:np) 361 ! 362 ! total grid points 363 ntot = igr%na*igr%np 364 igr%ntot = ntot 365 ALLOCATE (rco(3, ntot), wc(ntot)) 366 ig = 0 367 DO ir = 1, igr%np 368 DO ia = 1, igr%na 369 ig = ig + 1 370 rco(1:3, ig) = rang(1:3, ia)*rad(ir) 371 wc(ig) = wa(ia)*wr(ir) 372 END DO 373 END DO 374 ! grid for batches, odd number of cells 375 ng = NINT((REAL(ntot, dp)/32._dp)**0.33333_dp) 376 ng = ng + MOD(ng + 1, 2) 377 ! avarage number of points along radial grid 378 dco = 0.0_dp 379 ag = REAL(igr%np, dp)/ng 380 CPASSERT(SIZE(dco) >= (ng + 1)/2) 381 DO ig = 1, ng, 2 382 ir = MIN(NINT(ag)*ig, igr%np) 383 ia = (ig + 1)/2 384 dco(ia) = rad(ir) 385 END DO 386 ! batches 387 ALLOCATE (icell(ntot)) 388 icell = 0 389 nx = (ng - 1)/2 390 DO ig = 1, ntot 391 ix = grid_coord(rco(1, ig), dco, nx + 1) + nx 392 iy = grid_coord(rco(2, ig), dco, nx + 1) + nx 393 iz = grid_coord(rco(3, ig), dco, nx + 1) + nx 394 icell(ig) = iz*ng*ng + iy*ng + ix + 1 395 END DO 396 ! 397 igr%nbatch = ng*ng*ng 398 ALLOCATE (igr%batch(igr%nbatch)) 399 igr%batch(:)%np = 0 400 DO ig = 1, ntot 401 ia = icell(ig) 402 igr%batch(ia)%np = igr%batch(ia)%np + 1 403 END DO 404 DO ig = 1, igr%nbatch 405 np = igr%batch(ig)%np 406 ALLOCATE (igr%batch(ig)%rco(3, np), igr%batch(ig)%weight(np)) 407 igr%batch(ig)%np = 0 408 END DO 409 DO ig = 1, ntot 410 ia = icell(ig) 411 igr%batch(ia)%np = igr%batch(ia)%np + 1 412 np = igr%batch(ia)%np 413 igr%batch(ia)%rco(1:3, np) = rco(1:3, ig) 414 igr%batch(ia)%weight(np) = wc(ig) 415 END DO 416 ! 417 DEALLOCATE (rad, rad2, rang, wr, wa) 418 DEALLOCATE (rco, wc, icell) 419 ! 420 IF (ASSOCIATED(int_grid)) CALL deallocate_atom_int_grid(int_grid) 421 ALLOCATE (int_grid) 422 ALLOCATE (int_grid%rr(igr%nr), int_grid%wr(igr%nr), int_grid%wa(igr%na)) 423 int_grid%nr = igr%nr 424 int_grid%na = igr%na 425 int_grid%np = igr%np 426 int_grid%ntot = igr%ntot 427 int_grid%lebedev_grid = igr%lebedev_grid 428 int_grid%rr(:) = igr%rr(:) 429 int_grid%wr(:) = igr%wr(:) 430 int_grid%wa(:) = igr%wa(:) 431 ! 432 ! count batches 433 nbatch = 0 434 DO ig = 1, igr%nbatch 435 IF (igr%batch(ig)%np == 0) THEN 436 ! empty batch 437 ELSE IF (igr%batch(ig)%np <= 48) THEN 438 ! single batch 439 nbatch = nbatch + 1 440 ELSE 441 ! multiple batches 442 nbatch = nbatch + NINT(igr%batch(ig)%np/32._dp) 443 END IF 444 END DO 445 int_grid%nbatch = nbatch 446 ALLOCATE (int_grid%batch(nbatch)) 447 ! fill batches 448 n1 = 0 449 DO ig = 1, igr%nbatch 450 IF (igr%batch(ig)%np == 0) THEN 451 ! empty batch 452 ELSE IF (igr%batch(ig)%np <= 48) THEN 453 ! single batch 454 n1 = n1 + 1 455 np = igr%batch(ig)%np 456 ALLOCATE (int_grid%batch(n1)%rco(3, np), int_grid%batch(n1)%weight(np)) 457 int_grid%batch(n1)%np = np 458 int_grid%batch(n1)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, 1:np) 459 int_grid%batch(n1)%weight(1:np) = igr%batch(ig)%weight(1:np) 460 ELSE 461 ! multiple batches 462 n2 = NINT(igr%batch(ig)%np/32._dp) 463 n3 = igr%batch(ig)%np/n2 464 DO ia = n1 + 1, n1 + n2 465 nu = (ia - n1 - 1)*n3 + 1 466 no = nu + n3 - 1 467 IF (ia == n1 + n2) no = igr%batch(ig)%np 468 np = no - nu + 1 469 ALLOCATE (int_grid%batch(ia)%rco(3, np), int_grid%batch(ia)%weight(np)) 470 int_grid%batch(ia)%np = np 471 int_grid%batch(ia)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, nu:no) 472 int_grid%batch(ia)%weight(1:np) = igr%batch(ig)%weight(nu:no) 473 END DO 474 n1 = n1 + n2 475 END IF 476 END DO 477 CPASSERT(nbatch == n1) 478 ! batch center and radius 479 DO ig = 1, int_grid%nbatch 480 np = int_grid%batch(ig)%np 481 IF (np > 0) THEN 482 rm(1) = SUM(int_grid%batch(ig)%rco(1, 1:np)) 483 rm(2) = SUM(int_grid%batch(ig)%rco(2, 1:np)) 484 rm(3) = SUM(int_grid%batch(ig)%rco(3, 1:np)) 485 rm(1:3) = rm(1:3)/REAL(np, KIND=dp) 486 END IF 487 int_grid%batch(ig)%rcenter(1:3) = rm(1:3) 488 dmax = 0.0_dp 489 DO ia = 1, np 490 dd = SUM((int_grid%batch(ig)%rco(1:3, ia) - rm(1:3))**2) 491 dmax = MAX(dd, dmax) 492 END DO 493 int_grid%batch(ig)%rad = SQRT(dmax) 494 END DO 495 ! 496 CALL deallocate_atom_int_grid(igr) 497 ! 498 IF (PRESENT(iunit)) THEN 499 IF (iunit > 0) THEN 500 WRITE (iunit, "(/,A)") " Atomic Integration Grid Information" 501 WRITE (iunit, "(A,T51,3I10)") " Number of grid points [radial,angular,total]", & 502 int_grid%np, int_grid%na, int_grid%ntot 503 WRITE (iunit, "(A,T71,I10)") " Lebedev grid number", int_grid%lebedev_grid 504 WRITE (iunit, "(A,T61,F20.5)") " Maximum of radial grid [Bohr]", & 505 int_grid%rr(int_grid%np) 506 nbatch = int_grid%nbatch 507 WRITE (iunit, "(A,T71,I10)") " Total number of gridpoint batches", nbatch 508 n1 = int_grid%ntot 509 n2 = 0 510 n3 = NINT(REAL(int_grid%ntot, dp)/REAL(nbatch, dp)) 511 DO ig = 1, nbatch 512 n1 = MIN(n1, int_grid%batch(ig)%np) 513 n2 = MAX(n2, int_grid%batch(ig)%np) 514 END DO 515 WRITE (iunit, "(A,T51,3I10)") " Number of grid points/batch [min,max,ave]", n1, n2, n3 516 r1 = 1000._dp 517 r2 = 0.0_dp 518 r3 = 0.0_dp 519 DO ig = 1, int_grid%nbatch 520 r1 = MIN(r1, int_grid%batch(ig)%rad) 521 r2 = MAX(r2, int_grid%batch(ig)%rad) 522 r3 = r3 + int_grid%batch(ig)%rad 523 END DO 524 r3 = r3/REAL(ng*ng*ng, KIND=dp) 525 WRITE (iunit, "(A,T51,3F10.2)") " Batch radius (bohr) [min,max,ave]", r1, r2, r3 526 END IF 527 END IF 528 529 END SUBROUTINE initialize_atomic_grid 530 531! ************************************************************************************************** 532!> \brief ... 533!> \param x ... 534!> \param dco ... 535!> \param ng ... 536!> \return ... 537!> \retval igrid ... 538! ************************************************************************************************** 539 FUNCTION grid_coord(x, dco, ng) RESULT(igrid) 540 REAL(KIND=dp), INTENT(IN) :: x 541 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: dco 542 INTEGER, INTENT(IN) :: ng 543 INTEGER :: igrid 544 545 INTEGER :: ig 546 REAL(KIND=dp) :: xval 547 548 xval = ABS(x) 549 igrid = ng 550 DO ig = 1, ng 551 IF (xval <= dco(ig)) THEN 552 igrid = ig - 1 553 EXIT 554 END IF 555 END DO 556 IF (x < 0.0_dp) igrid = -igrid 557 CPASSERT(ABS(igrid) < ng) 558 END FUNCTION grid_coord 559 560! ************************************************************************************************** 561!> \brief ... 562!> \param int_grid ... 563! ************************************************************************************************** 564 SUBROUTINE deallocate_atom_int_grid(int_grid) 565 TYPE(atom_integration_grid_type), POINTER :: int_grid 566 567 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_atom_int_grid', & 568 routineP = moduleN//':'//routineN 569 570 INTEGER :: ib 571 572 IF (ASSOCIATED(int_grid)) THEN 573 IF (ALLOCATED(int_grid%rr)) DEALLOCATE (int_grid%rr) 574 IF (ALLOCATED(int_grid%wr)) DEALLOCATE (int_grid%wr) 575 IF (ALLOCATED(int_grid%wa)) DEALLOCATE (int_grid%wa) 576 ! batch 577 IF (ALLOCATED(int_grid%batch)) THEN 578 DO ib = 1, SIZE(int_grid%batch) 579 IF (ALLOCATED(int_grid%batch(ib)%rco)) DEALLOCATE (int_grid%batch(ib)%rco) 580 IF (ALLOCATED(int_grid%batch(ib)%weight)) DEALLOCATE (int_grid%batch(ib)%weight) 581 END DO 582 DEALLOCATE (int_grid%batch) 583 END IF 584 ! 585 DEALLOCATE (int_grid) 586 NULLIFY (int_grid) 587 END IF 588 589 END SUBROUTINE deallocate_atom_int_grid 590! ************************************************************************************************** 591!> \brief Generate a radial grid with n points by a quadrature rule. 592!> \param n ... 593!> \param r ... 594!> \param r2 ... 595!> \param wr ... 596!> \param radial_quadrature ... 597!> \date 20.09.1999 598!> \par Literature 599!> - A. D. Becke, J. Chem. Phys. 88, 2547 (1988) 600!> - J. M. Perez-Jorda, A. D. Becke and E. San-Fabian, 601!> J. Chem. Phys. 100, 6520 (1994) 602!> - M. Krack and A. M. Koester, J. Chem. Phys. 108, 3226 (1998) 603!> \author Matthias Krack 604!> \version 1.0 605! ************************************************************************************************** 606 SUBROUTINE radial_grid(n, r, r2, wr, radial_quadrature) 607 608 INTEGER, INTENT(IN) :: n 609 REAL(dp), DIMENSION(:), INTENT(INOUT) :: r, r2, wr 610 INTEGER, INTENT(IN) :: radial_quadrature 611 612 CHARACTER(len=*), PARAMETER :: routineN = 'radial_grid', routineP = moduleN//':'//routineN 613 614 INTEGER :: i 615 REAL(dp) :: cost, f, sint, sint2, t, w, x 616 617 f = pi/REAL(n + 1, dp) 618 619 SELECT CASE (radial_quadrature) 620 621 CASE (do_gapw_gcs) 622 623 ! *** Gauss-Chebyshev quadrature formula of the second kind *** 624 ! *** u [-1,+1] -> r [0,infinity] => r = (1 + u)/(1 - u) *** 625 626 DO i = 1, n 627 t = REAL(i, dp)*f 628 x = COS(t) 629 w = f*SIN(t)**2 630 r(i) = (1.0_dp + x)/(1.0_dp - x) 631 r2(i) = r(i)**2 632 wr(i) = w/SQRT(1.0_dp - x**2) 633 wr(i) = 2.0_dp*wr(i)*r2(i)/(1.0_dp - x)**2 634 END DO 635 636 CASE (do_gapw_gct) 637 638 ! *** Transformed Gauss-Chebyshev quadrature formula of the second kind *** 639 ! *** u [-1,+1] -> r [0,infinity] => r = (1 + u)/(1 - u) *** 640 641 DO i = 1, n 642 t = REAL(i, dp)*f 643 cost = COS(t) 644 sint = SIN(t) 645 sint2 = sint**2 646 x = REAL(2*i - n - 1, dp)/REAL(n + 1, dp) - & 647 2.0_dp*(1.0_dp + 2.0_dp*sint2/3.0_dp)*cost*sint/pi 648 w = 16.0_dp*sint2**2/REAL(3*(n + 1), dp) 649 r(n + 1 - i) = (1.0_dp + x)/(1.0_dp - x) 650 r2(n + 1 - i) = r(n + 1 - i)**2 651 wr(n + 1 - i) = 2.0_dp*w*r2(n + 1 - i)/(1.0_dp - x)**2 652 END DO 653 654 CASE (do_gapw_log) 655 656 ! *** Transformed Gauss-Chebyshev quadrature formula of the second kind *** 657 ! *** u [-1,+1] -> r [0,infinity] => r = ln(2/(1 - u))/ln(2) *** 658 659 DO i = 1, n 660 t = REAL(i, dp)*f 661 cost = COS(t) 662 sint = SIN(t) 663 sint2 = sint**2 664 x = REAL(2*i - n - 1, dp)/REAL(n + 1, dp) - & 665 2.0_dp*(1.0_dp + 2.0_dp*sint2/3.0_dp)*cost*sint/pi 666 w = 16.0_dp*sint2**2/REAL(3*(n + 1), dp) 667 r(n + 1 - i) = LOG(2.0_dp/(1.0_dp - x))/LOG(2.0_dp) 668 r2(n + 1 - i) = r(n + 1 - i)**2 669 wr(n + 1 - i) = w*r2(n + 1 - i)/(LOG(2.0_dp)*(1.0_dp - x)) 670 END DO 671 672 CASE DEFAULT 673 674 CPABORT("Invalid radial quadrature type specified") 675 676 END SELECT 677 678 END SUBROUTINE radial_grid 679 680END MODULE qs_grid_atom 681