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 Shared RI-RS grid I/O and contracted-Gaussian evaluation utilities.
10 : ! **************************************************************************************************
11 : MODULE rirs_grid_utils
12 : USE basis_set_types, ONLY: gto_basis_set_type
13 : USE cell_types, ONLY: cell_type,&
14 : pbc
15 : USE cp_files, ONLY: close_file,&
16 : open_file
17 : USE kinds, ONLY: default_path_length,&
18 : default_string_length,&
19 : dp
20 : USE orbital_pointers, ONLY: indco,&
21 : ncoset
22 : USE particle_types, ONLY: particle_type
23 : USE qs_kind_types, ONLY: get_qs_kind,&
24 : qs_kind_type
25 : #include "./base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rirs_grid_utils'
31 :
32 : PUBLIC :: evaluate_ao_basis_on_points, evaluate_ao_on_points, get_rirs_grid_filepath, &
33 : initialize_rirs_grid, read_rirs_grid_file, read_rirs_grid_npoints
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief Construct the path of an RI-RS grid file.
39 : !> \param element_symbol ...
40 : !> \param grid_select ...
41 : !> \param grid_file_suffix ...
42 : !> \param filepath ...
43 : ! **************************************************************************************************
44 74 : SUBROUTINE get_rirs_grid_filepath(element_symbol, grid_select, grid_file_suffix, filepath)
45 : CHARACTER(LEN=*), INTENT(IN) :: element_symbol
46 : INTEGER, INTENT(IN) :: grid_select
47 : CHARACTER(LEN=*), INTENT(IN) :: grid_file_suffix
48 : CHARACTER(LEN=*), INTENT(OUT) :: filepath
49 :
50 : CHARACTER(LEN=default_string_length) :: suffix
51 :
52 142 : SELECT CASE (grid_select)
53 : CASE (1)
54 68 : suffix = "_def2-tzvp-rs.ion"
55 : CASE (2)
56 4 : suffix = "_cc-pvtz-rs.ion"
57 : CASE (3)
58 2 : IF (LEN_TRIM(grid_file_suffix) > 0) THEN
59 2 : suffix = TRIM(grid_file_suffix)
60 : ELSE
61 0 : suffix = "_rirs.ion"
62 : END IF
63 : CASE DEFAULT
64 0 : CPABORT("Unknown grid_select (1=def2-TZVPP, 2=cc-pVTZ, 3=user-provided).")
65 : END SELECT
66 74 : filepath = "ri_rs_grid/"//TRIM(element_symbol)//TRIM(suffix)
67 74 : END SUBROUTINE get_rirs_grid_filepath
68 :
69 : ! **************************************************************************************************
70 : !> \brief Read the number of points from an RI-RS grid file.
71 : !> \param filename ...
72 : !> \return Number of grid points.
73 : ! **************************************************************************************************
74 24 : INTEGER FUNCTION read_rirs_grid_npoints(filename) RESULT(npoints)
75 : CHARACTER(LEN=*), INTENT(IN) :: filename
76 :
77 : INTEGER :: iunit
78 :
79 : CALL open_file(file_name=TRIM(filename), unit_number=iunit, &
80 12 : file_action='READ', file_status='OLD')
81 12 : CALL read_rirs_grid_header(iunit, filename, npoints)
82 12 : CALL close_file(unit_number=iunit)
83 12 : END FUNCTION read_rirs_grid_npoints
84 :
85 : ! **************************************************************************************************
86 : !> \brief Read Cartesian grid points from the existing CP2K RI-RS .ion format.
87 : !> \param filename Complete input filename; no suffix or element-name construction is performed.
88 : !> \param points Grid points in Bohr, indexed (alpha,l).
89 : ! **************************************************************************************************
90 62 : SUBROUTINE read_rirs_grid_file(filename, points)
91 : CHARACTER(LEN=*), INTENT(IN) :: filename
92 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
93 : INTENT(OUT) :: points
94 :
95 : CHARACTER(len=*), PARAMETER :: routineN = 'read_rirs_grid_file'
96 :
97 : CHARACTER(LEN=default_path_length) :: line
98 : INTEGER :: handle, ierr, iunit, l, npoints
99 : LOGICAL :: found_points
100 :
101 62 : CALL timeset(routineN, handle)
102 : CALL open_file(file_name=TRIM(filename), unit_number=iunit, &
103 62 : file_action='READ', file_status='OLD')
104 :
105 62 : CALL read_rirs_grid_header(iunit, filename, npoints)
106 :
107 62 : REWIND (iunit)
108 62 : found_points = .FALSE.
109 : DO
110 2224 : READ (iunit, '(A)', IOSTAT=ierr) line
111 2224 : IF (ierr /= 0) EXIT
112 2224 : IF (INDEX(line, '<grid_points>') > 0) THEN
113 : found_points = .TRUE.
114 : EXIT
115 : END IF
116 : END DO
117 62 : IF (.NOT. found_points) CPABORT('RI-RS .ion file has no <grid_points> block: '//TRIM(filename))
118 :
119 186 : ALLOCATE (points(3, npoints))
120 13018 : DO l = 1, npoints
121 12956 : READ (iunit, *, IOSTAT=ierr) points(:, l)
122 13018 : IF (ierr /= 0) CPABORT('Invalid grid point in RI-RS .ion file: '//TRIM(filename))
123 : END DO
124 62 : CALL close_file(unit_number=iunit)
125 62 : CALL timestop(handle)
126 124 : END SUBROUTINE read_rirs_grid_file
127 :
128 : ! **************************************************************************************************
129 : !> \brief Read and validate the point count in an open RI-RS grid file.
130 : !> \param iunit ...
131 : !> \param filename ...
132 : !> \param npoints ...
133 : ! **************************************************************************************************
134 74 : SUBROUTINE read_rirs_grid_header(iunit, filename, npoints)
135 : INTEGER, INTENT(IN) :: iunit
136 : CHARACTER(LEN=*), INTENT(IN) :: filename
137 : INTEGER, INTENT(OUT) :: npoints
138 :
139 : CHARACTER(LEN=default_path_length) :: line
140 : INTEGER :: colon, ierr
141 : LOGICAL :: found_size
142 :
143 74 : found_size = .FALSE.
144 74 : npoints = 0
145 : DO
146 2508 : READ (iunit, '(A)', IOSTAT=ierr) line
147 2508 : IF (ierr /= 0) EXIT
148 2508 : IF (INDEX(line, 'n points') > 0) THEN
149 74 : colon = INDEX(line, ':')
150 74 : IF (colon > 0) THEN
151 74 : READ (line(colon + 1:), *, IOSTAT=ierr) npoints
152 74 : found_size = ierr == 0 .AND. npoints > 0
153 : END IF
154 : EXIT
155 : END IF
156 : END DO
157 0 : IF (.NOT. found_size) CPABORT('RI-RS .ion file has no valid n points field: '//TRIM(filename))
158 74 : END SUBROUTINE read_rirs_grid_header
159 :
160 : ! **************************************************************************************************
161 : !> \brief Build one deterministic atom-specific grid from a tabulated source grid.
162 : !> \param points Source-grid coordinates on entry and selected coordinates on return.
163 : !> \param n_select Number of requested points.
164 : !> \param center_atom Atom on which the relative source grid is centred.
165 : !> \param particle_set Molecular atom positions.
166 : !> \param n_voronoi_candidates Number of candidates available in the atom's Voronoi volume.
167 : !> \author Jan Wilhelm
168 : ! **************************************************************************************************
169 18 : SUBROUTINE initialize_rirs_grid(points, n_select, center_atom, particle_set, n_voronoi_candidates)
170 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: points(:, :)
171 : INTEGER, INTENT(IN) :: n_select, center_atom
172 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
173 : INTEGER, INTENT(OUT) :: n_voronoi_candidates
174 :
175 18 : REAL(KIND=dp), ALLOCATABLE :: source_points(:, :)
176 :
177 18 : CPASSERT(SIZE(points, 1) == 3)
178 18 : CPASSERT(center_atom >= 1 .AND. center_atom <= SIZE(particle_set))
179 14664 : ALLOCATE (source_points, source=points)
180 18 : CALL filter_grid_to_voronoi(points, center_atom, particle_set)
181 18 : n_voronoi_candidates = SIZE(points, 2)
182 18 : IF (n_voronoi_candidates < n_select) CALL MOVE_ALLOC(source_points, points)
183 18 : CALL select_grid_points(points, n_select)
184 18 : END SUBROUTINE initialize_rirs_grid
185 :
186 : ! **************************************************************************************************
187 : !> \brief Retain source-grid points inside the Voronoi volume of one atom.
188 : !> \param points Source-grid coordinates on entry and Voronoi-filtered coordinates on return.
189 : !> \param center_atom Atom whose Voronoi volume is retained.
190 : !> \param particle_set Molecular atom positions.
191 : ! **************************************************************************************************
192 18 : SUBROUTINE filter_grid_to_voronoi(points, center_atom, particle_set)
193 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: points(:, :)
194 : INTEGER, INTENT(IN) :: center_atom
195 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
196 :
197 : INTEGER :: iatom, ipoint, n_keep
198 18 : LOGICAL, ALLOCATABLE :: keep(:)
199 : REAL(KIND=dp) :: other_distance_squared, &
200 : own_distance_squared, physical_point(3)
201 18 : REAL(KIND=dp), ALLOCATABLE :: filtered_points(:, :)
202 :
203 3702 : ALLOCATE (keep(SIZE(points, 2)), source=.TRUE.)
204 3666 : DO ipoint = 1, SIZE(points, 2)
205 14592 : physical_point = particle_set(center_atom)%r + points(:, ipoint)
206 14592 : own_distance_squared = SUM(points(:, ipoint)**2)
207 12600 : DO iatom = 1, SIZE(particle_set)
208 9924 : IF (iatom == center_atom) CYCLE
209 26568 : other_distance_squared = SUM((physical_point - particle_set(iatom)%r)**2)
210 9300 : IF (other_distance_squared < own_distance_squared) THEN
211 990 : keep(ipoint) = .FALSE.
212 990 : EXIT
213 : END IF
214 : END DO
215 : END DO
216 :
217 3666 : n_keep = COUNT(keep)
218 54 : ALLOCATE (filtered_points(3, n_keep))
219 54 : IF (n_keep > 0) filtered_points(:, :) = RESHAPE(PACK(points, SPREAD(keep, 1, 3)), [3, n_keep])
220 18 : CALL MOVE_ALLOC(filtered_points, points)
221 18 : END SUBROUTINE filter_grid_to_voronoi
222 :
223 : ! **************************************************************************************************
224 : !> \brief Select a deterministic maximin subset from a source grid.
225 : !> \param points Candidate coordinates on entry and selected coordinates on return.
226 : !> \param n_select Number of points to select.
227 : ! **************************************************************************************************
228 18 : SUBROUTINE select_grid_points(points, n_select)
229 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: points(:, :)
230 : INTEGER, INTENT(IN) :: n_select
231 :
232 : INTEGER :: candidate, i, isel, n_source
233 18 : LOGICAL, ALLOCATABLE :: available(:)
234 : REAL(KIND=dp) :: best_distance, distance_squared
235 18 : REAL(KIND=dp), ALLOCATABLE :: nearest_squared(:), selected_points(:, :)
236 :
237 18 : n_source = SIZE(points, 2)
238 18 : IF (n_select <= 0) CPABORT("GRID_SIZE point counts must be positive.")
239 18 : IF (n_select > n_source) THEN
240 0 : CPABORT("GRID_SIZE exceeds the number of available RI-RS source-grid points.")
241 : END IF
242 18 : IF (n_select == n_source) RETURN
243 :
244 126 : ALLOCATE (available(n_source), nearest_squared(n_source), selected_points(3, n_select))
245 2676 : available = .TRUE.
246 :
247 10704 : candidate = MINLOC(SUM(points**2, DIM=1), DIM=1)
248 72 : selected_points(:, 1) = points(:, candidate)
249 18 : available(candidate) = .FALSE.
250 10650 : nearest_squared(:) = SUM((points - SPREAD(selected_points(:, 1), 2, n_source))**2, DIM=1)
251 :
252 448 : DO isel = 2, n_select
253 : candidate = 0
254 : best_distance = -1.0_dp
255 81632 : DO i = 1, n_source
256 81632 : IF (available(i) .AND. nearest_squared(i) > best_distance) THEN
257 13704 : candidate = i
258 13704 : best_distance = nearest_squared(i)
259 : END IF
260 : END DO
261 430 : CPASSERT(candidate > 0)
262 1720 : selected_points(:, isel) = points(:, candidate)
263 430 : available(candidate) = .FALSE.
264 81650 : DO i = 1, n_source
265 81632 : IF (available(i)) THEN
266 287912 : distance_squared = SUM((points(:, i) - selected_points(:, isel))**2)
267 71978 : nearest_squared(i) = MIN(nearest_squared(i), distance_squared)
268 : END IF
269 : END DO
270 : END DO
271 :
272 18 : CALL MOVE_ALLOC(selected_points, points)
273 18 : END SUBROUTINE select_grid_points
274 :
275 : ! **************************************************************************************************
276 : !> \brief Evaluate contracted spherical AOs, and optionally their grid-coordinate derivatives.
277 : !>
278 : !> For atom A and displacement d(alpha)=r_l(alpha)-R_A(alpha), the value is
279 : !>
280 : !> phi_mu(r_l) = sum_(p,c) S_(pc,mu) d_x^lx d_y^ly d_z^lz exp(-zeta_p |d|^2).
281 : !>
282 : !> The optional derivative is evaluated analytically as
283 : !>
284 : !> d phi_mu(r_l)/d r_(l,alpha)
285 : !> = sum_(p,c) S_(pc,mu) exp(-zeta_p |d|^2)
286 : !> [d polynomial_c/d d_alpha - 2 zeta_p d_alpha polynomial_c].
287 : !>
288 : !> The routine is shared by production GW RI-RS and grid optimization. Production callers can
289 : !> pass a cutoff_squared; the optimizer omits it so moving points never cross a discontinuous
290 : !> AO-screening boundary.
291 : !>
292 : !> \param phi AO values, accumulated into phi(l,mu).
293 : !> \param grid_points Cartesian grid points, indexed (alpha,l).
294 : !> \param atom_index Source atom whose contracted AOs are evaluated.
295 : !> \param particle_set Molecular particles.
296 : !> \param qs_kind_set Quickstep atomic kinds containing the ORB bases.
297 : !> \param cell Simulation cell used for the minimum-image displacement.
298 : !> \param dphi Optional AO derivatives, accumulated into dphi(alpha,l,mu).
299 : !> \param cutoff_squared Optional squared AO cutoff radius.
300 : ! **************************************************************************************************
301 177 : SUBROUTINE evaluate_ao_on_points(phi, grid_points, atom_index, particle_set, qs_kind_set, cell, &
302 177 : dphi, cutoff_squared)
303 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: phi
304 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: grid_points
305 : INTEGER, INTENT(IN) :: atom_index
306 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
307 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
308 : TYPE(cell_type), POINTER :: cell
309 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
310 : OPTIONAL :: dphi
311 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff_squared
312 :
313 : CHARACTER(len=*), PARAMETER :: routineN = 'evaluate_ao_on_points'
314 :
315 : INTEGER :: handle, kind_index
316 : TYPE(gto_basis_set_type), POINTER :: basis
317 :
318 177 : CALL timeset(routineN, handle)
319 177 : CPASSERT(SIZE(grid_points, 1) == 3)
320 177 : CPASSERT(SIZE(phi, 1) == SIZE(grid_points, 2))
321 177 : IF (PRESENT(dphi)) THEN
322 0 : CPASSERT(SIZE(dphi, 1) == 3)
323 0 : CPASSERT(SIZE(dphi, 2) == SIZE(phi, 1))
324 0 : CPASSERT(SIZE(dphi, 3) == SIZE(phi, 2))
325 : END IF
326 :
327 177 : kind_index = particle_set(atom_index)%atomic_kind%kind_number
328 177 : CALL get_qs_kind(qs_kind_set(kind_index), basis_set=basis, basis_type='ORB')
329 177 : IF (.NOT. ASSOCIATED(basis)) THEN
330 0 : CALL timestop(handle)
331 0 : RETURN
332 : END IF
333 : CALL evaluate_ao_basis_on_points(phi, grid_points, basis, &
334 354 : particle_set(atom_index)%r, cell, dphi, cutoff_squared)
335 177 : CALL timestop(handle)
336 : END SUBROUTINE evaluate_ao_on_points
337 :
338 : ! **************************************************************************************************
339 : !> \brief Evaluate one explicitly supplied contracted Gaussian basis on Cartesian points.
340 : !> \param phi AO values, accumulated into phi(l,mu).
341 : !> \param grid_points Cartesian grid points, indexed (alpha,l).
342 : !> \param basis Contracted Gaussian basis to evaluate.
343 : !> \param source_position Centre of the basis.
344 : !> \param cell Simulation cell used for the minimum-image displacement.
345 : !> \param dphi Optional AO derivatives, accumulated into dphi(alpha,l,mu).
346 : !> \param cutoff_squared Optional squared AO cutoff radius.
347 : ! **************************************************************************************************
348 3921 : SUBROUTINE evaluate_ao_basis_on_points(phi, grid_points, basis, source_position, cell, &
349 3921 : dphi, cutoff_squared)
350 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: phi
351 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: grid_points
352 : TYPE(gto_basis_set_type), POINTER :: basis
353 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: source_position
354 : TYPE(cell_type), POINTER :: cell
355 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
356 : OPTIONAL :: dphi
357 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff_squared
358 :
359 : INTEGER :: alpha, first_sgf, ico, iend_co, ipgf, &
360 : iset, isgf, ishell, istart_co, l, &
361 : last_sgf, lx, ly, lz, n_cart_total, &
362 : point, row_index
363 : REAL(KIND=dp) :: exponent, exponential, polynomial, &
364 : polynomial_derivative(3), radius2, &
365 : relative(3), weight
366 :
367 3921 : CPASSERT(ASSOCIATED(basis))
368 3921 : CPASSERT(SIZE(grid_points, 1) == 3)
369 3921 : CPASSERT(SIZE(phi, 1) == SIZE(grid_points, 2))
370 3921 : IF (PRESENT(dphi)) THEN
371 3744 : CPASSERT(SIZE(dphi, 1) == 3)
372 3744 : CPASSERT(SIZE(dphi, 2) == SIZE(phi, 1))
373 3744 : CPASSERT(SIZE(dphi, 3) == SIZE(phi, 2))
374 : END IF
375 :
376 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
377 : !$OMP SHARED(phi, dphi, grid_points, cell, basis, source_position, cutoff_squared, ncoset, indco) &
378 : !$OMP PRIVATE(point, relative, radius2, iset, n_cart_total, ishell, l, istart_co, &
379 : !$OMP iend_co, first_sgf, last_sgf, ipgf, exponent, exponential, isgf, ico, &
380 3921 : !$OMP row_index, weight, lx, ly, lz, polynomial, polynomial_derivative, alpha)
381 : DO point = 1, SIZE(grid_points, 2)
382 : relative = pbc(grid_points(:, point) - source_position, cell)
383 : radius2 = DOT_PRODUCT(relative, relative)
384 : IF (PRESENT(cutoff_squared)) THEN
385 : IF (radius2 > cutoff_squared) CYCLE
386 : END IF
387 :
388 : DO iset = 1, basis%nset
389 : n_cart_total = ncoset(basis%lmax(iset))
390 : DO ishell = 1, basis%nshell(iset)
391 : l = basis%l(ishell, iset)
392 : istart_co = ncoset(l - 1) + 1
393 : iend_co = ncoset(l)
394 : first_sgf = basis%first_sgf(ishell, iset)
395 : last_sgf = basis%last_sgf(ishell, iset)
396 : DO ipgf = 1, basis%npgf(iset)
397 : exponent = basis%zet(ipgf, iset)
398 : exponential = EXP(-exponent*radius2)
399 : DO isgf = first_sgf, last_sgf
400 : DO ico = istart_co, iend_co
401 : row_index = (ipgf - 1)*n_cart_total + ico
402 : weight = basis%sphi(row_index, isgf)
403 : lx = indco(1, ico)
404 : ly = indco(2, ico)
405 : lz = indco(3, ico)
406 : polynomial = relative(1)**lx*relative(2)**ly*relative(3)**lz
407 : phi(point, isgf) = phi(point, isgf) + weight*polynomial*exponential
408 :
409 : IF (PRESENT(dphi)) THEN
410 : polynomial_derivative = 0.0_dp
411 : IF (lx > 0) polynomial_derivative(1) = REAL(lx, dp)*relative(1)**(lx - 1)* &
412 : relative(2)**ly*relative(3)**lz
413 : IF (ly > 0) polynomial_derivative(2) = REAL(ly, dp)*relative(1)**lx* &
414 : relative(2)**(ly - 1)*relative(3)**lz
415 : IF (lz > 0) polynomial_derivative(3) = REAL(lz, dp)*relative(1)**lx* &
416 : relative(2)**ly*relative(3)**(lz - 1)
417 : DO alpha = 1, 3
418 : dphi(alpha, point, isgf) = dphi(alpha, point, isgf) + weight*exponential* &
419 : (polynomial_derivative(alpha) - 2.0_dp*exponent*relative(alpha)*polynomial)
420 : END DO
421 : END IF
422 : END DO
423 : END DO
424 : END DO
425 : END DO
426 : END DO
427 : END DO
428 : !$OMP END PARALLEL DO
429 3921 : END SUBROUTINE evaluate_ao_basis_on_points
430 :
431 : END MODULE rirs_grid_utils
|
