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 Optimize automatically initialized atom-centred RI-RS grids.
10 : !> \par History
11 : !> 09.2026 created [Jan Wilhelm]
12 : ! **************************************************************************************************
13 : MODULE gw_ri_rs_grid_optimization
14 : USE basis_set_types, ONLY: gto_basis_set_type
15 : USE cell_types, ONLY: cell_type
16 : USE cp_lbfgs, ONLY: setulb
17 : USE gw_ri_rs_grid_initialization, ONLY: initialize_ri_rs_grid
18 : USE gw_ri_rs_utils, ONLY: evaluate_ao_basis_on_points,&
19 : filter_grid_to_voronoi,&
20 : get_rirs_cluster_atoms
21 : USE gw_utils_compute_integrals, ONLY: build_3c_integral_block_ctx,&
22 : gw_3c_ctx_create,&
23 : gw_3c_ctx_release,&
24 : gw_3c_ctx_type,&
25 : gw_3c_ws_create,&
26 : gw_3c_ws_release,&
27 : gw_3c_ws_type
28 : USE kinds, ONLY: dp
29 : USE machine, ONLY: m_flush,&
30 : m_walltime
31 : USE message_passing, ONLY: mp_para_env_type
32 : USE particle_types, ONLY: particle_type
33 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type,&
34 : rirs_grid_type
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 : PRIVATE
39 :
40 : PUBLIC :: optimize_ri_rs_grid
41 :
42 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_ri_rs_grid_optimization'
43 :
44 : TYPE, PRIVATE :: local_cluster_3c_integrals_type
45 : INTEGER, ALLOCATABLE :: atom_indices(:)
46 : REAL(KIND=dp), ALLOCATABLE :: Int_3c(:, :, :)
47 : END TYPE local_cluster_3c_integrals_type
48 :
49 : CONTAINS
50 :
51 : ! **************************************************************************************************
52 : !> \brief Initialize Lebedev grids and subsequently optimize their grid-point coordinates.
53 : !> \param bs_env ...
54 : ! **************************************************************************************************
55 10 : SUBROUTINE optimize_ri_rs_grid(bs_env)
56 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
57 :
58 : CHARACTER(len=*), PARAMETER :: routineN = 'optimize_ri_rs_grid'
59 :
60 : INTEGER :: handle
61 10 : TYPE(local_cluster_3c_integrals_type), ALLOCATABLE :: cluster_3c_int(:)
62 :
63 10 : CALL timeset(routineN, handle)
64 :
65 : ! Validate all state and user input before the initialization consumes them.
66 10 : CALL validate_ri_rs_grid_optimization_input(bs_env)
67 :
68 : ! setup Lebedev initial grid inside Voronoi volume, select points with Cholesky decomp.
69 10 : CALL initialize_ri_rs_grid(bs_env)
70 :
71 : ! Validate the initialized RI-RS state and report the start of the optimization.
72 10 : CALL prepare_ri_rs_grid_optimization(bs_env)
73 :
74 : ! C_A = {B: |R_A - R_B| < R_cut}; store (μν|P) for μ, ν, and P centered in C_A.
75 10 : CALL build_cluster_3c_int(bs_env, cluster_3c_int)
76 :
77 : ! E_A = Σ_{μνP∈C_A} [(μν|P) - Σ_{l∈G_A} ϕ_μ(r_l) ϕ_ν(r_l) Z_lP^(A)]².
78 : ! Minimize E_loc(norm) = N_atom^(-1) Σ_A [E_A / Σ_{μνP∈C_A} |(μν|P)|²].
79 10 : CALL optimize_grid_coordinates(bs_env, cluster_3c_int)
80 :
81 10 : CALL timestop(handle)
82 :
83 25 : END SUBROUTINE optimize_ri_rs_grid
84 :
85 : ! **************************************************************************************************
86 : !> \brief Validate state and user input needed to initialize and optimize RI-RS grids.
87 : !> \param bs_env ...
88 : ! **************************************************************************************************
89 10 : SUBROUTINE validate_ri_rs_grid_optimization_input(bs_env)
90 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
91 :
92 : INTEGER :: iatom, ikind
93 :
94 10 : CPASSERT(ASSOCIATED(bs_env%ri_rs%cell))
95 10 : CPASSERT(ASSOCIATED(bs_env%ri_rs%particle_set))
96 10 : CPASSERT(ASSOCIATED(bs_env%para_env))
97 10 : CPASSERT(ALLOCATED(bs_env%basis_set_AO))
98 10 : CPASSERT(ALLOCATED(bs_env%basis_set_RI))
99 10 : CPASSERT(ALLOCATED(bs_env%sizes_AO))
100 10 : CPASSERT(SIZE(bs_env%ri_rs%particle_set) == bs_env%n_atom)
101 10 : CPASSERT(SIZE(bs_env%sizes_AO) == bs_env%n_atom)
102 10 : CPASSERT(SIZE(bs_env%basis_set_AO) == SIZE(bs_env%basis_set_RI))
103 :
104 10 : IF (bs_env%ri_rs%grid_opt%rs_ao_ratio <= 0.0_dp) THEN
105 0 : CPABORT("GRID_OPTIMIZATION%RS_AO_RATIO must be positive")
106 : END IF
107 10 : IF (bs_env%ri_rs%grid_opt%cutoff_atomic_cluster <= 0.0_dp) THEN
108 0 : CPABORT("GRID_OPTIMIZATION%CUTOFF_ATOMIC_CLUSTER must be positive")
109 : END IF
110 10 : IF (bs_env%ri_rs%grid_opt%max_iter < 1) THEN
111 0 : CPABORT("GRID_OPTIMIZATION%MAX_ITER must be positive")
112 : END IF
113 10 : IF (bs_env%ri_rs%tikhonov < 0.0_dp) THEN
114 0 : CPABORT("RI_RS%TIKHONOV must not be negative")
115 : END IF
116 :
117 40 : DO iatom = 1, bs_env%n_atom
118 30 : ikind = bs_env%ri_rs%particle_set(iatom)%atomic_kind%kind_number
119 30 : CPASSERT(ikind >= 1 .AND. ikind <= SIZE(bs_env%basis_set_AO))
120 30 : CPASSERT(ASSOCIATED(bs_env%basis_set_AO(ikind)%gto_basis_set))
121 30 : CPASSERT(ASSOCIATED(bs_env%basis_set_RI(ikind)%gto_basis_set))
122 30 : IF (bs_env%sizes_AO(iatom) < 1 .OR. &
123 10 : bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf < 1) THEN
124 0 : CPABORT("Every atom in GRID_OPTIMIZATION needs ORB and RI_AUX functions")
125 : END IF
126 : END DO
127 10 : END SUBROUTINE validate_ri_rs_grid_optimization_input
128 :
129 : ! **************************************************************************************************
130 : !> \brief Validate initialized RI-RS state and announce the grid optimization.
131 : !> \param bs_env ...
132 : ! **************************************************************************************************
133 10 : SUBROUTINE prepare_ri_rs_grid_optimization(bs_env)
134 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
135 :
136 : INTEGER :: iatom, unit_nr
137 :
138 10 : CPASSERT(ALLOCATED(bs_env%ri_rs%atomic_grids))
139 10 : CPASSERT(SIZE(bs_env%ri_rs%atomic_grids) == bs_env%n_atom)
140 :
141 40 : DO iatom = 1, bs_env%n_atom
142 30 : CPASSERT(ALLOCATED(bs_env%ri_rs%atomic_grids(iatom)%raw_points))
143 30 : CPASSERT(SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 1) == 3)
144 40 : CPASSERT(SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 2) > 0)
145 : END DO
146 :
147 10 : unit_nr = bs_env%unit_nr
148 10 : IF (unit_nr > 0) THEN
149 5 : WRITE (unit_nr, '(T2,A)') 'Started RI-RS grid optimization'
150 5 : CALL m_flush(unit_nr)
151 : END IF
152 10 : END SUBROUTINE prepare_ri_rs_grid_optimization
153 :
154 : ! **************************************************************************************************
155 : !> \brief Evaluate the normalized local objective E_loc and its gradient for one L-BFGS trial vector.
156 : !> \param grid_coordinates Flattened atom-relative grid coordinates.
157 : !> \param atom_coordinate_offsets Starting coordinate offset for each atom.
158 : !> \param bs_env ...
159 : !> \param cluster_3c_int Rank-local cluster three-centre integrals.
160 : !> \param normalized_error Mean normalized squared three-centre-integral error.
161 : !> \param coordinate_gradient Derivative of normalized_error with respect to coordinates.
162 : !> \param max_abs_error Largest absolute three-centre-integral error.
163 : !> \param fit_successful Whether every local-cluster evaluation succeeded.
164 : ! **************************************************************************************************
165 256 : SUBROUTINE grid_objective(grid_coordinates, atom_coordinate_offsets, bs_env, cluster_3c_int, &
166 256 : normalized_error, coordinate_gradient, &
167 : max_abs_error, fit_successful)
168 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: grid_coordinates
169 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_coordinate_offsets
170 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
171 : TYPE(local_cluster_3c_integrals_type), &
172 : DIMENSION(:), INTENT(IN) :: cluster_3c_int
173 : REAL(KIND=dp), INTENT(OUT) :: normalized_error
174 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: coordinate_gradient
175 : REAL(KIND=dp), INTENT(OUT) :: max_abs_error
176 : LOGICAL, INTENT(OUT) :: fit_successful
177 :
178 : INTEGER :: evaluated_cluster_count, icluster, &
179 : successful_ranks
180 : LOGICAL :: cluster_fit_successful
181 : REAL(KIND=dp) :: cluster_max_abs_error, &
182 : cluster_normalized_error
183 256 : REAL(KIND=dp), ALLOCATABLE :: cluster_coordinate_gradient(:, :)
184 : TYPE(mp_para_env_type), POINTER :: para_env
185 :
186 256 : para_env => bs_env%para_env
187 256 : CPASSERT(SIZE(coordinate_gradient) == SIZE(grid_coordinates))
188 : CALL unpack_lbfgs_grid_coordinates(grid_coordinates, atom_coordinate_offsets, &
189 256 : bs_env%ri_rs%atomic_grids)
190 256 : normalized_error = 0.0_dp
191 20470 : coordinate_gradient = 0.0_dp
192 256 : evaluated_cluster_count = 0
193 256 : max_abs_error = 0.0_dp
194 256 : fit_successful = .TRUE.
195 640 : DO icluster = 1, SIZE(cluster_3c_int)
196 : CALL evaluate_local_cluster(cluster_3c_int(icluster), bs_env, &
197 : cluster_normalized_error, cluster_coordinate_gradient, &
198 384 : cluster_max_abs_error, cluster_fit_successful)
199 384 : IF (.NOT. cluster_fit_successful) THEN
200 0 : fit_successful = .FALSE.
201 0 : IF (ALLOCATED(cluster_coordinate_gradient)) DEALLOCATE (cluster_coordinate_gradient)
202 : EXIT
203 : END IF
204 384 : normalized_error = normalized_error + cluster_normalized_error
205 : CALL accumulate_atom_gradient(cluster_3c_int(icluster)%atom_indices, &
206 : bs_env%ri_rs%atomic_grids, &
207 : atom_coordinate_offsets, cluster_coordinate_gradient, &
208 384 : coordinate_gradient)
209 384 : evaluated_cluster_count = evaluated_cluster_count + 1
210 : max_abs_error = &
211 384 : MAX(max_abs_error, cluster_max_abs_error)
212 640 : DEALLOCATE (cluster_coordinate_gradient)
213 : END DO
214 256 : successful_ranks = MERGE(1, 0, fit_successful)
215 256 : CALL para_env%sum(normalized_error)
216 40684 : CALL para_env%sum(coordinate_gradient)
217 256 : CALL para_env%sum(evaluated_cluster_count)
218 256 : CALL para_env%sum(successful_ranks)
219 256 : CALL para_env%max(max_abs_error)
220 : fit_successful = evaluated_cluster_count == bs_env%n_atom .AND. &
221 256 : successful_ranks == para_env%num_pe
222 256 : IF (.NOT. fit_successful) THEN
223 0 : normalized_error = HUGE(normalized_error)
224 0 : coordinate_gradient = 0.0_dp
225 0 : RETURN
226 : END IF
227 : ! E_loc(norm) = (1/N_atom) Σ_A [E_A / Σ_{μνP ∈ C_A} |(μν|P)|²].
228 256 : normalized_error = normalized_error/REAL(bs_env%n_atom, dp)
229 20470 : coordinate_gradient = coordinate_gradient/REAL(bs_env%n_atom, dp)
230 256 : END SUBROUTINE grid_objective
231 :
232 : ! **************************************************************************************************
233 : !> \brief Drive L-BFGS, repeatedly evaluating E_loc and its gradient, and return the best trial.
234 : !> \param bs_env ...
235 : !> \param cluster_3c_int Rank-local cluster three-centre integrals.
236 : ! **************************************************************************************************
237 10 : SUBROUTINE optimize_grid_coordinates(bs_env, cluster_3c_int)
238 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
239 : TYPE(local_cluster_3c_integrals_type), &
240 : DIMENSION(:), INTENT(IN) :: cluster_3c_int
241 :
242 : INTEGER, PARAMETER :: lbfgs_history = 7
243 : REAL(KIND=dp), PARAMETER :: lbfgs_factr = 0.0_dp, &
244 : lbfgs_pgtol = 1.0E-9_dp
245 :
246 : CHARACTER(LEN=60) :: line_search_state, optimizer_task
247 : INTEGER :: evaluation_count, handle, iatom, &
248 : molecular_outside, npoints, unit_nr
249 10 : INTEGER, ALLOCATABLE :: atom_coordinate_offsets(:), &
250 10 : bound_types(:), integer_workspace(:)
251 : INTEGER, DIMENSION(44) :: integer_state
252 : LOGICAL :: evaluation_successful, fit_successful, &
253 : have_best
254 10 : LOGICAL, ALLOCATABLE :: inside(:)
255 : LOGICAL, DIMENSION(4) :: logical_state
256 : REAL(KIND=dp) :: best_normalized_error, max_abs_error, &
257 : normalized_error, start_time
258 10 : REAL(KIND=dp), ALLOCATABLE :: best_grid_coordinates(:), coordinate_gradient(:), &
259 10 : grid_coordinates(:), lower_bounds(:), points(:, :), upper_bounds(:), workspace(:)
260 : REAL(KIND=dp), DIMENSION(29) :: real_state
261 :
262 10 : CALL timeset("rirs_grid_LBFGS", handle)
263 10 : start_time = m_walltime()
264 :
265 : CALL pack_lbfgs_grid_coordinates(bs_env%ri_rs%atomic_grids, grid_coordinates, &
266 10 : atom_coordinate_offsets)
267 10 : CPASSERT(SIZE(grid_coordinates) > 0)
268 :
269 : ALLOCATE (best_grid_coordinates(SIZE(grid_coordinates)), &
270 : coordinate_gradient(SIZE(grid_coordinates)), &
271 : lower_bounds(SIZE(grid_coordinates)), &
272 : upper_bounds(SIZE(grid_coordinates)), &
273 : bound_types(SIZE(grid_coordinates)), &
274 : integer_workspace(3*SIZE(grid_coordinates)), &
275 : workspace(2*lbfgs_history*SIZE(grid_coordinates) + &
276 120 : 5*SIZE(grid_coordinates) + 11*lbfgs_history**2 + 8*lbfgs_history))
277 : ! L-BFGS-B ignores lower_bounds and upper_bounds when bound_types is zero.
278 10 : lower_bounds = 0.0_dp
279 10 : upper_bounds = 0.0_dp
280 10 : bound_types = 0
281 10 : optimizer_task = 'START'
282 10 : line_search_state = ''
283 10 : normalized_error = HUGE(normalized_error)
284 10 : coordinate_gradient = 0.0_dp
285 10 : workspace = 0.0_dp
286 10 : integer_workspace = 0
287 10 : logical_state = .FALSE.
288 10 : integer_state = 0
289 10 : real_state = 0.0_dp
290 10 : evaluation_count = 0
291 10 : have_best = .FALSE.
292 10 : best_normalized_error = HUGE(best_normalized_error)
293 2296 : best_grid_coordinates(:) = grid_coordinates
294 :
295 : DO
296 : CALL setulb(SIZE(grid_coordinates), lbfgs_history, grid_coordinates, &
297 : lower_bounds, upper_bounds, bound_types, &
298 : normalized_error, coordinate_gradient, lbfgs_factr, lbfgs_pgtol, &
299 : workspace, integer_workspace, optimizer_task, -1, line_search_state, &
300 444 : logical_state, integer_state, real_state, -1.0_dp)
301 454 : IF (optimizer_task(1:2) == 'FG') THEN
302 252 : IF (evaluation_count >= bs_env%ri_rs%grid_opt%max_iter) EXIT
303 : CALL grid_objective(grid_coordinates, atom_coordinate_offsets, bs_env, cluster_3c_int, &
304 : normalized_error, coordinate_gradient, &
305 246 : max_abs_error, evaluation_successful)
306 246 : evaluation_count = evaluation_count + 1
307 246 : IF (.NOT. evaluation_successful) EXIT
308 246 : IF (evaluation_count == 1 .AND. bs_env%unit_nr > 0) THEN
309 : WRITE (bs_env%unit_nr, '(T2,A,T61,ES20.12)') &
310 5 : 'RI-RS grid initial normalized error', normalized_error
311 5 : CALL m_flush(bs_env%unit_nr)
312 : END IF
313 246 : IF (.NOT. have_best .OR. normalized_error < best_normalized_error) THEN
314 214 : have_best = .TRUE.
315 214 : best_normalized_error = normalized_error
316 16182 : best_grid_coordinates(:) = grid_coordinates
317 : END IF
318 192 : ELSE IF (optimizer_task(1:5) == 'NEW_X') THEN
319 : CYCLE
320 : ELSE
321 188 : EXIT
322 : END IF
323 : END DO
324 2296 : IF (have_best) grid_coordinates(:) = best_grid_coordinates
325 : CALL grid_objective(grid_coordinates, atom_coordinate_offsets, bs_env, cluster_3c_int, &
326 : normalized_error, coordinate_gradient, &
327 10 : max_abs_error, fit_successful)
328 :
329 10 : IF (.NOT. fit_successful) CPABORT("RI-RS grid optimization produced no regularized fit")
330 :
331 40 : DO iatom = 1, SIZE(bs_env%ri_rs%atomic_grids)
332 : bs_env%ri_rs%atomic_grids(iatom)%npts = &
333 40 : SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 2)
334 : END DO
335 10 : bs_env%ri_rs%Z_lP_exists = .FALSE.
336 :
337 10 : unit_nr = bs_env%unit_nr
338 10 : IF (unit_nr > 0) THEN
339 5 : npoints = 0
340 5 : molecular_outside = 0
341 20 : DO iatom = 1, SIZE(bs_env%ri_rs%atomic_grids)
342 0 : ALLOCATE (points(3, bs_env%ri_rs%atomic_grids(iatom)%npts), &
343 75 : inside(bs_env%ri_rs%atomic_grids(iatom)%npts))
344 1539 : points(:, :) = bs_env%ri_rs%atomic_grids(iatom)%raw_points
345 15 : CALL filter_grid_to_voronoi(points, iatom, bs_env%ri_rs%particle_set, mask=inside)
346 396 : molecular_outside = molecular_outside + COUNT(.NOT. inside)
347 15 : npoints = npoints + SIZE(inside)
348 20 : DEALLOCATE (points, inside)
349 : END DO
350 : WRITE (unit_nr, '(T2,A,T69,F10.1,A)') &
351 5 : 'RI-RS grid optimization completed, execution time:', m_walltime() - start_time, ' s'
352 5 : WRITE (unit_nr, '(T2,A,T72,ES9.1)') 'Normalized 3C error:', normalized_error
353 5 : WRITE (unit_nr, '(T2,A,T72,ES9.1)') 'Maximum absolute 3C error:', max_abs_error
354 : WRITE (unit_nr, '(T2,A,T69,I12,A,I0)') &
355 5 : 'Optimized points outside molecular Voronoi cells:', molecular_outside, ' / ', npoints
356 5 : FLUSH (unit_nr)
357 : END IF
358 10 : CALL timestop(handle)
359 20 : END SUBROUTINE optimize_grid_coordinates
360 :
361 : ! **************************************************************************************************
362 : !> \brief Build each C_A and store its exact (μν|P); this routine owns the integral context.
363 : !> \param bs_env ...
364 : !> \param cluster_3c_int Rank-local cluster three-centre integrals.
365 : ! **************************************************************************************************
366 10 : SUBROUTINE build_cluster_3c_int(bs_env, cluster_3c_int)
367 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
368 : TYPE(local_cluster_3c_integrals_type), &
369 : ALLOCATABLE, INTENT(OUT) :: cluster_3c_int(:)
370 :
371 : INTEGER :: handle, i_cluster_atom_j, i_cluster_atom_k, i_cluster_atom_p, iatom, iatom_j, &
372 : iatom_k, iatom_p, icenter_atom, icluster_atom, ikind, ilocal_cluster, ip, &
373 : n_local_clusters, nAO_cluster, nRI_cluster, nRI_nonzero
374 20 : INTEGER, ALLOCATABLE :: AO_offset(:), RI_offset(:), &
375 10 : sizes_ref_RI(:)
376 : LOGICAL :: screened
377 : REAL(KIND=dp) :: cluster_radius
378 10 : REAL(KIND=dp), ALLOCATABLE :: Int_3c_nonzero(:, :, :)
379 : TYPE(cell_type), POINTER :: cell
380 130 : TYPE(gw_3c_ctx_type) :: integral_context
381 10 : TYPE(gw_3c_ws_type) :: workspace
382 : TYPE(mp_para_env_type), POINTER :: para_env
383 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
384 :
385 10 : CALL timeset("rirs_grid_build_clusters", handle)
386 :
387 10 : cell => bs_env%ri_rs%cell
388 10 : para_env => bs_env%para_env
389 10 : particle_set => bs_env%ri_rs%particle_set
390 10 : cluster_radius = bs_env%ri_rs%grid_opt%cutoff_atomic_cluster
391 30 : ALLOCATE (sizes_ref_RI(bs_env%n_atom))
392 40 : DO iatom = 1, bs_env%n_atom
393 30 : ikind = particle_set(iatom)%atomic_kind%kind_number
394 40 : sizes_ref_RI(iatom) = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
395 : END DO
396 : n_local_clusters = COUNT([(MOD(icenter_atom - 1, para_env%num_pe) == para_env%mepos, &
397 80 : icenter_atom=1, bs_env%n_atom)])
398 0 : ALLOCATE (cluster_3c_int(n_local_clusters), &
399 75 : AO_offset(bs_env%n_atom), RI_offset(bs_env%n_atom))
400 : CALL gw_3c_ctx_create(integral_context, bs_env, bs_env%ri_metric, &
401 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
402 10 : basis_i=bs_env%basis_set_RI)
403 10 : CALL gw_3c_ws_create(workspace, integral_context)
404 10 : ilocal_cluster = 0
405 40 : DO icenter_atom = 1, bs_env%n_atom
406 30 : IF (MOD(icenter_atom - 1, para_env%num_pe) /= para_env%mepos) CYCLE
407 15 : ilocal_cluster = ilocal_cluster + 1
408 : CALL get_rirs_cluster_atoms(particle_set, cell, icenter_atom, cluster_radius, &
409 15 : cluster_3c_int(ilocal_cluster)%atom_indices)
410 15 : AO_offset = 0
411 15 : RI_offset = 0
412 15 : nAO_cluster = 0
413 15 : nRI_cluster = 0
414 60 : DO icluster_atom = 1, SIZE(cluster_3c_int(ilocal_cluster)%atom_indices)
415 45 : iatom = cluster_3c_int(ilocal_cluster)%atom_indices(icluster_atom)
416 45 : AO_offset(iatom) = nAO_cluster
417 45 : RI_offset(iatom) = nRI_cluster
418 45 : nAO_cluster = nAO_cluster + bs_env%sizes_AO(iatom)
419 60 : nRI_cluster = nRI_cluster + sizes_ref_RI(iatom)
420 : END DO
421 0 : ALLOCATE (cluster_3c_int(ilocal_cluster)%Int_3c(nAO_cluster, nAO_cluster, nRI_cluster), &
422 92244 : source=0.0_dp)
423 60 : DO i_cluster_atom_p = 1, SIZE(cluster_3c_int(ilocal_cluster)%atom_indices)
424 45 : iatom_p = cluster_3c_int(ilocal_cluster)%atom_indices(i_cluster_atom_p)
425 195 : DO i_cluster_atom_k = 1, SIZE(cluster_3c_int(ilocal_cluster)%atom_indices)
426 135 : iatom_k = cluster_3c_int(ilocal_cluster)%atom_indices(i_cluster_atom_k)
427 585 : DO i_cluster_atom_j = 1, SIZE(cluster_3c_int(ilocal_cluster)%atom_indices)
428 405 : iatom_j = cluster_3c_int(ilocal_cluster)%atom_indices(i_cluster_atom_j)
429 : CALL build_3c_integral_block_ctx(cluster_3c_int(ilocal_cluster)%Int_3c, &
430 : integral_context, workspace, &
431 : atom_j=iatom_j, atom_k=iatom_k, atom_i=iatom_p, &
432 : j_offset=AO_offset(iatom_j), &
433 : k_offset=AO_offset(iatom_k), &
434 540 : i_offset=RI_offset(iatom_p), screened=screened)
435 : END DO
436 : END DO
437 : END DO
438 : ! Zero auxiliary columns contribute neither to the objective nor its derivatives.
439 : nRI_nonzero = COUNT([(ANY(cluster_3c_int(ilocal_cluster)%Int_3c(:, :, ip) /= 0.0_dp), &
440 3348 : ip=1, nRI_cluster)])
441 25 : IF (nRI_nonzero < nRI_cluster) THEN
442 0 : ALLOCATE (Int_3c_nonzero(nAO_cluster, nAO_cluster, nRI_nonzero))
443 0 : nRI_nonzero = 0
444 0 : DO ip = 1, nRI_cluster
445 0 : IF (.NOT. ANY(cluster_3c_int(ilocal_cluster)%Int_3c(:, :, ip) /= 0.0_dp)) CYCLE
446 0 : nRI_nonzero = nRI_nonzero + 1
447 0 : Int_3c_nonzero(:, :, nRI_nonzero) = cluster_3c_int(ilocal_cluster)%Int_3c(:, :, ip)
448 : END DO
449 0 : CALL MOVE_ALLOC(Int_3c_nonzero, cluster_3c_int(ilocal_cluster)%Int_3c)
450 : END IF
451 : END DO
452 10 : CALL gw_3c_ws_release(workspace)
453 10 : CALL gw_3c_ctx_release(integral_context)
454 10 : CALL timestop(handle)
455 20 : END SUBROUTINE build_cluster_3c_int
456 :
457 : ! **************************************************************************************************
458 : !> \brief Evaluate ϕ_μ(r_l) on G_A and compute the normalized local error E_A and its gradient.
459 : !> \param cluster_3c_int Local atom indices and exact three-centre integrals.
460 : !> \param bs_env ...
461 : !> \param normalized_error Normalized squared three-centre-integral error.
462 : !> \param coordinate_gradient Derivative with respect to the cluster's physical grid coordinates.
463 : !> \param max_abs_error Largest absolute three-centre-integral error.
464 : !> \param fit_successful Whether the regularized fitting equations were solved.
465 : ! **************************************************************************************************
466 384 : SUBROUTINE evaluate_local_cluster(cluster_3c_int, bs_env, normalized_error, coordinate_gradient, &
467 : max_abs_error, fit_successful)
468 : TYPE(local_cluster_3c_integrals_type), INTENT(IN) :: cluster_3c_int
469 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
470 : REAL(KIND=dp), INTENT(OUT) :: normalized_error
471 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: coordinate_gradient(:, :)
472 : REAL(KIND=dp), INTENT(OUT) :: max_abs_error
473 : LOGICAL, INTENT(OUT) :: fit_successful
474 :
475 : INTEGER :: AO_offset, grid_point_offset, iatom, &
476 : icluster_atom, ikind, n_grid_points, &
477 : nAO, nAO_atom
478 : REAL(KIND=dp), ALLOCATABLE :: dPhi_alpha_l_mu(:, :, :), &
479 : grid_points(:, :), Phi_l_mu(:, :)
480 : TYPE(cell_type), POINTER :: cell
481 : TYPE(gto_basis_set_type), POINTER :: AO_basis
482 384 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
483 :
484 384 : cell => bs_env%ri_rs%cell
485 384 : particle_set => bs_env%ri_rs%particle_set
486 384 : nAO = SIZE(cluster_3c_int%Int_3c, 1)
487 384 : n_grid_points = 0
488 1536 : DO icluster_atom = 1, SIZE(cluster_3c_int%atom_indices)
489 1152 : iatom = cluster_3c_int%atom_indices(icluster_atom)
490 1536 : n_grid_points = n_grid_points + SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 2)
491 : END DO
492 : ALLOCATE (grid_points(3, n_grid_points), Phi_l_mu(n_grid_points, nAO), &
493 3840 : dPhi_alpha_l_mu(3, n_grid_points, nAO), coordinate_gradient(3, n_grid_points))
494 384 : Phi_l_mu = 0.0_dp
495 384 : dPhi_alpha_l_mu = 0.0_dp
496 384 : grid_point_offset = 0
497 1536 : DO icluster_atom = 1, SIZE(cluster_3c_int%atom_indices)
498 1152 : iatom = cluster_3c_int%atom_indices(icluster_atom)
499 : grid_points(:, grid_point_offset + 1:grid_point_offset + &
500 : SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 2)) = &
501 : SPREAD(particle_set(iatom)%r, 2, &
502 : SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 2)) + &
503 41580 : bs_env%ri_rs%atomic_grids(iatom)%raw_points
504 1536 : grid_point_offset = grid_point_offset + SIZE(bs_env%ri_rs%atomic_grids(iatom)%raw_points, 2)
505 : END DO
506 : AO_offset = 0
507 1536 : DO icluster_atom = 1, SIZE(cluster_3c_int%atom_indices)
508 1152 : iatom = cluster_3c_int%atom_indices(icluster_atom)
509 1152 : ikind = particle_set(iatom)%atomic_kind%kind_number
510 1152 : AO_basis => bs_env%basis_set_AO(ikind)%gto_basis_set
511 1152 : nAO_atom = AO_basis%nsgf
512 : CALL evaluate_ao_basis_on_points(Phi_l_mu(:, AO_offset + 1:AO_offset + nAO_atom), &
513 : grid_points, AO_basis, particle_set(iatom)%r, cell, &
514 1152 : dphi=dPhi_alpha_l_mu(:, :, AO_offset + 1:AO_offset + nAO_atom))
515 1536 : AO_offset = AO_offset + nAO_atom
516 : END DO
517 : CALL evaluate_rirs_grid_cluster(Phi_l_mu, dPhi_alpha_l_mu, cluster_3c_int%Int_3c, &
518 : bs_env%ri_rs%tikhonov, normalized_error, coordinate_gradient, &
519 384 : max_abs_error, fit_successful)
520 768 : END SUBROUTINE evaluate_local_cluster
521 :
522 : ! **************************************************************************************************
523 : !> \brief Solve the cluster-local Z_lP for one C_A and evaluate its normalized error and exact gradient.
524 : !>
525 : !> The solve applies the same column Jacobi scaling and Tikhonov parameter as the production
526 : !> Z_lP construction. Z_prime_lP denotes the coefficients before undoing the Jacobi scaling.
527 : !> \param Phi_l_mu AO values ϕ_μ(r_l), indexed (l, μ).
528 : !> \param dPhi_alpha_l_mu Cartesian derivatives of ϕ_μ(r_l), indexed (α, l, μ).
529 : !> \param Int_3c Exact three-centre integrals, indexed (mu,nu,P).
530 : !> \param tikhonov Tikhonov parameter used by the production RI-RS solve.
531 : !> \param normalized_error Normalized squared residual for this cluster.
532 : !> \param coordinate_gradient Analytic derivative of normalized_error with respect to r(alpha,l).
533 : !> \param max_abs_error Largest absolute error in an unweighted three-centre integral.
534 : !> \param fit_successful False if the regularized normal equations cannot be solved.
535 : ! **************************************************************************************************
536 384 : SUBROUTINE evaluate_rirs_grid_cluster(Phi_l_mu, dPhi_alpha_l_mu, Int_3c, tikhonov, &
537 384 : normalized_error, coordinate_gradient, &
538 : max_abs_error, fit_successful)
539 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Phi_l_mu
540 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dPhi_alpha_l_mu, Int_3c
541 : REAL(KIND=dp), INTENT(IN) :: tikhonov
542 : REAL(KIND=dp), INTENT(OUT) :: normalized_error
543 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: coordinate_gradient
544 : REAL(KIND=dp), INTENT(OUT) :: max_abs_error
545 : LOGICAL, INTENT(OUT) :: fit_successful
546 :
547 : CHARACTER(len=*), PARAMETER :: routineN = 'evaluate_rirs_grid_cluster'
548 : REAL(KIND=dp), PARAMETER :: jacobi_floor = 1.0E-16_dp
549 :
550 : INTEGER :: handle, i_ao_pair, icartesian, &
551 : igrid_point, imu, inu, ip, &
552 : lapack_info, n_AO_pairs, &
553 : n_grid_points, nAO, nRI, phase_handle
554 384 : INTEGER, ALLOCATABLE :: mu_of_pair(:), nu_of_pair(:)
555 : REAL(KIND=dp) :: absolute_error, dPhi_munu_l, &
556 : Int_3c_norm2, jacobi_scale, &
557 : Phi_munu_norm2, scaling_projection, &
558 : symmetry_factor
559 384 : REAL(KIND=dp), ALLOCATABLE :: D_inverse_Z_prime_lP(:, :), &
560 384 : D_inverse_Z_prime_times_Z_prime_transpose_ll(:, :), D_ll(:, :), d_lP(:, :), &
561 384 : dE_dPhi_munu_l(:, :), dPhi_munu_l_scaled(:), Int_3c_munu_P(:, :), jacobi_scaling(:), &
562 384 : Phi_munu_l(:, :), Phi_munu_l_scaled(:, :), R_munu_P(:, :), Z_prime_lP(:, :)
563 :
564 384 : CALL timeset(routineN, handle)
565 :
566 384 : n_grid_points = SIZE(Phi_l_mu, 1)
567 384 : nAO = SIZE(Phi_l_mu, 2)
568 384 : nRI = SIZE(Int_3c, 3)
569 384 : n_AO_pairs = nAO*(nAO + 1)/2
570 :
571 384 : CPASSERT(SIZE(dPhi_alpha_l_mu, 1) == 3)
572 384 : CPASSERT(SIZE(dPhi_alpha_l_mu, 2) == n_grid_points)
573 384 : CPASSERT(SIZE(dPhi_alpha_l_mu, 3) == nAO)
574 384 : CPASSERT(SIZE(Int_3c, 1) == nAO)
575 384 : CPASSERT(SIZE(Int_3c, 2) == nAO)
576 384 : CPASSERT(SIZE(coordinate_gradient, 1) == 3)
577 384 : CPASSERT(SIZE(coordinate_gradient, 2) == n_grid_points)
578 384 : CPASSERT(tikhonov >= 0.0_dp)
579 :
580 384 : normalized_error = HUGE(normalized_error)
581 40812 : coordinate_gradient = 0.0_dp
582 384 : max_abs_error = HUGE(max_abs_error)
583 384 : fit_successful = .FALSE.
584 384 : IF (n_grid_points < 1 .OR. nRI < 1) THEN
585 0 : CALL timestop(handle)
586 0 : RETURN
587 : END IF
588 :
589 : ALLOCATE (Phi_munu_l(n_AO_pairs, n_grid_points), Phi_munu_l_scaled(n_AO_pairs, n_grid_points), &
590 : Int_3c_munu_P(n_AO_pairs, nRI), mu_of_pair(n_AO_pairs), nu_of_pair(n_AO_pairs), &
591 5376 : jacobi_scaling(n_grid_points))
592 384 : CALL timeset("rirs_cluster_AO_products", phase_handle)
593 384 : i_ao_pair = 0
594 3072 : DO inu = 1, nAO
595 13824 : DO imu = 1, inu
596 10752 : i_ao_pair = i_ao_pair + 1
597 10752 : mu_of_pair(i_ao_pair) = imu
598 10752 : nu_of_pair(i_ao_pair) = inu
599 10752 : symmetry_factor = SQRT(REAL(2 - MERGE(1, 0, imu == inu), dp))
600 293748 : DO igrid_point = 1, n_grid_points
601 : Phi_munu_l(i_ao_pair, igrid_point) = &
602 293748 : symmetry_factor*Phi_l_mu(igrid_point, imu)*Phi_l_mu(igrid_point, inu)
603 : END DO
604 2280768 : DO ip = 1, nRI
605 2278080 : Int_3c_munu_P(i_ao_pair, ip) = symmetry_factor*Int_3c(imu, inu, ip)
606 : END DO
607 : END DO
608 : END DO
609 :
610 2348688 : Int_3c_norm2 = SUM(Int_3c_munu_P*Int_3c_munu_P)
611 384 : CALL timestop(phase_handle)
612 384 : IF (Int_3c_norm2 <= TINY(1.0_dp)) THEN
613 0 : DEALLOCATE (Phi_munu_l, Phi_munu_l_scaled, Int_3c_munu_P, mu_of_pair, nu_of_pair, jacobi_scaling)
614 0 : CALL timestop(handle)
615 0 : RETURN
616 : END IF
617 :
618 : ! D'_ll' = d_l D_ll' d_l' + λδ_ll', with d_l = 1/sqrt(D_ll).
619 10491 : DO igrid_point = 1, n_grid_points
620 293103 : Phi_munu_norm2 = SUM(Phi_munu_l(:, igrid_point)*Phi_munu_l(:, igrid_point))
621 10107 : jacobi_scaling(igrid_point) = 1.0_dp/SQRT(MAX(Phi_munu_norm2, jacobi_floor))
622 : Phi_munu_l_scaled(:, igrid_point) = &
623 293487 : jacobi_scaling(igrid_point)*Phi_munu_l(:, igrid_point)
624 : END DO
625 3456 : ALLOCATE (D_ll(n_grid_points, n_grid_points), d_lP(n_grid_points, nRI), Z_prime_lP(n_grid_points, nRI))
626 384 : CALL timeset("rirs_cluster_D_ll", phase_handle)
627 : CALL dgemm('T', 'N', n_grid_points, n_grid_points, n_AO_pairs, 1.0_dp, Phi_munu_l_scaled, n_AO_pairs, &
628 384 : Phi_munu_l_scaled, n_AO_pairs, 0.0_dp, D_ll, n_grid_points)
629 384 : CALL timestop(phase_handle)
630 10491 : DO igrid_point = 1, n_grid_points
631 10491 : D_ll(igrid_point, igrid_point) = D_ll(igrid_point, igrid_point) + tikhonov
632 : END DO
633 384 : CALL timeset("rirs_cluster_d_lP", phase_handle)
634 : CALL dgemm('T', 'N', n_grid_points, nRI, n_AO_pairs, 1.0_dp, Phi_munu_l_scaled, n_AO_pairs, &
635 384 : Int_3c_munu_P, n_AO_pairs, 0.0_dp, d_lP, n_grid_points)
636 384 : CALL timestop(phase_handle)
637 384 : CALL dpotrf('L', n_grid_points, D_ll, n_grid_points, lapack_info)
638 384 : IF (lapack_info /= 0) THEN
639 0 : DEALLOCATE (Phi_munu_l, Phi_munu_l_scaled, Int_3c_munu_P, D_ll, mu_of_pair, nu_of_pair, d_lP, &
640 0 : jacobi_scaling, Z_prime_lP)
641 0 : CALL timestop(handle)
642 0 : RETURN
643 : END IF
644 1975140 : Z_prime_lP(:, :) = d_lP
645 384 : CALL dpotrs('L', n_grid_points, nRI, D_ll, n_grid_points, Z_prime_lP, n_grid_points, lapack_info)
646 384 : DEALLOCATE (d_lP)
647 384 : IF (lapack_info /= 0) THEN
648 0 : DEALLOCATE (Phi_munu_l, Phi_munu_l_scaled, Int_3c_munu_P, D_ll, mu_of_pair, nu_of_pair, &
649 0 : jacobi_scaling, Z_prime_lP)
650 0 : CALL timestop(handle)
651 0 : RETURN
652 : END IF
653 :
654 384 : CALL timeset("rirs_cluster_residual", phase_handle)
655 1536 : ALLOCATE (R_munu_P(n_AO_pairs, nRI))
656 2348688 : R_munu_P(:, :) = Int_3c_munu_P
657 : ! R_μνP = (μν|P) - Σ_l ϕ_μ(r_l) ϕ_ν(r_l) d_l Z'_lP.
658 : CALL dgemm('N', 'N', n_AO_pairs, nRI, n_grid_points, -1.0_dp, Phi_munu_l_scaled, n_AO_pairs, &
659 384 : Z_prime_lP, n_grid_points, 1.0_dp, R_munu_P, n_AO_pairs)
660 2348688 : normalized_error = SUM(R_munu_P*R_munu_P)/Int_3c_norm2
661 384 : CALL timestop(phase_handle)
662 :
663 : ! For λ>0, the residual derivative includes the response of the regularized coefficients:
664 : ! dE/dϕ = [-2 R Z^T - 2λ(R Y^T - ϕ Y Z^T)] / ||Int_3c||^2.
665 2688 : ALLOCATE (D_inverse_Z_prime_lP(n_grid_points, nRI), dE_dPhi_munu_l(n_AO_pairs, n_grid_points))
666 1975140 : D_inverse_Z_prime_lP(:, :) = Z_prime_lP
667 384 : CALL dpotrs('L', n_grid_points, nRI, D_ll, n_grid_points, D_inverse_Z_prime_lP, n_grid_points, lapack_info)
668 384 : IF (lapack_info /= 0) THEN
669 0 : DEALLOCATE (Phi_munu_l, Phi_munu_l_scaled, Int_3c_munu_P, dE_dPhi_munu_l, D_ll, mu_of_pair, nu_of_pair, &
670 0 : R_munu_P, jacobi_scaling, D_inverse_Z_prime_lP, Z_prime_lP)
671 0 : CALL timestop(handle)
672 0 : RETURN
673 : END IF
674 384 : CALL timeset("rirs_cluster_gradient_products", phase_handle)
675 : CALL dgemm('N', 'T', n_AO_pairs, n_grid_points, nRI, -2.0_dp/Int_3c_norm2, R_munu_P, n_AO_pairs, &
676 384 : Z_prime_lP, n_grid_points, 0.0_dp, dE_dPhi_munu_l, n_AO_pairs)
677 384 : IF (tikhonov > 0.0_dp) THEN
678 : CALL dgemm('N', 'T', n_AO_pairs, n_grid_points, nRI, -2.0_dp*tikhonov/Int_3c_norm2, R_munu_P, n_AO_pairs, &
679 384 : D_inverse_Z_prime_lP, n_grid_points, 1.0_dp, dE_dPhi_munu_l, n_AO_pairs)
680 1536 : ALLOCATE (D_inverse_Z_prime_times_Z_prime_transpose_ll(n_grid_points, n_grid_points))
681 : CALL dgemm('N', 'T', n_grid_points, n_grid_points, nRI, 1.0_dp, D_inverse_Z_prime_lP, n_grid_points, &
682 384 : Z_prime_lP, n_grid_points, 0.0_dp, D_inverse_Z_prime_times_Z_prime_transpose_ll, n_grid_points)
683 : CALL dgemm('N', 'N', n_AO_pairs, n_grid_points, n_grid_points, &
684 : 2.0_dp*tikhonov/Int_3c_norm2, Phi_munu_l_scaled, n_AO_pairs, &
685 384 : D_inverse_Z_prime_times_Z_prime_transpose_ll, n_grid_points, 1.0_dp, dE_dPhi_munu_l, n_AO_pairs)
686 384 : DEALLOCATE (D_inverse_Z_prime_times_Z_prime_transpose_ll)
687 : END IF
688 :
689 384 : CALL timestop(phase_handle)
690 384 : CALL timeset("rirs_cluster_gradient_coordinates", phase_handle)
691 1152 : ALLOCATE (dPhi_munu_l_scaled(n_AO_pairs))
692 10491 : DO igrid_point = 1, n_grid_points
693 10107 : jacobi_scale = jacobi_scaling(igrid_point)
694 293103 : Phi_munu_norm2 = SUM(Phi_munu_l(:, igrid_point)*Phi_munu_l(:, igrid_point))
695 40812 : DO icartesian = 1, 3
696 879309 : DO i_ao_pair = 1, n_AO_pairs
697 848988 : imu = mu_of_pair(i_ao_pair)
698 848988 : inu = nu_of_pair(i_ao_pair)
699 848988 : symmetry_factor = SQRT(REAL(2 - MERGE(1, 0, imu == inu), dp))
700 : dPhi_munu_l = symmetry_factor*( &
701 : dPhi_alpha_l_mu(icartesian, igrid_point, imu)* &
702 : Phi_l_mu(igrid_point, inu) + &
703 : Phi_l_mu(igrid_point, imu)* &
704 848988 : dPhi_alpha_l_mu(icartesian, igrid_point, inu))
705 879309 : dPhi_munu_l_scaled(i_ao_pair) = jacobi_scale*dPhi_munu_l
706 : END DO
707 30321 : IF (Phi_munu_norm2 > jacobi_floor) THEN
708 : scaling_projection = &
709 829719 : DOT_PRODUCT(Phi_munu_l(:, igrid_point), dPhi_munu_l_scaled)/jacobi_scale
710 : dPhi_munu_l_scaled(:) = dPhi_munu_l_scaled - &
711 829719 : jacobi_scale**3*Phi_munu_l(:, igrid_point)*scaling_projection
712 : END IF
713 : coordinate_gradient(icartesian, igrid_point) = &
714 889416 : DOT_PRODUCT(dE_dPhi_munu_l(:, igrid_point), dPhi_munu_l_scaled)
715 : END DO
716 : END DO
717 :
718 384 : CALL timestop(phase_handle)
719 384 : max_abs_error = 0.0_dp
720 81360 : DO ip = 1, nRI
721 2348688 : DO i_ao_pair = 1, n_AO_pairs
722 : symmetry_factor = &
723 2267328 : SQRT(REAL(2 - MERGE(1, 0, mu_of_pair(i_ao_pair) == nu_of_pair(i_ao_pair)), dp))
724 2267328 : absolute_error = ABS(R_munu_P(i_ao_pair, ip))/symmetry_factor
725 2348304 : max_abs_error = MAX(max_abs_error, absolute_error)
726 : END DO
727 : END DO
728 :
729 384 : fit_successful = .TRUE.
730 0 : DEALLOCATE (Phi_munu_l, Phi_munu_l_scaled, Int_3c_munu_P, dE_dPhi_munu_l, dPhi_munu_l_scaled, D_ll, &
731 384 : mu_of_pair, nu_of_pair, R_munu_P, jacobi_scaling, D_inverse_Z_prime_lP, Z_prime_lP)
732 384 : CALL timestop(handle)
733 2688 : END SUBROUTINE evaluate_rirs_grid_cluster
734 :
735 : ! **************************************************************************************************
736 : !> \brief Sum gradient contributions from overlapping C_A into each atom-centred grid coordinate.
737 : !> \param cluster_atoms Atom indices in the local cluster.
738 : !> \param grids Atom-centred RI-RS grids.
739 : !> \param atom_coordinate_offsets Starting coordinate offset for each atom.
740 : !> \param cluster_coordinate_gradient Gradient for the cluster's physical grid points.
741 : !> \param coordinate_gradient Global flattened atom-relative gradient to update.
742 : ! **************************************************************************************************
743 384 : SUBROUTINE accumulate_atom_gradient(cluster_atoms, grids, atom_coordinate_offsets, &
744 384 : cluster_coordinate_gradient, coordinate_gradient)
745 : INTEGER, DIMENSION(:), INTENT(IN) :: cluster_atoms
746 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN) :: grids
747 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_coordinate_offsets
748 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: cluster_coordinate_gradient
749 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: coordinate_gradient
750 :
751 : INTEGER :: grid_point_offset, iatom, icartesian, &
752 : icluster_atom, igrid_point
753 :
754 384 : grid_point_offset = 0
755 1536 : DO icluster_atom = 1, SIZE(cluster_atoms)
756 1152 : iatom = cluster_atoms(icluster_atom)
757 11259 : DO igrid_point = 1, SIZE(grids(iatom)%raw_points, 2)
758 41580 : DO icartesian = 1, 3
759 : coordinate_gradient(atom_coordinate_offsets(iatom) + 3*(igrid_point - 1) + icartesian) = &
760 : coordinate_gradient(atom_coordinate_offsets(iatom) + 3*(igrid_point - 1) + icartesian) + &
761 40428 : cluster_coordinate_gradient(icartesian, grid_point_offset + igrid_point)
762 : END DO
763 : END DO
764 1536 : grid_point_offset = grid_point_offset + SIZE(grids(iatom)%raw_points, 2)
765 : END DO
766 384 : END SUBROUTINE accumulate_atom_gradient
767 :
768 : ! **************************************************************************************************
769 : !> \brief Pack all atom-centred grids into the single Cartesian vector required by L-BFGS.
770 : !> \param grids Atom-centred RI-RS grids.
771 : !> \param grid_coordinates Flattened atom-relative grid coordinates.
772 : !> \param atom_coordinate_offsets Starting coordinate offset for each atom.
773 : !>
774 : !> SIZE(grid_coordinates) = 3 Σ_A N_grid,A. Coordinates are ordered Cartesian component first,
775 : !> then grid point, then atom, matching the column-major layout of raw_points(3,N_grid,A).
776 : ! **************************************************************************************************
777 10 : SUBROUTINE pack_lbfgs_grid_coordinates(grids, grid_coordinates, atom_coordinate_offsets)
778 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN) :: grids
779 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: grid_coordinates(:)
780 : INTEGER, ALLOCATABLE, INTENT(OUT) :: atom_coordinate_offsets(:)
781 :
782 : INTEGER :: atom_coordinate_count, iatom, &
783 : n_grid_coordinates
784 :
785 30 : ALLOCATE (atom_coordinate_offsets(SIZE(grids)))
786 10 : n_grid_coordinates = 0
787 40 : DO iatom = 1, SIZE(grids)
788 30 : atom_coordinate_offsets(iatom) = n_grid_coordinates
789 100 : n_grid_coordinates = n_grid_coordinates + SIZE(grids(iatom)%raw_points)
790 : END DO
791 30 : ALLOCATE (grid_coordinates(n_grid_coordinates))
792 :
793 40 : DO iatom = 1, SIZE(grids)
794 90 : atom_coordinate_count = SIZE(grids(iatom)%raw_points)
795 : grid_coordinates(atom_coordinate_offsets(iatom) + 1: &
796 : atom_coordinate_offsets(iatom) + atom_coordinate_count) = &
797 70 : RESHAPE(grids(iatom)%raw_points, [atom_coordinate_count])
798 : END DO
799 10 : END SUBROUTINE pack_lbfgs_grid_coordinates
800 :
801 : ! **************************************************************************************************
802 : !> \brief Restore the optimized Cartesian vector to the persistent atom-centred RI-RS grids.
803 : !> \param grid_coordinates Flattened atom-relative grid coordinates.
804 : !> \param atom_coordinate_offsets Starting coordinate offset for each atom.
805 : !> \param grids Atom-centred RI-RS grids to update.
806 : ! **************************************************************************************************
807 256 : SUBROUTINE unpack_lbfgs_grid_coordinates(grid_coordinates, atom_coordinate_offsets, grids)
808 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: grid_coordinates
809 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_coordinate_offsets
810 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(INOUT) :: grids
811 :
812 : INTEGER :: atom_coordinate_count, iatom
813 :
814 256 : CPASSERT(SIZE(atom_coordinate_offsets) == SIZE(grids))
815 1024 : DO iatom = 1, SIZE(grids)
816 2304 : atom_coordinate_count = SIZE(grids(iatom)%raw_points)
817 : grids(iatom)%raw_points(:, :) = RESHAPE( &
818 : grid_coordinates(atom_coordinate_offsets(iatom) + 1: &
819 : atom_coordinate_offsets(iatom) + &
820 : atom_coordinate_count), &
821 29512 : SHAPE(grids(iatom)%raw_points))
822 : END DO
823 256 : END SUBROUTINE unpack_lbfgs_grid_coordinates
824 :
825 : END MODULE gw_ri_rs_grid_optimization
|
