Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 3038 : 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) &
142 0 : CPABORT("Not all KINDs have the same number of PAO_DESCRIPTOR sections")
143 98 : weight = pao_descriptors(idesc)%weight
144 98 : beta = pao_descriptors(idesc)%beta
145 392 : CALL pao_calc_gaussian(basis_set, block_V=block_V, Rab=Rab, lpot=0, beta=beta, weight=weight)
146 : END DO
147 :
148 : ! diagonalize block_V
149 3038 : V_evecs(:, :) = block_V(:, :)
150 98 : CALL diamat_all(V_evecs, V_evals)
151 :
152 : ! use eigenvalues of V_block as descriptor
153 98 : IF (PRESENT(descriptor)) &
154 528 : descriptor((idesc - 1)*N + 1:idesc*N) = V_evals(:)
155 :
156 : ! FORCES ----------------------------------------------------------------------------------
157 152 : IF (PRESENT(forces)) THEN
158 10 : CPASSERT(PRESENT(descr_grad))
159 310 : block_M = 0.0_dp
160 60 : DO k = 1, N
161 50 : w = descr_grad((idesc - 1)*N + k)
162 5060 : block_M(:, :) = block_M(:, :) + w*MATMUL(V_evecs(:, k:k), TRANSPOSE(V_evecs(:, k:k)))
163 : END DO
164 30 : DO jatom = 1, natoms
165 20 : IF (jatom == iatom) CYCLE
166 40 : Ra = particle_set(iatom)%r
167 40 : Rb = particle_set(jatom)%r
168 10 : Rab = pbc(ra, rb, cell)
169 10 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
170 10 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=pao_descriptors)
171 10 : weight = pao_descriptors(idesc)%weight
172 10 : beta = pao_descriptors(idesc)%beta
173 940 : block_D = 0.0_dp
174 10 : CALL pao_calc_gaussian(basis_set, block_D=block_D, Rab=Rab, lpot=0, beta=beta, weight=weight)
175 60 : DO i = 1, 3
176 930 : forces(iatom, i) = forces(iatom, i) - SUM(block_M*block_D(:, :, i))
177 950 : forces(jatom, i) = forces(jatom, i) + SUM(block_M*block_D(:, :, i))
178 : END DO
179 : END DO
180 : END IF
181 :
182 : END DO
183 :
184 54 : CALL timestop(handle)
185 108 : END SUBROUTINE calc_descriptor_pot
186 :
187 : ! **************************************************************************************************
188 : !> \brief Calculates a descriptor based on the eigenvalues of local overlap matrix
189 : !> \param particle_set ...
190 : !> \param qs_kind_set ...
191 : !> \param cell ...
192 : !> \param iatom ...
193 : !> \param descriptor ...
194 : !> \param descr_grad ...
195 : !> \param forces ...
196 : ! **************************************************************************************************
197 42 : SUBROUTINE calc_descriptor_overlap(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
198 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
199 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
200 : TYPE(cell_type), POINTER :: cell
201 : INTEGER, INTENT(IN) :: iatom
202 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: descriptor
203 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: descr_grad
204 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
205 :
206 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_descriptor_overlap'
207 :
208 : INTEGER :: handle, idesc, ikind, j, jatom, jkind, &
209 : k, katom, kkind, N, natoms, ndesc
210 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: neighbor_order
211 : REAL(dp) :: beta_sum, deriv, exponent, integral, jbeta, jweight, kbeta, kweight, &
212 : normalization, Rij2, Rik2, Rjk2, sbeta, screening_radius, screening_volume, w
213 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: S_evals
214 42 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: block_M, block_S, S_evecs
215 : REAL(dp), DIMENSION(3) :: Ri, Rij, Rik, Rj, Rjk, Rk
216 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: neighbor_dist
217 42 : TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: ipao_descriptors, jpao_descriptors, &
218 42 : kpao_descriptors
219 :
220 42 : CALL timeset(routineN, handle)
221 :
222 42 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
223 42 : CALL get_qs_kind(qs_kind_set(ikind), pao_descriptors=ipao_descriptors)
224 :
225 42 : natoms = SIZE(particle_set)
226 42 : ndesc = SIZE(ipao_descriptors)
227 42 : IF (ndesc == 0) CPABORT("No PAO_DESCRIPTOR section found")
228 :
229 : ! determine largest screening radius
230 42 : screening_radius = 0.0_dp
231 118 : DO idesc = 1, ndesc
232 118 : screening_radius = MAX(screening_radius, ipao_descriptors(idesc)%screening_radius)
233 : END DO
234 :
235 : ! estimate maximum number of neighbors within screening
236 42 : screening_volume = fourpi/3.0_dp*screening_radius**3
237 42 : N = INT(screening_volume/35.0_dp) ! rule of thumb
238 :
239 336 : ALLOCATE (block_S(N, N), S_evals(N), S_evecs(N, N))
240 118 : IF (PRESENT(descriptor)) ALLOCATE (descriptor(N*ndesc))
241 50 : IF (PRESENT(forces)) ALLOCATE (block_M(N, N))
242 :
243 : !find neighbors
244 : !TODO: this is a quadratic algorithm, use a neighbor-list instead
245 210 : ALLOCATE (neighbor_dist(natoms), neighbor_order(natoms))
246 168 : Ri = particle_set(iatom)%r
247 132 : DO jatom = 1, natoms
248 360 : Rj = particle_set(jatom)%r
249 90 : Rij = pbc(Ri, Rj, cell)
250 402 : neighbor_dist(jatom) = SQRT(SUM(Rij**2))
251 : END DO
252 42 : CALL sort(neighbor_dist, natoms, neighbor_order)
253 42 : CPASSERT(neighbor_order(1) == iatom) !central atom should be closesd to itself
254 :
255 : ! check if N was chosen large enough
256 42 : IF (natoms > N) THEN
257 0 : IF (neighbor_dist(N + 1) < screening_radius) &
258 0 : CPABORT("PAO heuristic for descriptor size broke down")
259 : END IF
260 :
261 118 : DO idesc = 1, ndesc
262 76 : sbeta = ipao_descriptors(idesc)%screening
263 :
264 : ! construct matrix S_block from neighboring atoms
265 1653532 : block_S = 0.0_dp
266 234 : DO j = 1, MIN(natoms, N)
267 568 : DO k = 1, MIN(natoms, N)
268 334 : jatom = neighbor_order(j)
269 334 : katom = neighbor_order(k)
270 :
271 : ! get weigths and betas
272 334 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
273 334 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=jpao_descriptors)
274 334 : CALL get_atomic_kind(particle_set(katom)%atomic_kind, kind_number=kkind)
275 334 : CALL get_qs_kind(qs_kind_set(kkind), pao_descriptors=kpao_descriptors)
276 334 : IF (SIZE(jpao_descriptors) /= ndesc .OR. SIZE(kpao_descriptors) /= ndesc) &
277 0 : CPABORT("Not all KINDs have the same number of PAO_DESCRIPTOR sections")
278 334 : jweight = jpao_descriptors(idesc)%weight
279 334 : jbeta = jpao_descriptors(idesc)%beta
280 334 : kweight = kpao_descriptors(idesc)%weight
281 334 : kbeta = kpao_descriptors(idesc)%beta
282 334 : beta_sum = sbeta + jbeta + kbeta
283 :
284 : ! get distances
285 1336 : Rj = particle_set(jatom)%r
286 1336 : Rk = particle_set(katom)%r
287 334 : Rij = pbc(Ri, Rj, cell)
288 334 : Rik = pbc(Ri, Rk, cell)
289 334 : Rjk = pbc(Rj, Rk, cell)
290 1336 : Rij2 = SUM(Rij**2)
291 1336 : Rik2 = SUM(Rik**2)
292 1336 : Rjk2 = SUM(Rjk**2)
293 :
294 : ! calculate integral over three Gaussians
295 334 : exponent = -(sbeta*jbeta*Rij2 + sbeta*kbeta*Rik2 + jbeta*kbeta*Rjk2)/beta_sum
296 334 : integral = EXP(exponent)*rootpi/SQRT(beta_sum)
297 334 : normalization = SQRT(jbeta*kbeta)/rootpi**2
298 1160 : block_S(j, k) = jweight*kweight*normalization*integral
299 : END DO
300 : END DO
301 :
302 : ! diagonalize V_block
303 1653532 : S_evecs(:, :) = block_S(:, :)
304 76 : CALL diamat_all(S_evecs, S_evals)
305 :
306 : ! use eigenvalues of S_block as descriptor
307 76 : IF (PRESENT(descriptor)) &
308 10360 : descriptor((idesc - 1)*N + 1:idesc*N) = S_evals(:)
309 :
310 : ! FORCES ----------------------------------------------------------------------------------
311 118 : IF (PRESENT(forces)) THEN
312 6 : CPASSERT(PRESENT(descr_grad))
313 130542 : block_M = 0.0_dp
314 888 : DO k = 1, N
315 882 : w = descr_grad((idesc - 1)*N + k)
316 57699564 : block_M(:, :) = block_M(:, :) + w*MATMUL(S_evecs(:, k:k), TRANSPOSE(S_evecs(:, k:k)))
317 : END DO
318 :
319 20 : DO j = 1, MIN(natoms, N)
320 54 : DO k = 1, MIN(natoms, N)
321 34 : jatom = neighbor_order(j)
322 34 : katom = neighbor_order(k)
323 :
324 : ! get weigths and betas
325 34 : CALL get_atomic_kind(particle_set(jatom)%atomic_kind, kind_number=jkind)
326 34 : CALL get_qs_kind(qs_kind_set(jkind), pao_descriptors=jpao_descriptors)
327 34 : CALL get_atomic_kind(particle_set(katom)%atomic_kind, kind_number=kkind)
328 34 : CALL get_qs_kind(qs_kind_set(kkind), pao_descriptors=kpao_descriptors)
329 34 : jweight = jpao_descriptors(idesc)%weight
330 34 : jbeta = jpao_descriptors(idesc)%beta
331 34 : kweight = kpao_descriptors(idesc)%weight
332 34 : kbeta = kpao_descriptors(idesc)%beta
333 34 : beta_sum = sbeta + jbeta + kbeta
334 :
335 : ! get distances
336 136 : Rj = particle_set(jatom)%r
337 136 : Rk = particle_set(katom)%r
338 34 : Rij = pbc(Ri, Rj, cell)
339 34 : Rik = pbc(Ri, Rk, cell)
340 34 : Rjk = pbc(Rj, Rk, cell)
341 136 : Rij2 = SUM(Rij**2)
342 136 : Rik2 = SUM(Rik**2)
343 136 : Rjk2 = SUM(Rjk**2)
344 :
345 : ! calculate integral over three Gaussians
346 34 : exponent = -(sbeta*jbeta*Rij2 + sbeta*kbeta*Rik2 + jbeta*kbeta*Rjk2)/beta_sum
347 34 : integral = EXP(exponent)*rootpi/SQRT(beta_sum)
348 34 : normalization = SQRT(jbeta*kbeta)/rootpi**2
349 34 : deriv = 2.0_dp/beta_sum*block_M(j, k)
350 34 : w = jweight*kweight*normalization*integral*deriv
351 136 : forces(iatom, :) = forces(iatom, :) - sbeta*jbeta*Rij*w
352 136 : forces(jatom, :) = forces(jatom, :) + sbeta*jbeta*Rij*w
353 136 : forces(iatom, :) = forces(iatom, :) - sbeta*kbeta*Rik*w
354 136 : forces(katom, :) = forces(katom, :) + sbeta*kbeta*Rik*w
355 136 : forces(jatom, :) = forces(jatom, :) - jbeta*kbeta*Rjk*w
356 218 : forces(katom, :) = forces(katom, :) + jbeta*kbeta*Rjk*w
357 : END DO
358 : END DO
359 : END IF
360 : END DO
361 :
362 42 : CALL timestop(handle)
363 84 : END SUBROUTINE calc_descriptor_overlap
364 :
365 : ! **************************************************************************************************
366 : !> \brief Calculates a descriptor based on distance between two atoms
367 : !> \param particle_set ...
368 : !> \param qs_kind_set ...
369 : !> \param cell ...
370 : !> \param iatom ...
371 : !> \param descriptor ...
372 : !> \param descr_grad ...
373 : !> \param forces ...
374 : ! **************************************************************************************************
375 112 : SUBROUTINE calc_descriptor_r12(particle_set, qs_kind_set, cell, iatom, descriptor, descr_grad, forces)
376 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
377 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
378 : TYPE(cell_type), POINTER :: cell
379 : INTEGER, INTENT(IN) :: iatom
380 : REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: descriptor
381 : REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: descr_grad
382 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: forces
383 :
384 : REAL(dp), DIMENSION(3) :: G, R1, R12, R2
385 :
386 112 : CPASSERT(SIZE(particle_set) == 2)
387 :
388 : MARK_USED(qs_kind_set)
389 : MARK_USED(iatom)
390 : MARK_USED(cell)
391 :
392 448 : R1 = particle_set(1)%r
393 448 : R2 = particle_set(2)%r
394 112 : R12 = pbc(R1, R2, cell)
395 :
396 112 : IF (PRESENT(descriptor)) THEN
397 104 : ALLOCATE (descriptor(1))
398 416 : descriptor(1) = SQRT(SUM(R12**2))
399 : END IF
400 :
401 112 : IF (PRESENT(forces)) THEN
402 8 : CPASSERT(PRESENT(descr_grad))
403 56 : G = R12/SQRT(SUM(R12**2))*descr_grad(1)
404 32 : forces(1, :) = forces(1, :) + G
405 32 : forces(2, :) = forces(2, :) - G
406 : END IF
407 112 : END SUBROUTINE calc_descriptor_r12
408 :
409 932 : END MODULE pao_ml_descriptor
|
