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 : ! **************************************************************************************************
9 : !> \brief Initialize atom-owned RI-RS grids by Cholesky selection from Lebedev grids.
10 : ! **************************************************************************************************
11 : MODULE gw_ri_rs_grid_initialization
12 : USE basis_set_types, ONLY: gto_basis_set_type
13 : USE gw_ri_rs_utils, ONLY: evaluate_ao_basis_on_points,&
14 : filter_grid_to_voronoi,&
15 : get_rirs_cluster_atoms
16 : USE input_constants, ONLY: do_gapw_log
17 : USE kinds, ONLY: default_string_length,&
18 : dp
19 : USE lebedev, ONLY: get_number_of_lebedev_grid,&
20 : lebedev_grid
21 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type,&
22 : rirs_grid_type
23 : USE qs_grid_atom, ONLY: allocate_grid_atom,&
24 : create_grid_atom,&
25 : deallocate_grid_atom,&
26 : grid_atom_type
27 : USE util, ONLY: locate
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_ri_rs_grid_initialization'
33 :
34 : PUBLIC :: initialize_ri_rs_grid
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Construct the initial RI-RS grid of every atom.
40 : !>
41 : !> For every atom B, the Lebedev grid points are
42 : !>
43 : !> (1) r_l = R_B + ρ_s Ω_a,
44 : !>
45 : !> where ρ_s is a radial quadrature point and Ω_a is a Lebedev direction. Lebedev grids are
46 : !> tabulated according to the angular degree L up to which they integrate exactly. To integrate a
47 : !> three-centre integral (μν|P), the smallest available Lebedev grid satisfying
48 : !>
49 : !> (2) L >= 2 l_AO,max + l_RI,max + ΔL
50 : !>
51 : !> is used. Here, l_AO,max is the maximum angular momentum of the atomic AO basis functions ϕ_μ,
52 : !> l_RI,max is the maximum angular momentum of the auxiliary basis functions φ_P, and the internal
53 : !> angular buffer is ΔL. Enough radial points ρ_s are used that the Lebedev grid of every atom B
54 : !> contains at least
55 : !>
56 : !> (3) N_initial^B >= α_initial N_AO^B,
57 : !>
58 : !> points, where N_AO^B is the number of atomic orbitals on atom B and
59 : !> α_initial = max(30, 2 RS_AO_RATIO). The factor two provides candidates that can be discarded by
60 : !> the molecular Voronoi filter. The requested initial RI-RS grid of atom A contains
61 : !>
62 : !> (4) N_R^A = ceil(RS_AO_RATIO * N_AO^A)
63 : !>
64 : !> points. Around every atom A, a cluster is defined as
65 : !>
66 : !> (5) C_A = {B : |R_B - R_A| <= CUTOFF_ATOMIC_CLUSTER}.
67 : !>
68 : !> The Cholesky selection uses all Lebedev grid points of every atom B in C_A. The molecular
69 : !> Voronoi cell of atom A is
70 : !>
71 : !> (6) V_A = {r_l : |r_l - R_A| <= |r_l - R_B| for every atom B}.
72 : !>
73 : !> For cluster C_A, Cholesky selection is performed on
74 : !>
75 : !> (7) D_ll' = [Σ_(μ in C_A) ϕ_μ(r_l) ϕ_μ(r_l')]^2.
76 : !>
77 : !> The first grid point is the point with the largest diagonal element,
78 : !>
79 : !> (8) d_l^(0) = D_ll, q_1 = arg max_l d_l^(0).
80 : !>
81 : !> For every selected point q_k, the Cholesky column and diagonal are updated according to
82 : !>
83 : !> (9) L_lk = [D_lq_k - Σ_(j<k) L_lj L_q_kj]/sqrt(d_q_k^(k-1)),
84 : !>
85 : !> (10) d_l^(k) = max(0, d_l^(k-1) - L_lk^2),
86 : !>
87 : !> and the next point q_(k+1) is the point with the largest d_l^(k). Points outside V_A take part
88 : !> in Eqs. (7)-(10), but only selected points inside V_A are placed into the RI-RS grid of atom A.
89 : !> If D_ll' reaches its numerical rank before N_R^A points have been retained, let S_A contain the
90 : !> retained points. Every unused point r_l in V_A is assigned the distance
91 : !>
92 : !> (11) ρ_l = min_(q in S_A) |r_l - q|,
93 : !>
94 : !> and the point
95 : !>
96 : !> (12) q_new = arg max_(r_l in V_A and r_l not in S_A) ρ_l
97 : !>
98 : !> is appended to S_A. Equations (11)-(12) are repeated until Eq. (4) is satisfied. This maximin
99 : !> distance is only a geometric selection criterion, not the three-centre-integral fitting error.
100 : !> Each cluster C_A is processed independently on one MPI rank; only the completed atom grids are
101 : !> communicated.
102 : !> \param bs_env Band-structure environment containing GW parameters.
103 : ! **************************************************************************************************
104 10 : SUBROUTINE initialize_ri_rs_grid(bs_env)
105 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
106 :
107 : CHARACTER(LEN=*), PARAMETER :: routineN = 'initialize_ri_rs_grid'
108 : INTEGER, PARAMETER :: l_additional = 2, radial_quadrature = do_gapw_log
109 : REAL(KIND=dp), PARAMETER :: minimum_candidate_ratio = 30.0_dp
110 :
111 : INTEGER :: handle, ikind
112 : REAL(KIND=dp) :: candidate_ratio
113 : TYPE(rirs_grid_type), ALLOCATABLE :: radial_lebedev_grids(:)
114 :
115 10 : CALL timeset(routineN, handle)
116 :
117 10 : candidate_ratio = MAX(minimum_candidate_ratio, 2.0_dp*bs_env%ri_rs%grid_opt%rs_ao_ratio)
118 :
119 : ! Build the Lebedev grids of Eqs. (1)-(3) for all elements in the calculation
120 50 : ALLOCATE (radial_lebedev_grids(SIZE(bs_env%basis_set_AO)))
121 30 : DO ikind = 1, SIZE(bs_env%basis_set_AO)
122 : CALL build_lebedev_grid(bs_env%basis_set_AO(ikind)%gto_basis_set, &
123 : bs_env%basis_set_RI(ikind)%gto_basis_set, &
124 : candidate_ratio, l_additional, radial_quadrature, &
125 30 : radial_lebedev_grids(ikind))
126 : END DO
127 :
128 : ! Apply Eqs. (5)-(12) for every atom and keep points inside Voronoi volume, Eq. (6)
129 10 : CALL cholesky_selection_and_voronoi_filtering(bs_env, radial_lebedev_grids)
130 :
131 10 : CALL broadcast_ri_rs_grids(bs_env)
132 :
133 10 : CALL timestop(handle)
134 :
135 40 : END SUBROUTINE initialize_ri_rs_grid
136 :
137 : ! **************************************************************************************************
138 : !> \brief Build one fixed-orientation Lebedev grid according to Eqs. (1)-(3).
139 : !> \param ao ...
140 : !> \param ri ...
141 : !> \param ratio α_initial in Eq. (3).
142 : !> \param l_additional ΔL in Eq. (2).
143 : !> \param radial_quadrature ...
144 : !> \param grid ...
145 : ! **************************************************************************************************
146 20 : SUBROUTINE build_lebedev_grid(ao, ri, ratio, l_additional, radial_quadrature, grid)
147 : TYPE(gto_basis_set_type), POINTER :: ao, ri
148 : REAL(KIND=dp), INTENT(IN) :: ratio
149 : INTEGER, INTENT(IN) :: l_additional, radial_quadrature
150 : TYPE(rirs_grid_type), INTENT(OUT) :: grid
151 :
152 : INTEGER :: degree, ir, nang, nrad, offset, rule
153 : TYPE(grid_atom_type), POINTER :: radial_grid
154 :
155 20 : CPASSERT(ASSOCIATED(ao) .AND. ASSOCIATED(ri))
156 270 : degree = 2*MAXVAL(ao%lmax) + MAXVAL(ri%lmax) + l_additional
157 20 : rule = get_number_of_lebedev_grid(l=degree)
158 20 : nang = lebedev_grid(rule)%n
159 20 : nrad = MAX(2, CEILING(ratio*REAL(ao%nsgf, dp)/REAL(nang, dp)))
160 :
161 20 : NULLIFY (radial_grid)
162 20 : CALL allocate_grid_atom(radial_grid)
163 20 : CALL create_grid_atom(radial_grid, nrad, nang, 0, rule, radial_quadrature)
164 20 : grid%npts = nrad*nang
165 60 : ALLOCATE (grid%raw_points(3, grid%npts))
166 148 : DO ir = 1, nrad
167 128 : offset = (ir - 1)*nang
168 : grid%raw_points(:, offset + 1:offset + nang) = &
169 10388 : radial_grid%rad(ir)*lebedev_grid(rule)%r
170 : END DO
171 20 : CALL deallocate_grid_atom(radial_grid)
172 20 : END SUBROUTINE build_lebedev_grid
173 :
174 : ! **************************************************************************************************
175 : !> \brief Apply Eqs. (5)-(12) independently for every atom and keep points inside Eq. (6).
176 : !> \param bs_env ...
177 : !> \param radial_lebedev_grids ...
178 : ! **************************************************************************************************
179 10 : SUBROUTINE cholesky_selection_and_voronoi_filtering(bs_env, radial_lebedev_grids)
180 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
181 : TYPE(rirs_grid_type), INTENT(IN) :: radial_lebedev_grids(:)
182 :
183 : INTEGER :: atom_mpi_rank, iatom, ikind, npoints
184 :
185 10 : CPASSERT(.NOT. ALLOCATED(bs_env%ri_rs%atomic_grids))
186 60 : ALLOCATE (bs_env%ri_rs%atomic_grids(bs_env%n_atom))
187 40 : DO iatom = 1, bs_env%n_atom
188 30 : ikind = bs_env%ri_rs%particle_set(iatom)%atomic_kind%kind_number
189 : npoints = CEILING(bs_env%ri_rs%grid_opt%rs_ao_ratio* &
190 30 : bs_env%basis_set_AO(ikind)%gto_basis_set%nsgf)
191 30 : bs_env%ri_rs%atomic_grids(iatom)%npts = npoints
192 90 : ALLOCATE (bs_env%ri_rs%atomic_grids(iatom)%raw_points(3, npoints))
193 :
194 30 : atom_mpi_rank = MOD(iatom - 1, bs_env%para_env%num_pe)
195 40 : IF (atom_mpi_rank == bs_env%para_env%mepos) THEN
196 15 : CALL select_cholesky_grid_iatom(bs_env, radial_lebedev_grids, iatom)
197 : END IF
198 : END DO
199 10 : END SUBROUTINE cholesky_selection_and_voronoi_filtering
200 :
201 : ! **************************************************************************************************
202 : !> \brief Select one atom grid with Eqs. (5)-(12) and the molecular Voronoi cell in Eq. (6).
203 : !> \param bs_env ...
204 : !> \param radial_lebedev_grids ...
205 : !> \param iatom ...
206 : ! **************************************************************************************************
207 15 : SUBROUTINE select_cholesky_grid_iatom(bs_env, radial_lebedev_grids, iatom)
208 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
209 : TYPE(rirs_grid_type), INTENT(IN) :: radial_lebedev_grids(:)
210 : INTEGER, INTENT(IN) :: iatom
211 :
212 : CHARACTER(LEN=3*default_string_length) :: failure
213 : INTEGER :: igrid_point, npoints, nretained
214 15 : INTEGER, ALLOCATABLE :: ao_point_indices(:), atom_n_ao(:), atom_point_offsets(:), &
215 15 : atom_value_offsets(:), cluster_atoms(:), selected_indices(:)
216 : LOGICAL, ALLOCATABLE :: inside_voronoi(:)
217 15 : REAL(KIND=dp), ALLOCATABLE :: ao_values(:), cluster_points(:, :)
218 :
219 : CALL get_rirs_cluster_atoms(bs_env%ri_rs%particle_set, bs_env%ri_rs%cell, iatom, &
220 15 : bs_env%ri_rs%grid_opt%cutoff_atomic_cluster, cluster_atoms)
221 : CALL collect_cluster_lebedev_points(bs_env, radial_lebedev_grids, iatom, &
222 15 : cluster_atoms, cluster_points)
223 15 : npoints = SIZE(cluster_points, 2)
224 45 : ALLOCATE (inside_voronoi(npoints))
225 : CALL filter_grid_to_voronoi(cluster_points, iatom, bs_env%ri_rs%particle_set, &
226 15 : mask=inside_voronoi)
227 4599 : IF (COUNT(inside_voronoi) < bs_env%ri_rs%atomic_grids(iatom)%npts) THEN
228 0 : WRITE (failure, '(A,I0,A,I0,A,I0,A)') 'Atom ', iatom, ': only ', COUNT(inside_voronoi), &
229 0 : ' Voronoi grid points for ', bs_env%ri_rs%atomic_grids(iatom)%npts, &
230 0 : ' points; reduce RS_AO_RATIO.'
231 0 : CPABORT(TRIM(failure))
232 : END IF
233 :
234 : CALL evaluate_cluster_ao_values(bs_env, iatom, cluster_atoms, cluster_points, &
235 : ao_point_indices, atom_point_offsets, atom_value_offsets, &
236 15 : atom_n_ao, ao_values)
237 45 : ALLOCATE (selected_indices(bs_env%ri_rs%atomic_grids(iatom)%npts))
238 : CALL select_cholesky_grid_points(ao_point_indices, atom_point_offsets, atom_value_offsets, &
239 : atom_n_ao, ao_values, npoints, &
240 15 : inside_voronoi, selected_indices)
241 396 : nretained = COUNT(selected_indices > 0)
242 15 : IF (nretained == 0) THEN
243 0 : WRITE (failure, '(A,I0,A)') 'Atom ', iatom, &
244 0 : ': Cholesky selection found no numerically independent point.'
245 0 : CPABORT(TRIM(failure))
246 15 : ELSE IF (nretained < SIZE(selected_indices)) THEN
247 5 : CALL complete_grid_by_maximin_distance(cluster_points, inside_voronoi, selected_indices)
248 : END IF
249 396 : DO igrid_point = 1, SIZE(selected_indices)
250 : bs_env%ri_rs%atomic_grids(iatom)%raw_points(:, igrid_point) = &
251 1539 : cluster_points(:, selected_indices(igrid_point))
252 : END DO
253 15 : END SUBROUTINE select_cholesky_grid_iatom
254 :
255 : ! **************************************************************************************************
256 : !> \brief Complete a rank-saturated atom grid with the geometric maximin rule in Eqs. (11)-(12).
257 : !> \param cluster_points Coordinates r_l of all existing cluster-grid points.
258 : !> \param inside_voronoi True for points inside the molecular Voronoi cell V_A.
259 : !> \param selected_indices Cholesky indices on entry and the completed indices on exit.
260 : ! **************************************************************************************************
261 5 : SUBROUTINE complete_grid_by_maximin_distance(cluster_points, inside_voronoi, selected_indices)
262 : REAL(KIND=dp), INTENT(IN) :: cluster_points(:, :)
263 : LOGICAL, INTENT(IN) :: inside_voronoi(:)
264 : INTEGER, INTENT(INOUT) :: selected_indices(:)
265 :
266 : INTEGER :: igrid_point, inew_point, iselected, &
267 : nretained
268 5 : LOGICAL, ALLOCATABLE :: unselected_inside(:)
269 : REAL(KIND=dp) :: distance_sq, largest_distance_sq
270 5 : REAL(KIND=dp), ALLOCATABLE :: nearest_distance_sq(:)
271 :
272 5 : CPASSERT(SIZE(cluster_points, 2) == SIZE(inside_voronoi))
273 2269 : CPASSERT(COUNT(inside_voronoi) >= SIZE(selected_indices))
274 :
275 321 : nretained = COUNT(selected_indices > 0)
276 5 : CPASSERT(nretained > 0)
277 : ALLOCATE (unselected_inside(SIZE(inside_voronoi)), &
278 25 : nearest_distance_sq(SIZE(inside_voronoi)))
279 2269 : unselected_inside(:) = inside_voronoi
280 2269 : nearest_distance_sq(:) = HUGE(1.0_dp)
281 :
282 : ! Initialize ρ_l^2 in Eq. (11) from all retained Cholesky points.
283 59 : DO iselected = 1, nretained
284 54 : unselected_inside(selected_indices(iselected)) = .FALSE.
285 23723 : DO igrid_point = 1, SIZE(inside_voronoi)
286 23664 : IF (.NOT. unselected_inside(igrid_point)) CYCLE
287 : distance_sq = SUM((cluster_points(:, igrid_point) - &
288 40372 : cluster_points(:, selected_indices(iselected)))**2)
289 : nearest_distance_sq(igrid_point) = &
290 23718 : MIN(nearest_distance_sq(igrid_point), distance_sq)
291 : END DO
292 : END DO
293 :
294 267 : DO WHILE (nretained < SIZE(selected_indices))
295 149334 : inew_point = 0
296 149334 : largest_distance_sq = -1.0_dp
297 149334 : DO igrid_point = 1, SIZE(inside_voronoi)
298 149072 : IF (.NOT. unselected_inside(igrid_point)) CYCLE
299 : ! Strict comparison makes the lowest cluster-point index win an exact tie.
300 44771 : IF (nearest_distance_sq(igrid_point) > largest_distance_sq) THEN
301 701 : largest_distance_sq = nearest_distance_sq(igrid_point)
302 701 : inew_point = igrid_point
303 : END IF
304 : END DO
305 262 : CPASSERT(inew_point > 0)
306 :
307 262 : nretained = nretained + 1
308 262 : selected_indices(nretained) = inew_point
309 262 : unselected_inside(inew_point) = .FALSE.
310 149334 : DO igrid_point = 1, SIZE(inside_voronoi)
311 149072 : IF (.NOT. unselected_inside(igrid_point)) CYCLE
312 : distance_sq = SUM((cluster_points(:, igrid_point) - &
313 176988 : cluster_points(:, inew_point))**2)
314 : nearest_distance_sq(igrid_point) = &
315 149334 : MIN(nearest_distance_sq(igrid_point), distance_sq)
316 : END DO
317 : END DO
318 5 : END SUBROUTINE complete_grid_by_maximin_distance
319 :
320 : ! **************************************************************************************************
321 : !> \brief Collect all cluster Lebedev points used in Eqs. (5) and (7).
322 : !> \param bs_env ...
323 : !> \param radial_lebedev_grids ...
324 : !> \param iatom ...
325 : !> \param cluster_atoms ...
326 : !> \param cluster_points ...
327 : ! **************************************************************************************************
328 45 : SUBROUTINE collect_cluster_lebedev_points(bs_env, radial_lebedev_grids, iatom, &
329 15 : cluster_atoms, cluster_points)
330 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
331 : TYPE(rirs_grid_type), INTENT(IN) :: radial_lebedev_grids(:)
332 : INTEGER, INTENT(IN) :: iatom, cluster_atoms(:)
333 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: cluster_points(:, :)
334 :
335 : INTEGER :: cluster_iatom, i_cluster_atom, ikind, n, &
336 : npoints, offset
337 :
338 15 : npoints = 0
339 60 : DO i_cluster_atom = 1, SIZE(cluster_atoms)
340 45 : ikind = bs_env%ri_rs%particle_set(cluster_atoms(i_cluster_atom))%atomic_kind%kind_number
341 60 : npoints = npoints + radial_lebedev_grids(ikind)%npts
342 : END DO
343 45 : ALLOCATE (cluster_points(3, npoints))
344 :
345 15 : offset = 0
346 60 : DO i_cluster_atom = 1, SIZE(cluster_atoms)
347 45 : cluster_iatom = cluster_atoms(i_cluster_atom)
348 45 : ikind = bs_env%ri_rs%particle_set(cluster_iatom)%atomic_kind%kind_number
349 45 : n = radial_lebedev_grids(ikind)%npts
350 : cluster_points(:, offset + 1:offset + n) = radial_lebedev_grids(ikind)%raw_points + &
351 : SPREAD(bs_env%ri_rs%particle_set(cluster_iatom)%r - &
352 18516 : bs_env%ri_rs%particle_set(iatom)%r, 2, n)
353 60 : offset = offset + n
354 : END DO
355 15 : END SUBROUTINE collect_cluster_lebedev_points
356 :
357 : ! **************************************************************************************************
358 : !> \brief Evaluate ϕ_μ(r_l) needed for D_ll' in Eq. (7), using AO locality. Atom b uses the
359 : !> point rows atom_point_offsets(b):atom_point_offsets(b+1)-1. Its column-major block
360 : !> ϕ_μ(r_l) has shape (number of point rows, atom_n_ao(b)) and starts at
361 : !> atom_value_offsets(b) in ao_values.
362 : !> \param bs_env ...
363 : !> \param iatom ...
364 : !> \param cluster_atoms ...
365 : !> \param cluster_points ...
366 : !> \param ao_point_indices Cluster-point index for every stored point row.
367 : !> \param atom_point_offsets First stored point row for each atom, followed by the final bound.
368 : !> \param atom_value_offsets First AO value for each atom, followed by the final bound.
369 : !> \param atom_n_ao Number of AO functions for each atom.
370 : !> \param ao_values Compact atom-local values ϕ_μ(r_l).
371 : ! **************************************************************************************************
372 15 : SUBROUTINE evaluate_cluster_ao_values(bs_env, iatom, cluster_atoms, cluster_points, &
373 : ao_point_indices, atom_point_offsets, atom_value_offsets, &
374 : atom_n_ao, ao_values)
375 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
376 : INTEGER, INTENT(IN) :: iatom, cluster_atoms(:)
377 : REAL(KIND=dp), INTENT(IN) :: cluster_points(:, :)
378 : INTEGER, ALLOCATABLE, INTENT(OUT) :: ao_point_indices(:), &
379 : atom_point_offsets(:), &
380 : atom_value_offsets(:), atom_n_ao(:)
381 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: ao_values(:)
382 :
383 : INTEGER :: cluster_iatom, first_point, first_value, i_cluster_atom, igrid_point, ikind, &
384 : last_point, last_value, nactive, npoints
385 : INTEGER, ALLOCATABLE :: point_indices(:)
386 : LOGICAL, ALLOCATABLE :: active(:)
387 : REAL(KIND=dp), ALLOCATABLE :: absolute_points(:, :), &
388 15 : atom_points(:, :), atom_values(:, :)
389 : TYPE(gto_basis_set_type), POINTER :: ao
390 :
391 15 : npoints = SIZE(cluster_points, 2)
392 45 : ALLOCATE (absolute_points(3, npoints))
393 : absolute_points(:, :) = cluster_points + &
394 18366 : SPREAD(bs_env%ri_rs%particle_set(iatom)%r, 2, npoints)
395 45 : ALLOCATE (point_indices(npoints))
396 9198 : point_indices(:) = [(igrid_point, igrid_point=1, npoints)]
397 0 : ALLOCATE (active(npoints), atom_n_ao(SIZE(cluster_atoms)), &
398 0 : atom_point_offsets(SIZE(cluster_atoms) + 1), &
399 105 : atom_value_offsets(SIZE(cluster_atoms) + 1))
400 :
401 15 : atom_point_offsets(1) = 1
402 15 : atom_value_offsets(1) = 1
403 60 : DO i_cluster_atom = 1, SIZE(cluster_atoms)
404 45 : cluster_iatom = cluster_atoms(i_cluster_atom)
405 45 : ikind = bs_env%ri_rs%particle_set(cluster_iatom)%atomic_kind%kind_number
406 45 : ao => bs_env%basis_set_AO(ikind)%gto_basis_set
407 45 : CPASSERT(ao%kind_radius > 0.0_dp)
408 : active(:) = SUM((absolute_points - &
409 : SPREAD(bs_env%ri_rs%particle_set(cluster_iatom)%r, 2, npoints))**2, DIM=1) &
410 55053 : <= ao%kind_radius**2
411 13797 : nactive = COUNT(active)
412 45 : atom_n_ao(i_cluster_atom) = ao%nsgf
413 45 : atom_point_offsets(i_cluster_atom + 1) = atom_point_offsets(i_cluster_atom) + nactive
414 : atom_value_offsets(i_cluster_atom + 1) = atom_value_offsets(i_cluster_atom) + &
415 60 : nactive*atom_n_ao(i_cluster_atom)
416 : END DO
417 :
418 0 : ALLOCATE (ao_point_indices(atom_point_offsets(SIZE(cluster_atoms) + 1) - 1), &
419 75 : ao_values(atom_value_offsets(SIZE(cluster_atoms) + 1) - 1))
420 60 : DO i_cluster_atom = 1, SIZE(cluster_atoms)
421 45 : cluster_iatom = cluster_atoms(i_cluster_atom)
422 45 : ikind = bs_env%ri_rs%particle_set(cluster_iatom)%atomic_kind%kind_number
423 45 : ao => bs_env%basis_set_AO(ikind)%gto_basis_set
424 : active(:) = SUM((absolute_points - &
425 : SPREAD(bs_env%ri_rs%particle_set(cluster_iatom)%r, 2, npoints))**2, DIM=1) &
426 55053 : <= ao%kind_radius**2
427 45 : first_point = atom_point_offsets(i_cluster_atom)
428 45 : last_point = atom_point_offsets(i_cluster_atom + 1) - 1
429 45 : first_value = atom_value_offsets(i_cluster_atom)
430 45 : last_value = atom_value_offsets(i_cluster_atom + 1) - 1
431 45 : nactive = last_point - first_point + 1
432 45 : ao_point_indices(first_point:last_point) = PACK(point_indices, active)
433 270 : ALLOCATE (atom_points(3, nactive), atom_values(nactive, ao%nsgf))
434 135 : atom_points(:, :) = RESHAPE(PACK(absolute_points, SPREAD(active, 1, 3)), [3, nactive])
435 45 : atom_values(:, :) = 0.0_dp
436 : CALL evaluate_ao_basis_on_points(atom_values, atom_points, ao, &
437 45 : bs_env%ri_rs%particle_set(cluster_iatom)%r, bs_env%ri_rs%cell)
438 90 : ao_values(first_value:last_value) = RESHAPE(atom_values, [nactive*ao%nsgf])
439 60 : DEALLOCATE (atom_points, atom_values)
440 : END DO
441 15 : END SUBROUTINE evaluate_cluster_ao_values
442 :
443 : ! **************************************************************************************************
444 : !> \brief Apply the Cholesky selection of Eqs. (7)-(10) and retain selected points in V_A.
445 : !> \param ao_point_indices Cluster-point index for every stored point row.
446 : !> \param atom_point_offsets First stored point row for each atom, followed by the final bound.
447 : !> \param atom_value_offsets First AO value for each atom, followed by the final bound.
448 : !> \param atom_n_ao Number of AO functions for each atom.
449 : !> \param ao_values Compact atom-local values ϕ_μ(r_l).
450 : !> \param npoints Number of points in the cluster Lebedev grids.
451 : !> \param inside_voronoi True for points inside the molecular Voronoi cell V_A.
452 : !> \param selected_indices Selected indices in V_A; zero denotes rank exhaustion.
453 : ! **************************************************************************************************
454 15 : SUBROUTINE select_cholesky_grid_points(ao_point_indices, atom_point_offsets, atom_value_offsets, &
455 30 : atom_n_ao, ao_values, npoints, inside_voronoi, &
456 15 : selected_indices)
457 : INTEGER, INTENT(IN) :: ao_point_indices(:), &
458 : atom_point_offsets(:), &
459 : atom_value_offsets(:), atom_n_ao(:)
460 : REAL(KIND=dp), CONTIGUOUS, INTENT(IN), TARGET :: ao_values(:)
461 : INTEGER, INTENT(IN) :: npoints
462 : LOGICAL, INTENT(IN) :: inside_voronoi(npoints)
463 : INTEGER, INTENT(OUT) :: selected_indices(:)
464 :
465 : INTEGER :: capacity, first_point, first_value, iblock, igrid_point, irow, last_point, &
466 : last_value, n_ao, new_capacity, nretained, nrows, nselected, selected_point
467 : REAL(KIND=dp) :: selected_diagonal
468 15 : REAL(KIND=dp), ALLOCATABLE :: diagonal(:), factor(:, :), grown(:, :), initial_diagonal(:), &
469 15 : selected_ao_values(:), selected_factor_values(:), selection_matrix_column(:), values(:)
470 15 : REAL(KIND=dp), POINTER :: atom_ao_values(:, :)
471 :
472 : ALLOCATE (diagonal(npoints), initial_diagonal(npoints), &
473 90 : selection_matrix_column(npoints), values(npoints))
474 :
475 : ! d_l^(0) = D_ll = [Σ_μ ϕ_μ(r_l)^2]^2.
476 15 : diagonal(:) = 0.0_dp
477 60 : DO iblock = 1, SIZE(atom_n_ao)
478 45 : first_point = atom_point_offsets(iblock)
479 45 : last_point = atom_point_offsets(iblock + 1) - 1
480 45 : first_value = atom_value_offsets(iblock)
481 45 : last_value = atom_value_offsets(iblock + 1) - 1
482 45 : nrows = last_point - first_point + 1
483 45 : atom_ao_values(1:nrows, 1:atom_n_ao(iblock)) => ao_values(first_value:last_value)
484 12894 : DO irow = 1, nrows
485 12834 : igrid_point = ao_point_indices(first_point + irow - 1)
486 41961 : diagonal(igrid_point) = diagonal(igrid_point) + SUM(atom_ao_values(irow, :)**2)
487 : END DO
488 : END DO
489 4599 : diagonal(:) = diagonal**2
490 4599 : initial_diagonal(:) = diagonal
491 :
492 15 : capacity = MIN(32, npoints)
493 : ALLOCATE (factor(npoints, capacity), selected_ao_values(MAXVAL(atom_n_ao)), &
494 150 : selected_factor_values(npoints))
495 396 : selected_indices(:) = 0
496 15 : nselected = 0
497 15 : nretained = 0
498 373 : DO WHILE (nretained < SIZE(selected_indices) .AND. nselected < npoints)
499 119275 : selected_point = MAXLOC(diagonal, DIM=1)
500 363 : selected_diagonal = diagonal(selected_point)
501 363 : IF (selected_diagonal <= 0.0_dp) EXIT
502 :
503 358 : IF (nselected == capacity) THEN
504 0 : new_capacity = MIN(npoints, capacity + MAX(32, capacity/2))
505 0 : ALLOCATE (grown(npoints, new_capacity))
506 0 : grown(:, :capacity) = factor(:, :capacity)
507 0 : CALL MOVE_ALLOC(grown, factor)
508 0 : capacity = new_capacity
509 : END IF
510 :
511 : ! Form D_lq from atom-local AO products without storing the dense D_ll' matrix.
512 358 : selection_matrix_column(:) = 0.0_dp
513 1432 : DO iblock = 1, SIZE(atom_n_ao)
514 1074 : first_point = atom_point_offsets(iblock)
515 1074 : last_point = atom_point_offsets(iblock + 1) - 1
516 1074 : nrows = last_point - first_point + 1
517 1074 : IF (nrows == 0) CYCLE
518 1074 : irow = locate(ao_point_indices(first_point:last_point), selected_point)
519 1074 : IF (irow == 0) CYCLE
520 1074 : first_value = atom_value_offsets(iblock)
521 1074 : last_value = atom_value_offsets(iblock + 1) - 1
522 1074 : n_ao = atom_n_ao(iblock)
523 1074 : atom_ao_values(1:nrows, 1:n_ao) => ao_values(first_value:last_value)
524 3580 : selected_ao_values(:n_ao) = atom_ao_values(irow, :)
525 : CALL dgemv('N', nrows, n_ao, 1.0_dp, atom_ao_values, nrows, &
526 1074 : selected_ao_values, 1, 0.0_dp, values, 1)
527 326784 : DO irow = 1, nrows
528 325352 : igrid_point = ao_point_indices(first_point + irow - 1)
529 : selection_matrix_column(igrid_point) = &
530 326426 : selection_matrix_column(igrid_point) + values(irow)
531 : END DO
532 : END DO
533 117006 : selection_matrix_column(:) = selection_matrix_column**2
534 :
535 : ! L_lk = [D_lq_k - sum_(j<k) L_lj L_q_kj]/sqrt(d_q_k^(k-1)).
536 358 : IF (nselected > 0) THEN
537 4710 : selected_factor_values(:nselected) = factor(selected_point, :nselected)
538 : CALL dgemv('N', npoints, nselected, -1.0_dp, factor, npoints, &
539 343 : selected_factor_values, 1, 1.0_dp, selection_matrix_column, 1)
540 : END IF
541 358 : nselected = nselected + 1
542 117006 : factor(:, nselected) = selection_matrix_column/SQRT(selected_diagonal)
543 117006 : diagonal(:) = MAX(0.0_dp, diagonal - factor(:, nselected)**2)
544 117006 : WHERE (diagonal <= 64.0_dp*EPSILON(1.0_dp)*initial_diagonal) diagonal = 0.0_dp
545 358 : diagonal(selected_point) = 0.0_dp
546 :
547 363 : IF (inside_voronoi(selected_point)) THEN
548 119 : nretained = nretained + 1
549 119 : selected_indices(nretained) = selected_point
550 : END IF
551 : END DO
552 15 : END SUBROUTINE select_cholesky_grid_points
553 :
554 : ! **************************************************************************************************
555 : !> \brief Communicate the independently constructed atom grids after Eqs. (5)-(12).
556 : !> \param bs_env ...
557 : ! **************************************************************************************************
558 10 : SUBROUTINE broadcast_ri_rs_grids(bs_env)
559 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
560 :
561 : INTEGER :: atom_mpi_rank, iatom
562 :
563 40 : DO iatom = 1, bs_env%n_atom
564 30 : atom_mpi_rank = MOD(iatom - 1, bs_env%para_env%num_pe)
565 40 : CALL bs_env%para_env%bcast(bs_env%ri_rs%atomic_grids(iatom)%raw_points, atom_mpi_rank)
566 : END DO
567 10 : END SUBROUTINE broadcast_ri_rs_grids
568 :
569 : END MODULE gw_ri_rs_grid_initialization
|
