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 Feature vectors for describing chemical environments in a rotationally invariant fashion.
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE pao_ml_descriptor
13 : USE atomic_kind_types, ONLY: get_atomic_kind
14 : USE basis_set_types, ONLY: gto_basis_set_type
15 : USE cell_types, ONLY: cell_type,&
16 : pbc
17 : USE kinds, ONLY: dp
18 : USE mathconstants, ONLY: fourpi,&
19 : rootpi
20 : USE mathlib, ONLY: diamat_all
21 : USE pao_input, ONLY: pao_ml_desc_overlap,&
22 : pao_ml_desc_pot,&
23 : pao_ml_desc_r12
24 : USE pao_potentials, ONLY: pao_calc_gaussian
25 : USE pao_types, ONLY: pao_env_type
26 : USE particle_types, ONLY: particle_type
27 : USE qs_kind_types, ONLY: get_qs_kind,&
28 : pao_descriptor_type,&
29 : qs_kind_type
30 : USE util, ONLY: sort
31 : #include "./base/base_uses.f90"
32 :
33 : IMPLICIT NONE
34 :
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_ml_descriptor'
38 :
39 : PUBLIC :: pao_ml_calc_descriptor
40 :
41 : CONTAINS
42 :
43 : ! **************************************************************************************************
44 : !> \brief Calculates a descriptor for chemical environment of given atom
45 : !> \param pao ...
46 : !> \param particle_set ...
47 : !> \param qs_kind_set ...
48 : !> \param cell ...
49 : !> \param iatom ...
50 : !> \param descriptor ...
51 : !> \param descr_grad ...
52 : !> \param forces ...
53 : ! **************************************************************************************************
54 208 : SUBROUTINE pao_ml_calc_descriptor(pao, particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
55 : TYPE(pao_env_type), POINTER :: pao
56 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
57 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
58 : TYPE(cell_type), POINTER :: cell
59 : INTEGER, INTENT(IN) :: iatom
60 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: descriptor
61 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: descr_grad
62 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
63 :
64 : CHARACTER(len=*), PARAMETER :: routineN = 'pao_ml_calc_descriptor'
65 :
66 : INTEGER :: handle
67 :
68 208 : CALL timeset(routineN, handle)
69 :
70 208 : CPASSERT(PRESENT(forces) .EQV. PRESENT(descr_grad))
71 :
72 208 : SELECT CASE (pao%ml_descriptor)
73 : CASE (pao_ml_desc_pot)
74 150 : CALL calc_descriptor_pot(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
75 : CASE (pao_ml_desc_overlap)
76 118 : CALL calc_descriptor_overlap(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
77 : CASE (pao_ml_desc_r12)
78 320 : CALL calc_descriptor_r12(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
79 : CASE DEFAULT
80 208 : CPABORT("PAO: unknown descriptor")
81 : END SELECT
82 :
83 208 : CALL timestop(handle)
84 208 : END SUBROUTINE pao_ml_calc_descriptor
85 :
86 : ! **************************************************************************************************
87 : !> \brief Calculates a descriptor based on the eigenvalues of V_neighbors
88 : !> \param particle_set ...
89 : !> \param qs_kind_set ...
90 : !> \param cell ...
91 : !> \param iatom ...
92 : !> \param descriptor ...
93 : !> \param descr_grad ...
94 : !> \param forces ...
95 : ! **************************************************************************************************
96 54 : SUBROUTINE calc_descriptor_pot(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
97 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
98 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
99 : TYPE(cell_type), POINTER :: cell
100 : INTEGER, INTENT(IN) :: iatom
101 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: descriptor
102 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: descr_grad
103 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
104 :
105 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_descriptor_pot'
106 :
107 : INTEGER :: handle, i, idesc, ikind, jatom, jkind, &
108 : k, N, natoms, ndesc
109 : REAL(dp) :: beta, w, weight
110 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: V_evals
111 54 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: block_M, block_V, V_evecs
112 54 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_D
113 : REAL(dp), DIMENSION(3) :: Ra, Rab, Rb
114 : TYPE(gto_basis_set_type), POINTER :: basis_set
115 54 : TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: pao_descriptors
116 :
117 54 : CALL timeset(routineN, handle)
118 :
119 54 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
120 54 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set, pao_descriptors=pao_descriptors)
121 54 : N = basis_set%nsgf
122 54 : natoms = SIZE(particle_set)
123 54 : ndesc = SIZE(pao_descriptors)
124 54 : IF (ndesc == 0) CPABORT("No PAO_DESCRIPTOR section found")
125 :
126 432 : ALLOCATE (block_V(N, N), V_evecs(N, N), V_evals(N))
127 150 : IF (PRESENT(descriptor)) ALLOCATE (descriptor(N*ndesc))
128 90 : IF (PRESENT(forces)) ALLOCATE (block_D(N, N, 3), block_M(N, N))
129 :
130 152 : DO idesc = 1, ndesc
131 :
132 : ! construct matrix V_block from neighboring atoms
133 98 : block_V = 0.0_dp
134 294 : DO jatom = 1, natoms
135 196 : IF (jatom == iatom) CYCLE
136 392 : Ra = particle_set(iatom)%r
137 392 : Rb = particle_set(jatom)%r
138 98 : Rab = pbc(ra, rb, cell)
139 98 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
140 98 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=pao_descriptors)
141 98 : IF (SIZE(pao_descriptors) /= ndesc) THEN
142 0 : CPABORT("Not all KINDs have the same number of PAO_DESCRIPTOR sections")
143 : END IF
144 98 : weight = pao_descriptors(idesc)%weight
145 98 : beta = pao_descriptors(idesc)%beta
146 392 : CALL pao_calc_gaussian(basis_set, block_V=block_V, Rab=Rab, lpot=0, beta=beta, weight=weight)
147 : END DO
148 :
149 : ! diagonalize block_V
150 3038 : V_evecs(:, :) = block_V(:, :)
151 98 : CALL diamat_all(V_evecs, V_evals)
152 :
153 : ! use eigenvalues of V_block as descriptor
154 98 : IF (PRESENT(descriptor)) THEN
155 528 : descriptor((idesc - 1)*N + 1:idesc*N) = V_evals(:)
156 : END IF
157 :
158 : ! FORCES ----------------------------------------------------------------------------------
159 152 : IF (PRESENT(forces)) THEN
160 10 : CPASSERT(PRESENT(descr_grad))
161 10 : block_M = 0.0_dp
162 60 : DO k = 1, N
163 50 : w = descr_grad((idesc - 1)*N + k)
164 3460 : block_M(:, :) = block_M(:, :) + w*MATMUL(V_evecs(:, k:k), TRANSPOSE(V_evecs(:, k:k)))
165 : END DO
166 30 : DO jatom = 1, natoms
167 20 : IF (jatom == iatom) CYCLE
168 40 : Ra = particle_set(iatom)%r
169 40 : Rb = particle_set(jatom)%r
170 10 : Rab = pbc(ra, rb, cell)
171 10 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
172 10 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=pao_descriptors)
173 10 : weight = pao_descriptors(idesc)%weight
174 10 : beta = pao_descriptors(idesc)%beta
175 10 : block_D = 0.0_dp
176 10 : CALL pao_calc_gaussian(basis_set, block_D=block_D, Rab=Rab, lpot=0, beta=beta, weight=weight)
177 60 : DO i = 1, 3
178 930 : forces(iatom, i) = forces(iatom, i) - SUM(block_M*block_D(:, :, i))
179 950 : forces(jatom, i) = forces(jatom, i) + SUM(block_M*block_D(:, :, i))
180 : END DO
181 : END DO
182 : END IF
183 :
184 : END DO
185 :
186 54 : CALL timestop(handle)
187 108 : END SUBROUTINE calc_descriptor_pot
188 :
189 : ! **************************************************************************************************
190 : !> \brief Calculates a descriptor based on the eigenvalues of local overlap matrix
191 : !> \param particle_set ...
192 : !> \param qs_kind_set ...
193 : !> \param cell ...
194 : !> \param iatom ...
195 : !> \param descriptor ...
196 : !> \param descr_grad ...
197 : !> \param forces ...
198 : ! **************************************************************************************************
199 42 : SUBROUTINE calc_descriptor_overlap(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
200 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
201 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
202 : TYPE(cell_type), POINTER :: cell
203 : INTEGER, INTENT(IN) :: iatom
204 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: descriptor
205 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: descr_grad
206 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
207 :
208 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_descriptor_overlap'
209 :
210 : INTEGER :: handle, idesc, ikind, j, jatom, jkind, &
211 : k, katom, kkind, N, natoms, ndesc
212 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: neighbor_order
213 : REAL(dp) :: beta_sum, deriv, exponent, integral, jbeta, jweight, kbeta, kweight, &
214 : normalization, Rij2, Rik2, Rjk2, sbeta, screening_radius, screening_volume, w
215 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
216 42 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: block_M, block_S, S_evecs
217 : REAL(dp), DIMENSION(3) :: Ri, Rij, Rik, Rj, Rjk, Rk
218 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: neighbor_dist
219 42 : TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: ipao_descriptors, jpao_descriptors, &
220 42 : kpao_descriptors
221 :
222 42 : CALL timeset(routineN, handle)
223 :
224 42 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
225 42 : CALL get_qs_kind(qs_kind_set(ikind), pao_descriptors=ipao_descriptors)
226 :
227 42 : natoms = SIZE(particle_set)
228 42 : ndesc = SIZE(ipao_descriptors)
229 42 : IF (ndesc == 0) CPABORT("No PAO_DESCRIPTOR section found")
230 :
231 : ! determine largest screening radius
232 42 : screening_radius = 0.0_dp
233 118 : DO idesc = 1, ndesc
234 118 : screening_radius = MAX(screening_radius, ipao_descriptors(idesc)%screening_radius)
235 : END DO
236 :
237 : ! estimate maximum number of neighbors within screening
238 42 : screening_volume = fourpi/3.0_dp*screening_radius**3
239 42 : N = INT(screening_volume/35.0_dp) ! rule of thumb
240 :
241 336 : ALLOCATE (block_S(N, N), S_evals(N), S_evecs(N, N))
242 118 : IF (PRESENT(descriptor)) ALLOCATE (descriptor(N*ndesc))
243 50 : IF (PRESENT(forces)) ALLOCATE (block_M(N, N))
244 :
245 : !find neighbors
246 : !TODO: this is a quadratic algorithm, use a neighbor-list instead
247 210 : ALLOCATE (neighbor_dist(natoms), neighbor_order(natoms))
248 168 : Ri = particle_set(iatom)%r
249 132 : DO jatom = 1, natoms
250 360 : Rj = particle_set(jatom)%r
251 90 : Rij = pbc(Ri, Rj, cell)
252 402 : neighbor_dist(jatom) = SQRT(SUM(Rij**2))
253 : END DO
254 42 : CALL sort(neighbor_dist, natoms, neighbor_order)
255 42 : CPASSERT(neighbor_order(1) == iatom) !central atom should be closesd to itself
256 :
257 : ! check if N was chosen large enough
258 42 : IF (natoms > N) THEN
259 0 : IF (neighbor_dist(N + 1) < screening_radius) THEN
260 0 : CPABORT("PAO heuristic for descriptor size broke down")
261 : END IF
262 : END IF
263 :
264 118 : DO idesc = 1, ndesc
265 76 : sbeta = ipao_descriptors(idesc)%screening
266 :
267 : ! construct matrix S_block from neighboring atoms
268 76 : block_S = 0.0_dp
269 234 : DO j = 1, MIN(natoms, N)
270 568 : DO k = 1, MIN(natoms, N)
271 334 : jatom = neighbor_order(j)
272 334 : katom = neighbor_order(k)
273 :
274 : ! get weigths and betas
275 334 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
276 334 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=jpao_descriptors)
277 334 : CALL get_atomic_kind(particle_set(katom)%atomic_kind, kind_number=kkind)
278 334 : CALL get_qs_kind(qs_kind_set(kkind), pao_descriptors=kpao_descriptors)
279 334 : IF (SIZE(jpao_descriptors) /= ndesc .OR. SIZE(kpao_descriptors) /= ndesc) THEN
280 0 : CPABORT("Not all KINDs have the same number of PAO_DESCRIPTOR sections")
281 : END IF
282 334 : jweight = jpao_descriptors(idesc)%weight
283 334 : jbeta = jpao_descriptors(idesc)%beta
284 334 : kweight = kpao_descriptors(idesc)%weight
285 334 : kbeta = kpao_descriptors(idesc)%beta
286 334 : beta_sum = sbeta + jbeta + kbeta
287 :
288 : ! get distances
289 1336 : Rj = particle_set(jatom)%r
290 1336 : Rk = particle_set(katom)%r
291 334 : Rij = pbc(Ri, Rj, cell)
292 334 : Rik = pbc(Ri, Rk, cell)
293 334 : Rjk = pbc(Rj, Rk, cell)
294 1336 : Rij2 = SUM(Rij**2)
295 1336 : Rik2 = SUM(Rik**2)
296 1336 : Rjk2 = SUM(Rjk**2)
297 :
298 : ! calculate integral over three Gaussians
299 334 : exponent = -(sbeta*jbeta*Rij2 + sbeta*kbeta*Rik2 + jbeta*kbeta*Rjk2)/beta_sum
300 334 : integral = EXP(exponent)*rootpi/SQRT(beta_sum)
301 334 : normalization = SQRT(jbeta*kbeta)/rootpi**2
302 1160 : block_S(j, k) = jweight*kweight*normalization*integral
303 : END DO
304 : END DO
305 :
306 : ! diagonalize V_block
307 1653532 : S_evecs(:, :) = block_S(:, :)
308 76 : CALL diamat_all(S_evecs, S_evals)
309 :
310 : ! use eigenvalues of S_block as descriptor
311 76 : IF (PRESENT(descriptor)) THEN
312 10360 : descriptor((idesc - 1)*N + 1:idesc*N) = S_evals(:)
313 : END IF
314 :
315 : ! FORCES ----------------------------------------------------------------------------------
316 118 : IF (PRESENT(forces)) THEN
317 6 : CPASSERT(PRESENT(descr_grad))
318 6 : block_M = 0.0_dp
319 888 : DO k = 1, N
320 882 : w = descr_grad((idesc - 1)*N + k)
321 38510772 : block_M(:, :) = block_M(:, :) + w*MATMUL(S_evecs(:, k:k), TRANSPOSE(S_evecs(:, k:k)))
322 : END DO
323 :
324 20 : DO j = 1, MIN(natoms, N)
325 54 : DO k = 1, MIN(natoms, N)
326 34 : jatom = neighbor_order(j)
327 34 : katom = neighbor_order(k)
328 :
329 : ! get weigths and betas
330 34 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
331 34 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=jpao_descriptors)
332 34 : CALL get_atomic_kind(particle_set(katom)%atomic_kind, kind_number=kkind)
333 34 : CALL get_qs_kind(qs_kind_set(kkind), pao_descriptors=kpao_descriptors)
334 34 : jweight = jpao_descriptors(idesc)%weight
335 34 : jbeta = jpao_descriptors(idesc)%beta
336 34 : kweight = kpao_descriptors(idesc)%weight
337 34 : kbeta = kpao_descriptors(idesc)%beta
338 34 : beta_sum = sbeta + jbeta + kbeta
339 :
340 : ! get distances
341 136 : Rj = particle_set(jatom)%r
342 136 : Rk = particle_set(katom)%r
343 34 : Rij = pbc(Ri, Rj, cell)
344 34 : Rik = pbc(Ri, Rk, cell)
345 34 : Rjk = pbc(Rj, Rk, cell)
346 136 : Rij2 = SUM(Rij**2)
347 136 : Rik2 = SUM(Rik**2)
348 136 : Rjk2 = SUM(Rjk**2)
349 :
350 : ! calculate integral over three Gaussians
351 34 : exponent = -(sbeta*jbeta*Rij2 + sbeta*kbeta*Rik2 + jbeta*kbeta*Rjk2)/beta_sum
352 34 : integral = EXP(exponent)*rootpi/SQRT(beta_sum)
353 34 : normalization = SQRT(jbeta*kbeta)/rootpi**2
354 34 : deriv = 2.0_dp/beta_sum*block_M(j, k)
355 34 : w = jweight*kweight*normalization*integral*deriv
356 136 : forces(iatom, :) = forces(iatom, :) - sbeta*jbeta*Rij*w
357 136 : forces(jatom, :) = forces(jatom, :) + sbeta*jbeta*Rij*w
358 136 : forces(iatom, :) = forces(iatom, :) - sbeta*kbeta*Rik*w
359 136 : forces(katom, :) = forces(katom, :) + sbeta*kbeta*Rik*w
360 136 : forces(jatom, :) = forces(jatom, :) - jbeta*kbeta*Rjk*w
361 218 : forces(katom, :) = forces(katom, :) + jbeta*kbeta*Rjk*w
362 : END DO
363 : END DO
364 : END IF
365 : END DO
366 :
367 42 : CALL timestop(handle)
368 84 : END SUBROUTINE calc_descriptor_overlap
369 :
370 : ! **************************************************************************************************
371 : !> \brief Calculates a descriptor based on distance between two atoms
372 : !> \param particle_set ...
373 : !> \param qs_kind_set ...
374 : !> \param cell ...
375 : !> \param iatom ...
376 : !> \param descriptor ...
377 : !> \param descr_grad ...
378 : !> \param forces ...
379 : ! **************************************************************************************************
380 112 : SUBROUTINE calc_descriptor_r12(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
381 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
382 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
383 : TYPE(cell_type), POINTER :: cell
384 : INTEGER, INTENT(IN) :: iatom
385 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: descriptor
386 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: descr_grad
387 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
388 :
389 : REAL(dp), DIMENSION(3) :: G, R1, R12, R2
390 :
391 112 : CPASSERT(SIZE(particle_set) == 2)
392 :
393 : MARK_USED(qs_kind_set)
394 : MARK_USED(iatom)
395 : MARK_USED(cell)
396 :
397 448 : R1 = particle_set(1)%r
398 448 : R2 = particle_set(2)%r
399 112 : R12 = pbc(R1, R2, cell)
400 :
401 112 : IF (PRESENT(descriptor)) THEN
402 104 : ALLOCATE (descriptor(1))
403 416 : descriptor(1) = SQRT(SUM(R12**2))
404 : END IF
405 :
406 112 : IF (PRESENT(forces)) THEN
407 8 : CPASSERT(PRESENT(descr_grad))
408 56 : G = R12/SQRT(SUM(R12**2))*descr_grad(1)
409 32 : forces(1, :) = forces(1, :) + G
410 32 : forces(2, :) = forces(2, :) - G
411 : END IF
412 112 : END SUBROUTINE calc_descriptor_r12
413 :
414 932 : END MODULE pao_ml_descriptor
|
