(git:e7e05ae)
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 ! **************************************************************************************************
8 MODULE 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 
62 CONTAINS
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 
668 END MODULE qs_grid_atom
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