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