Line data Source code
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 : ! **************************************************************************************************
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 = -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 :
63 : CONTAINS
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 12943 : SUBROUTINE allocate_grid_atom(grid_atom)
74 :
75 : TYPE(grid_atom_type), POINTER :: grid_atom
76 :
77 12943 : IF (ASSOCIATED(grid_atom)) CALL deallocate_grid_atom(grid_atom)
78 :
79 12943 : 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 12943 : 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 12943 : SUBROUTINE deallocate_grid_atom(grid_atom)
107 : TYPE(grid_atom_type), POINTER :: grid_atom
108 :
109 12943 : IF (ASSOCIATED(grid_atom)) THEN
110 :
111 12943 : IF (ASSOCIATED(grid_atom%rad)) THEN
112 12935 : DEALLOCATE (grid_atom%rad)
113 : END IF
114 :
115 12943 : IF (ASSOCIATED(grid_atom%rad2)) THEN
116 12935 : DEALLOCATE (grid_atom%rad2)
117 : END IF
118 :
119 12943 : IF (ASSOCIATED(grid_atom%wr)) THEN
120 12935 : DEALLOCATE (grid_atom%wr)
121 : END IF
122 :
123 12943 : IF (ASSOCIATED(grid_atom%wa)) THEN
124 12935 : DEALLOCATE (grid_atom%wa)
125 : END IF
126 :
127 12943 : IF (ASSOCIATED(grid_atom%weight)) THEN
128 12935 : DEALLOCATE (grid_atom%weight)
129 : END IF
130 :
131 12943 : IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
132 472 : DEALLOCATE (grid_atom%gapw_weight_s)
133 : END IF
134 :
135 12943 : IF (ASSOCIATED(grid_atom%azi)) THEN
136 12935 : DEALLOCATE (grid_atom%azi)
137 : END IF
138 :
139 12943 : IF (ASSOCIATED(grid_atom%cos_azi)) THEN
140 12935 : DEALLOCATE (grid_atom%cos_azi)
141 : END IF
142 :
143 12943 : IF (ASSOCIATED(grid_atom%sin_azi)) THEN
144 12935 : DEALLOCATE (grid_atom%sin_azi)
145 : END IF
146 :
147 12943 : IF (ASSOCIATED(grid_atom%pol)) THEN
148 12935 : DEALLOCATE (grid_atom%pol)
149 : END IF
150 :
151 12943 : IF (ASSOCIATED(grid_atom%cos_pol)) THEN
152 12935 : DEALLOCATE (grid_atom%cos_pol)
153 : END IF
154 :
155 12943 : IF (ASSOCIATED(grid_atom%sin_pol)) THEN
156 12935 : DEALLOCATE (grid_atom%sin_pol)
157 : END IF
158 :
159 12943 : IF (ASSOCIATED(grid_atom%usin_azi)) THEN
160 12935 : DEALLOCATE (grid_atom%usin_azi)
161 : END IF
162 :
163 12943 : IF (ASSOCIATED(grid_atom%rad2l)) THEN
164 12935 : DEALLOCATE (grid_atom%rad2l)
165 : END IF
166 :
167 12943 : IF (ASSOCIATED(grid_atom%oorad2l)) THEN
168 12935 : DEALLOCATE (grid_atom%oorad2l)
169 : END IF
170 :
171 12943 : DEALLOCATE (grid_atom)
172 : ELSE
173 : CALL cp_abort(__LOCATION__, &
174 : "The pointer grid_atom is not associated and "// &
175 0 : "cannot be deallocated")
176 : END IF
177 12943 : 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 12935 : 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 12935 : REAL(dp), DIMENSION(:), POINTER :: rad, rad2, wr
198 :
199 12935 : CALL timeset(routineN, handle)
200 :
201 12935 : NULLIFY (rad, rad2, wr)
202 :
203 12935 : IF (ASSOCIATED(grid_atom)) THEN
204 :
205 : ! Allocate the radial grid arrays
206 12935 : CALL reallocate(grid_atom%rad, 1, nr)
207 12935 : CALL reallocate(grid_atom%rad2, 1, nr)
208 12935 : CALL reallocate(grid_atom%wr, 1, nr)
209 12935 : CALL reallocate(grid_atom%wa, 1, na)
210 12935 : CALL reallocate(grid_atom%weight, 1, na, 1, nr)
211 12935 : CALL reallocate(grid_atom%azi, 1, na)
212 12935 : CALL reallocate(grid_atom%cos_azi, 1, na)
213 12935 : CALL reallocate(grid_atom%sin_azi, 1, na)
214 12935 : CALL reallocate(grid_atom%pol, 1, na)
215 12935 : CALL reallocate(grid_atom%cos_pol, 1, na)
216 12935 : CALL reallocate(grid_atom%sin_pol, 1, na)
217 12935 : CALL reallocate(grid_atom%usin_azi, 1, na)
218 12935 : CALL reallocate(grid_atom%rad2l, 1, nr, 0, llmax + 1)
219 12935 : CALL reallocate(grid_atom%oorad2l, 1, nr, 0, llmax + 1)
220 :
221 : ! Calculate the radial grid for this kind
222 12935 : rad => grid_atom%rad
223 12935 : rad2 => grid_atom%rad2
224 12935 : wr => grid_atom%wr
225 :
226 12935 : grid_atom%quadrature = quadrature
227 12935 : CALL radial_grid(nr, rad, rad2, wr, quadrature)
228 :
229 4269331 : grid_atom%rad2l(:, 0) = 1._dp
230 4269331 : grid_atom%oorad2l(:, 0) = 1._dp
231 45829 : DO l = 1, llmax + 1
232 17886758 : grid_atom%rad2l(:, l) = grid_atom%rad2l(:, l - 1)*rad(:)
233 17899693 : grid_atom%oorad2l(:, l) = grid_atom%oorad2l(:, l - 1)/rad(:)
234 : END DO
235 :
236 12935 : IF (ll > 0) THEN
237 149076 : grid_atom%wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:na)
238 155096 : DO ir = 1, nr
239 9508976 : DO ia = 1, na
240 9506220 : grid_atom%weight(ia, ir) = grid_atom%wr(ir)*grid_atom%wa(ia)
241 : END DO
242 : END DO
243 :
244 149076 : DO ia = 1, na
245 : ! polar angle: pol = acos(r(3))
246 146320 : cosia = lebedev_grid(ll)%r(3, ia)
247 146320 : grid_atom%cos_pol(ia) = cosia
248 : ! azimuthal angle: pol = atan(r(2)/r(1))
249 146320 : 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 5512 : grid_atom%azi(ia) = 0.0_dp
252 : ELSE
253 140808 : grid_atom%azi(ia) = ATAN2(lebedev_grid(ll)%r(2, ia), lebedev_grid(ll)%r(1, ia))
254 : END IF
255 146320 : grid_atom%cos_azi(ia) = COS(grid_atom%azi(ia))
256 146320 : pol = ACOS(cosia)
257 146320 : grid_atom%pol(ia) = pol
258 146320 : grid_atom%sin_pol(ia) = SIN(grid_atom%pol(ia))
259 :
260 146320 : grid_atom%sin_azi(ia) = SIN(grid_atom%azi(ia))
261 149076 : IF (ABS(grid_atom%sin_azi(ia)) > EPSILON(1.0_dp)) THEN
262 123768 : grid_atom%usin_azi(ia) = 1.0_dp/grid_atom%sin_azi(ia)
263 : ELSE
264 22552 : grid_atom%usin_azi(ia) = 1.0_dp
265 : END IF
266 :
267 : END DO
268 :
269 : END IF
270 :
271 : ELSE
272 0 : CPABORT("The pointer grid_atom is not associated")
273 : END IF
274 :
275 12935 : CALL timestop(handle)
276 :
277 12935 : 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 10 : 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 10 : INTEGER, ALLOCATABLE, DIMENSION(:) :: icell
301 : REAL(KIND=dp) :: ag, dd, dmax, r1, r2, r3
302 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rad, rad2, wa, wc, wr
303 10 : 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 10 : ALLOCATE (igr)
309 :
310 : ! type of quadrature grid
311 10 : IF (PRESENT(quadrature)) THEN
312 0 : my_quad = quadrature
313 : ELSE
314 10 : my_quad = do_gapw_log
315 : END IF
316 :
317 : ! radial grid
318 10 : CPASSERT(nr > 1)
319 50 : ALLOCATE (rad(nr), rad2(nr), wr(nr))
320 10 : CALL radial_grid(nr, rad, rad2, wr, my_quad)
321 : !
322 10 : igr%nr = nr
323 20 : ALLOCATE (igr%rr(nr))
324 20 : ALLOCATE (igr%wr(nr))
325 : ! store grid points always in ascending order
326 10 : IF (rad(1) > rad(nr)) THEN
327 510 : DO ir = nr, 1, -1
328 500 : igr%rr(nr - ir + 1) = rad(ir)
329 510 : igr%wr(nr - ir + 1) = wr(ir)
330 : END DO
331 : ELSE
332 0 : igr%rr(1:nr) = rad(1:nr)
333 0 : igr%wr(1:nr) = wr(1:nr)
334 : END IF
335 : ! only include grid points smaller than rmax
336 : np = 0
337 510 : DO ir = 1, nr
338 510 : IF (igr%rr(ir) < rmax) THEN
339 500 : np = np + 1
340 500 : rad(np) = igr%rr(ir)
341 500 : wr(np) = igr%wr(ir)
342 : END IF
343 : END DO
344 10 : igr%np = np
345 : !
346 : ! angular grid
347 10 : CPASSERT(na > 1)
348 10 : ll = get_number_of_lebedev_grid(n=na)
349 10 : np = lebedev_grid(ll)%n
350 10 : la = lebedev_grid(ll)%l
351 50 : ALLOCATE (rang(3, np), wa(np))
352 390 : wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:np)
353 1530 : rang(1:3, 1:np) = lebedev_grid(ll)%r(1:3, 1:np)
354 10 : igr%lebedev_grid = ll
355 20 : ALLOCATE (igr%wa(np))
356 10 : igr%na = np
357 390 : igr%wa(1:np) = wa(1:np)
358 : !
359 : ! total grid points
360 10 : ntot = igr%na*igr%np
361 10 : igr%ntot = ntot
362 50 : ALLOCATE (rco(3, ntot), wc(ntot))
363 510 : ig = 0
364 510 : DO ir = 1, igr%np
365 19510 : DO ia = 1, igr%na
366 19000 : ig = ig + 1
367 76000 : rco(1:3, ig) = rang(1:3, ia)*rad(ir)
368 19500 : wc(ig) = wa(ia)*wr(ir)
369 : END DO
370 : END DO
371 : ! grid for batches, odd number of cells
372 10 : ng = NINT((REAL(ntot, dp)/32._dp)**0.33333_dp)
373 10 : ng = ng + MOD(ng + 1, 2)
374 : ! avarage number of points along radial grid
375 10 : dco = 0.0_dp
376 10 : ag = REAL(igr%np, dp)/ng
377 10 : CPASSERT(SIZE(dco) >= (ng + 1)/2)
378 10 : DO ig = 1, ng, 2
379 30 : ir = MIN(NINT(ag)*ig, igr%np)
380 30 : ia = (ig + 1)/2
381 30 : dco(ia) = rad(ir)
382 : END DO
383 : ! batches
384 30 : ALLOCATE (icell(ntot))
385 19010 : icell = 0
386 10 : nx = (ng - 1)/2
387 19010 : DO ig = 1, ntot
388 19000 : ix = grid_coord(rco(1, ig), dco, nx + 1) + nx
389 19000 : iy = grid_coord(rco(2, ig), dco, nx + 1) + nx
390 19000 : iz = grid_coord(rco(3, ig), dco, nx + 1) + nx
391 19010 : icell(ig) = iz*ng*ng + iy*ng + ix + 1
392 : END DO
393 : !
394 10 : igr%nbatch = ng*ng*ng
395 1310 : ALLOCATE (igr%batch(igr%nbatch))
396 1260 : igr%batch(:)%np = 0
397 19010 : DO ig = 1, ntot
398 19000 : ia = icell(ig)
399 19010 : igr%batch(ia)%np = igr%batch(ia)%np + 1
400 : END DO
401 1260 : DO ig = 1, igr%nbatch
402 1250 : np = igr%batch(ig)%np
403 5290 : ALLOCATE (igr%batch(ig)%rco(3, np), igr%batch(ig)%weight(np))
404 1260 : igr%batch(ig)%np = 0
405 : END DO
406 19010 : DO ig = 1, ntot
407 19000 : ia = icell(ig)
408 19000 : igr%batch(ia)%np = igr%batch(ia)%np + 1
409 19000 : np = igr%batch(ia)%np
410 76000 : igr%batch(ia)%rco(1:3, np) = rco(1:3, ig)
411 19010 : igr%batch(ia)%weight(np) = wc(ig)
412 : END DO
413 : !
414 10 : DEALLOCATE (rad, rad2, rang, wr, wa)
415 10 : DEALLOCATE (rco, wc, icell)
416 : !
417 10 : IF (ASSOCIATED(int_grid)) CALL deallocate_atom_int_grid(int_grid)
418 10 : ALLOCATE (int_grid)
419 70 : ALLOCATE (int_grid%rr(igr%nr), int_grid%wr(igr%nr), int_grid%wa(igr%na))
420 10 : int_grid%nr = igr%nr
421 10 : int_grid%na = igr%na
422 10 : int_grid%np = igr%np
423 10 : int_grid%ntot = igr%ntot
424 10 : int_grid%lebedev_grid = igr%lebedev_grid
425 1010 : int_grid%rr(:) = igr%rr(:)
426 1010 : int_grid%wr(:) = igr%wr(:)
427 770 : int_grid%wa(:) = igr%wa(:)
428 : !
429 : ! count batches
430 10 : nbatch = 0
431 1260 : DO ig = 1, igr%nbatch
432 1260 : IF (igr%batch(ig)%np == 0) THEN
433 : ! empty batch
434 770 : ELSE IF (igr%batch(ig)%np <= 48) THEN
435 : ! single batch
436 760 : nbatch = nbatch + 1
437 : ELSE
438 : ! multiple batches
439 10 : nbatch = nbatch + NINT(igr%batch(ig)%np/32._dp)
440 : END IF
441 : END DO
442 10 : int_grid%nbatch = nbatch
443 940 : ALLOCATE (int_grid%batch(nbatch))
444 : ! fill batches
445 10 : n1 = 0
446 1260 : DO ig = 1, igr%nbatch
447 1260 : IF (igr%batch(ig)%np == 0) THEN
448 : ! empty batch
449 770 : ELSE IF (igr%batch(ig)%np <= 48) THEN
450 : ! single batch
451 760 : n1 = n1 + 1
452 760 : np = igr%batch(ig)%np
453 3800 : ALLOCATE (int_grid%batch(n1)%rco(3, np), int_grid%batch(n1)%weight(np))
454 760 : int_grid%batch(n1)%np = np
455 121720 : int_grid%batch(n1)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, 1:np)
456 31000 : int_grid%batch(n1)%weight(1:np) = igr%batch(ig)%weight(1:np)
457 : ELSE
458 : ! multiple batches
459 10 : n2 = NINT(igr%batch(ig)%np/32._dp)
460 10 : n3 = igr%batch(ig)%np/n2
461 130 : DO ia = n1 + 1, n1 + n2
462 120 : nu = (ia - n1 - 1)*n3 + 1
463 120 : no = nu + n3 - 1
464 120 : IF (ia == n1 + n2) no = igr%batch(ig)%np
465 120 : np = no - nu + 1
466 600 : ALLOCATE (int_grid%batch(ia)%rco(3, np), int_grid%batch(ia)%weight(np))
467 120 : int_grid%batch(ia)%np = np
468 31160 : int_grid%batch(ia)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, nu:no)
469 7890 : int_grid%batch(ia)%weight(1:np) = igr%batch(ig)%weight(nu:no)
470 : END DO
471 10 : n1 = n1 + n2
472 : END IF
473 : END DO
474 10 : CPASSERT(nbatch == n1)
475 : ! batch center and radius
476 890 : DO ig = 1, int_grid%nbatch
477 880 : np = int_grid%batch(ig)%np
478 880 : IF (np > 0) THEN
479 19880 : rm(1) = SUM(int_grid%batch(ig)%rco(1, 1:np))
480 19880 : rm(2) = SUM(int_grid%batch(ig)%rco(2, 1:np))
481 19880 : rm(3) = SUM(int_grid%batch(ig)%rco(3, 1:np))
482 3520 : rm(1:3) = rm(1:3)/REAL(np, KIND=dp)
483 : ELSE
484 0 : rm(:) = 0.0_dp
485 : END IF
486 3520 : int_grid%batch(ig)%rcenter(1:3) = rm(1:3)
487 : dmax = 0.0_dp
488 19880 : DO ia = 1, np
489 76000 : dd = SUM((int_grid%batch(ig)%rco(1:3, ia) - rm(1:3))**2)
490 19880 : dmax = MAX(dd, dmax)
491 : END DO
492 890 : int_grid%batch(ig)%rad = SQRT(dmax)
493 : END DO
494 : !
495 10 : CALL deallocate_atom_int_grid(igr)
496 : !
497 10 : IF (PRESENT(iunit)) THEN
498 10 : IF (iunit > 0) THEN
499 5 : WRITE (iunit, "(/,A)") " Atomic Integration Grid Information"
500 5 : WRITE (iunit, "(A,T51,3I10)") " Number of grid points [radial,angular,total]", &
501 10 : int_grid%np, int_grid%na, int_grid%ntot
502 5 : WRITE (iunit, "(A,T71,I10)") " Lebedev grid number", int_grid%lebedev_grid
503 5 : WRITE (iunit, "(A,T61,F20.5)") " Maximum of radial grid [Bohr]", &
504 10 : int_grid%rr(int_grid%np)
505 5 : nbatch = int_grid%nbatch
506 5 : WRITE (iunit, "(A,T71,I10)") " Total number of gridpoint batches", nbatch
507 5 : n1 = int_grid%ntot
508 5 : n2 = 0
509 5 : n3 = NINT(REAL(int_grid%ntot, dp)/REAL(nbatch, dp))
510 445 : DO ig = 1, nbatch
511 440 : n1 = MIN(n1, int_grid%batch(ig)%np)
512 445 : n2 = MAX(n2, int_grid%batch(ig)%np)
513 : END DO
514 5 : WRITE (iunit, "(A,T51,3I10)") " Number of grid points/batch [min,max,ave]", n1, n2, n3
515 5 : r1 = 1000._dp
516 5 : r2 = 0.0_dp
517 5 : r3 = 0.0_dp
518 445 : DO ig = 1, int_grid%nbatch
519 440 : r1 = MIN(r1, int_grid%batch(ig)%rad)
520 440 : r2 = MAX(r2, int_grid%batch(ig)%rad)
521 445 : r3 = r3 + int_grid%batch(ig)%rad
522 : END DO
523 5 : r3 = r3/REAL(ng*ng*ng, KIND=dp)
524 5 : WRITE (iunit, "(A,T51,3F10.2)") " Batch radius (bohr) [min,max,ave]", r1, r2, r3
525 : END IF
526 : END IF
527 :
528 10 : 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 57000 : 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 57000 : xval = ABS(x)
548 57000 : igrid = ng
549 100920 : DO ig = 1, ng
550 100920 : IF (xval <= dco(ig)) THEN
551 57000 : igrid = ig - 1
552 57000 : EXIT
553 : END IF
554 : END DO
555 57000 : IF (x < 0.0_dp) igrid = -igrid
556 57000 : CPASSERT(ABS(igrid) < ng)
557 57000 : END FUNCTION grid_coord
558 :
559 : ! **************************************************************************************************
560 : !> \brief ...
561 : !> \param int_grid ...
562 : ! **************************************************************************************************
563 20 : SUBROUTINE deallocate_atom_int_grid(int_grid)
564 : TYPE(atom_integration_grid_type), POINTER :: int_grid
565 :
566 : INTEGER :: ib
567 :
568 20 : IF (ASSOCIATED(int_grid)) THEN
569 20 : IF (ALLOCATED(int_grid%rr)) DEALLOCATE (int_grid%rr)
570 20 : IF (ALLOCATED(int_grid%wr)) DEALLOCATE (int_grid%wr)
571 20 : IF (ALLOCATED(int_grid%wa)) DEALLOCATE (int_grid%wa)
572 : ! batch
573 20 : IF (ALLOCATED(int_grid%batch)) THEN
574 2150 : DO ib = 1, SIZE(int_grid%batch)
575 2130 : IF (ALLOCATED(int_grid%batch(ib)%rco)) DEALLOCATE (int_grid%batch(ib)%rco)
576 2150 : IF (ALLOCATED(int_grid%batch(ib)%weight)) DEALLOCATE (int_grid%batch(ib)%weight)
577 : END DO
578 2150 : DEALLOCATE (int_grid%batch)
579 : END IF
580 : !
581 20 : DEALLOCATE (int_grid)
582 : NULLIFY (int_grid)
583 : END IF
584 :
585 20 : 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 12945 : 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 12945 : f = pi/REAL(n + 1, KIND=dp)
612 :
613 12945 : 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 2406 : DO i = 1, n
618 2400 : t = REAL(i, dp)*f
619 2400 : x = COS(t)
620 2400 : w = f*SIN(t)**2
621 2400 : r(i) = (1.0_dp + x)/(1.0_dp - x)
622 2400 : r2(i) = r(i)**2
623 2400 : wr(i) = w/SQRT(1.0_dp - x**2)
624 2406 : 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 1604 : DO i = 1, n
630 1600 : t = REAL(i, dp)*f
631 1600 : cost = COS(t)
632 1600 : sint = SIN(t)
633 1600 : sint2 = sint**2
634 : x = REAL(2*i - n - 1, dp)/REAL(n + 1, dp) - &
635 1600 : 2.0_dp*(1.0_dp + 2.0_dp*sint2/3.0_dp)*cost*sint/pi
636 1600 : w = 16.0_dp*sint2**2/REAL(3*(n + 1), dp)
637 1600 : r(n + 1 - i) = (1.0_dp + x)/(1.0_dp - x)
638 1600 : r2(n + 1 - i) = r(n + 1 - i)**2
639 1604 : 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 4265831 : DO i = 1, n
645 4252896 : t = REAL(i, dp)*f
646 4252896 : cost = COS(t)
647 4252896 : sint = SIN(t)
648 4252896 : sint2 = sint**2
649 : x = REAL(2*i - n - 1, dp)/REAL(n + 1, dp) - &
650 4252896 : 2.0_dp*(1.0_dp + 2.0_dp*sint2/3.0_dp)*cost*sint/pi
651 4252896 : w = 16.0_dp*sint2**2/REAL(3*(n + 1), dp)
652 4252896 : r(n + 1 - i) = LOG(2.0_dp/(1.0_dp - x))/LOG(2.0_dp)
653 4252896 : r2(n + 1 - i) = r(n + 1 - i)**2
654 4265831 : wr(n + 1 - i) = w*r2(n + 1 - i)/(LOG(2.0_dp)*(1.0_dp - x))
655 : END DO
656 : CASE DEFAULT
657 12945 : CPABORT("Invalid radial quadrature type specified")
658 : END SELECT
659 :
660 12945 : END SUBROUTINE radial_grid
661 :
662 0 : END MODULE qs_grid_atom
|
