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 Geometry and AO evaluation for real-space resolution of the identity.
10 : ! **************************************************************************************************
11 : MODULE ri_rs_utils
12 : USE atomic_kind_types, ONLY: get_atomic_kind_set
13 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
14 : gto_basis_set_type
15 : USE cell_types, ONLY: cell_type,&
16 : get_cell,&
17 : pbc
18 : USE kinds, ONLY: dp
19 : USE orbital_pointers, ONLY: indco,&
20 : ncoset
21 : USE particle_types, ONLY: particle_type
22 : USE physcon, ONLY: angstrom
23 : USE qs_kind_types, ONLY: get_qs_kind,&
24 : qs_kind_type
25 : USE ri_rs_types, ONLY: ri_rs_env
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 : PRIVATE
30 :
31 : PUBLIC :: evaluate_ao_basis_on_points, evaluate_ao_on_points, filter_grid_to_voronoi, &
32 : get_rirs_cluster_atoms, precompute_ri_rs_radii
33 :
34 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ri_rs_utils'
35 :
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief Compute per-atom AO and RI basis radii from the most diffuse Gaussian
40 : !> primitive in the AO ("ORB") and RI auxiliary ("RI_AUX") basis sets.
41 : !> Stores per-atom screening radii in ri_rs and prints a per-kind table.
42 : !> Radius: r_kind = sqrt(-log(eps)/alpha_min_kind), with eps = eps_filter.
43 : !> \param ri_rs RI-RS geometry references and screening radii
44 : !> \param basis_set_AO AO basis sets, one per atomic kind
45 : !> \param basis_set_RI auxiliary basis sets, one per atomic kind
46 : !> \param eps screening threshold
47 : !> \param unit_nr output unit
48 : ! **************************************************************************************************
49 48 : SUBROUTINE precompute_ri_rs_radii(ri_rs, basis_set_AO, basis_set_RI, eps, unit_nr)
50 : TYPE(ri_rs_env), INTENT(INOUT) :: ri_rs
51 : TYPE(gto_basis_set_p_type), INTENT(IN) :: basis_set_AO(:), basis_set_RI(:)
52 : REAL(KIND=dp), INTENT(IN) :: eps
53 : INTEGER, INTENT(IN) :: unit_nr
54 :
55 : CHARACTER(LEN=*), PARAMETER :: routineN = 'precompute_ri_rs_radii'
56 : REAL(KIND=dp), PARAMETER :: min_exponent_for_radius = 1.0E-3_dp
57 :
58 : INTEGER :: handle, i, iatom, ikind, j, natom, nkind
59 48 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
60 48 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha_min_ao_kind, alpha_min_ri_kind
61 48 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet_ao, zet_ri
62 :
63 48 : CALL timeset(routineN, handle)
64 :
65 48 : CPASSERT(ASSOCIATED(ri_rs%atomic_kind_set))
66 48 : CPASSERT(ASSOCIATED(ri_rs%particle_set))
67 48 : CPASSERT(0.0_dp < eps .AND. eps < 1.0_dp)
68 48 : nkind = SIZE(ri_rs%atomic_kind_set)
69 48 : natom = SIZE(ri_rs%particle_set)
70 :
71 192 : ALLOCATE (alpha_min_ao_kind(nkind), alpha_min_ri_kind(nkind))
72 118 : alpha_min_ao_kind = HUGE(1.0_dp)
73 118 : alpha_min_ri_kind = HUGE(1.0_dp)
74 :
75 118 : DO ikind = 1, nkind
76 70 : zet_ao => basis_set_AO(ikind)%gto_basis_set%zet
77 70 : zet_ri => basis_set_RI(ikind)%gto_basis_set%zet
78 236 : DO i = 1, SIZE(zet_ao, 1)
79 546 : DO j = 1, SIZE(zet_ao, 2)
80 476 : IF (zet_ao(i, j) > min_exponent_for_radius) THEN
81 310 : alpha_min_ao_kind(ikind) = MIN(alpha_min_ao_kind(ikind), zet_ao(i, j))
82 : END IF
83 : END DO
84 : END DO
85 188 : DO i = 1, SIZE(zet_ri, 1)
86 674 : DO j = 1, SIZE(zet_ri, 2)
87 604 : IF (zet_ri(i, j) > min_exponent_for_radius) THEN
88 534 : alpha_min_ri_kind(ikind) = MIN(alpha_min_ri_kind(ikind), zet_ri(i, j))
89 : END IF
90 : END DO
91 : END DO
92 : END DO
93 :
94 48 : CALL get_atomic_kind_set(atomic_kind_set=ri_rs%atomic_kind_set, kind_of=kind_of)
95 :
96 144 : ALLOCATE (ri_rs%radius_ao_per_atom(natom))
97 96 : ALLOCATE (ri_rs%radius_ri_per_atom(natom))
98 166 : DO iatom = 1, natom
99 118 : ikind = kind_of(iatom)
100 118 : ri_rs%radius_ao_per_atom(iatom) = SQRT(-LOG(eps)/alpha_min_ao_kind(ikind))
101 166 : ri_rs%radius_ri_per_atom(iatom) = SQRT(-LOG(eps)/alpha_min_ri_kind(ikind))
102 : END DO
103 :
104 48 : IF (unit_nr > 0) THEN
105 24 : WRITE (unit_nr, '(T2,A)') 'RI-RS basis radii (Å):'
106 24 : WRITE (unit_nr, '(T4,A6,2X,A4,2A14)') 'Kind', 'Elem', 'r_AO (Å)', 'r_RI (Å)'
107 59 : DO ikind = 1, nkind
108 : WRITE (unit_nr, '(T4,I6,2X,A4,2F14.4)') &
109 35 : ikind, &
110 35 : ri_rs%atomic_kind_set(ikind)%element_symbol, &
111 35 : SQRT(-LOG(eps)/alpha_min_ao_kind(ikind))*angstrom, &
112 94 : SQRT(-LOG(eps)/alpha_min_ri_kind(ikind))*angstrom
113 : END DO
114 24 : WRITE (unit_nr, '(A)') ' '
115 : END IF
116 :
117 48 : DEALLOCATE (alpha_min_ao_kind, alpha_min_ri_kind, kind_of)
118 :
119 48 : CALL timestop(handle)
120 :
121 48 : END SUBROUTINE precompute_ri_rs_radii
122 :
123 : !> \brief Retain source-grid points inside the Voronoi volume of one atom.
124 : !> \param points Source-grid coordinates on entry and Voronoi-filtered coordinates on return.
125 : !> \param icenter_atom Atom whose Voronoi volume is retained.
126 : !> \param particle_set Molecular atom positions.
127 : !> \param mask Optional membership mask; if present, leave points unchanged.
128 : !> \param atom_indices Optional subset of nuclei defining the Voronoi partition.
129 : ! **************************************************************************************************
130 30 : SUBROUTINE filter_grid_to_voronoi(points, icenter_atom, particle_set, mask, atom_indices)
131 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: points(:, :)
132 : INTEGER, INTENT(IN) :: icenter_atom
133 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
134 : LOGICAL, INTENT(OUT), OPTIONAL :: mask(:)
135 : INTEGER, INTENT(IN), OPTIONAL :: atom_indices(:)
136 :
137 : INTEGER :: iatom, iatom_index, ipoint, n_atoms, &
138 : n_keep
139 30 : LOGICAL, ALLOCATABLE :: keep(:)
140 : REAL(KIND=dp) :: displacement(3), distance2, max_radius2, &
141 : tolerance
142 30 : REAL(KIND=dp), ALLOCATABLE :: filtered_points(:, :)
143 :
144 5055 : ALLOCATE (keep(SIZE(points, 2)), SOURCE=.TRUE.)
145 19890 : max_radius2 = MAXVAL(SUM(points**2, DIM=1))
146 30 : n_atoms = SIZE(particle_set)
147 30 : IF (PRESENT(atom_indices)) n_atoms = SIZE(atom_indices)
148 120 : DO iatom_index = 1, n_atoms
149 90 : IF (PRESENT(atom_indices)) THEN
150 0 : iatom = atom_indices(iatom_index)
151 0 : CPASSERT(iatom >= 1 .AND. iatom <= SIZE(particle_set))
152 : ELSE
153 : iatom = iatom_index
154 : END IF
155 90 : IF (iatom == icenter_atom) CYCLE
156 240 : displacement(:) = particle_set(iatom)%r - particle_set(icenter_atom)%r
157 240 : distance2 = SUM(displacement**2)
158 : ! |R_B-R_A| > 2 max_l|r_l-R_A| cannot cut this finite point set.
159 60 : IF (distance2 > 4.0_dp*max_radius2) CYCLE
160 48 : tolerance = 32.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, distance2)
161 9929 : DO ipoint = 1, SIZE(points, 2)
162 9851 : IF (.NOT. keep(ipoint)) CYCLE
163 : ! |r_l-R_B|² - |r_l-R_A|² = |R_B-R_A|² - 2(r_l-R_A)·(R_B-R_A).
164 28934 : IF (2.0_dp*DOT_PRODUCT(points(:, ipoint), displacement) > distance2 + tolerance) THEN
165 3009 : keep(ipoint) = .FALSE.
166 4202 : ELSE IF (iatom < icenter_atom) THEN
167 5516 : IF (ABS(2.0_dp*DOT_PRODUCT(points(:, ipoint), displacement) - distance2) <= tolerance) THEN
168 49 : keep(ipoint) = .FALSE.
169 : END IF
170 : END IF
171 : END DO
172 : END DO
173 30 : IF (PRESENT(mask)) THEN
174 30 : CPASSERT(SIZE(mask) == SIZE(keep))
175 4995 : mask(:) = keep
176 30 : RETURN
177 : END IF
178 0 : n_keep = COUNT(keep)
179 0 : ALLOCATE (filtered_points(3, n_keep))
180 0 : IF (n_keep > 0) filtered_points(:, :) = RESHAPE(PACK(points, SPREAD(keep, 1, 3)), [3, n_keep])
181 0 : CALL MOVE_ALLOC(filtered_points, points)
182 60 : END SUBROUTINE filter_grid_to_voronoi
183 :
184 : ! **************************************************************************************************
185 : !> \brief Form C_A from nuclei within the specified radius, in global atom order.
186 : !> \param particle_set Molecular nuclei.
187 : !> \param cell Simulation cell.
188 : !> \param icenter_atom Central atom A.
189 : !> \param radius Cluster radius.
190 : !> \param atom_indices Cluster atom indices including A.
191 : ! **************************************************************************************************
192 30 : SUBROUTINE get_rirs_cluster_atoms(particle_set, cell, icenter_atom, radius, atom_indices)
193 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
194 : TYPE(cell_type), POINTER :: cell
195 : INTEGER, INTENT(IN) :: icenter_atom
196 : REAL(KIND=dp), INTENT(IN) :: radius
197 : INTEGER, ALLOCATABLE, INTENT(OUT) :: atom_indices(:)
198 :
199 : INTEGER :: iatom, ncluster
200 30 : INTEGER, ALLOCATABLE :: work(:)
201 :
202 90 : ALLOCATE (work(SIZE(particle_set)))
203 120 : ncluster = 0
204 120 : DO iatom = 1, SIZE(particle_set)
205 750 : IF (SUM(pbc(particle_set(iatom)%r - particle_set(icenter_atom)%r, cell)**2) <= radius**2) THEN
206 90 : ncluster = ncluster + 1
207 180 : work(ncluster) = iatom
208 : END IF
209 : END DO
210 180 : ALLOCATE (atom_indices(ncluster), SOURCE=work(:ncluster))
211 30 : END SUBROUTINE get_rirs_cluster_atoms
212 :
213 : ! **************************************************************************************************
214 : !> \brief Evaluate contracted spherical AOs, and optionally their grid-coordinate derivatives.
215 : !>
216 : !> For atom A and displacement d(alpha)=r_l(alpha)-R_A(alpha), the value is
217 : !>
218 : !> ϕ_μ(r_l) = Σ_(p,c) S_(pc,μ) d_x^lx d_y^ly d_z^lz exp(-ζ_p |d|^2).
219 : !>
220 : !> The optional derivative is evaluated analytically as
221 : !>
222 : !> d ϕ_μ(r_l)/d r_(l,α)
223 : !> = sum_(p,c) S_(pc,mu) exp(-zeta_p |d|^2)
224 : !> [d polynomial_c/d d_alpha - 2 zeta_p d_alpha polynomial_c].
225 : !>
226 : !> \param phi AO values, accumulated into phi(l,mu).
227 : !> \param grid_points Cartesian grid points, indexed (alpha,l).
228 : !> \param iatom Source atom whose contracted AOs are evaluated.
229 : !> \param particle_set Molecular particles.
230 : !> \param qs_kind_set Quickstep atomic kinds containing the ORB bases.
231 : !> \param cell Simulation cell used for the minimum-image displacement.
232 : !> \param dphi Optional AO derivatives, accumulated into dphi(alpha,l,mu).
233 : !> \param cutoff_squared Optional squared AO cutoff radius.
234 : ! **************************************************************************************************
235 136 : SUBROUTINE evaluate_ao_on_points(phi, grid_points, iatom, particle_set, qs_kind_set, cell, &
236 136 : dphi, cutoff_squared)
237 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: phi
238 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: grid_points
239 : INTEGER, INTENT(IN) :: iatom
240 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
241 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
242 : TYPE(cell_type), POINTER :: cell
243 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
244 : OPTIONAL :: dphi
245 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff_squared
246 :
247 : CHARACTER(len=*), PARAMETER :: routineN = 'evaluate_ao_on_points'
248 :
249 : INTEGER :: handle, ikind
250 : TYPE(gto_basis_set_type), POINTER :: basis
251 :
252 136 : CALL timeset(routineN, handle)
253 136 : CPASSERT(SIZE(grid_points, 1) == 3)
254 136 : CPASSERT(SIZE(phi, 1) == SIZE(grid_points, 2))
255 136 : IF (PRESENT(dphi)) THEN
256 0 : CPASSERT(SIZE(dphi, 1) == 3)
257 0 : CPASSERT(SIZE(dphi, 2) == SIZE(phi, 1))
258 0 : CPASSERT(SIZE(dphi, 3) == SIZE(phi, 2))
259 : END IF
260 :
261 136 : ikind = particle_set(iatom)%atomic_kind%kind_number
262 136 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis, basis_type='ORB')
263 136 : IF (.NOT. ASSOCIATED(basis)) THEN
264 0 : CALL timestop(handle)
265 0 : RETURN
266 : END IF
267 : CALL evaluate_ao_basis_on_points(phi, grid_points, basis, &
268 272 : particle_set(iatom)%r, cell, dphi, cutoff_squared)
269 136 : CALL timestop(handle)
270 : END SUBROUTINE evaluate_ao_on_points
271 :
272 : ! **************************************************************************************************
273 : !> \brief Evaluate one explicitly supplied contracted Gaussian basis on Cartesian points.
274 : !> \param phi AO values, accumulated into phi(l,mu).
275 : !> \param grid_points Cartesian grid points, indexed (alpha,l).
276 : !> \param basis Contracted Gaussian basis to evaluate.
277 : !> \param source_position Centre of the basis.
278 : !> \param cell Simulation cell defining minimum-image or image-summed displacements.
279 : !> \param dphi Optional AO derivatives, accumulated into dphi(alpha,l,mu).
280 : !> \param cutoff_squared Optional squared AO cutoff radius.
281 : !> \param periodic_images Sum images -1:1 along periodic axes and only image 0 along other axes.
282 : ! **************************************************************************************************
283 1392 : SUBROUTINE evaluate_ao_basis_on_points(phi, grid_points, basis, source_position, cell, &
284 1392 : dphi, cutoff_squared, periodic_images)
285 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: phi
286 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: grid_points
287 : TYPE(gto_basis_set_type), POINTER :: basis
288 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: source_position
289 : TYPE(cell_type), POINTER :: cell
290 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
291 : OPTIONAL :: dphi
292 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: cutoff_squared
293 : LOGICAL, INTENT(IN), OPTIONAL :: periodic_images
294 :
295 : INTEGER :: first_sgf, ialpha, ico, iend_co, image_extent(3), ipgf, ipoint, irow, iset, isgf, &
296 : ishell, istart_co, ix, iy, iz, l, last_sgf, lx, ly, lz, n_cart_total
297 : LOGICAL :: sum_images
298 : REAL(KIND=dp) :: exponent, exponential, hmat(3, 3), polynomial, polynomial_derivative(3), &
299 : radius2, relative(3), relative_raw(3), weight
300 :
301 1392 : CPASSERT(ASSOCIATED(basis))
302 1392 : CPASSERT(SIZE(grid_points, 1) == 3)
303 1392 : CPASSERT(SIZE(phi, 1) == SIZE(grid_points, 2))
304 1392 : IF (PRESENT(dphi)) THEN
305 1107 : CPASSERT(SIZE(dphi, 1) == 3)
306 1107 : CPASSERT(SIZE(dphi, 2) == SIZE(phi, 1))
307 1107 : CPASSERT(SIZE(dphi, 3) == SIZE(phi, 2))
308 : END IF
309 :
310 1392 : sum_images = .FALSE.
311 1392 : IF (PRESENT(periodic_images)) sum_images = periodic_images
312 1392 : image_extent = 0
313 1392 : hmat = 0.0_dp
314 1392 : IF (sum_images) THEN
315 0 : image_extent = cell%perd
316 0 : CALL get_cell(cell, h=hmat)
317 : END IF
318 :
319 : !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
320 : !$OMP SHARED(phi, dphi, grid_points, cell, basis, source_position, cutoff_squared, ncoset, indco) &
321 : !$OMP SHARED(sum_images, image_extent, hmat) &
322 : !$OMP PRIVATE(ipoint, relative, radius2, iset, n_cart_total, ishell, l, istart_co, &
323 : !$OMP iend_co, first_sgf, last_sgf, ipgf, exponent, exponential, isgf, ico, &
324 : !$OMP irow, weight, lx, ly, lz, polynomial, polynomial_derivative, ialpha, &
325 1392 : !$OMP relative_raw, ix, iy, iz)
326 : DO ipoint = 1, SIZE(grid_points, 2)
327 : relative_raw = grid_points(:, ipoint) - source_position
328 : IF (.NOT. sum_images) relative_raw = pbc(relative_raw, cell)
329 : DO ix = -image_extent(1), image_extent(1)
330 : DO iy = -image_extent(2), image_extent(2)
331 : DO iz = -image_extent(3), image_extent(3)
332 : relative = relative_raw
333 : IF (sum_images) relative = relative_raw - MATMUL(hmat, REAL([ix, iy, iz], dp))
334 : radius2 = DOT_PRODUCT(relative, relative)
335 : IF (PRESENT(cutoff_squared)) THEN
336 : IF (radius2 > cutoff_squared) CYCLE
337 : END IF
338 :
339 : DO iset = 1, basis%nset
340 : n_cart_total = ncoset(basis%lmax(iset))
341 : DO ishell = 1, basis%nshell(iset)
342 : l = basis%l(ishell, iset)
343 : istart_co = ncoset(l - 1) + 1
344 : iend_co = ncoset(l)
345 : first_sgf = basis%first_sgf(ishell, iset)
346 : last_sgf = basis%last_sgf(ishell, iset)
347 : DO ipgf = 1, basis%npgf(iset)
348 : exponent = basis%zet(ipgf, iset)
349 : exponential = EXP(-exponent*radius2)
350 : DO isgf = first_sgf, last_sgf
351 : DO ico = istart_co, iend_co
352 : irow = (ipgf - 1)*n_cart_total + ico
353 : weight = basis%sphi(irow, isgf)
354 : lx = indco(1, ico)
355 : ly = indco(2, ico)
356 : lz = indco(3, ico)
357 : polynomial = relative(1)**lx*relative(2)**ly*relative(3)**lz
358 : phi(ipoint, isgf) = phi(ipoint, isgf) + weight*polynomial*exponential
359 :
360 : IF (PRESENT(dphi)) THEN
361 : polynomial_derivative = 0.0_dp
362 : IF (lx > 0) polynomial_derivative(1) = REAL(lx, dp)*relative(1)**(lx - 1)* &
363 : relative(2)**ly*relative(3)**lz
364 : IF (ly > 0) polynomial_derivative(2) = REAL(ly, dp)*relative(1)**lx* &
365 : relative(2)**(ly - 1)*relative(3)**lz
366 : IF (lz > 0) polynomial_derivative(3) = REAL(lz, dp)*relative(1)**lx* &
367 : relative(2)**ly*relative(3)**(lz - 1)
368 : DO ialpha = 1, 3
369 : dphi(ialpha, ipoint, isgf) = dphi(ialpha, ipoint, isgf) + weight*exponential* &
370 : (polynomial_derivative(ialpha) - &
371 : 2.0_dp*exponent*relative(ialpha)*polynomial)
372 : END DO
373 : END IF
374 : END DO
375 : END DO
376 : END DO
377 : END DO
378 : END DO
379 : END DO
380 : END DO
381 : END DO
382 : END DO
383 : !$OMP END PARALLEL DO
384 1392 : END SUBROUTINE evaluate_ao_basis_on_points
385 :
386 : END MODULE ri_rs_utils
|
