(git:374b731)
Loading...
Searching...
No Matches
nnp_acsf.F
Go to the documentation of this file.
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 Functionality for atom centered symmetry functions
10!> for neural network potentials
11!> \author Christoph Schran (christoph.schran@rub.de)
12!> \date 2020-10-10
13! **************************************************************************************************
14MODULE nnp_acsf
15 USE cell_types, ONLY: cell_type,&
16 pbc
17 USE cp_log_handling, ONLY: cp_get_default_logger,&
18 cp_logger_get_default_unit_nr,&
19 cp_logger_type
20 USE kinds, ONLY: default_string_length,&
21 dp
22 USE mathconstants, ONLY: pi
23 USE message_passing, ONLY: mp_para_env_type
24 USE nnp_environment_types, ONLY: nnp_acsf_ang_type,&
25 nnp_acsf_rad_type,&
26 nnp_cut_cos,&
27 nnp_cut_tanh,&
28 nnp_env_get,&
29 nnp_neighbor_type,&
30 nnp_type
31 USE periodic_table, ONLY: get_ptable_info
32#include "./base/base_uses.f90"
33
34 IMPLICIT NONE
35
36 PRIVATE
37
38 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
39 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_acsf'
40 ! Public subroutines ***
41 PUBLIC :: nnp_calc_acsf, &
42 nnp_init_acsf_groups, &
43 nnp_sort_acsf, &
44 nnp_sort_ele, &
45 nnp_write_acsf
46
47CONTAINS
48
49! **************************************************************************************************
50!> \brief Calculate atom centered symmetry functions for given atom i
51!> \param nnp ...
52!> \param i ...
53!> \param dsymdxyz ...
54!> \param stress ...
55!> \date 2020-10-10
56!> \author Christoph Schran (christoph.schran@rub.de)
57! **************************************************************************************************
58 SUBROUTINE nnp_calc_acsf(nnp, i, dsymdxyz, stress)
59
60 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
61 INTEGER, INTENT(IN) :: i
62 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT), &
63 OPTIONAL :: dsymdxyz, stress
64
65 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_calc_acsf'
66
67 INTEGER :: handle, handle_sf, ind, j, jj, k, kk, l, &
68 m, off, s, sf
69 REAL(kind=dp) :: r1, r2, r3
70 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: symtmp
71 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: forcetmp
72 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: force3tmp
73 REAL(kind=dp), DIMENSION(3) :: rvect1, rvect2, rvect3
74 TYPE(nnp_neighbor_type) :: neighbor
75
76 CALL timeset(routinen, handle)
77
78 !determine index of atom type
79 ind = nnp%ele_ind(i)
80
81 ! compute neighbors of atom i
82 CALL nnp_neighbor_create(nnp, ind, nnp%num_atoms, neighbor)
83 CALL nnp_compute_neighbors(nnp, neighbor, i)
84
85 ! Reset y:
86 nnp%rad(ind)%y = 0.0_dp
87 nnp%ang(ind)%y = 0.0_dp
88
89 !calc forces
90 IF (PRESENT(dsymdxyz)) THEN
91 !loop over radial sym fnct grps
92 CALL timeset('nnp_acsf_radial', handle_sf)
93 DO s = 1, nnp%rad(ind)%n_symfgrp
94 ALLOCATE (symtmp(nnp%rad(ind)%symfgrp(s)%n_symf))
95 ALLOCATE (forcetmp(3, nnp%rad(ind)%symfgrp(s)%n_symf))
96 !loop over associated neighbors
97 DO j = 1, neighbor%n_rad(s)
98 rvect1 = neighbor%dist_rad(1:3, j, s)
99 r1 = neighbor%dist_rad(4, j, s)
100 CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, symtmp, forcetmp)
101 jj = neighbor%ind_rad(j, s)
102 DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
103 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
104 ! update forces into dsymdxyz
105 DO l = 1, 3
106 dsymdxyz(l, m, i) = dsymdxyz(l, m, i) + forcetmp(l, sf)
107 dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) - forcetmp(l, sf)
108 END DO
109 IF (PRESENT(stress)) THEN
110 DO l = 1, 3
111 stress(:, l, m) = stress(:, l, m) + rvect1(:)*forcetmp(l, sf)
112 END DO
113 END IF
114 nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
115 END DO
116 END DO
117 DEALLOCATE (symtmp)
118 DEALLOCATE (forcetmp)
119 END DO
120 CALL timestop(handle_sf)
121 !loop over angular sym fnct grps
122 CALL timeset('nnp_acsf_angular', handle_sf)
123 off = nnp%n_rad(ind)
124 DO s = 1, nnp%ang(ind)%n_symfgrp
125 ALLOCATE (symtmp(nnp%ang(ind)%symfgrp(s)%n_symf))
126 ALLOCATE (force3tmp(3, 3, nnp%ang(ind)%symfgrp(s)%n_symf))
127 !loop over associated neighbors
128 IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
129 DO j = 1, neighbor%n_ang1(s)
130 rvect1 = neighbor%dist_ang1(1:3, j, s)
131 r1 = neighbor%dist_ang1(4, j, s)
132 DO k = j + 1, neighbor%n_ang1(s)
133 rvect2 = neighbor%dist_ang1(1:3, k, s)
134 r2 = neighbor%dist_ang1(4, k, s)
135 rvect3 = rvect2 - rvect1
136 r3 = norm2(rvect3(:))
137 IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
138 jj = neighbor%ind_ang1(j, s)
139 kk = neighbor%ind_ang1(k, s)
140 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
141 r1, r2, r3, symtmp, force3tmp)
142 DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
143 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
144 ! update forces into dsymdxy
145 DO l = 1, 3
146 dsymdxyz(l, m, i) = dsymdxyz(l, m, i) &
147 + force3tmp(l, 1, sf)
148 dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) &
149 + force3tmp(l, 2, sf)
150 dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
151 + force3tmp(l, 3, sf)
152 END DO
153 IF (PRESENT(stress)) THEN
154 DO l = 1, 3
155 stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
156 stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
157 END DO
158 END IF
159 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
160 END DO
161 END IF
162 END DO
163 END DO
164 ELSE
165 DO j = 1, neighbor%n_ang1(s)
166 rvect1 = neighbor%dist_ang1(1:3, j, s)
167 r1 = neighbor%dist_ang1(4, j, s)
168 DO k = 1, neighbor%n_ang2(s)
169 rvect2 = neighbor%dist_ang2(1:3, k, s)
170 r2 = neighbor%dist_ang2(4, k, s)
171 rvect3 = rvect2 - rvect1
172 r3 = norm2(rvect3(:))
173 IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
174 jj = neighbor%ind_ang1(j, s)
175 kk = neighbor%ind_ang1(k, s)
176 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, &
177 r1, r2, r3, symtmp, force3tmp)
178 !loop over associated sym fncts
179 DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
180 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
181 jj = neighbor%ind_ang1(j, s)
182 kk = neighbor%ind_ang2(k, s)
183 ! update forces into dsymdxy
184 DO l = 1, 3
185 dsymdxyz(l, m, i) = dsymdxyz(l, m, i) &
186 + force3tmp(l, 1, sf)
187 dsymdxyz(l, m, jj) = dsymdxyz(l, m, jj) &
188 + force3tmp(l, 2, sf)
189 dsymdxyz(l, m, kk) = dsymdxyz(l, m, kk) &
190 + force3tmp(l, 3, sf)
191 END DO
192 IF (PRESENT(stress)) THEN
193 DO l = 1, 3
194 stress(:, l, m) = stress(:, l, m) - rvect1(:)*force3tmp(l, 2, sf)
195 stress(:, l, m) = stress(:, l, m) - rvect2(:)*force3tmp(l, 3, sf)
196 END DO
197 END IF
198 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
199 END DO
200 END IF
201 END DO
202 END DO
203 END IF
204 DEALLOCATE (symtmp)
205 DEALLOCATE (force3tmp)
206 END DO
207 CALL timestop(handle_sf)
208 !no forces
209 ELSE
210 !loop over radial sym fnct grps
211 CALL timeset('nnp_acsf_radial', handle_sf)
212 DO s = 1, nnp%rad(ind)%n_symfgrp
213 ALLOCATE (symtmp(nnp%rad(ind)%symfgrp(s)%n_symf))
214 !loop over associated neighbors
215 DO j = 1, neighbor%n_rad(s)
216 rvect1 = neighbor%dist_rad(1:3, j, s)
217 r1 = neighbor%dist_rad(4, j, s)
218 CALL nnp_calc_rad(nnp, ind, s, rvect1, r1, symtmp)
219 DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
220 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
221 nnp%rad(ind)%y(m) = nnp%rad(ind)%y(m) + symtmp(sf)
222 END DO
223 END DO
224 DEALLOCATE (symtmp)
225 END DO
226 CALL timestop(handle_sf)
227 !loop over angular sym fnct grps
228 CALL timeset('nnp_acsf_angular', handle_sf)
229 off = nnp%n_rad(ind)
230 DO s = 1, nnp%ang(ind)%n_symfgrp
231 ALLOCATE (symtmp(nnp%ang(ind)%symfgrp(s)%n_symf))
232 !loop over associated neighbors
233 IF (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2)) THEN
234 DO j = 1, neighbor%n_ang1(s)
235 rvect1 = neighbor%dist_ang1(1:3, j, s)
236 r1 = neighbor%dist_ang1(4, j, s)
237 DO k = j + 1, neighbor%n_ang1(s)
238 rvect2 = neighbor%dist_ang1(1:3, k, s)
239 r2 = neighbor%dist_ang1(4, k, s)
240 rvect3 = rvect2 - rvect1
241 r3 = norm2(rvect3(:))
242 IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
243 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, symtmp)
244 DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
245 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
246 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
247 END DO
248 END IF
249 END DO
250 END DO
251 ELSE
252 DO j = 1, neighbor%n_ang1(s)
253 rvect1 = neighbor%dist_ang1(1:3, j, s)
254 r1 = neighbor%dist_ang1(4, j, s)
255 DO k = 1, neighbor%n_ang2(s)
256 rvect2 = neighbor%dist_ang2(1:3, k, s)
257 r2 = neighbor%dist_ang2(4, k, s)
258 rvect3 = rvect2 - rvect1
259 r3 = norm2(rvect3(:))
260 IF (r3 < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
261 CALL nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, symtmp)
262 !loop over associated sym fncts
263 DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
264 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
265 nnp%ang(ind)%y(m - off) = nnp%ang(ind)%y(m - off) + symtmp(sf)
266 END DO
267 END IF
268 END DO
269 END DO
270 END IF
271 DEALLOCATE (symtmp)
272 END DO
273 CALL timestop(handle_sf)
274 END IF
275
276 !check extrapolation
277 CALL nnp_check_extrapolation(nnp, ind)
278
279 IF (PRESENT(dsymdxyz)) THEN
280 IF (PRESENT(stress)) THEN
281 CALL nnp_scale_acsf(nnp, ind, dsymdxyz, stress)
282 ELSE
283 CALL nnp_scale_acsf(nnp, ind, dsymdxyz)
284 END IF
285 ELSE
286 CALL nnp_scale_acsf(nnp, ind)
287 END IF
288
289 CALL nnp_neighbor_release(neighbor)
290 CALL timestop(handle)
291
292 END SUBROUTINE nnp_calc_acsf
293
294! **************************************************************************************************
295!> \brief Check if the nnp is extrapolating
296!> \param nnp ...
297!> \param ind ...
298!> \date 2020-10-10
299!> \author Christoph Schran (christoph.schran@rub.de)
300! **************************************************************************************************
301 SUBROUTINE nnp_check_extrapolation(nnp, ind)
302
303 TYPE(nnp_type), INTENT(INOUT) :: nnp
304 INTEGER, INTENT(IN) :: ind
305
306 REAL(kind=dp), PARAMETER :: threshold = 0.0001_dp
307
308 INTEGER :: j
309 LOGICAL :: extrapolate
310
311 extrapolate = nnp%output_expol
312
313 DO j = 1, nnp%n_rad(ind)
314 IF (nnp%rad(ind)%y(j) - &
315 nnp%rad(ind)%loc_max(j) > threshold) THEN
316 extrapolate = .true.
317 ELSE IF (-nnp%rad(ind)%y(j) + &
318 nnp%rad(ind)%loc_min(j) > threshold) THEN
319 extrapolate = .true.
320 END IF
321 END DO
322 DO j = 1, nnp%n_ang(ind)
323 IF (nnp%ang(ind)%y(j) - &
324 nnp%ang(ind)%loc_max(j) > threshold) THEN
325 extrapolate = .true.
326 ELSE IF (-nnp%ang(ind)%y(j) + &
327 nnp%ang(ind)%loc_min(j) > threshold) THEN
328 extrapolate = .true.
329 END IF
330 END DO
331
332 nnp%output_expol = extrapolate
333
334 END SUBROUTINE nnp_check_extrapolation
335
336! **************************************************************************************************
337!> \brief Scale and center symetry functions (and gradients)
338!> \param nnp ...
339!> \param ind ...
340!> \param dsymdxyz ...
341!> \param stress ...
342!> \date 2020-10-10
343!> \author Christoph Schran (christoph.schran@rub.de)
344! **************************************************************************************************
345 SUBROUTINE nnp_scale_acsf(nnp, ind, dsymdxyz, stress)
346
347 TYPE(nnp_type), INTENT(INOUT) :: nnp
348 INTEGER, INTENT(IN) :: ind
349 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT), &
350 OPTIONAL :: dsymdxyz, stress
351
352 INTEGER :: j, k, off
353
354 IF (nnp%center_acsf) THEN
355 DO j = 1, nnp%n_rad(ind)
356 nnp%arc(ind)%layer(1)%node(j) = &
357 (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j))
358 END DO
359 off = nnp%n_rad(ind)
360 DO j = 1, nnp%n_ang(ind)
361 nnp%arc(ind)%layer(1)%node(j + off) = &
362 (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j))
363 END DO
364
365 IF (nnp%scale_acsf) THEN
366 DO j = 1, nnp%n_rad(ind)
367 nnp%arc(ind)%layer(1)%node(j) = &
368 nnp%arc(ind)%layer(1)%node(j)/ &
369 (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))* &
370 (nnp%scmax - nnp%scmin) + nnp%scmin
371 END DO
372 off = nnp%n_rad(ind)
373 DO j = 1, nnp%n_ang(ind)
374 nnp%arc(ind)%layer(1)%node(j + off) = &
375 nnp%arc(ind)%layer(1)%node(j + off)/ &
376 (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))* &
377 (nnp%scmax - nnp%scmin) + nnp%scmin
378 END DO
379 END IF
380 ELSE IF (nnp%scale_acsf) THEN
381 DO j = 1, nnp%n_rad(ind)
382 nnp%arc(ind)%layer(1)%node(j) = &
383 (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_min(j))/ &
384 (nnp%rad(ind)%loc_max(j) - nnp%rad(ind)%loc_min(j))* &
385 (nnp%scmax - nnp%scmin) + nnp%scmin
386 END DO
387 off = nnp%n_rad(ind)
388 DO j = 1, nnp%n_ang(ind)
389 nnp%arc(ind)%layer(1)%node(j + off) = &
390 (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_min(j))/ &
391 (nnp%ang(ind)%loc_max(j) - nnp%ang(ind)%loc_min(j))* &
392 (nnp%scmax - nnp%scmin) + nnp%scmin
393 END DO
394 ELSE IF (nnp%scale_sigma_acsf) THEN
395 DO j = 1, nnp%n_rad(ind)
396 nnp%arc(ind)%layer(1)%node(j) = &
397 (nnp%rad(ind)%y(j) - nnp%rad(ind)%loc_av(j))/ &
398 nnp%rad(ind)%sigma(j)* &
399 (nnp%scmax - nnp%scmin) + nnp%scmin
400 END DO
401 off = nnp%n_rad(ind)
402 DO j = 1, nnp%n_ang(ind)
403 nnp%arc(ind)%layer(1)%node(j + off) = &
404 (nnp%ang(ind)%y(j) - nnp%ang(ind)%loc_av(j))/ &
405 nnp%ang(ind)%sigma(j)* &
406 (nnp%scmax - nnp%scmin) + nnp%scmin
407 END DO
408 ELSE
409 DO j = 1, nnp%n_rad(ind)
410 nnp%arc(ind)%layer(1)%node(j) = nnp%rad(ind)%y(j)
411 END DO
412 off = nnp%n_rad(ind)
413 DO j = 1, nnp%n_ang(ind)
414 nnp%arc(ind)%layer(1)%node(j + off) = nnp%ang(ind)%y(j)
415 END DO
416 END IF
417
418 IF (PRESENT(dsymdxyz)) THEN
419 IF (nnp%scale_acsf) THEN
420 DO k = 1, nnp%num_atoms
421 DO j = 1, nnp%n_rad(ind)
422 dsymdxyz(:, j, k) = dsymdxyz(:, j, k)/ &
423 (nnp%rad(ind)%loc_max(j) - &
424 nnp%rad(ind)%loc_min(j))* &
425 (nnp%scmax - nnp%scmin)
426 END DO
427 END DO
428 off = nnp%n_rad(ind)
429 DO k = 1, nnp%num_atoms
430 DO j = 1, nnp%n_ang(ind)
431 dsymdxyz(:, j + off, k) = dsymdxyz(:, j + off, k)/ &
432 (nnp%ang(ind)%loc_max(j) - &
433 nnp%ang(ind)%loc_min(j))* &
434 (nnp%scmax - nnp%scmin)
435 END DO
436 END DO
437 ELSE IF (nnp%scale_sigma_acsf) THEN
438 DO k = 1, nnp%num_atoms
439 DO j = 1, nnp%n_rad(ind)
440 dsymdxyz(:, j, k) = dsymdxyz(:, j, k)/ &
441 nnp%rad(ind)%sigma(j)* &
442 (nnp%scmax - nnp%scmin)
443 END DO
444 END DO
445 off = nnp%n_rad(ind)
446 DO k = 1, nnp%num_atoms
447 DO j = 1, nnp%n_ang(ind)
448 dsymdxyz(:, j + off, k) = dsymdxyz(:, j + off, k)/ &
449 nnp%ang(ind)%sigma(j)* &
450 (nnp%scmax - nnp%scmin)
451 END DO
452 END DO
453 END IF
454 END IF
455
456 IF (PRESENT(stress)) THEN
457 IF (nnp%scale_acsf) THEN
458 DO j = 1, nnp%n_rad(ind)
459 stress(:, :, j) = stress(:, :, j)/ &
460 (nnp%rad(ind)%loc_max(j) - &
461 nnp%rad(ind)%loc_min(j))* &
462 (nnp%scmax - nnp%scmin)
463 END DO
464 off = nnp%n_rad(ind)
465 DO j = 1, nnp%n_ang(ind)
466 stress(:, :, j + off) = stress(:, :, j + off)/ &
467 (nnp%ang(ind)%loc_max(j) - &
468 nnp%ang(ind)%loc_min(j))* &
469 (nnp%scmax - nnp%scmin)
470 END DO
471 ELSE IF (nnp%scale_sigma_acsf) THEN
472 DO j = 1, nnp%n_rad(ind)
473 stress(:, :, j) = stress(:, :, j)/ &
474 nnp%rad(ind)%sigma(j)* &
475 (nnp%scmax - nnp%scmin)
476 END DO
477 off = nnp%n_rad(ind)
478 DO j = 1, nnp%n_ang(ind)
479 stress(:, :, j + off) = stress(:, :, j + off)/ &
480 nnp%ang(ind)%sigma(j)* &
481 (nnp%scmax - nnp%scmin)
482 END DO
483 END IF
484 END IF
485
486 END SUBROUTINE nnp_scale_acsf
487
488! **************************************************************************************************
489!> \brief Calculate radial symmetry function and gradient (optinal)
490!> for given displacment vecotr rvect of atom i and j
491!> \param nnp ...
492!> \param ind ...
493!> \param s ...
494!> \param rvect ...
495!> \param r ...
496!> \param sym ...
497!> \param force ...
498!> \date 2020-10-10
499!> \author Christoph Schran (christoph.schran@rub.de)
500! **************************************************************************************************
501 SUBROUTINE nnp_calc_rad(nnp, ind, s, rvect, r, sym, force)
502
503 TYPE(nnp_type), INTENT(IN) :: nnp
504 INTEGER, INTENT(IN) :: ind, s
505 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rvect
506 REAL(kind=dp), INTENT(IN) :: r
507 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: sym
508 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
509 OPTIONAL :: force
510
511 INTEGER :: k, sf
512 REAL(kind=dp) :: dfcutdr, dsymdr, eta, fcut, funccut, rs, &
513 tmp
514 REAL(kind=dp), DIMENSION(3) :: drdx
515
516 !init
517 drdx = 0.0_dp
518 fcut = 0.0_dp
519 dfcutdr = 0.0_dp
520
521 !Calculate cutoff function and partial derivative
522 funccut = nnp%rad(ind)%symfgrp(s)%cutoff !cutoff
523 SELECT CASE (nnp%cut_type)
524 CASE (nnp_cut_cos)
525 tmp = pi*r/funccut
526 fcut = 0.5_dp*(cos(tmp) + 1.0_dp)
527 IF (PRESENT(force)) THEN
528 dfcutdr = 0.5_dp*(-sin(tmp))*(pi/funccut)
529 END IF
530 CASE (nnp_cut_tanh)
531 tmp = tanh(1.0_dp - r/funccut)
532 fcut = tmp**3
533 IF (PRESENT(force)) THEN
534 dfcutdr = (-3.0_dp/funccut)*(tmp**2 - tmp**4)
535 END IF
536 CASE DEFAULT
537 cpabort("NNP| Cutoff function unknown")
538 END SELECT
539
540 IF (PRESENT(force)) drdx(:) = rvect(:)/r
541
542 !loop over sym fncts of sym fnct group s
543 DO sf = 1, nnp%rad(ind)%symfgrp(s)%n_symf
544 k = nnp%rad(ind)%symfgrp(s)%symf(sf) !symf indice
545 eta = nnp%rad(ind)%eta(k) !eta
546 rs = nnp%rad(ind)%rs(k) !rshift
547
548 ! Calculate radial symmetry function
549 sym(sf) = exp(-eta*(r - rs)**2)
550
551 ! Calculate partial derivatives of symmetry function and distance
552 ! and combine them to obtain force
553 IF (PRESENT(force)) THEN
554 dsymdr = sym(sf)*(-2.0_dp*eta*(r - rs))
555 force(:, sf) = fcut*dsymdr*drdx(:) + sym(sf)*dfcutdr*drdx(:)
556 END IF
557
558 ! combine radial symmetry function and cutoff function
559 sym(sf) = sym(sf)*fcut
560 END DO
561
562 END SUBROUTINE nnp_calc_rad
563
564! **************************************************************************************************
565!> \brief Calculate angular symmetry function and gradient (optinal)
566!> for given displacment vectors rvect1,2,3 of atom i,j and k
567!> \param nnp ...
568!> \param ind ...
569!> \param s ...
570!> \param rvect1 ...
571!> \param rvect2 ...
572!> \param rvect3 ...
573!> \param r1 ...
574!> \param r2 ...
575!> \param r3 ...
576!> \param sym ...
577!> \param force ...
578!> \date 2020-10-10
579!> \author Christoph Schran (christoph.schran@rub.de)
580! **************************************************************************************************
581 SUBROUTINE nnp_calc_ang(nnp, ind, s, rvect1, rvect2, rvect3, r1, r2, r3, sym, force)
582
583 TYPE(nnp_type), INTENT(IN) :: nnp
584 INTEGER, INTENT(IN) :: ind, s
585 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rvect1, rvect2, rvect3
586 REAL(kind=dp), INTENT(IN) :: r1, r2, r3
587 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: sym
588 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT), &
589 OPTIONAL :: force
590
591 INTEGER :: i, m, sf
592 REAL(kind=dp) :: angular, costheta, dfcutdr1, dfcutdr2, dfcutdr3, dsymdr1, dsymdr2, dsymdr3, &
593 eta, f, fcut, fcut1, fcut2, fcut3, ftot, g, lam, pref, prefzeta, rsqr1, rsqr2, rsqr3, &
594 symtmp, tmp, tmp1, tmp2, tmp3, tmpzeta, zeta
595 REAL(kind=dp), DIMENSION(3) :: dangulardx1, dangulardx2, dangulardx3, dcosthetadx1, &
596 dcosthetadx2, dcosthetadx3, dfdx1, dfdx2, dfdx3, dgdx1, dgdx2, dgdx3, dr1dx, dr2dx, dr3dx
597
598 rsqr1 = r1**2
599 rsqr2 = r2**2
600 rsqr3 = r3**2
601
602 !init
603 dangulardx1 = 0.0_dp
604 dangulardx2 = 0.0_dp
605 dangulardx3 = 0.0_dp
606 dr1dx = 0.0_dp
607 dr2dx = 0.0_dp
608 dr3dx = 0.0_dp
609 ftot = 0.0_dp
610 dfcutdr1 = 0.0_dp
611 dfcutdr2 = 0.0_dp
612 dfcutdr3 = 0.0_dp
613
614 ! Calculate cos(theta)
615 ! use law of cosine for theta
616 ! cos(ang(r1,r2)) = (r3**2 - r1**2 - r2**2)/(-2*r1*r2)
617 ! | f | g |
618 f = (rsqr3 - rsqr1 - rsqr2)
619 g = -2.0_dp*r1*r2
620 costheta = f/g
621
622 ! Calculate cutoff function and partial derivatives
623 fcut = nnp%ang(ind)%symfgrp(s)%cutoff !cutoff
624 SELECT CASE (nnp%cut_type)
625 CASE (nnp_cut_cos)
626 tmp1 = pi*r1/fcut
627 tmp2 = pi*r2/fcut
628 tmp3 = pi*r3/fcut
629 fcut1 = 0.5_dp*(cos(tmp1) + 1.0_dp)
630 fcut2 = 0.5_dp*(cos(tmp2) + 1.0_dp)
631 fcut3 = 0.5_dp*(cos(tmp3) + 1.0_dp)
632 ftot = fcut1*fcut2*fcut3
633 IF (PRESENT(force)) THEN
634 pref = 0.5_dp*(pi/fcut)
635 dfcutdr1 = pref*(-sin(tmp1))*fcut2*fcut3
636 dfcutdr2 = pref*(-sin(tmp2))*fcut1*fcut3
637 dfcutdr3 = pref*(-sin(tmp3))*fcut1*fcut2
638 END IF
639 CASE (nnp_cut_tanh)
640 tmp1 = tanh(1.0_dp - r1/fcut)
641 tmp2 = tanh(1.0_dp - r2/fcut)
642 tmp3 = tanh(1.0_dp - r3/fcut)
643 fcut1 = tmp1**3
644 fcut2 = tmp2**3
645 fcut3 = tmp3**3
646 ftot = fcut1*fcut2*fcut3
647 IF (PRESENT(force)) THEN
648 pref = -3.0_dp/fcut
649 dfcutdr1 = pref*(tmp1**2 - tmp1**4)*fcut2*fcut3
650 dfcutdr2 = pref*(tmp2**2 - tmp2**4)*fcut1*fcut3
651 dfcutdr3 = pref*(tmp3**2 - tmp3**4)*fcut1*fcut2
652 END IF
653 CASE DEFAULT
654 cpabort("NNP| Cutoff function unknown")
655 END SELECT
656
657 IF (PRESENT(force)) THEN
658 dr1dx(:) = rvect1(:)/r1
659 dr2dx(:) = rvect2(:)/r2
660 dr3dx(:) = rvect3(:)/r3
661 END IF
662
663 !loop over associated sym fncts
664 DO sf = 1, nnp%ang(ind)%symfgrp(s)%n_symf
665 m = nnp%ang(ind)%symfgrp(s)%symf(sf)
666 lam = nnp%ang(ind)%lam(m) !lambda
667 zeta = nnp%ang(ind)%zeta(m) !zeta
668 prefzeta = nnp%ang(ind)%prefzeta(m) ! 2**(1-zeta)
669 eta = nnp%ang(ind)%eta(m) !eta
670
671 tmp = (1.0_dp + lam*costheta)
672
673 IF (tmp <= 0.0_dp) THEN
674 sym(sf) = 0.0_dp
675 IF (PRESENT(force)) force(:, :, sf) = 0.0_dp
676 cycle
677 END IF
678
679 ! Calculate symmetry function
680 ! Calculate angular symmetry function
681 ! ang = (1+lam*cos(theta_ijk))**zeta
682 i = nint(zeta)
683 IF (1.0_dp*i == zeta) THEN
684 tmpzeta = tmp**(i - 1) ! integer power is a LOT faster
685 ELSE
686 tmpzeta = tmp**(zeta - 1.0_dp)
687 END IF
688 angular = tmpzeta*tmp
689 ! exponential part:
690 ! exp(-eta*(r1**2+r2**2+r3**2))
691 symtmp = exp(-eta*(rsqr1 + rsqr2 + rsqr3))
692
693 ! Partial derivatives (f/g)' = (f'g - fg')/g^2
694 IF (PRESENT(force)) THEN
695 pref = zeta*(tmpzeta)/g**2
696 DO i = 1, 3
697 dfdx1(i) = -2.0_dp*lam*(rvect1(i) + rvect2(i))
698 dfdx2(i) = 2.0_dp*lam*(rvect3(i) + rvect1(i))
699 dfdx3(i) = 2.0_dp*lam*(rvect2(i) - rvect3(i))
700
701 tmp1 = 2.0_dp*r2*dr1dx(i)
702 tmp2 = 2.0_dp*r1*dr2dx(i)
703 dgdx1(i) = -(tmp1 + tmp2)
704 dgdx2(i) = tmp1
705 dgdx3(i) = tmp2
706
707 dcosthetadx1(i) = dfdx1(i)*g - lam*f*dgdx1(i)
708 dcosthetadx2(i) = dfdx2(i)*g - lam*f*dgdx2(i)
709 dcosthetadx3(i) = dfdx3(i)*g - lam*f*dgdx3(i)
710
711 dangulardx1(i) = pref*dcosthetadx1(i)
712 dangulardx2(i) = pref*dcosthetadx2(i)
713 dangulardx3(i) = pref*dcosthetadx3(i)
714 END DO
715
716 ! Calculate partial derivatives of exponential part and distance
717 ! and combine partial derivatives to obtain force
718 pref = prefzeta
719 tmp = -2.0_dp*symtmp*eta
720 dsymdr1 = tmp*r1
721 dsymdr2 = tmp*r2
722 dsymdr3 = tmp*r3
723
724 ! G(r1,r2,r3) = pref*angular(r1,r2,r3)*exp(r1,r2,r3)*fcut(r1,r2,r3)
725 ! dG/dx1 = (dangular/dx1* exp * fcut +
726 ! angular * dexp/dx1* fcut +
727 ! angular * exp * dfcut/dx1)*pref
728 ! dr1/dx1 = -dr1/dx2
729 tmp = pref*symtmp*ftot
730 tmp1 = pref*angular*(ftot*dsymdr1 + dfcutdr1*symtmp)
731 tmp2 = pref*angular*(ftot*dsymdr2 + dfcutdr2*symtmp)
732 tmp3 = pref*angular*(ftot*dsymdr3 + dfcutdr3*symtmp)
733 DO i = 1, 3
734 force(i, 1, sf) = tmp*dangulardx1(i) + tmp1*dr1dx(i) + tmp2*dr2dx(i)
735 force(i, 2, sf) = tmp*dangulardx2(i) - tmp1*dr1dx(i) + tmp3*dr3dx(i)
736 force(i, 3, sf) = tmp*dangulardx3(i) - tmp2*dr2dx(i) - tmp3*dr3dx(i)
737 END DO
738 END IF
739
740 ! Don't forget prefactor: 2**(1-ang%zeta)
741 pref = prefzeta
742 ! combine angular and exponential part with cutoff function
743 sym(sf) = pref*angular*symtmp*ftot
744 END DO
745
746 END SUBROUTINE nnp_calc_ang
747
748! **************************************************************************************************
749!> \brief Sort element array according to atomic number
750!> \param ele ...
751!> \param nuc_ele ...
752!> \date 2020-10-10
753!> \author Christoph Schran (christoph.schran@rub.de)
754! **************************************************************************************************
755 SUBROUTINE nnp_sort_ele(ele, nuc_ele)
756 CHARACTER(len=2), DIMENSION(:), INTENT(INOUT) :: ele
757 INTEGER, DIMENSION(:), INTENT(INOUT) :: nuc_ele
758
759 CHARACTER(len=2) :: tmp_ele
760 INTEGER :: i, j, loc, minimum, tmp_nuc_ele
761
762 ! Determine atomic number
763 DO i = 1, SIZE(ele)
764 CALL get_ptable_info(ele(i), number=nuc_ele(i))
765 END DO
766
767 ! Sort both arrays
768 DO i = 1, SIZE(ele) - 1
769 minimum = nuc_ele(i)
770 loc = i
771 DO j = i + 1, SIZE(ele)
772 IF (nuc_ele(j) .LT. minimum) THEN
773 loc = j
774 minimum = nuc_ele(j)
775 END IF
776 END DO
777 tmp_nuc_ele = nuc_ele(i)
778 nuc_ele(i) = nuc_ele(loc)
779 nuc_ele(loc) = tmp_nuc_ele
780
781 tmp_ele = ele(i)
782 ele(i) = ele(loc)
783 ele(loc) = tmp_ele
784
785 END DO
786
787 END SUBROUTINE nnp_sort_ele
788
789! **************************************************************************************************
790!> \brief Sort symmetry functions according to different criteria
791!> \param nnp ...
792!> \date 2020-10-10
793!> \author Christoph Schran (christoph.schran@rub.de)
794! **************************************************************************************************
795 SUBROUTINE nnp_sort_acsf(nnp)
796 TYPE(nnp_type), INTENT(INOUT) :: nnp
797
798 INTEGER :: i, j, k, loc
799
800 ! First sort is according to symmetry function type
801 ! This is done manually, since data structures are separate
802 ! Note: Bubble sort is OK here, since small sort + special
803 DO i = 1, nnp%n_ele
804 ! sort by cutoff
805 ! rad
806 DO j = 1, nnp%n_rad(i) - 1
807 loc = j
808 DO k = j + 1, nnp%n_rad(i)
809 IF (nnp%rad(i)%funccut(loc) .GT. nnp%rad(i)%funccut(k)) THEN
810 loc = k
811 END IF
812 END DO
813 ! swap symfnct
814 CALL nnp_swaprad(nnp%rad(i), j, loc)
815 END DO
816
817 ! sort by eta
818 ! rad
819 DO j = 1, nnp%n_rad(i) - 1
820 loc = j
821 DO k = j + 1, nnp%n_rad(i)
822 IF (nnp%rad(i)%funccut(loc) .EQ. nnp%rad(i)%funccut(k) .AND. &
823 nnp%rad(i)%eta(loc) .GT. nnp%rad(i)%eta(k)) THEN
824 loc = k
825 END IF
826 END DO
827 ! swap symfnct
828 CALL nnp_swaprad(nnp%rad(i), j, loc)
829 END DO
830
831 ! sort by rshift
832 ! rad
833 DO j = 1, nnp%n_rad(i) - 1
834 loc = j
835 DO k = j + 1, nnp%n_rad(i)
836 IF (nnp%rad(i)%funccut(loc) .EQ. nnp%rad(i)%funccut(k) .AND. &
837 nnp%rad(i)%eta(loc) .EQ. nnp%rad(i)%eta(k) .AND. &
838 nnp%rad(i)%rs(loc) .GT. nnp%rad(i)%rs(k)) THEN
839 loc = k
840 END IF
841 END DO
842 ! swap symfnct
843 CALL nnp_swaprad(nnp%rad(i), j, loc)
844 END DO
845
846 ! sort by ele
847 ! rad
848 DO j = 1, nnp%n_rad(i) - 1
849 loc = j
850 DO k = j + 1, nnp%n_rad(i)
851 IF (nnp%rad(i)%funccut(loc) .EQ. nnp%rad(i)%funccut(k) .AND. &
852 nnp%rad(i)%eta(loc) .EQ. nnp%rad(i)%eta(k) .AND. &
853 nnp%rad(i)%rs(loc) .EQ. nnp%rad(i)%rs(k) .AND. &
854 nnp%rad(i)%nuc_ele(loc) .GT. nnp%rad(i)%nuc_ele(k)) THEN
855 loc = k
856 END IF
857 END DO
858 ! swap symfnct
859 CALL nnp_swaprad(nnp%rad(i), j, loc)
860 END DO
861
862 ! ang
863 ! sort by cutoff
864 DO j = 1, nnp%n_ang(i) - 1
865 loc = j
866 DO k = j + 1, nnp%n_ang(i)
867 IF (nnp%ang(i)%funccut(loc) .GT. nnp%ang(i)%funccut(k)) THEN
868 loc = k
869 END IF
870 END DO
871 ! swap symfnct
872 CALL nnp_swapang(nnp%ang(i), j, loc)
873 END DO
874
875 ! sort by eta
876 ! ang
877 DO j = 1, nnp%n_ang(i) - 1
878 loc = j
879 DO k = j + 1, nnp%n_ang(i)
880 IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
881 nnp%ang(i)%eta(loc) .GT. nnp%ang(i)%eta(k)) THEN
882 loc = k
883 END IF
884 END DO
885 ! swap symfnct
886 CALL nnp_swapang(nnp%ang(i), j, loc)
887 END DO
888
889 ! sort by zeta
890 ! ang
891 DO j = 1, nnp%n_ang(i) - 1
892 loc = j
893 DO k = j + 1, nnp%n_ang(i)
894 IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
895 nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
896 nnp%ang(i)%zeta(loc) .GT. nnp%ang(i)%zeta(k)) THEN
897 loc = k
898 END IF
899 END DO
900 ! swap symfnct
901 CALL nnp_swapang(nnp%ang(i), j, loc)
902 END DO
903
904 ! sort by lambda
905 ! ang
906 DO j = 1, nnp%n_ang(i) - 1
907 loc = j
908 DO k = j + 1, nnp%n_ang(i)
909 IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
910 nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
911 nnp%ang(i)%zeta(loc) .EQ. nnp%ang(i)%zeta(k) .AND. &
912 nnp%ang(i)%lam(loc) .GT. nnp%ang(i)%lam(k)) THEN
913 loc = k
914 END IF
915 END DO
916 ! swap symfnct
917 CALL nnp_swapang(nnp%ang(i), j, loc)
918 END DO
919
920 ! sort by ele
921 ! ang, ele1
922 DO j = 1, nnp%n_ang(i) - 1
923 loc = j
924 DO k = j + 1, nnp%n_ang(i)
925 IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
926 nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
927 nnp%ang(i)%zeta(loc) .EQ. nnp%ang(i)%zeta(k) .AND. &
928 nnp%ang(i)%lam(loc) .EQ. nnp%ang(i)%lam(k) .AND. &
929 nnp%ang(i)%nuc_ele1(loc) .GT. nnp%ang(i)%nuc_ele1(k)) THEN
930 loc = k
931 END IF
932 END DO
933 ! swap symfnct
934 CALL nnp_swapang(nnp%ang(i), j, loc)
935 END DO
936 ! ang, ele2
937 DO j = 1, nnp%n_ang(i) - 1
938 loc = j
939 DO k = j + 1, nnp%n_ang(i)
940 IF (nnp%ang(i)%funccut(loc) .EQ. nnp%ang(i)%funccut(k) .AND. &
941 nnp%ang(i)%eta(loc) .EQ. nnp%ang(i)%eta(k) .AND. &
942 nnp%ang(i)%zeta(loc) .EQ. nnp%ang(i)%zeta(k) .AND. &
943 nnp%ang(i)%lam(loc) .EQ. nnp%ang(i)%lam(k) .AND. &
944 nnp%ang(i)%nuc_ele1(loc) .EQ. nnp%ang(i)%nuc_ele1(k) .AND. &
945 nnp%ang(i)%nuc_ele2(loc) .GT. nnp%ang(i)%nuc_ele2(k)) THEN
946 loc = k
947 END IF
948 END DO
949 ! swap symfnct
950 CALL nnp_swapang(nnp%ang(i), j, loc)
951 END DO
952 END DO
953
954 END SUBROUTINE nnp_sort_acsf
955
956! **************************************************************************************************
957!> \brief Swap two radial symmetry functions
958!> \param rad ...
959!> \param i ...
960!> \param j ...
961!> \date 2020-10-10
962!> \author Christoph Schran (christoph.schran@rub.de)
963! **************************************************************************************************
964 SUBROUTINE nnp_swaprad(rad, i, j)
965 TYPE(nnp_acsf_rad_type), INTENT(INOUT) :: rad
966 INTEGER, INTENT(IN) :: i, j
967
968 CHARACTER(len=2) :: tmpc
969 INTEGER :: tmpi
970 REAL(kind=dp) :: tmpr
971
972 tmpr = rad%funccut(i)
973 rad%funccut(i) = rad%funccut(j)
974 rad%funccut(j) = tmpr
975
976 tmpr = rad%eta(i)
977 rad%eta(i) = rad%eta(j)
978 rad%eta(j) = tmpr
979
980 tmpr = rad%rs(i)
981 rad%rs(i) = rad%rs(j)
982 rad%rs(j) = tmpr
983
984 tmpc = rad%ele(i)
985 rad%ele(i) = rad%ele(j)
986 rad%ele(j) = tmpc
987
988 tmpi = rad%nuc_ele(i)
989 rad%nuc_ele(i) = rad%nuc_ele(j)
990 rad%nuc_ele(j) = tmpi
991
992 END SUBROUTINE nnp_swaprad
993
994! **************************************************************************************************
995!> \brief Swap two angular symmetry functions
996!> \param ang ...
997!> \param i ...
998!> \param j ...
999!> \date 2020-10-10
1000!> \author Christoph Schran (christoph.schran@rub.de)
1001! **************************************************************************************************
1002 SUBROUTINE nnp_swapang(ang, i, j)
1003 TYPE(nnp_acsf_ang_type), INTENT(INOUT) :: ang
1004 INTEGER, INTENT(IN) :: i, j
1005
1006 CHARACTER(len=2) :: tmpc
1007 INTEGER :: tmpi
1008 REAL(kind=dp) :: tmpr
1009
1010 tmpr = ang%funccut(i)
1011 ang%funccut(i) = ang%funccut(j)
1012 ang%funccut(j) = tmpr
1013
1014 tmpr = ang%eta(i)
1015 ang%eta(i) = ang%eta(j)
1016 ang%eta(j) = tmpr
1017
1018 tmpr = ang%zeta(i)
1019 ang%zeta(i) = ang%zeta(j)
1020 ang%zeta(j) = tmpr
1021
1022 tmpr = ang%prefzeta(i)
1023 ang%prefzeta(i) = ang%prefzeta(j)
1024 ang%prefzeta(j) = tmpr
1025
1026 tmpr = ang%lam(i)
1027 ang%lam(i) = ang%lam(j)
1028 ang%lam(j) = tmpr
1029
1030 tmpc = ang%ele1(i)
1031 ang%ele1(i) = ang%ele1(j)
1032 ang%ele1(j) = tmpc
1033
1034 tmpi = ang%nuc_ele1(i)
1035 ang%nuc_ele1(i) = ang%nuc_ele1(j)
1036 ang%nuc_ele1(j) = tmpi
1037
1038 tmpc = ang%ele2(i)
1039 ang%ele2(i) = ang%ele2(j)
1040 ang%ele2(j) = tmpc
1041
1042 tmpi = ang%nuc_ele2(i)
1043 ang%nuc_ele2(i) = ang%nuc_ele2(j)
1044 ang%nuc_ele2(j) = tmpi
1045
1046 END SUBROUTINE nnp_swapang
1047
1048! **************************************************************************************************
1049!> \brief Initialize symmetry function groups
1050!> \param nnp ...
1051!> \date 2020-10-10
1052!> \author Christoph Schran (christoph.schran@rub.de)
1053! **************************************************************************************************
1054 SUBROUTINE nnp_init_acsf_groups(nnp)
1055
1056 TYPE(nnp_type), INTENT(INOUT) :: nnp
1057
1058 INTEGER :: ang, i, j, k, rad, s
1059 REAL(kind=dp) :: funccut
1060
1061 !find out how many symmetry functions groups are needed
1062 DO i = 1, nnp%n_ele
1063 nnp%rad(i)%n_symfgrp = 0
1064 nnp%ang(i)%n_symfgrp = 0
1065 !search radial symmetry functions
1066 DO j = 1, nnp%n_ele
1067 funccut = -1.0_dp
1068 DO s = 1, nnp%n_rad(i)
1069 IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1070 IF (abs(nnp%rad(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1071 nnp%rad(i)%n_symfgrp = nnp%rad(i)%n_symfgrp + 1
1072 funccut = nnp%rad(i)%funccut(s)
1073 END IF
1074 END IF
1075 END DO
1076 END DO
1077 !search angular symmetry functions
1078 DO j = 1, nnp%n_ele
1079 DO k = j, nnp%n_ele
1080 funccut = -1.0_dp
1081 DO s = 1, nnp%n_ang(i)
1082 IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1083 nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1084 (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1085 nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1086 IF (abs(nnp%ang(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1087 nnp%ang(i)%n_symfgrp = nnp%ang(i)%n_symfgrp + 1
1088 funccut = nnp%ang(i)%funccut(s)
1089 END IF
1090 END IF
1091 END DO
1092 END DO
1093 END DO
1094 END DO
1095
1096 !allocate memory for symmetry functions group
1097 DO i = 1, nnp%n_ele
1098 ALLOCATE (nnp%rad(i)%symfgrp(nnp%rad(i)%n_symfgrp))
1099 ALLOCATE (nnp%ang(i)%symfgrp(nnp%ang(i)%n_symfgrp))
1100 DO j = 1, nnp%rad(i)%n_symfgrp
1101 nnp%rad(i)%symfgrp(j)%n_symf = 0
1102 END DO
1103 DO j = 1, nnp%ang(i)%n_symfgrp
1104 nnp%ang(i)%symfgrp(j)%n_symf = 0
1105 END DO
1106 END DO
1107
1108 !init symmetry functions group
1109 DO i = 1, nnp%n_ele
1110 rad = 0
1111 ang = 0
1112 DO j = 1, nnp%n_ele
1113 funccut = -1.0_dp
1114 DO s = 1, nnp%n_rad(i)
1115 IF (nnp%rad(i)%ele(s) == nnp%ele(j)) THEN
1116 IF (abs(nnp%rad(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1117 rad = rad + 1
1118 funccut = nnp%rad(i)%funccut(s)
1119 nnp%rad(i)%symfgrp(rad)%cutoff = funccut
1120 ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele(1))
1121 ALLOCATE (nnp%rad(i)%symfgrp(rad)%ele_ind(1))
1122 nnp%rad(i)%symfgrp(rad)%ele(1) = nnp%ele(j)
1123 nnp%rad(i)%symfgrp(rad)%ele_ind(1) = j
1124 END IF
1125 nnp%rad(i)%symfgrp(rad)%n_symf = nnp%rad(i)%symfgrp(rad)%n_symf + 1
1126 END IF
1127 END DO
1128 END DO
1129 DO j = 1, nnp%n_ele
1130 DO k = j, nnp%n_ele
1131 funccut = -1.0_dp
1132 DO s = 1, nnp%n_ang(i)
1133 IF ((nnp%ang(i)%ele1(s) == nnp%ele(j) .AND. &
1134 nnp%ang(i)%ele2(s) == nnp%ele(k)) .OR. &
1135 (nnp%ang(i)%ele1(s) == nnp%ele(k) .AND. &
1136 nnp%ang(i)%ele2(s) == nnp%ele(j))) THEN
1137 IF (abs(nnp%ang(i)%funccut(s) - funccut) > 1.0e-5_dp) THEN
1138 ang = ang + 1
1139 funccut = nnp%ang(i)%funccut(s)
1140 nnp%ang(i)%symfgrp(ang)%cutoff = funccut
1141 ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele(2))
1142 ALLOCATE (nnp%ang(i)%symfgrp(ang)%ele_ind(2))
1143 nnp%ang(i)%symfgrp(ang)%ele(1) = nnp%ele(j)
1144 nnp%ang(i)%symfgrp(ang)%ele(2) = nnp%ele(k)
1145 nnp%ang(i)%symfgrp(ang)%ele_ind(1) = j
1146 nnp%ang(i)%symfgrp(ang)%ele_ind(2) = k
1147 END IF
1148 nnp%ang(i)%symfgrp(ang)%n_symf = nnp%ang(i)%symfgrp(ang)%n_symf + 1
1149 END IF
1150 END DO
1151 END DO
1152 END DO
1153 END DO
1154
1155 !add symmetry functions to associated groups
1156 DO i = 1, nnp%n_ele
1157 DO j = 1, nnp%rad(i)%n_symfgrp
1158 ALLOCATE (nnp%rad(i)%symfgrp(j)%symf(nnp%rad(i)%symfgrp(j)%n_symf))
1159 rad = 0
1160 DO s = 1, nnp%n_rad(i)
1161 IF (nnp%rad(i)%ele(s) == nnp%rad(i)%symfgrp(j)%ele(1)) THEN
1162 IF (abs(nnp%rad(i)%funccut(s) - nnp%rad(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
1163 rad = rad + 1
1164 nnp%rad(i)%symfgrp(j)%symf(rad) = s
1165 END IF
1166 END IF
1167 END DO
1168 END DO
1169 DO j = 1, nnp%ang(i)%n_symfgrp
1170 ALLOCATE (nnp%ang(i)%symfgrp(j)%symf(nnp%ang(i)%symfgrp(j)%n_symf))
1171 ang = 0
1172 DO s = 1, nnp%n_ang(i)
1173 IF ((nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(1) .AND. &
1174 nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(2)) .OR. &
1175 (nnp%ang(i)%ele1(s) == nnp%ang(i)%symfgrp(j)%ele(2) .AND. &
1176 nnp%ang(i)%ele2(s) == nnp%ang(i)%symfgrp(j)%ele(1))) THEN
1177 IF (abs(nnp%ang(i)%funccut(s) - nnp%ang(i)%symfgrp(j)%cutoff) <= 1.0e-5_dp) THEN
1178 ang = ang + 1
1179 nnp%ang(i)%symfgrp(j)%symf(ang) = s
1180 END IF
1181 END IF
1182 END DO
1183 END DO
1184 END DO
1185
1186 END SUBROUTINE nnp_init_acsf_groups
1187
1188! **************************************************************************************************
1189!> \brief Write symmetry function information
1190!> \param nnp ...
1191!> \param para_env ...
1192!> \param printtag ...
1193!> \date 2020-10-10
1194!> \author Christoph Schran (christoph.schran@rub.de)
1195! **************************************************************************************************
1196 SUBROUTINE nnp_write_acsf(nnp, para_env, printtag)
1197 TYPE(nnp_type), INTENT(INOUT) :: nnp
1198 TYPE(mp_para_env_type), POINTER :: para_env
1199 CHARACTER(LEN=*), INTENT(IN) :: printtag
1200
1201 CHARACTER(len=default_string_length) :: my_label
1202 INTEGER :: i, j, unit_nr
1203 TYPE(cp_logger_type), POINTER :: logger
1204
1205 NULLIFY (logger)
1206 logger => cp_get_default_logger()
1207
1208 my_label = trim(printtag)//"| "
1209 IF (para_env%is_source()) THEN
1210 unit_nr = cp_logger_get_default_unit_nr(logger)
1211 WRITE (unit_nr, '(1X,A,1X,10(I2,1X))') trim(my_label)//" Activation functions:", nnp%actfnct(:)
1212 DO i = 1, nnp%n_ele
1213 WRITE (unit_nr, *) trim(my_label)//" short range atomic symmetry functions element "// &
1214 nnp%ele(i)//":"
1215 DO j = 1, nnp%n_rad(i)
1216 WRITE (unit_nr, '(1X,A,1X,I3,1X,A2,1X,I2,1X,A2,11X,3(F6.3,1X))') trim(my_label), j, nnp%ele(i), 2, &
1217 nnp%rad(i)%ele(j), nnp%rad(i)%eta(j), &
1218 nnp%rad(i)%rs(j), nnp%rad(i)%funccut(j)
1219 END DO
1220 DO j = 1, nnp%n_ang(i)
1221 WRITE (unit_nr, '(1X,A,1X,I3,1X,A2,1X,I2,2(1X,A2),1X,4(F6.3,1X))') &
1222 trim(my_label), j, nnp%ele(i), 3, &
1223 nnp%ang(i)%ele1(j), nnp%ang(i)%ele2(j), &
1224 nnp%ang(i)%eta(j), nnp%ang(i)%lam(j), &
1225 nnp%ang(i)%zeta(j), nnp%ang(i)%funccut(j)
1226 END DO
1227 END DO
1228 END IF
1229
1230 RETURN
1231
1232 END SUBROUTINE nnp_write_acsf
1233
1234! **************************************************************************************************
1235!> \brief Create neighbor object
1236!> \param nnp ...
1237!> \param ind ...
1238!> \param nat ...
1239!> \param neighbor ...
1240!> \date 2020-10-10
1241!> \author Christoph Schran (christoph.schran@rub.de)
1242! **************************************************************************************************
1243 SUBROUTINE nnp_neighbor_create(nnp, ind, nat, neighbor)
1244
1245 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
1246 INTEGER, INTENT(IN) :: ind, nat
1247 TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1248
1249 INTEGER :: n
1250 TYPE(cell_type), POINTER :: cell
1251
1252 NULLIFY (cell)
1253 CALL nnp_env_get(nnp_env=nnp, cell=cell)
1254
1255 CALL nnp_compute_pbc_copies(neighbor%pbc_copies, cell, nnp%max_cut)
1256 n = (sum(neighbor%pbc_copies) + 1)*nat
1257 ALLOCATE (neighbor%dist_rad(4, n, nnp%rad(ind)%n_symfgrp))
1258 ALLOCATE (neighbor%dist_ang1(4, n, nnp%ang(ind)%n_symfgrp))
1259 ALLOCATE (neighbor%dist_ang2(4, n, nnp%ang(ind)%n_symfgrp))
1260 ALLOCATE (neighbor%ind_rad(n, nnp%rad(ind)%n_symfgrp))
1261 ALLOCATE (neighbor%ind_ang1(n, nnp%ang(ind)%n_symfgrp))
1262 ALLOCATE (neighbor%ind_ang2(n, nnp%ang(ind)%n_symfgrp))
1263 ALLOCATE (neighbor%n_rad(nnp%rad(ind)%n_symfgrp))
1264 ALLOCATE (neighbor%n_ang1(nnp%ang(ind)%n_symfgrp))
1265 ALLOCATE (neighbor%n_ang2(nnp%ang(ind)%n_symfgrp))
1266 neighbor%n_rad(:) = 0
1267 neighbor%n_ang1(:) = 0
1268 neighbor%n_ang2(:) = 0
1269
1270 END SUBROUTINE nnp_neighbor_create
1271
1272! **************************************************************************************************
1273!> \brief Destroy neighbor object
1274!> \param neighbor ...
1275!> \date 2020-10-10
1276!> \author Christoph Schran (christoph.schran@rub.de)
1277! **************************************************************************************************
1278 SUBROUTINE nnp_neighbor_release(neighbor)
1279
1280 TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1281
1282 DEALLOCATE (neighbor%dist_rad)
1283 DEALLOCATE (neighbor%dist_ang1)
1284 DEALLOCATE (neighbor%dist_ang2)
1285 DEALLOCATE (neighbor%ind_rad)
1286 DEALLOCATE (neighbor%ind_ang1)
1287 DEALLOCATE (neighbor%ind_ang2)
1288 DEALLOCATE (neighbor%n_rad)
1289 DEALLOCATE (neighbor%n_ang1)
1290 DEALLOCATE (neighbor%n_ang2)
1291
1292 END SUBROUTINE nnp_neighbor_release
1293
1294! **************************************************************************************************
1295!> \brief Generate neighboring list for an atomic configuration
1296!> \param nnp ...
1297!> \param neighbor ...
1298!> \param i ...
1299!> \date 2020-10-10
1300!> \author Christoph Schran (christoph.schran@rub.de)
1301! **************************************************************************************************
1302 SUBROUTINE nnp_compute_neighbors(nnp, neighbor, i)
1303
1304 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
1305 TYPE(nnp_neighbor_type), INTENT(INOUT) :: neighbor
1306 INTEGER, INTENT(IN) :: i
1307
1308 INTEGER :: c1, c2, c3, ind, j, s
1309 INTEGER, DIMENSION(3) :: nl
1310 REAL(kind=dp) :: norm
1311 REAL(kind=dp), DIMENSION(3) :: dr
1312 TYPE(cell_type), POINTER :: cell
1313
1314 NULLIFY (cell)
1315 CALL nnp_env_get(nnp_env=nnp, cell=cell)
1316
1317 ind = nnp%ele_ind(i)
1318
1319 DO j = 1, nnp%num_atoms
1320 DO c1 = 1, 2*neighbor%pbc_copies(1) + 1
1321 nl(1) = -neighbor%pbc_copies(1) + c1 - 1
1322 DO c2 = 1, 2*neighbor%pbc_copies(2) + 1
1323 nl(2) = -neighbor%pbc_copies(2) + c2 - 1
1324 DO c3 = 1, 2*neighbor%pbc_copies(3) + 1
1325 nl(3) = -neighbor%pbc_copies(3) + c3 - 1
1326 IF (j == i .AND. nl(1) == 0 .AND. nl(2) == 0 .AND. nl(3) == 0) cycle
1327 dr(:) = nnp%coord(:, i) - nnp%coord(:, j)
1328 !Apply pbc, but subtract nl boxes from periodic image
1329 dr = pbc(dr, cell, nl)
1330 norm = norm2(dr(:))
1331 DO s = 1, nnp%rad(ind)%n_symfgrp
1332 IF (nnp%ele_ind(j) == nnp%rad(ind)%symfgrp(s)%ele_ind(1)) THEN
1333 IF (norm < nnp%rad(ind)%symfgrp(s)%cutoff) THEN
1334 neighbor%n_rad(s) = neighbor%n_rad(s) + 1
1335 neighbor%ind_rad(neighbor%n_rad(s), s) = j
1336 neighbor%dist_rad(1:3, neighbor%n_rad(s), s) = dr(:)
1337 neighbor%dist_rad(4, neighbor%n_rad(s), s) = norm
1338 END IF
1339 END IF
1340 END DO
1341 DO s = 1, nnp%ang(ind)%n_symfgrp
1342 IF (norm < nnp%ang(ind)%symfgrp(s)%cutoff) THEN
1343 IF (nnp%ele_ind(j) == nnp%ang(ind)%symfgrp(s)%ele_ind(1)) THEN
1344 neighbor%n_ang1(s) = neighbor%n_ang1(s) + 1
1345 neighbor%ind_ang1(neighbor%n_ang1(s), s) = j
1346 neighbor%dist_ang1(1:3, neighbor%n_ang1(s), s) = dr(:)
1347 neighbor%dist_ang1(4, neighbor%n_ang1(s), s) = norm
1348 END IF
1349 IF (nnp%ele_ind(j) == nnp%ang(ind)%symfgrp(s)%ele_ind(2)) THEN
1350 neighbor%n_ang2(s) = neighbor%n_ang2(s) + 1
1351 neighbor%ind_ang2(neighbor%n_ang2(s), s) = j
1352 neighbor%dist_ang2(1:3, neighbor%n_ang2(s), s) = dr(:)
1353 neighbor%dist_ang2(4, neighbor%n_ang2(s), s) = norm
1354 END IF
1355 END IF
1356 END DO
1357 END DO
1358 END DO
1359 END DO
1360 END DO
1361
1362 END SUBROUTINE nnp_compute_neighbors
1363
1364! **************************************************************************************************
1365!> \brief Determine required pbc copies for small cells
1366!> \param pbc_copies ...
1367!> \param cell ...
1368!> \param cutoff ...
1369!> \date 2020-10-10
1370!> \author Christoph Schran (christoph.schran@rub.de)
1371! **************************************************************************************************
1372 SUBROUTINE nnp_compute_pbc_copies(pbc_copies, cell, cutoff)
1373 INTEGER, DIMENSION(3), INTENT(INOUT) :: pbc_copies
1374 TYPE(cell_type), INTENT(IN), POINTER :: cell
1375 REAL(kind=dp), INTENT(IN) :: cutoff
1376
1377 REAL(kind=dp) :: proja, projb, projc
1378 REAL(kind=dp), DIMENSION(3) :: axb, axc, bxc
1379
1380 axb(1) = cell%hmat(2, 1)*cell%hmat(3, 2) - cell%hmat(3, 1)*cell%hmat(2, 2)
1381 axb(2) = cell%hmat(3, 1)*cell%hmat(1, 2) - cell%hmat(1, 1)*cell%hmat(3, 2)
1382 axb(3) = cell%hmat(1, 1)*cell%hmat(2, 2) - cell%hmat(2, 1)*cell%hmat(1, 2)
1383 axb(:) = axb(:)/norm2(axb(:))
1384
1385 axc(1) = cell%hmat(2, 1)*cell%hmat(3, 3) - cell%hmat(3, 1)*cell%hmat(2, 3)
1386 axc(2) = cell%hmat(3, 1)*cell%hmat(1, 3) - cell%hmat(1, 1)*cell%hmat(3, 3)
1387 axc(3) = cell%hmat(1, 1)*cell%hmat(2, 3) - cell%hmat(2, 1)*cell%hmat(1, 3)
1388 axc(:) = axc(:)/norm2(axc(:))
1389
1390 bxc(1) = cell%hmat(2, 2)*cell%hmat(3, 3) - cell%hmat(3, 2)*cell%hmat(2, 3)
1391 bxc(2) = cell%hmat(3, 2)*cell%hmat(1, 3) - cell%hmat(1, 2)*cell%hmat(3, 3)
1392 bxc(3) = cell%hmat(1, 2)*cell%hmat(2, 3) - cell%hmat(2, 2)*cell%hmat(1, 3)
1393 bxc(:) = bxc(:)/norm2(bxc(:))
1394
1395 proja = abs(sum(cell%hmat(:, 1)*bxc(:)))*0.5_dp
1396 projb = abs(sum(cell%hmat(:, 2)*axc(:)))*0.5_dp
1397 projc = abs(sum(cell%hmat(:, 3)*axb(:)))*0.5_dp
1398
1399 pbc_copies(:) = 0
1400 DO WHILE ((pbc_copies(1) + 1)*proja <= cutoff)
1401 pbc_copies(1) = pbc_copies(1) + 1
1402 END DO
1403 DO WHILE ((pbc_copies(2) + 1)*projb <= cutoff)
1404 pbc_copies(2) = pbc_copies(2) + 1
1405 END DO
1406 DO WHILE ((pbc_copies(3) + 1)*projc <= cutoff)
1407 pbc_copies(3) = pbc_copies(3) + 1
1408 END DO
1409 ! Apply non periodic setting
1410 pbc_copies(:) = pbc_copies(:)*cell%perd(:)
1411
1412 END SUBROUTINE nnp_compute_pbc_copies
1413
1414END MODULE nnp_acsf
cell_types::pbc
Definition cell_types.F:92
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
nnp_acsf
Functionality for atom centered symmetry functions for neural network potentials.
Definition nnp_acsf.F:14
nnp_acsf::nnp_calc_acsf
subroutine, public nnp_calc_acsf(nnp, i, dsymdxyz, stress)
Calculate atom centered symmetry functions for given atom i.
Definition nnp_acsf.F:59
nnp_acsf::nnp_write_acsf
subroutine, public nnp_write_acsf(nnp, para_env, printtag)
Write symmetry function information.
Definition nnp_acsf.F:1197
nnp_acsf::nnp_init_acsf_groups
subroutine, public nnp_init_acsf_groups(nnp)
Initialize symmetry function groups.
Definition nnp_acsf.F:1055
nnp_acsf::nnp_sort_acsf
subroutine, public nnp_sort_acsf(nnp)
Sort symmetry functions according to different criteria.
Definition nnp_acsf.F:796
nnp_acsf::nnp_sort_ele
subroutine, public nnp_sort_ele(ele, nuc_ele)
Sort element array according to atomic number.
Definition nnp_acsf.F:756
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:13
nnp_environment_types::nnp_cut_tanh
integer, parameter, public nnp_cut_tanh
Definition nnp_environment_types.F:63
nnp_environment_types::nnp_cut_cos
integer, parameter, public nnp_cut_cos
derived data types
Definition nnp_environment_types.F:63
nnp_environment_types::nnp_env_get
subroutine, public nnp_env_get(nnp_env, nnp_forces, subsys, atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, nnp_input, force_env_input, cell, cell_ref, use_ref_cell, nnp_potential_energy, virial)
Returns various attributes of the nnp environment.
Definition nnp_environment_types.F:418
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
nnp_environment_types::nnp_acsf_ang_type
Set of angular symmetry function type.
Definition nnp_environment_types.F:204
nnp_environment_types::nnp_acsf_rad_type
Set of radial symmetry function type.
Definition nnp_environment_types.F:173
nnp_environment_types::nnp_neighbor_type
Contains neighbors list of an atom.
Definition nnp_environment_types.F:230
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:83