Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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), DIMENSION(:), POINTER :: rad => NULL(), rad2 => NULL(), &
48 : wr => NULL(), wa => NULL(), &
49 : azi => NULL(), cos_azi => NULL(), sin_azi => NULL(), &
50 : pol => NULL(), cos_pol => NULL(), sin_pol => NULL(), usin_azi => NULL()
51 : REAL(dp), DIMENSION(:, :), &
52 : POINTER :: rad2l => NULL(), oorad2l => NULL(), weight => NULL()
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 11459 : SUBROUTINE allocate_grid_atom(grid_atom)
73 :
74 : TYPE(grid_atom_type), POINTER :: grid_atom
75 :
76 11459 : IF (ASSOCIATED(grid_atom)) CALL deallocate_grid_atom(grid_atom)
77 :
78 11459 : 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 11459 : 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 11459 : SUBROUTINE deallocate_grid_atom(grid_atom)
105 : TYPE(grid_atom_type), POINTER :: grid_atom
106 :
107 11459 : IF (ASSOCIATED(grid_atom)) THEN
108 :
109 11459 : IF (ASSOCIATED(grid_atom%rad)) THEN
110 11459 : DEALLOCATE (grid_atom%rad)
111 : END IF
112 :
113 11459 : IF (ASSOCIATED(grid_atom%rad2)) THEN
114 11459 : DEALLOCATE (grid_atom%rad2)
115 : END IF
116 :
117 11459 : IF (ASSOCIATED(grid_atom%wr)) THEN
118 11459 : DEALLOCATE (grid_atom%wr)
119 : END IF
120 :
121 11459 : IF (ASSOCIATED(grid_atom%wa)) THEN
122 11459 : DEALLOCATE (grid_atom%wa)
123 : END IF
124 :
125 11459 : IF (ASSOCIATED(grid_atom%weight)) THEN
126 11459 : DEALLOCATE (grid_atom%weight)
127 : END IF
128 :
129 11459 : IF (ASSOCIATED(grid_atom%azi)) THEN
130 11459 : DEALLOCATE (grid_atom%azi)
131 : END IF
132 :
133 11459 : IF (ASSOCIATED(grid_atom%cos_azi)) THEN
134 11459 : DEALLOCATE (grid_atom%cos_azi)
135 : END IF
136 :
137 11459 : IF (ASSOCIATED(grid_atom%sin_azi)) THEN
138 11459 : DEALLOCATE (grid_atom%sin_azi)
139 : END IF
140 :
141 11459 : IF (ASSOCIATED(grid_atom%pol)) THEN
142 11459 : DEALLOCATE (grid_atom%pol)
143 : END IF
144 :
145 11459 : IF (ASSOCIATED(grid_atom%cos_pol)) THEN
146 11459 : DEALLOCATE (grid_atom%cos_pol)
147 : END IF
148 :
149 11459 : IF (ASSOCIATED(grid_atom%sin_pol)) THEN
150 11459 : DEALLOCATE (grid_atom%sin_pol)
151 : END IF
152 :
153 11459 : IF (ASSOCIATED(grid_atom%usin_azi)) THEN
154 11459 : DEALLOCATE (grid_atom%usin_azi)
155 : END IF
156 :
157 11459 : IF (ASSOCIATED(grid_atom%rad2l)) THEN
158 11459 : DEALLOCATE (grid_atom%rad2l)
159 : END IF
160 :
161 11459 : IF (ASSOCIATED(grid_atom%oorad2l)) THEN
162 11459 : DEALLOCATE (grid_atom%oorad2l)
163 : END IF
164 :
165 11459 : DEALLOCATE (grid_atom)
166 : ELSE
167 : CALL cp_abort(__LOCATION__, &
168 : "The pointer grid_atom is not associated and "// &
169 0 : "cannot be deallocated")
170 : END IF
171 11459 : 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 11459 : 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 11459 : REAL(dp), DIMENSION(:), POINTER :: rad, rad2, wr
192 :
193 11459 : CALL timeset(routineN, handle)
194 :
195 11459 : NULLIFY (rad, rad2, wr)
196 :
197 11459 : IF (ASSOCIATED(grid_atom)) THEN
198 :
199 : ! Allocate the radial grid arrays
200 11459 : CALL reallocate(grid_atom%rad, 1, nr)
201 11459 : CALL reallocate(grid_atom%rad2, 1, nr)
202 11459 : CALL reallocate(grid_atom%wr, 1, nr)
203 11459 : CALL reallocate(grid_atom%wa, 1, na)
204 11459 : CALL reallocate(grid_atom%weight, 1, na, 1, nr)
205 11459 : CALL reallocate(grid_atom%azi, 1, na)
206 11459 : CALL reallocate(grid_atom%cos_azi, 1, na)
207 11459 : CALL reallocate(grid_atom%sin_azi, 1, na)
208 11459 : CALL reallocate(grid_atom%pol, 1, na)
209 11459 : CALL reallocate(grid_atom%cos_pol, 1, na)
210 11459 : CALL reallocate(grid_atom%sin_pol, 1, na)
211 11459 : CALL reallocate(grid_atom%usin_azi, 1, na)
212 11459 : CALL reallocate(grid_atom%rad2l, 1, nr, 0, llmax + 1)
213 11459 : CALL reallocate(grid_atom%oorad2l, 1, nr, 0, llmax + 1)
214 :
215 : ! Calculate the radial grid for this kind
216 11459 : rad => grid_atom%rad
217 11459 : rad2 => grid_atom%rad2
218 11459 : wr => grid_atom%wr
219 :
220 11459 : grid_atom%quadrature = quadrature
221 11459 : CALL radial_grid(nr, rad, rad2, wr, quadrature)
222 :
223 3909443 : grid_atom%rad2l(:, 0) = 1._dp
224 3909443 : grid_atom%oorad2l(:, 0) = 1._dp
225 39577 : DO l = 1, llmax + 1
226 16237254 : grid_atom%rad2l(:, l) = grid_atom%rad2l(:, l - 1)*rad(:)
227 16248713 : grid_atom%oorad2l(:, l) = grid_atom%oorad2l(:, l - 1)/rad(:)
228 : END DO
229 :
230 11459 : IF (ll > 0) THEN
231 112740 : grid_atom%wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:na)
232 118636 : DO ir = 1, nr
233 7212636 : DO ia = 1, na
234 7210560 : grid_atom%weight(ia, ir) = grid_atom%wr(ir)*grid_atom%wa(ia)
235 : END DO
236 : END DO
237 :
238 112740 : DO ia = 1, na
239 : ! polar angle: pol = acos(r(3))
240 110664 : cosia = lebedev_grid(ll)%r(3, ia)
241 110664 : grid_atom%cos_pol(ia) = cosia
242 : ! azimuthal angle: pol = atan(r(2)/r(1))
243 110664 : 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 4152 : grid_atom%azi(ia) = 0.0_dp
246 : ELSE
247 106512 : grid_atom%azi(ia) = ATAN2(lebedev_grid(ll)%r(2, ia), lebedev_grid(ll)%r(1, ia))
248 : END IF
249 110664 : grid_atom%cos_azi(ia) = COS(grid_atom%azi(ia))
250 110664 : pol = ACOS(cosia)
251 110664 : grid_atom%pol(ia) = pol
252 110664 : grid_atom%sin_pol(ia) = SIN(grid_atom%pol(ia))
253 :
254 110664 : grid_atom%sin_azi(ia) = SIN(grid_atom%azi(ia))
255 112740 : IF (ABS(grid_atom%sin_azi(ia)) > EPSILON(1.0_dp)) THEN
256 93640 : grid_atom%usin_azi(ia) = 1.0_dp/grid_atom%sin_azi(ia)
257 : ELSE
258 17024 : grid_atom%usin_azi(ia) = 1.0_dp
259 : END IF
260 :
261 : END DO
262 :
263 : END IF
264 :
265 : ELSE
266 0 : CPABORT("The pointer grid_atom is not associated")
267 : END IF
268 :
269 11459 : CALL timestop(handle)
270 :
271 11459 : 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 2 : 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 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: icell
295 : REAL(KIND=dp) :: ag, dd, dmax, r1, r2, r3
296 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rad, rad2, wa, wc, wr
297 2 : 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 2 : ALLOCATE (igr)
303 :
304 : ! type of quadrature grid
305 2 : IF (PRESENT(quadrature)) THEN
306 0 : my_quad = quadrature
307 : ELSE
308 2 : my_quad = do_gapw_log
309 : END IF
310 :
311 : ! radial grid
312 2 : CPASSERT(nr > 1)
313 10 : ALLOCATE (rad(nr), rad2(nr), wr(nr))
314 2 : CALL radial_grid(nr, rad, rad2, wr, my_quad)
315 : !
316 2 : igr%nr = nr
317 4 : ALLOCATE (igr%rr(nr))
318 4 : ALLOCATE (igr%wr(nr))
319 : ! store grid points always in ascending order
320 2 : IF (rad(1) > rad(nr)) THEN
321 102 : DO ir = nr, 1, -1
322 100 : igr%rr(nr - ir + 1) = rad(ir)
323 102 : igr%wr(nr - ir + 1) = wr(ir)
324 : END DO
325 : ELSE
326 0 : igr%rr(1:nr) = rad(1:nr)
327 0 : igr%wr(1:nr) = wr(1:nr)
328 : END IF
329 : ! only include grid points smaller than rmax
330 : np = 0
331 102 : DO ir = 1, nr
332 102 : IF (igr%rr(ir) < rmax) THEN
333 100 : np = np + 1
334 100 : rad(np) = igr%rr(ir)
335 100 : wr(np) = igr%wr(ir)
336 : END IF
337 : END DO
338 2 : igr%np = np
339 : !
340 : ! angular grid
341 2 : CPASSERT(na > 1)
342 2 : ll = get_number_of_lebedev_grid(n=na)
343 2 : np = lebedev_grid(ll)%n
344 2 : la = lebedev_grid(ll)%l
345 10 : ALLOCATE (rang(3, np), wa(np))
346 78 : wa(1:na) = 4._dp*pi*lebedev_grid(ll)%w(1:np)
347 306 : rang(1:3, 1:np) = lebedev_grid(ll)%r(1:3, 1:np)
348 2 : igr%lebedev_grid = ll
349 4 : ALLOCATE (igr%wa(np))
350 2 : igr%na = np
351 78 : igr%wa(1:np) = wa(1:np)
352 : !
353 : ! total grid points
354 2 : ntot = igr%na*igr%np
355 2 : igr%ntot = ntot
356 10 : ALLOCATE (rco(3, ntot), wc(ntot))
357 102 : ig = 0
358 102 : DO ir = 1, igr%np
359 3902 : DO ia = 1, igr%na
360 3800 : ig = ig + 1
361 15200 : rco(1:3, ig) = rang(1:3, ia)*rad(ir)
362 3900 : wc(ig) = wa(ia)*wr(ir)
363 : END DO
364 : END DO
365 : ! grid for batches, odd number of cells
366 2 : ng = NINT((REAL(ntot, dp)/32._dp)**0.33333_dp)
367 2 : ng = ng + MOD(ng + 1, 2)
368 : ! avarage number of points along radial grid
369 2 : dco = 0.0_dp
370 2 : ag = REAL(igr%np, dp)/ng
371 2 : CPASSERT(SIZE(dco) >= (ng + 1)/2)
372 2 : DO ig = 1, ng, 2
373 6 : ir = MIN(NINT(ag)*ig, igr%np)
374 6 : ia = (ig + 1)/2
375 6 : dco(ia) = rad(ir)
376 : END DO
377 : ! batches
378 6 : ALLOCATE (icell(ntot))
379 3802 : icell = 0
380 2 : nx = (ng - 1)/2
381 3802 : DO ig = 1, ntot
382 3800 : ix = grid_coord(rco(1, ig), dco, nx + 1) + nx
383 3800 : iy = grid_coord(rco(2, ig), dco, nx + 1) + nx
384 3800 : iz = grid_coord(rco(3, ig), dco, nx + 1) + nx
385 3802 : icell(ig) = iz*ng*ng + iy*ng + ix + 1
386 : END DO
387 : !
388 2 : igr%nbatch = ng*ng*ng
389 262 : ALLOCATE (igr%batch(igr%nbatch))
390 252 : igr%batch(:)%np = 0
391 3802 : DO ig = 1, ntot
392 3800 : ia = icell(ig)
393 3802 : igr%batch(ia)%np = igr%batch(ia)%np + 1
394 : END DO
395 252 : DO ig = 1, igr%nbatch
396 250 : np = igr%batch(ig)%np
397 1058 : ALLOCATE (igr%batch(ig)%rco(3, np), igr%batch(ig)%weight(np))
398 252 : igr%batch(ig)%np = 0
399 : END DO
400 3802 : DO ig = 1, ntot
401 3800 : ia = icell(ig)
402 3800 : igr%batch(ia)%np = igr%batch(ia)%np + 1
403 3800 : np = igr%batch(ia)%np
404 15200 : igr%batch(ia)%rco(1:3, np) = rco(1:3, ig)
405 3802 : igr%batch(ia)%weight(np) = wc(ig)
406 : END DO
407 : !
408 2 : DEALLOCATE (rad, rad2, rang, wr, wa)
409 2 : DEALLOCATE (rco, wc, icell)
410 : !
411 2 : IF (ASSOCIATED(int_grid)) CALL deallocate_atom_int_grid(int_grid)
412 2 : ALLOCATE (int_grid)
413 14 : ALLOCATE (int_grid%rr(igr%nr), int_grid%wr(igr%nr), int_grid%wa(igr%na))
414 2 : int_grid%nr = igr%nr
415 2 : int_grid%na = igr%na
416 2 : int_grid%np = igr%np
417 2 : int_grid%ntot = igr%ntot
418 2 : int_grid%lebedev_grid = igr%lebedev_grid
419 202 : int_grid%rr(:) = igr%rr(:)
420 202 : int_grid%wr(:) = igr%wr(:)
421 154 : int_grid%wa(:) = igr%wa(:)
422 : !
423 : ! count batches
424 2 : nbatch = 0
425 252 : DO ig = 1, igr%nbatch
426 252 : IF (igr%batch(ig)%np == 0) THEN
427 : ! empty batch
428 154 : ELSE IF (igr%batch(ig)%np <= 48) THEN
429 : ! single batch
430 152 : nbatch = nbatch + 1
431 : ELSE
432 : ! multiple batches
433 2 : nbatch = nbatch + NINT(igr%batch(ig)%np/32._dp)
434 : END IF
435 : END DO
436 2 : int_grid%nbatch = nbatch
437 188 : ALLOCATE (int_grid%batch(nbatch))
438 : ! fill batches
439 2 : n1 = 0
440 252 : DO ig = 1, igr%nbatch
441 252 : IF (igr%batch(ig)%np == 0) THEN
442 : ! empty batch
443 154 : ELSE IF (igr%batch(ig)%np <= 48) THEN
444 : ! single batch
445 152 : n1 = n1 + 1
446 152 : np = igr%batch(ig)%np
447 760 : ALLOCATE (int_grid%batch(n1)%rco(3, np), int_grid%batch(n1)%weight(np))
448 152 : int_grid%batch(n1)%np = np
449 24344 : int_grid%batch(n1)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, 1:np)
450 6200 : int_grid%batch(n1)%weight(1:np) = igr%batch(ig)%weight(1:np)
451 : ELSE
452 : ! multiple batches
453 2 : n2 = NINT(igr%batch(ig)%np/32._dp)
454 2 : n3 = igr%batch(ig)%np/n2
455 26 : DO ia = n1 + 1, n1 + n2
456 24 : nu = (ia - n1 - 1)*n3 + 1
457 24 : no = nu + n3 - 1
458 24 : IF (ia == n1 + n2) no = igr%batch(ig)%np
459 24 : np = no - nu + 1
460 120 : ALLOCATE (int_grid%batch(ia)%rco(3, np), int_grid%batch(ia)%weight(np))
461 24 : int_grid%batch(ia)%np = np
462 6232 : int_grid%batch(ia)%rco(1:3, 1:np) = igr%batch(ig)%rco(1:3, nu:no)
463 1578 : int_grid%batch(ia)%weight(1:np) = igr%batch(ig)%weight(nu:no)
464 : END DO
465 2 : n1 = n1 + n2
466 : END IF
467 : END DO
468 2 : CPASSERT(nbatch == n1)
469 : ! batch center and radius
470 178 : DO ig = 1, int_grid%nbatch
471 176 : np = int_grid%batch(ig)%np
472 176 : IF (np > 0) THEN
473 3976 : rm(1) = SUM(int_grid%batch(ig)%rco(1, 1:np))
474 3976 : rm(2) = SUM(int_grid%batch(ig)%rco(2, 1:np))
475 3976 : rm(3) = SUM(int_grid%batch(ig)%rco(3, 1:np))
476 704 : rm(1:3) = rm(1:3)/REAL(np, KIND=dp)
477 : ELSE
478 0 : rm(:) = 0.0_dp
479 : END IF
480 704 : int_grid%batch(ig)%rcenter(1:3) = rm(1:3)
481 : dmax = 0.0_dp
482 3976 : DO ia = 1, np
483 15200 : dd = SUM((int_grid%batch(ig)%rco(1:3, ia) - rm(1:3))**2)
484 3976 : dmax = MAX(dd, dmax)
485 : END DO
486 178 : int_grid%batch(ig)%rad = SQRT(dmax)
487 : END DO
488 : !
489 2 : CALL deallocate_atom_int_grid(igr)
490 : !
491 2 : IF (PRESENT(iunit)) THEN
492 2 : IF (iunit > 0) THEN
493 1 : WRITE (iunit, "(/,A)") " Atomic Integration Grid Information"
494 1 : WRITE (iunit, "(A,T51,3I10)") " Number of grid points [radial,angular,total]", &
495 2 : int_grid%np, int_grid%na, int_grid%ntot
496 1 : WRITE (iunit, "(A,T71,I10)") " Lebedev grid number", int_grid%lebedev_grid
497 1 : WRITE (iunit, "(A,T61,F20.5)") " Maximum of radial grid [Bohr]", &
498 2 : int_grid%rr(int_grid%np)
499 1 : nbatch = int_grid%nbatch
500 1 : WRITE (iunit, "(A,T71,I10)") " Total number of gridpoint batches", nbatch
501 1 : n1 = int_grid%ntot
502 1 : n2 = 0
503 1 : n3 = NINT(REAL(int_grid%ntot, dp)/REAL(nbatch, dp))
504 89 : DO ig = 1, nbatch
505 88 : n1 = MIN(n1, int_grid%batch(ig)%np)
506 89 : n2 = MAX(n2, int_grid%batch(ig)%np)
507 : END DO
508 1 : WRITE (iunit, "(A,T51,3I10)") " Number of grid points/batch [min,max,ave]", n1, n2, n3
509 1 : r1 = 1000._dp
510 1 : r2 = 0.0_dp
511 1 : r3 = 0.0_dp
512 89 : DO ig = 1, int_grid%nbatch
513 88 : r1 = MIN(r1, int_grid%batch(ig)%rad)
514 88 : r2 = MAX(r2, int_grid%batch(ig)%rad)
515 89 : r3 = r3 + int_grid%batch(ig)%rad
516 : END DO
517 1 : r3 = r3/REAL(ng*ng*ng, KIND=dp)
518 1 : WRITE (iunit, "(A,T51,3F10.2)") " Batch radius (bohr) [min,max,ave]", r1, r2, r3
519 : END IF
520 : END IF
521 :
522 2 : 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 11400 : 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 11400 : xval = ABS(x)
542 11400 : igrid = ng
543 20184 : DO ig = 1, ng
544 20184 : IF (xval <= dco(ig)) THEN
545 11400 : igrid = ig - 1
546 11400 : EXIT
547 : END IF
548 : END DO
549 11400 : IF (x < 0.0_dp) igrid = -igrid
550 11400 : CPASSERT(ABS(igrid) < ng)
551 11400 : END FUNCTION grid_coord
552 :
553 : ! **************************************************************************************************
554 : !> \brief ...
555 : !> \param int_grid ...
556 : ! **************************************************************************************************
557 4 : SUBROUTINE deallocate_atom_int_grid(int_grid)
558 : TYPE(atom_integration_grid_type), POINTER :: int_grid
559 :
560 : INTEGER :: ib
561 :
562 4 : IF (ASSOCIATED(int_grid)) THEN
563 4 : IF (ALLOCATED(int_grid%rr)) DEALLOCATE (int_grid%rr)
564 4 : IF (ALLOCATED(int_grid%wr)) DEALLOCATE (int_grid%wr)
565 4 : IF (ALLOCATED(int_grid%wa)) DEALLOCATE (int_grid%wa)
566 : ! batch
567 4 : IF (ALLOCATED(int_grid%batch)) THEN
568 430 : DO ib = 1, SIZE(int_grid%batch)
569 426 : IF (ALLOCATED(int_grid%batch(ib)%rco)) DEALLOCATE (int_grid%batch(ib)%rco)
570 430 : IF (ALLOCATED(int_grid%batch(ib)%weight)) DEALLOCATE (int_grid%batch(ib)%weight)
571 : END DO
572 430 : DEALLOCATE (int_grid%batch)
573 : END IF
574 : !
575 4 : DEALLOCATE (int_grid)
576 : NULLIFY (int_grid)
577 : END IF
578 :
579 4 : 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 11461 : 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 11461 : f = pi/REAL(n + 1, dp)
606 :
607 11461 : 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 2406 : DO i = 1, n
615 2400 : t = REAL(i, dp)*f
616 2400 : x = COS(t)
617 2400 : w = f*SIN(t)**2
618 2400 : r(i) = (1.0_dp + x)/(1.0_dp - x)
619 2400 : r2(i) = r(i)**2
620 2400 : wr(i) = w/SQRT(1.0_dp - x**2)
621 2406 : 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 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 :
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 3905535 : DO i = 1, n
648 3894084 : t = REAL(i, dp)*f
649 3894084 : cost = COS(t)
650 3894084 : sint = SIN(t)
651 3894084 : sint2 = sint**2
652 : x = REAL(2*i - n - 1, dp)/REAL(n + 1, dp) - &
653 3894084 : 2.0_dp*(1.0_dp + 2.0_dp*sint2/3.0_dp)*cost*sint/pi
654 3894084 : w = 16.0_dp*sint2**2/REAL(3*(n + 1), dp)
655 3894084 : r(n + 1 - i) = LOG(2.0_dp/(1.0_dp - x))/LOG(2.0_dp)
656 3894084 : r2(n + 1 - i) = r(n + 1 - i)**2
657 3905535 : 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 11461 : CPABORT("Invalid radial quadrature type specified")
663 :
664 : END SELECT
665 :
666 11461 : END SUBROUTINE radial_grid
667 :
668 0 : END MODULE qs_grid_atom
|
