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