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