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 Local-environment optimizer for atom-centred real-space RI grids.
10 : !> \author Jan Wilhelm
11 : ! **************************************************************************************************
12 : MODULE gw_optimize_ri_rs_grid
13 : USE basis_set_types, ONLY: gto_basis_set_type
14 : USE cell_types, ONLY: cell_type,&
15 : get_cell,&
16 : pbc,&
17 : use_perd_none
18 : USE cp_lbfgs, ONLY: setulb
19 : USE gw_integrals, ONLY: build_3c_integral_block_ctx,&
20 : gw_3c_ctx_create,&
21 : gw_3c_ctx_release,&
22 : gw_3c_ctx_type,&
23 : gw_3c_ws_create,&
24 : gw_3c_ws_release,&
25 : gw_3c_ws_type
26 : USE input_constants, ONLY: do_potential_truncated
27 : USE kinds, ONLY: dp
28 : USE libint_2c_3c, ONLY: libint_potential_type
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE particle_types, ONLY: particle_type
31 : USE physcon, ONLY: angstrom
32 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type,&
33 : rirs_grid_type
34 : USE qs_environment_types, ONLY: qs_environment_type
35 : USE rirs_grid_utils, ONLY: evaluate_ao_basis_on_points
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 : PRIVATE
40 :
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_optimize_ri_rs_grid'
42 : REAL(KIND=dp), PARAMETER, PRIVATE :: optimizer_accuracy = 1.0E-9_dp
43 : REAL(KIND=dp), PARAMETER, PRIVATE :: truncated_coulomb_cutoff = 3.0_dp/angstrom
44 :
45 : TYPE :: local_cluster_type
46 : INTEGER, ALLOCATABLE :: atoms(:)
47 : REAL(KIND=dp), ALLOCATABLE :: three_center(:, :, :)
48 : END TYPE local_cluster_type
49 :
50 : PUBLIC :: optimize_ri_rs_grid
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief Optimize the selected RI-RS grids stored in a GW band-structure environment.
56 : !> \param qs_env Quickstep environment used by the upstream three-center-integral context.
57 : !> \param bs_env GW environment containing all configuration, basis, geometry, parallel, and grid data.
58 : ! **************************************************************************************************
59 6 : SUBROUTINE optimize_ri_rs_grid(qs_env, bs_env)
60 : TYPE(qs_environment_type), POINTER :: qs_env
61 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'optimize_ri_rs_grid'
64 :
65 : INTEGER :: handle, iatom, ikind, n_variable, &
66 : periodic(3), unit_nr
67 6 : INTEGER, ALLOCATABLE :: ao_size(:), grid_offsets(:), ri_size(:)
68 : LOGICAL :: successful
69 : REAL(KIND=dp) :: f, maximum_absolute_error
70 12 : REAL(KIND=dp), ALLOCATABLE :: g(:), lower(:), upper(:), x(:)
71 : TYPE(cell_type), POINTER :: cell
72 78 : TYPE(gw_3c_ctx_type) :: integral_context
73 6 : TYPE(local_cluster_type), ALLOCATABLE :: clusters(:)
74 : TYPE(mp_para_env_type), POINTER :: para_env
75 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
76 :
77 6 : CALL timeset(routineN, handle)
78 : NULLIFY (cell, para_env, particle_set)
79 6 : cell => bs_env%ri_rs%cell
80 6 : particle_set => bs_env%ri_rs%particle_set
81 6 : para_env => bs_env%para_env
82 6 : CPASSERT(ASSOCIATED(cell))
83 6 : CPASSERT(ASSOCIATED(particle_set))
84 6 : CPASSERT(ASSOCIATED(para_env))
85 6 : CPASSERT(ALLOCATED(bs_env%basis_set_AO))
86 6 : CPASSERT(ALLOCATED(bs_env%basis_set_RI))
87 6 : CPASSERT(ALLOCATED(bs_env%ri_rs%grid_cache))
88 :
89 6 : IF (bs_env%ri_rs%grid_opt%max_iter < 1) THEN
90 0 : CPABORT("GRID_OPTIMIZATION%MAX_ITER must be positive")
91 : END IF
92 6 : IF (bs_env%ri_rs%grid_opt%cutoff_atomic_cluster <= 0.0_dp) THEN
93 0 : CPABORT("GRID_OPTIMIZATION%CUTOFF_ATOMIC_CLUSTER must be positive")
94 : END IF
95 6 : IF (bs_env%ri_rs%tikhonov < 0.0_dp) THEN
96 0 : CPABORT("RI_RS%TIKHONOV must not be negative")
97 : END IF
98 :
99 6 : CALL get_cell(cell, periodic=periodic)
100 24 : IF (ANY(periodic /= use_perd_none)) THEN
101 0 : CPABORT("GRID_OPTIMIZATION currently supports nonperiodic local environments only")
102 : END IF
103 :
104 6 : unit_nr = bs_env%unit_nr
105 24 : ALLOCATE (ao_size(SIZE(particle_set)), ri_size(SIZE(particle_set)))
106 24 : DO iatom = 1, SIZE(particle_set)
107 18 : ikind = particle_set(iatom)%atomic_kind%kind_number
108 18 : CPASSERT(ASSOCIATED(bs_env%basis_set_AO(ikind)%gto_basis_set))
109 18 : CPASSERT(ASSOCIATED(bs_env%basis_set_RI(ikind)%gto_basis_set))
110 18 : ao_size(iatom) = bs_env%basis_set_AO(ikind)%gto_basis_set%nsgf
111 24 : ri_size(iatom) = bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf
112 : END DO
113 48 : IF (ANY(ao_size < 1) .OR. ANY(ri_size < 1)) THEN
114 0 : CPABORT("Every atom in GRID_OPTIMIZATION needs ORB and RI_AUX functions")
115 : END IF
116 :
117 6 : CPASSERT(SIZE(bs_env%ri_rs%grid_cache) == SIZE(particle_set))
118 18 : ALLOCATE (grid_offsets(SIZE(bs_env%ri_rs%grid_cache)))
119 6 : n_variable = 0
120 24 : DO iatom = 1, SIZE(bs_env%ri_rs%grid_cache)
121 18 : grid_offsets(iatom) = n_variable
122 24 : n_variable = n_variable + 3*SIZE(bs_env%ri_rs%grid_cache(iatom)%raw_points, 2)
123 : END DO
124 : BLOCK
125 : TYPE(libint_potential_type) :: potential
126 6 : potential%potential_type = do_potential_truncated
127 6 : potential%cutoff_radius = truncated_coulomb_cutoff
128 : potential%omega = 0.0_dp
129 6 : potential%filename = "t_c_g.dat"
130 : CALL gw_3c_ctx_create(integral_context, qs_env, potential, &
131 6 : bs_env%basis_set_AO, bs_env%basis_set_AO, bs_env%basis_set_RI)
132 : END BLOCK
133 : CALL build_local_clusters(integral_context, particle_set, cell, ao_size, ri_size, &
134 6 : bs_env%ri_rs%grid_opt%cutoff_atomic_cluster, para_env, clusters)
135 :
136 36 : ALLOCATE (x(n_variable), g(n_variable), lower(n_variable), upper(n_variable))
137 6 : CALL pack_atom_grids(bs_env%ri_rs%grid_cache, grid_offsets, x)
138 1350 : lower(:) = x - 100.0_dp/angstrom
139 1350 : upper(:) = x + 100.0_dp/angstrom
140 : CALL optimize_grid_coordinates(x, lower, upper, optimizer_accuracy, &
141 6 : bs_env%ri_rs%grid_opt%max_iter, bs_env, clusters, grid_offsets)
142 : CALL grid_objective(x, bs_env, clusters, grid_offsets, f, g, &
143 6 : maximum_absolute_error, successful)
144 6 : IF (.NOT. successful) CPABORT("RI-RS grid optimization produced no regularized fit")
145 :
146 6 : IF (unit_nr > 0) THEN
147 3 : WRITE (unit_nr, '(/,T2,A)') 'RI-RS grid optimization terminated'
148 3 : WRITE (unit_nr, '(T2,A,T72,ES9.1)') 'Normalized 3C error:', f
149 : WRITE (unit_nr, '(T2,A,T72,ES9.1)') &
150 3 : 'Maximum absolute 3C error:', maximum_absolute_error
151 : END IF
152 :
153 24 : DO iatom = 1, SIZE(bs_env%ri_rs%grid_cache)
154 : bs_env%ri_rs%grid_cache(iatom)%npts = &
155 24 : SIZE(bs_env%ri_rs%grid_cache(iatom)%raw_points, 2)
156 : END DO
157 6 : IF (unit_nr > 0) FLUSH (unit_nr)
158 :
159 6 : bs_env%ri_rs%Z_lP_exists = .FALSE.
160 6 : CALL gw_3c_ctx_release(integral_context)
161 6 : CALL timestop(handle)
162 :
163 27 : END SUBROUTINE optimize_ri_rs_grid
164 :
165 : ! **************************************************************************************************
166 : !> \brief Evaluate the regularized three-centre-integral fitting error.
167 : !> \param coordinates Flattened atom-relative grid coordinates.
168 : !> \param bs_env GW environment containing RI-RS configuration and grid data.
169 : !> \param clusters Rank-local clusters and their exact three-centre integrals.
170 : !> \param grid_offsets Starting coordinate offset for each atom.
171 : !> \param value Mean normalized squared three-centre-integral error.
172 : !> \param gradient Derivative of value with respect to coordinates.
173 : !> \param maximum_absolute_error Largest absolute three-centre-integral error.
174 : !> \param valid Whether every local-cluster evaluation succeeded.
175 : ! **************************************************************************************************
176 832 : SUBROUTINE grid_objective(coordinates, bs_env, clusters, grid_offsets, value, gradient, &
177 : maximum_absolute_error, valid)
178 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: coordinates
179 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
180 : TYPE(local_cluster_type), DIMENSION(:), INTENT(IN) :: clusters
181 : INTEGER, DIMENSION(:), INTENT(IN) :: grid_offsets
182 : REAL(KIND=dp), INTENT(OUT) :: value
183 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: gradient
184 : REAL(KIND=dp), INTENT(OUT) :: maximum_absolute_error
185 : LOGICAL, INTENT(OUT) :: valid
186 :
187 : INTEGER :: all_valid, icluster, successful_clusters
188 : REAL(KIND=dp) :: cluster_maximum_absolute_error, &
189 : cluster_value
190 832 : REAL(KIND=dp), ALLOCATABLE :: cluster_gradient(:, :)
191 : TYPE(cell_type), POINTER :: cell
192 : TYPE(mp_para_env_type), POINTER :: para_env
193 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
194 :
195 832 : cell => bs_env%ri_rs%cell
196 832 : particle_set => bs_env%ri_rs%particle_set
197 832 : para_env => bs_env%para_env
198 832 : CALL unpack_atom_grids(coordinates, bs_env%ri_rs%grid_cache, grid_offsets)
199 832 : value = 0.0_dp
200 135214 : gradient = 0.0_dp
201 832 : successful_clusters = 0
202 832 : maximum_absolute_error = 0.0_dp
203 832 : valid = .TRUE.
204 2080 : DO icluster = 1, SIZE(clusters)
205 : CALL evaluate_local_cluster(clusters(icluster), bs_env%ri_rs%grid_cache, &
206 : bs_env, particle_set, cell, &
207 : cluster_value, cluster_gradient, &
208 1248 : cluster_maximum_absolute_error, valid)
209 1248 : IF (.NOT. valid) THEN
210 0 : IF (ALLOCATED(cluster_gradient)) DEALLOCATE (cluster_gradient)
211 : EXIT
212 : END IF
213 1248 : value = value + cluster_value
214 : CALL accumulate_atom_gradient(clusters(icluster)%atoms, bs_env%ri_rs%grid_cache, &
215 1248 : grid_offsets, cluster_gradient, gradient)
216 1248 : successful_clusters = successful_clusters + 1
217 : maximum_absolute_error = &
218 1248 : MAX(maximum_absolute_error, cluster_maximum_absolute_error)
219 2080 : DEALLOCATE (cluster_gradient)
220 : END DO
221 832 : all_valid = MERGE(1, 0, valid)
222 832 : CALL para_env%sum(value)
223 269596 : CALL para_env%sum(gradient)
224 832 : CALL para_env%sum(successful_clusters)
225 832 : CALL para_env%sum(all_valid)
226 832 : CALL para_env%max(maximum_absolute_error)
227 832 : valid = successful_clusters == SIZE(particle_set) .AND. all_valid == para_env%num_pe
228 832 : IF (.NOT. valid) THEN
229 0 : value = HUGE(value)
230 0 : gradient = 0.0_dp
231 0 : RETURN
232 : END IF
233 832 : value = value/REAL(SIZE(particle_set), dp)
234 135214 : gradient = gradient/REAL(SIZE(particle_set), dp)
235 832 : END SUBROUTINE grid_objective
236 :
237 : ! **************************************************************************************************
238 : !> \brief Minimize the RI-RS grid objective with CP2K's bound-constrained L-BFGS implementation.
239 : !> \param x Coordinates on entry and best coordinates found on return.
240 : !> \param lower Lower bound for each coordinate.
241 : !> \param upper Upper bound for each coordinate.
242 : !> \param accuracy Projected-gradient convergence threshold.
243 : !> \param max_evaluations Maximum number of objective evaluations.
244 : !> \param bs_env GW environment containing RI-RS configuration and grid data.
245 : !> \param clusters Rank-local clusters and their exact three-centre integrals.
246 : !> \param grid_offsets Starting coordinate offset for each atom.
247 : ! **************************************************************************************************
248 6 : SUBROUTINE optimize_grid_coordinates(x, lower, upper, accuracy, max_evaluations, &
249 6 : bs_env, clusters, grid_offsets)
250 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: x
251 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: lower, upper
252 : REAL(KIND=dp), INTENT(IN) :: accuracy
253 : INTEGER, INTENT(IN) :: max_evaluations
254 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
255 : TYPE(local_cluster_type), DIMENSION(:), INTENT(IN) :: clusters
256 : INTEGER, DIMENSION(:), INTENT(IN) :: grid_offsets
257 :
258 : INTEGER, PARAMETER :: memory = 7
259 :
260 : CHARACTER(LEN=60) :: csave, task
261 : INTEGER :: evaluations
262 : REAL(KIND=dp), DIMENSION(29) :: dsave
263 : REAL(KIND=dp), DIMENSION(2*memory*SIZE(x)+5*SIZE(x&
264 12 : )+11*memory**2+8*memory) :: wa
265 : REAL(KIND=dp) :: best_f, f, maximum_absolute_error
266 : LOGICAL, DIMENSION(4) :: lsave
267 : LOGICAL :: evaluation_ok, have_best
268 12 : INTEGER, DIMENSION(SIZE(x)) :: bound_type
269 : INTEGER, DIMENSION(44) :: isave
270 12 : INTEGER, DIMENSION(3*SIZE(x)) :: iwa
271 12 : REAL(KIND=dp), DIMENSION(SIZE(x)) :: best_x, g
272 :
273 6 : CPASSERT(SIZE(x) > 0)
274 30 : CPASSERT(ALL(SHAPE(lower) == SHAPE(x)) .AND. ALL(SHAPE(upper) == SHAPE(x)))
275 1350 : CPASSERT(ALL(lower <= upper))
276 6 : CPASSERT(accuracy > 0.0_dp .AND. max_evaluations > 0)
277 :
278 1350 : bound_type = 2
279 6 : task = 'START'
280 6 : csave = ''
281 6 : f = HUGE(f)
282 1350 : g = 0.0_dp
283 29112 : wa = 0.0_dp
284 4038 : iwa = 0
285 6 : lsave = .FALSE.
286 6 : isave = 0
287 6 : dsave = 0.0_dp
288 6 : evaluations = 0
289 6 : have_best = .FALSE.
290 6 : best_f = HUGE(best_f)
291 1350 : best_x = x
292 :
293 : DO
294 : CALL setulb(SIZE(x), memory, x, lower, upper, bound_type, f, g, &
295 1450 : 0.0_dp, accuracy, wa, iwa, task, -1, csave, lsave, isave, dsave, -1.0_dp)
296 1456 : IF (task(1:2) == 'FG') THEN
297 830 : IF (evaluations >= max_evaluations) EXIT
298 : CALL grid_objective(x, bs_env, clusters, grid_offsets, f, g, &
299 826 : maximum_absolute_error, evaluation_ok)
300 826 : evaluations = evaluations + 1
301 826 : IF (.NOT. evaluation_ok) EXIT
302 826 : IF (.NOT. have_best .OR. f < best_f) THEN
303 682 : have_best = .TRUE.
304 682 : best_f = f
305 112372 : best_x = x
306 : END IF
307 620 : ELSE IF (task(1:5) == 'NEW_X') THEN
308 : CYCLE
309 : ELSE
310 618 : EXIT
311 : END IF
312 : END DO
313 1350 : IF (have_best) x = best_x
314 6 : END SUBROUTINE optimize_grid_coordinates
315 :
316 : ! **************************************************************************************************
317 : !> \brief Build local clusters and precompute exact (mu nu|P) values on their owning rank.
318 : !> \param context Three-centre-integral evaluation context.
319 : !> \param particle_set Atomic positions and kinds.
320 : !> \param cell Simulation cell.
321 : !> \param ao_size Number of orbital basis functions on each atom.
322 : !> \param ri_size Number of auxiliary basis functions on each atom.
323 : !> \param cutoff Radius of each atom-centred local cluster.
324 : !> \param para_env MPI environment used to distribute cluster ownership.
325 : !> \param clusters Rank-local clusters and their exact three-centre integrals.
326 : ! **************************************************************************************************
327 6 : SUBROUTINE build_local_clusters(context, particle_set, cell, ao_size, ri_size, cutoff, &
328 : para_env, clusters)
329 : TYPE(gw_3c_ctx_type), INTENT(IN) :: context
330 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
331 : TYPE(cell_type), POINTER :: cell
332 : INTEGER, DIMENSION(:), INTENT(IN) :: ao_size, ri_size
333 : REAL(KIND=dp), INTENT(IN) :: cutoff
334 : TYPE(mp_para_env_type), POINTER :: para_env
335 : TYPE(local_cluster_type), ALLOCATABLE, INTENT(OUT) :: clusters(:)
336 :
337 : INTEGER :: center, ia, iatom, ja, jatom, ka, katom, &
338 : n_ao, n_cluster, n_ri
339 6 : INTEGER, ALLOCATABLE :: ao_offset(:), atom_buffer(:), &
340 6 : ri_offset(:)
341 : LOGICAL :: screened
342 6 : TYPE(gw_3c_ws_type) :: workspace
343 :
344 : n_cluster = COUNT([(MOD(center - 1, para_env%num_pe) == para_env%mepos, &
345 48 : center=1, SIZE(particle_set))])
346 0 : ALLOCATE (clusters(n_cluster), atom_buffer(SIZE(particle_set)), &
347 51 : ao_offset(SIZE(particle_set)), ri_offset(SIZE(particle_set)))
348 6 : CALL gw_3c_ws_create(workspace, context)
349 6 : n_cluster = 0
350 24 : DO center = 1, SIZE(particle_set)
351 18 : IF (MOD(center - 1, para_env%num_pe) /= para_env%mepos) CYCLE
352 9 : n_cluster = n_cluster + 1
353 9 : ia = 0
354 36 : DO iatom = 1, SIZE(particle_set)
355 225 : IF (SUM(pbc(particle_set(iatom)%r - particle_set(center)%r, cell)**2) <= cutoff**2) THEN
356 27 : ia = ia + 1
357 54 : atom_buffer(ia) = iatom
358 : END IF
359 : END DO
360 54 : ALLOCATE (clusters(n_cluster)%atoms(ia), source=atom_buffer(1:ia))
361 9 : ao_offset = 0
362 9 : ri_offset = 0
363 9 : n_ao = 0
364 9 : n_ri = 0
365 36 : DO ia = 1, SIZE(clusters(n_cluster)%atoms)
366 27 : iatom = clusters(n_cluster)%atoms(ia)
367 27 : ao_offset(iatom) = n_ao
368 27 : ri_offset(iatom) = n_ri
369 27 : n_ao = n_ao + ao_size(iatom)
370 36 : n_ri = n_ri + ri_size(iatom)
371 : END DO
372 10818 : ALLOCATE (clusters(n_cluster)%three_center(n_ao, n_ao, n_ri), source=0.0_dp)
373 42 : DO ia = 1, SIZE(clusters(n_cluster)%atoms)
374 27 : iatom = clusters(n_cluster)%atoms(ia)
375 117 : DO ka = 1, SIZE(clusters(n_cluster)%atoms)
376 81 : katom = clusters(n_cluster)%atoms(ka)
377 351 : DO ja = 1, SIZE(clusters(n_cluster)%atoms)
378 243 : jatom = clusters(n_cluster)%atoms(ja)
379 : CALL build_3c_integral_block_ctx( &
380 : clusters(n_cluster)%three_center, context, workspace, &
381 : atom_j=jatom, atom_k=katom, atom_i=iatom, &
382 : j_offset=ao_offset(jatom), k_offset=ao_offset(katom), &
383 324 : i_offset=ri_offset(iatom), screened=screened)
384 : END DO
385 : END DO
386 : END DO
387 : END DO
388 6 : CALL gw_3c_ws_release(workspace)
389 6 : END SUBROUTINE build_local_clusters
390 :
391 : ! **************************************************************************************************
392 : !> \brief Evaluate one complete cluster and return derivatives for its physical grid points.
393 : !> \param cluster Local atoms and exact three-centre integrals.
394 : !> \param grids Atom-centred RI-RS grids.
395 : !> \param bs_env GW environment containing orbital basis sets and RI-RS parameters.
396 : !> \param particle_set Atomic positions and kinds.
397 : !> \param cell Simulation cell.
398 : !> \param value Normalized squared three-centre-integral error.
399 : !> \param gradient Derivative with respect to the cluster's physical grid coordinates.
400 : !> \param maximum_absolute_error Largest absolute three-centre-integral error.
401 : !> \param successful Whether the regularized fitting equations were solved.
402 : ! **************************************************************************************************
403 1248 : SUBROUTINE evaluate_local_cluster(cluster, grids, bs_env, particle_set, cell, &
404 : value, gradient, maximum_absolute_error, successful)
405 : TYPE(local_cluster_type), INTENT(IN) :: cluster
406 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN) :: grids
407 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
408 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
409 : TYPE(cell_type), POINTER :: cell
410 : REAL(KIND=dp), INTENT(OUT) :: value
411 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: gradient(:, :)
412 : REAL(KIND=dp), INTENT(OUT) :: maximum_absolute_error
413 : LOGICAL, INTENT(OUT) :: successful
414 :
415 : INTEGER :: ao_offset, ia, iatom, ikind, n_ao, &
416 : n_atom_ao, n_grid, point_offset
417 1248 : REAL(KIND=dp), ALLOCATABLE :: dphi(:, :, :), phi(:, :), points(:, :)
418 : TYPE(gto_basis_set_type), POINTER :: basis
419 :
420 1248 : n_ao = SIZE(cluster%three_center, 1)
421 1248 : n_grid = 0
422 4992 : DO ia = 1, SIZE(cluster%atoms)
423 3744 : iatom = cluster%atoms(ia)
424 4992 : n_grid = n_grid + SIZE(grids(iatom)%raw_points, 2)
425 : END DO
426 : ALLOCATE (points(3, n_grid), phi(n_grid, n_ao), &
427 12480 : dphi(3, n_grid, n_ao), gradient(3, n_grid))
428 1248 : phi = 0.0_dp
429 1248 : dphi = 0.0_dp
430 1248 : point_offset = 0
431 4992 : DO ia = 1, SIZE(cluster%atoms)
432 3744 : iatom = cluster%atoms(ia)
433 : points(:, point_offset + 1:point_offset + SIZE(grids(iatom)%raw_points, 2)) = &
434 : SPREAD(particle_set(iatom)%r, 2, SIZE(grids(iatom)%raw_points, 2)) + &
435 272508 : grids(iatom)%raw_points
436 4992 : point_offset = point_offset + SIZE(grids(iatom)%raw_points, 2)
437 : END DO
438 : ao_offset = 0
439 4992 : DO ia = 1, SIZE(cluster%atoms)
440 3744 : iatom = cluster%atoms(ia)
441 3744 : ikind = particle_set(iatom)%atomic_kind%kind_number
442 3744 : basis => bs_env%basis_set_AO(ikind)%gto_basis_set
443 3744 : n_atom_ao = basis%nsgf
444 : CALL evaluate_ao_basis_on_points( &
445 : phi(:, ao_offset + 1:ao_offset + n_atom_ao), points, basis, &
446 : particle_set(iatom)%r, cell, &
447 3744 : dphi=dphi(:, :, ao_offset + 1:ao_offset + n_atom_ao))
448 4992 : ao_offset = ao_offset + n_atom_ao
449 : END DO
450 : CALL evaluate_rirs_grid_cluster(phi, dphi, cluster%three_center, bs_env%ri_rs%tikhonov, &
451 1248 : value, gradient, maximum_absolute_error, successful)
452 1248 : END SUBROUTINE evaluate_local_cluster
453 :
454 : ! **************************************************************************************************
455 : !> \brief Evaluate the regularized local-cluster three-centre fit and its coordinate gradient.
456 : !>
457 : !> The solve applies the same column Jacobi scaling and Tikhonov parameter as the production
458 : !> Z_lP construction. The objective contains only the normalized three-centre-integral residual.
459 : !> \param phi AO values, indexed (l,mu).
460 : !> \param dphi Cartesian derivatives of AO values, indexed (alpha,l,mu).
461 : !> \param three_center Exact three-centre integrals, indexed (mu,nu,P).
462 : !> \param tikhonov Tikhonov parameter used by the production RI-RS solve.
463 : !> \param value Normalized squared residual for this cluster.
464 : !> \param gradient Analytic derivative d value/d r(alpha,l).
465 : !> \param maximum_absolute_error Largest absolute error in an unweighted three-centre integral.
466 : !> \param successful False if the regularized normal equations cannot be solved.
467 : ! **************************************************************************************************
468 1248 : SUBROUTINE evaluate_rirs_grid_cluster(phi, dphi, three_center, tikhonov, value, gradient, &
469 : maximum_absolute_error, successful)
470 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi
471 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dphi, three_center
472 : REAL(KIND=dp), INTENT(IN) :: tikhonov
473 : REAL(KIND=dp), INTENT(OUT) :: value
474 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: gradient
475 : REAL(KIND=dp), INTENT(OUT) :: maximum_absolute_error
476 : LOGICAL, INTENT(OUT) :: successful
477 :
478 : CHARACTER(len=*), PARAMETER :: routineN = 'evaluate_rirs_grid_cluster'
479 : REAL(KIND=dp), PARAMETER :: jacobi_floor = 1.0E-16_dp
480 :
481 : INTEGER :: alpha, handle, info, ipair, l, mu, n_ao, &
482 : n_grid, n_pair, n_ri, nu, p
483 1248 : INTEGER, ALLOCATABLE :: pair_mu(:), pair_nu(:)
484 : REAL(KIND=dp) :: absolute_error, column_norm2, d_a, &
485 : denom, factor, projection, scale
486 1248 : REAL(KIND=dp), ALLOCATABLE :: a_matrix(:, :), a_scaled(:, :), b_matrix(:, :), &
487 1248 : derivative_a(:, :), derivative_scaled_a(:), gram(:, :), residual(:, :), rhs(:, :), &
488 1248 : scale_columns(:), y_matrix(:, :), yz_matrix(:, :), z_matrix(:, :)
489 :
490 1248 : CALL timeset(routineN, handle)
491 :
492 1248 : n_grid = SIZE(phi, 1)
493 1248 : n_ao = SIZE(phi, 2)
494 1248 : n_ri = SIZE(three_center, 3)
495 1248 : n_pair = n_ao*(n_ao + 1)/2
496 :
497 1248 : CPASSERT(SIZE(dphi, 1) == 3)
498 1248 : CPASSERT(SIZE(dphi, 2) == n_grid)
499 1248 : CPASSERT(SIZE(dphi, 3) == n_ao)
500 1248 : CPASSERT(SIZE(three_center, 1) == n_ao)
501 1248 : CPASSERT(SIZE(three_center, 2) == n_ao)
502 1248 : CPASSERT(SIZE(gradient, 1) == 3)
503 1248 : CPASSERT(SIZE(gradient, 2) == n_grid)
504 1248 : CPASSERT(tikhonov >= 0.0_dp)
505 :
506 1248 : value = HUGE(value)
507 270012 : gradient = 0.0_dp
508 1248 : maximum_absolute_error = HUGE(maximum_absolute_error)
509 1248 : successful = .FALSE.
510 1248 : IF (n_grid < 1 .OR. n_ri < 1) THEN
511 0 : CALL timestop(handle)
512 0 : RETURN
513 : END IF
514 :
515 : ALLOCATE (a_matrix(n_pair, n_grid), a_scaled(n_pair, n_grid), &
516 : b_matrix(n_pair, n_ri), pair_mu(n_pair), pair_nu(n_pair), &
517 17472 : scale_columns(n_grid))
518 9984 : ipair = 0
519 9984 : DO nu = 1, n_ao
520 44928 : DO mu = 1, nu
521 34944 : ipair = ipair + 1
522 34944 : pair_mu(ipair) = mu
523 34944 : pair_nu(ipair) = nu
524 34944 : factor = SQRT(REAL(2 - MERGE(1, 0, mu == nu), dp))
525 1916292 : DO l = 1, n_grid
526 1916292 : a_matrix(ipair, l) = factor*phi(l, mu)*phi(l, nu)
527 : END DO
528 777504 : DO p = 1, n_ri
529 768768 : b_matrix(ipair, p) = factor*three_center(mu, nu, p)
530 : END DO
531 : END DO
532 : END DO
533 :
534 761280 : denom = SUM(b_matrix*b_matrix)
535 1248 : IF (denom <= TINY(1.0_dp)) THEN
536 0 : DEALLOCATE (a_matrix, a_scaled, b_matrix, pair_mu, pair_nu, scale_columns)
537 0 : CALL timestop(handle)
538 0 : RETURN
539 : END IF
540 :
541 : ! Production Z_lP uses D'=d(A^T A)d+lambda I with d_l=1/sqrt((A^T A)_ll).
542 68439 : DO l = 1, n_grid
543 1948539 : column_norm2 = SUM(a_matrix(:, l)*a_matrix(:, l))
544 67191 : scale_columns(l) = 1.0_dp/SQRT(MAX(column_norm2, jacobi_floor))
545 1949787 : a_scaled(:, l) = scale_columns(l)*a_matrix(:, l)
546 : END DO
547 11232 : ALLOCATE (gram(n_grid, n_grid), rhs(n_grid, n_ri), z_matrix(n_grid, n_ri))
548 : CALL dgemm('T', 'N', n_grid, n_grid, n_pair, 1.0_dp, a_scaled, n_pair, &
549 1248 : a_scaled, n_pair, 0.0_dp, gram, n_grid)
550 68439 : DO l = 1, n_grid
551 68439 : gram(l, l) = gram(l, l) + tikhonov
552 : END DO
553 : CALL dgemm('T', 'N', n_grid, n_ri, n_pair, 1.0_dp, a_scaled, n_pair, &
554 1248 : b_matrix, n_pair, 0.0_dp, rhs, n_grid)
555 1248 : CALL dpotrf('L', n_grid, gram, n_grid, info)
556 1248 : IF (info /= 0) THEN
557 0 : DEALLOCATE (a_matrix, a_scaled, b_matrix, gram, pair_mu, pair_nu, rhs, &
558 0 : scale_columns, z_matrix)
559 0 : CALL timestop(handle)
560 0 : RETURN
561 : END IF
562 1438467 : z_matrix(:, :) = rhs
563 1248 : CALL dpotrs('L', n_grid, n_ri, gram, n_grid, z_matrix, n_grid, info)
564 1248 : DEALLOCATE (rhs)
565 1248 : IF (info /= 0) THEN
566 0 : DEALLOCATE (a_matrix, a_scaled, b_matrix, gram, pair_mu, pair_nu, &
567 0 : scale_columns, z_matrix)
568 0 : CALL timestop(handle)
569 0 : RETURN
570 : END IF
571 :
572 4992 : ALLOCATE (residual(n_pair, n_ri))
573 761280 : residual(:, :) = b_matrix
574 : CALL dgemm('N', 'N', n_pair, n_ri, n_grid, -1.0_dp, a_scaled, n_pair, &
575 1248 : z_matrix, n_grid, 1.0_dp, residual, n_pair)
576 761280 : value = SUM(residual*residual)/denom
577 :
578 : ! For lambda > 0, the residual-only derivative includes the response of the regularized
579 : ! coefficients. With Y=(A^T A+lambda I)^(-1) Z:
580 : ! dE/dA = [-2 R Z^T - 2 lambda (R Y^T - A Y Z^T)] / ||B||^2.
581 8736 : ALLOCATE (y_matrix(n_grid, n_ri), derivative_a(n_pair, n_grid))
582 1438467 : y_matrix(:, :) = z_matrix
583 1248 : CALL dpotrs('L', n_grid, n_ri, gram, n_grid, y_matrix, n_grid, info)
584 1248 : IF (info /= 0) THEN
585 0 : DEALLOCATE (a_matrix, a_scaled, b_matrix, derivative_a, gram, pair_mu, pair_nu, &
586 0 : residual, scale_columns, y_matrix, z_matrix)
587 0 : CALL timestop(handle)
588 0 : RETURN
589 : END IF
590 : CALL dgemm('N', 'T', n_pair, n_grid, n_ri, -2.0_dp/denom, residual, n_pair, &
591 1248 : z_matrix, n_grid, 0.0_dp, derivative_a, n_pair)
592 1248 : IF (tikhonov > 0.0_dp) THEN
593 : CALL dgemm('N', 'T', n_pair, n_grid, n_ri, -2.0_dp*tikhonov/denom, residual, n_pair, &
594 1248 : y_matrix, n_grid, 1.0_dp, derivative_a, n_pair)
595 4992 : ALLOCATE (yz_matrix(n_grid, n_grid))
596 : CALL dgemm('N', 'T', n_grid, n_grid, n_ri, 1.0_dp, y_matrix, n_grid, &
597 1248 : z_matrix, n_grid, 0.0_dp, yz_matrix, n_grid)
598 : CALL dgemm('N', 'N', n_pair, n_grid, n_grid, 2.0_dp*tikhonov/denom, a_scaled, n_pair, &
599 1248 : yz_matrix, n_grid, 1.0_dp, derivative_a, n_pair)
600 1248 : DEALLOCATE (yz_matrix)
601 : END IF
602 :
603 3744 : ALLOCATE (derivative_scaled_a(n_pair))
604 68439 : DO l = 1, n_grid
605 67191 : scale = scale_columns(l)
606 1948539 : column_norm2 = SUM(a_matrix(:, l)*a_matrix(:, l))
607 270012 : DO alpha = 1, 3
608 5845617 : DO ipair = 1, n_pair
609 5644044 : mu = pair_mu(ipair)
610 5644044 : nu = pair_nu(ipair)
611 5644044 : factor = SQRT(REAL(2 - MERGE(1, 0, mu == nu), dp))
612 : d_a = factor*(dphi(alpha, l, mu)*phi(l, nu) + &
613 5644044 : phi(l, mu)*dphi(alpha, l, nu))
614 5845617 : derivative_scaled_a(ipair) = scale*d_a
615 : END DO
616 201573 : IF (column_norm2 > jacobi_floor) THEN
617 4951431 : projection = DOT_PRODUCT(a_matrix(:, l), derivative_scaled_a)/scale
618 : derivative_scaled_a(:) = derivative_scaled_a - &
619 4951431 : scale**3*a_matrix(:, l)*projection
620 : END IF
621 5912808 : gradient(alpha, l) = DOT_PRODUCT(derivative_a(:, l), derivative_scaled_a)
622 : END DO
623 : END DO
624 :
625 1248 : maximum_absolute_error = 0.0_dp
626 27456 : DO p = 1, n_ri
627 761280 : DO ipair = 1, n_pair
628 733824 : factor = SQRT(REAL(2 - MERGE(1, 0, pair_mu(ipair) == pair_nu(ipair)), dp))
629 733824 : absolute_error = ABS(residual(ipair, p))/factor
630 760032 : maximum_absolute_error = MAX(maximum_absolute_error, absolute_error)
631 : END DO
632 : END DO
633 :
634 1248 : successful = .TRUE.
635 0 : DEALLOCATE (a_matrix, a_scaled, b_matrix, derivative_a, derivative_scaled_a, gram, &
636 1248 : pair_mu, pair_nu, residual, scale_columns, y_matrix, z_matrix)
637 1248 : CALL timestop(handle)
638 1248 : END SUBROUTINE evaluate_rirs_grid_cluster
639 :
640 : ! **************************************************************************************************
641 : !> \brief Sum each physical occurrence into its atom-relative coordinate.
642 : !> \param atoms Atom indices in the local cluster.
643 : !> \param grids Atom-centred RI-RS grids.
644 : !> \param grid_offsets Starting coordinate offset for each atom.
645 : !> \param physical_gradient Gradient for the cluster's physical grid points.
646 : !> \param gradient Global flattened atom-relative gradient to update.
647 : ! **************************************************************************************************
648 1248 : SUBROUTINE accumulate_atom_gradient(atoms, grids, grid_offsets, physical_gradient, gradient)
649 : INTEGER, DIMENSION(:), INTENT(IN) :: atoms
650 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN) :: grids
651 : INTEGER, DIMENSION(:), INTENT(IN) :: grid_offsets
652 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: physical_gradient
653 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: gradient
654 :
655 : INTEGER :: alpha, ia, iatom, l, point_offset
656 :
657 1248 : point_offset = 0
658 4992 : DO ia = 1, SIZE(atoms)
659 3744 : iatom = atoms(ia)
660 70935 : DO l = 1, SIZE(grids(iatom)%raw_points, 2)
661 272508 : DO alpha = 1, 3
662 : gradient(grid_offsets(iatom) + 3*(l - 1) + alpha) = &
663 : gradient(grid_offsets(iatom) + 3*(l - 1) + alpha) + &
664 268764 : physical_gradient(alpha, point_offset + l)
665 : END DO
666 : END DO
667 4992 : point_offset = point_offset + SIZE(grids(iatom)%raw_points, 2)
668 : END DO
669 1248 : END SUBROUTINE accumulate_atom_gradient
670 :
671 : ! **************************************************************************************************
672 : !> \brief Flatten atom-relative Cartesian coordinates into the optimizer vector.
673 : !> \param grids Atom-centred RI-RS grids.
674 : !> \param grid_offsets Starting coordinate offset for each atom.
675 : !> \param coordinates Flattened atom-relative grid coordinates.
676 : ! **************************************************************************************************
677 6 : SUBROUTINE pack_atom_grids(grids, grid_offsets, coordinates)
678 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(IN) :: grids
679 : INTEGER, DIMENSION(:), INTENT(IN) :: grid_offsets
680 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: coordinates
681 :
682 : INTEGER :: iatom, n_coordinate
683 :
684 6 : CPASSERT(SIZE(grid_offsets) == SIZE(grids))
685 24 : DO iatom = 1, SIZE(grids)
686 54 : n_coordinate = SIZE(grids(iatom)%raw_points)
687 : coordinates(grid_offsets(iatom) + 1:grid_offsets(iatom) + n_coordinate) = &
688 42 : RESHAPE(grids(iatom)%raw_points, [n_coordinate])
689 : END DO
690 6 : END SUBROUTINE pack_atom_grids
691 :
692 : ! **************************************************************************************************
693 : !> \brief Restore atom-relative Cartesian coordinates from the optimizer vector.
694 : !> \param coordinates Flattened atom-relative grid coordinates.
695 : !> \param grids Atom-centred RI-RS grids to update.
696 : !> \param grid_offsets Starting coordinate offset for each atom.
697 : ! **************************************************************************************************
698 832 : SUBROUTINE unpack_atom_grids(coordinates, grids, grid_offsets)
699 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: coordinates
700 : TYPE(rirs_grid_type), DIMENSION(:), INTENT(INOUT) :: grids
701 : INTEGER, DIMENSION(:), INTENT(IN) :: grid_offsets
702 :
703 : INTEGER :: iatom, n_coordinate
704 :
705 832 : CPASSERT(SIZE(grid_offsets) == SIZE(grids))
706 3328 : DO iatom = 1, SIZE(grids)
707 7488 : n_coordinate = SIZE(grids(iatom)%raw_points)
708 : grids(iatom)%raw_points(:, :) = RESHAPE( &
709 : coordinates(grid_offsets(iatom) + 1:grid_offsets(iatom) + n_coordinate), &
710 187496 : SHAPE(grids(iatom)%raw_points))
711 : END DO
712 832 : END SUBROUTINE unpack_atom_grids
713 :
714 0 : END MODULE gw_optimize_ri_rs_grid
|
