(git:fc9bb57)
Loading...
Searching...
No Matches
nnp_environment.F
Go to the documentation of this file.
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 Methods dealing with Neural Network potentials
10!> \author Christoph Schran (christoph.schran@rub.de)
11!> \author Dhruv Sharma (ds2173@cam.ac.uk)
12!> \date 2020-10-10
13! **************************************************************************************************
15
17 USE bibliography, ONLY: behler2007,&
21 cite_reference
22 USE cell_types, ONLY: cell_type
41 USE kinds, ONLY: default_path_length,&
42 dp
51 USE nnp_environment_types, ONLY: &
55 USE nnp_model, ONLY: nnp_write_arc
61#include "./base/base_uses.f90"
62
63 IMPLICIT NONE
64
65 PRIVATE
66
67 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_environment'
69
70 PUBLIC :: nnp_init
71 PUBLIC :: nnp_init_model
72
73CONTAINS
74
75! **************************************************************************************************
76!> \brief Read and initialize all the information for neural network potentials
77!> \param nnp_env ...
78!> \param root_section ...
79!> \param para_env ...
80!> \param force_env_section ...
81!> \param subsys_section ...
82!> \param use_motion_section ...
83!> \date 2020-10-10
84!> \author Christoph Schran (christoph.schran@rub.de)
85! **************************************************************************************************
86 SUBROUTINE nnp_init(nnp_env, root_section, para_env, force_env_section, subsys_section, &
87 use_motion_section)
88 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
89 TYPE(section_vals_type), INTENT(IN), POINTER :: root_section
90 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
91 TYPE(section_vals_type), INTENT(INOUT), POINTER :: force_env_section, subsys_section
92 LOGICAL, INTENT(IN) :: use_motion_section
93
94 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_init'
95
96 INTEGER :: handle
97 LOGICAL :: explicit
98 TYPE(cp_subsys_type), POINTER :: subsys
99 TYPE(section_vals_type), POINTER :: nnp_section
100
101 CALL timeset(routinen, handle)
102 CALL cite_reference(behler2007)
103 CALL cite_reference(behler2011)
104 CALL cite_reference(schran2020a)
105 CALL cite_reference(schran2020b)
106
107 cpassert(ASSOCIATED(nnp_env))
108
109 NULLIFY (nnp_section, subsys)
110
111 IF (.NOT. ASSOCIATED(subsys_section)) THEN
112 subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
113 END IF
114 nnp_section => section_vals_get_subs_vals(force_env_section, "NNP")
115 CALL section_vals_get(nnp_section, explicit=explicit)
116 IF (.NOT. explicit) THEN
117 cpwarn("NNP section not explicitly stated. Using default file names.")
118 END IF
119
120 CALL nnp_env_set(nnp_env=nnp_env, nnp_input=nnp_section, &
121 force_env_input=force_env_section)
122
123 CALL cp_subsys_create(subsys, para_env, root_section, &
124 force_env_section=force_env_section, subsys_section=subsys_section, &
125 use_motion_section=use_motion_section)
126
127 CALL nnp_init_subsys(nnp_env=nnp_env, subsys=subsys, &
128 subsys_section=subsys_section)
129
130 CALL timestop(handle)
131
132 END SUBROUTINE nnp_init
133
134! **************************************************************************************************
135!> \brief Read and initialize all the information for neural network potentials
136!> \param nnp_env ...
137!> \param subsys ...
138!> \param subsys_section ...
139!> \date 2020-10-10
140!> \author Christoph Schran (christoph.schran@rub.de)
141! **************************************************************************************************
142 SUBROUTINE nnp_init_subsys(nnp_env, subsys, subsys_section)
143 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
144 TYPE(cp_subsys_type), INTENT(IN), POINTER :: subsys
145 TYPE(section_vals_type), INTENT(IN), POINTER :: subsys_section
146
147 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_init_subsys'
148
149 INTEGER :: handle, natom
150 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
151 TYPE(cell_type), POINTER :: my_cell, my_cell_ref
152 TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
153 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
154 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
155 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
156
157 CALL timeset(routinen, handle)
158
159 NULLIFY (atomic_kind_set, molecule_kind_set, my_cell, my_cell_ref, &
160 particle_set, molecule_set, local_molecules, local_particles)
161
162 particle_set => subsys%particles%els
163 atomic_kind_set => subsys%atomic_kinds%els
164 molecule_kind_set => subsys%molecule_kinds%els
165 molecule_set => subsys%molecules%els
166 my_cell => subsys%cell
167 my_cell_ref => subsys%cell_ref
168
169 !Print the molecule kind set
170 CALL write_molecule_kind_set(molecule_kind_set, subsys_section)
171
172 !Set up cell
173 CALL nnp_env_set(nnp_env=nnp_env, subsys=subsys, &
174 cell=my_cell, cell_ref=my_cell_ref, &
175 use_ref_cell=subsys%use_ref_cell)
176
177 !Print the atomic coordinates
178 CALL write_fist_particle_coordinates(particle_set, subsys_section)
179 CALL write_particle_distances(particle_set, cell=my_cell, &
180 subsys_section=subsys_section)
181 CALL write_structure_data(particle_set, cell=my_cell, &
182 input_section=subsys_section)
183
184 !Distribute molecules and atoms using the new data structures
185 CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
186 particle_set=particle_set, &
187 local_particles=local_particles, &
188 molecule_kind_set=molecule_kind_set, &
189 molecule_set=molecule_set, &
190 local_molecules=local_molecules, &
191 force_env_section=nnp_env%force_env_input)
192
193 natom = SIZE(particle_set)
194
195 ALLOCATE (nnp_env%nnp_forces(3, natom))
196
197 nnp_env%nnp_forces(:, :) = 0.0_dp
198
199 nnp_env%nnp_potential_energy = 0.0_dp
200
201 ! Set up arrays for calculation:
202 nnp_env%num_atoms = natom
203 ALLOCATE (nnp_env%ele_ind(natom))
204 ALLOCATE (nnp_env%nuc_atoms(natom))
205 ALLOCATE (nnp_env%coord(3, natom))
206 ALLOCATE (nnp_env%atoms(natom))
207 ALLOCATE (nnp_env%sort(natom))
208 ALLOCATE (nnp_env%sort_inv(natom))
209
210 CALL nnp_env_set(nnp_env=nnp_env, &
211 local_molecules=local_molecules, &
212 local_particles=local_particles)
213
214 CALL distribution_1d_release(local_particles)
215 CALL distribution_1d_release(local_molecules)
216
217 CALL nnp_init_model(nnp_env=nnp_env, printtag="NNP")
218
219 CALL timestop(handle)
220
221 END SUBROUTINE nnp_init_subsys
222
223! **************************************************************************************************
224!> \brief Initialize the Neural Network Potential
225!> \param nnp_env ...
226!> \param printtag ...
227!> \date 2020-10-10
228!> \author Christoph Schran (christoph.schran@rub.de)
229! **************************************************************************************************
230 SUBROUTINE nnp_init_model(nnp_env, printtag)
231 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
232 CHARACTER(LEN=*), INTENT(IN) :: printtag
233
234 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_init_model'
235 INTEGER, PARAMETER :: def_str_len = 256
236
237 CHARACTER(len=1), ALLOCATABLE, DIMENSION(:) :: cactfnct
238 CHARACTER(len=2) :: ele
239 CHARACTER(len=def_str_len) :: dummy, line
240 CHARACTER(len=default_path_length) :: file_name
241 INTEGER :: handle, i, io, j, k, nele, nuc_ele, &
242 symfnct_type, unit_nr
243 LOGICAL :: atom_e_found, explicit, first, found
244 REAL(kind=dp) :: energy
245 REAL(kind=dp), DIMENSION(7) :: test_array
246 REAL(kind=dp), DIMENSION(:), POINTER :: work
247 TYPE(cp_logger_type), POINTER :: logger
248 TYPE(cp_parser_type) :: parser
249 TYPE(section_vals_type), POINTER :: bias_section, model_section
250
251 CALL timeset(routinen, handle)
252
253 NULLIFY (logger)
254
255 logger => cp_get_default_logger()
256
257 unit_nr = cp_logger_get_default_unit_nr(logger)
258 IF (unit_nr > 0) THEN
259 WRITE (unit_nr, *) ""
260 WRITE (unit_nr, *) trim(printtag)//"| Neural Network Potential Force Environment"
261 END IF
262
263 model_section => section_vals_get_subs_vals(nnp_env%nnp_input, "MODEL")
264 CALL section_vals_get(model_section, n_repetition=nnp_env%n_committee)
265 ALLOCATE (nnp_env%atomic_energy(nnp_env%num_atoms, nnp_env%n_committee))
266 ALLOCATE (nnp_env%committee_energy(nnp_env%n_committee))
267 ALLOCATE (nnp_env%myforce(3, nnp_env%num_atoms, nnp_env%n_committee))
268 ALLOCATE (nnp_env%committee_forces(3, nnp_env%num_atoms, nnp_env%n_committee))
269 ALLOCATE (nnp_env%committee_stress(3, 3, nnp_env%n_committee))
270
271 CALL section_vals_val_get(nnp_env%nnp_input, "NNP_INPUT_FILE_NAME", c_val=file_name)
272 CALL section_vals_val_get(nnp_env%nnp_input, "RAD_SPLINE_N", i_val=nnp_env%rad_spline_n)
273 IF (nnp_env%rad_spline_n < 2) THEN
274 cpabort("NNP| RAD_SPLINE_N must be >= 2.")
275 END IF
276 CALL section_vals_val_get(nnp_env%nnp_input, "VERLET_SKIN", r_val=nnp_env%verlet_skin)
277 CALL parser_create(parser, file_name, para_env=logger%para_env)
278
279 ! read number of elements and cut_type and check for scale and center
280 nnp_env%scale_acsf = .false.
281 nnp_env%scale_sigma_acsf = .false.
282 ! Defaults for scale min and max:
283 nnp_env%scmin = 0.0_dp
284 nnp_env%scmax = 1.0_dp
285 nnp_env%center_acsf = .false.
286 nnp_env%normnodes = .false.
287 nnp_env%n_hlayer = 0
288
289 IF (unit_nr > 0) THEN
290 WRITE (unit_nr, *) trim(printtag)//"| Reading NNP input from file: ", trim(file_name)
291 END IF
292
293 CALL parser_search_string(parser, "number_of_elements", .true., found, line, &
294 search_from_begin_of_file=.true.)
295 IF (found) THEN
296 READ (line, *) dummy, nnp_env%n_ele
297 ELSE
298 CALL cp_abort(__location__, trim(printtag)// &
299 "| number of elements missing in NNP_INPUT_FILE")
300 END IF
301
302 CALL parser_search_string(parser, "scale_symmetry_functions_sigma", .true., found, &
303 search_from_begin_of_file=.true.)
304 nnp_env%scale_sigma_acsf = found
305
306 CALL parser_search_string(parser, "scale_symmetry_functions", .true., found, &
307 search_from_begin_of_file=.true.)
308 nnp_env%scale_acsf = found
309
310 ! parser_search_string matches substrings, so a search for
311 ! "scale_symmetry_functions" also fires on "..._sigma". Re-search from the
312 ! current parser position to confirm the bare keyword appears on its own line.
313 CALL parser_search_string(parser, "scale_symmetry_functions", .true., found)
314 IF (found .AND. nnp_env%scale_sigma_acsf) THEN
315 cpwarn("Ignoring sigma ACSF scaling; both keywords were set.")
316 nnp_env%scale_sigma_acsf = .false.
317 ELSE IF (.NOT. found .AND. nnp_env%scale_sigma_acsf) THEN
318 nnp_env%scale_acsf = .false.
319 END IF
320
321 CALL parser_search_string(parser, "scale_min_short_atomic", .true., found, line, &
322 search_from_begin_of_file=.true.)
323 IF (found) READ (line, *) dummy, nnp_env%scmin
324
325 CALL parser_search_string(parser, "scale_max_short_atomic", .true., found, line, &
326 search_from_begin_of_file=.true.)
327 IF (found) READ (line, *) dummy, nnp_env%scmax
328
329 CALL parser_search_string(parser, "center_symmetry_functions", .true., found, &
330 search_from_begin_of_file=.true.)
331 nnp_env%center_acsf = found
332 ! n2p2 overwrites sigma scaling, if centering is requested:
333 IF (nnp_env%scale_sigma_acsf .AND. nnp_env%center_acsf) THEN
334 nnp_env%scale_sigma_acsf = .false.
335 END IF
336 ! Print warning if centering and scaling is requested:
337 IF (nnp_env%center_acsf .AND. nnp_env%scale_acsf) THEN
338 IF ((abs(nnp_env%scmin) > epsilon(0.0_dp)*1.0e+4_dp) .OR. (abs(nnp_env%scmax - 1.0_dp) > epsilon(0.0_dp)*1.0e+4_dp)) THEN
339 CALL cp_warn(__location__, &
340 "Centering and scaling of symmetry functions requested while scale_min_short_atomic != 0 and/or "// &
341 "scale_max_short_atomic != 1. Make sure that scaling and centering of symmetry functions in CP2K "// &
342 "is consistent with your training code. "// &
343 "In CP2K: G* = (G - ave(G)) / (max(G) - min(G)) * (Smax - Smin) + Smin")
344 END IF
345 END IF
346
347 CALL parser_search_string(parser, "normalize_nodes", .true., found, &
348 search_from_begin_of_file=.true.)
349 nnp_env%normnodes = found
350
351 CALL parser_search_string(parser, "cutoff_type", .true., found, line, &
352 search_from_begin_of_file=.true.)
353 IF (found) THEN
354 READ (line, *) dummy, nnp_env%cut_type
355 ELSE
356 CALL cp_abort(__location__, trim(printtag)// &
357 "| no cutoff type specified in NNP_INPUT_FILE")
358 END IF
359
360 CALL parser_search_string(parser, "global_hidden_layers_short", .true., found, line, &
361 search_from_begin_of_file=.true.)
362 IF (found) THEN
363 READ (line, *) dummy, nnp_env%n_hlayer
364 ELSE
365 CALL cp_abort(__location__, trim(printtag)// &
366 "| number of hidden layers missing in NNP_INPUT_FILE")
367 END IF
368 nnp_env%n_layer = nnp_env%n_hlayer + 2
369
370 nele = nnp_env%n_ele
371 ALLOCATE (nnp_env%rad(nele))
372 ALLOCATE (nnp_env%ang(nele))
373 ALLOCATE (nnp_env%n_rad(nele))
374 ALLOCATE (nnp_env%n_ang(nele))
375 ALLOCATE (nnp_env%actfnct(nnp_env%n_hlayer + 1))
376 ALLOCATE (cactfnct(nnp_env%n_hlayer + 1))
377 ALLOCATE (nnp_env%ele(nele))
378 ALLOCATE (nnp_env%nuc_ele(nele))
379 ALLOCATE (nnp_env%arc(nele))
380 DO i = 1, nele
381 ALLOCATE (nnp_env%arc(i)%layer(nnp_env%n_layer))
382 ALLOCATE (nnp_env%arc(i)%n_nodes(nnp_env%n_layer))
383 END DO
384 ALLOCATE (nnp_env%n_hnodes(nnp_env%n_hlayer))
385 ALLOCATE (nnp_env%atom_energies(nele))
386 nnp_env%atom_energies = 0.0_dp
387
388 ! read elements, broadcast and sort
389 CALL parser_reset(parser)
390 DO
391 CALL parser_search_string(parser, "elements", .true., found, line)
392 IF (found) THEN
393 READ (line, *) dummy
394 IF (trim(adjustl(dummy)) == "elements") THEN
395 READ (line, *) dummy, nnp_env%ele(:)
396 CALL nnp_sort_ele(nnp_env%ele, nnp_env%nuc_ele)
397 EXIT
398 END IF
399 ELSE
400 CALL cp_abort(__location__, trim(printtag)// &
401 "| elements not specified in NNP_INPUT_FILE")
402 END IF
403 END DO
404
405 CALL parser_search_string(parser, "remove_atom_energies", .true., atom_e_found, &
406 search_from_begin_of_file=.true.)
407
408 IF (atom_e_found) THEN
409 CALL parser_reset(parser)
410 i = 0
411 DO
412 CALL parser_search_string(parser, "atom_energy", .true., found, line)
413 IF (found) THEN
414 READ (line, *) dummy, ele, energy
415 DO j = 1, nele
416 IF (nnp_env%ele(j) == trim(ele)) THEN
417 i = i + 1
418 nnp_env%atom_energies(j) = energy
419 END IF
420 END DO
421 IF (i == nele) EXIT
422 ELSE
423 CALL cp_abort(__location__, trim(printtag)// &
424 "| atom energies are not specified")
425 END IF
426 END DO
427 END IF
428
429 CALL parser_search_string(parser, "global_nodes_short", .true., found, line, &
430 search_from_begin_of_file=.true.)
431 IF (found) THEN
432 READ (line, *) dummy, nnp_env%n_hnodes(:)
433 ELSE
434 CALL cp_abort(__location__, trim(printtag)// &
435 "NNP| global_nodes_short not specified in NNP_INPUT_FILE")
436 END IF
437
438 CALL parser_search_string(parser, "global_activation_short", .true., found, line, &
439 search_from_begin_of_file=.true.)
440 IF (found) THEN
441 READ (line, *) dummy, cactfnct(:)
442 ELSE
443 CALL cp_abort(__location__, trim(printtag)// &
444 "| global_activation_short not specified in NNP_INPUT_FILE")
445 END IF
446
447 DO i = 1, nnp_env%n_hlayer + 1
448 SELECT CASE (cactfnct(i))
449 CASE ("t")
450 nnp_env%actfnct(i) = nnp_actfnct_tanh
451 CASE ("g")
452 nnp_env%actfnct(i) = nnp_actfnct_gaus
453 CASE ("l")
454 nnp_env%actfnct(i) = nnp_actfnct_lin
455 CASE ("c")
456 nnp_env%actfnct(i) = nnp_actfnct_cos
457 CASE ("s")
458 nnp_env%actfnct(i) = nnp_actfnct_sig
459 CASE ("S")
460 nnp_env%actfnct(i) = nnp_actfnct_invsig
461 CASE ("e")
462 nnp_env%actfnct(i) = nnp_actfnct_exp
463 CASE ("p")
464 nnp_env%actfnct(i) = nnp_actfnct_softplus
465 CASE ("h")
466 nnp_env%actfnct(i) = nnp_actfnct_quad
467 CASE DEFAULT
468 CALL cp_abort(__location__, trim(printtag)// &
469 "| Activation function unkown")
470 END SELECT
471 END DO
472
473 ! determine n_rad and n_ang
474 DO i = 1, nele
475 nnp_env%n_rad(i) = 0
476 nnp_env%n_ang(i) = 0
477 END DO
478
479 ! count symfunctions
480 CALL parser_reset(parser)
481 first = .true.
482 DO
483 CALL parser_search_string(parser, "symfunction_short", .true., found, line)
484 IF (found) THEN
485 READ (line, *) dummy, ele, symfnct_type
486 DO i = 1, nele
487 IF (trim(ele) == nnp_env%ele(i)) THEN
488 IF (symfnct_type == 2) THEN
489 nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
490 ELSE IF (symfnct_type == 3) THEN
491 nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
492 ELSE
493 CALL cp_abort(__location__, trim(printtag)// &
494 "| Symmetry function type not supported")
495 END IF
496 END IF
497 END DO
498 first = .false.
499 ELSE
500 IF (first) CALL cp_abort(__location__, trim(printtag)// &
501 "| no symfunction_short specified in NNP_INPUT_FILE")
502 ! no additional symfnct found
503 EXIT
504 END IF
505 END DO
506
507 DO i = 1, nele
508 ALLOCATE (nnp_env%rad(i)%y(nnp_env%n_rad(i)))
509 ALLOCATE (nnp_env%rad(i)%funccut(nnp_env%n_rad(i)))
510 ALLOCATE (nnp_env%rad(i)%eta(nnp_env%n_rad(i)))
511 ALLOCATE (nnp_env%rad(i)%rs(nnp_env%n_rad(i)))
512 ALLOCATE (nnp_env%rad(i)%loc_min(nnp_env%n_rad(i)))
513 ALLOCATE (nnp_env%rad(i)%loc_max(nnp_env%n_rad(i)))
514 ALLOCATE (nnp_env%rad(i)%loc_av(nnp_env%n_rad(i)))
515 ALLOCATE (nnp_env%rad(i)%sigma(nnp_env%n_rad(i)))
516 ALLOCATE (nnp_env%rad(i)%ele(nnp_env%n_rad(i)))
517 ALLOCATE (nnp_env%rad(i)%nuc_ele(nnp_env%n_rad(i)))
518 nnp_env%rad(i)%funccut = 0.0_dp
519 nnp_env%rad(i)%eta = 0.0_dp
520 nnp_env%rad(i)%rs = 0.0_dp
521 nnp_env%rad(i)%ele = 'X'
522 nnp_env%rad(i)%nuc_ele = 0
523
524 ALLOCATE (nnp_env%ang(i)%y(nnp_env%n_ang(i)))
525 ALLOCATE (nnp_env%ang(i)%funccut(nnp_env%n_ang(i)))
526 ALLOCATE (nnp_env%ang(i)%eta(nnp_env%n_ang(i)))
527 ALLOCATE (nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
528 ALLOCATE (nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)))
529 ALLOCATE (nnp_env%ang(i)%lam(nnp_env%n_ang(i)))
530 ALLOCATE (nnp_env%ang(i)%loc_min(nnp_env%n_ang(i)))
531 ALLOCATE (nnp_env%ang(i)%loc_max(nnp_env%n_ang(i)))
532 ALLOCATE (nnp_env%ang(i)%loc_av(nnp_env%n_ang(i)))
533 ALLOCATE (nnp_env%ang(i)%sigma(nnp_env%n_ang(i)))
534 ALLOCATE (nnp_env%ang(i)%ele1(nnp_env%n_ang(i)))
535 ALLOCATE (nnp_env%ang(i)%ele2(nnp_env%n_ang(i)))
536 ALLOCATE (nnp_env%ang(i)%nuc_ele1(nnp_env%n_ang(i)))
537 ALLOCATE (nnp_env%ang(i)%nuc_ele2(nnp_env%n_ang(i)))
538 nnp_env%ang(i)%funccut = 0.0_dp
539 nnp_env%ang(i)%eta = 0.0_dp
540 nnp_env%ang(i)%zeta = 0.0_dp
541 nnp_env%ang(i)%prefzeta = 1.0_dp
542 nnp_env%ang(i)%lam = 0.0_dp
543 nnp_env%ang(i)%ele1 = 'X'
544 nnp_env%ang(i)%ele2 = 'X'
545 nnp_env%ang(i)%nuc_ele1 = 0
546 nnp_env%ang(i)%nuc_ele2 = 0
547
548 ! set number of nodes
549 nnp_env%arc(i)%n_nodes(1) = nnp_env%n_rad(i) + nnp_env%n_ang(i)
550 nnp_env%arc(i)%n_nodes(2:nnp_env%n_layer - 1) = nnp_env%n_hnodes
551 nnp_env%arc(i)%n_nodes(nnp_env%n_layer) = 1
552 DO j = 1, nnp_env%n_layer
553 ALLOCATE (nnp_env%arc(i)%layer(j)%node(nnp_env%arc(i)%n_nodes(j)))
554 ALLOCATE (nnp_env%arc(i)%layer(j)%node_grad(nnp_env%arc(i)%n_nodes(j)))
555 ALLOCATE (nnp_env%arc(i)%layer(j)%tmp_der(nnp_env%arc(i)%n_nodes(1), nnp_env%arc(i)%n_nodes(j)))
556 END DO
557 END DO
558
559 ! read, bcast and sort symfnct parameters
560 DO i = 1, nele
561 nnp_env%n_rad(i) = 0
562 nnp_env%n_ang(i) = 0
563 END DO
564 CALL parser_reset(parser)
565 first = .true.
566 nnp_env%max_cut = 0.0_dp
567 DO
568 CALL parser_search_string(parser, "symfunction_short", .true., found, line)
569 IF (found) THEN
570 READ (line, *) dummy, ele, symfnct_type
571 DO i = 1, nele
572 IF (trim(ele) == nnp_env%ele(i)) THEN
573 IF (symfnct_type == 2) THEN
574 nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
575 READ (line, *) dummy, ele, symfnct_type, &
576 nnp_env%rad(i)%ele(nnp_env%n_rad(i)), &
577 nnp_env%rad(i)%eta(nnp_env%n_rad(i)), &
578 nnp_env%rad(i)%rs(nnp_env%n_rad(i)), &
579 nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
580 IF (nnp_env%max_cut < nnp_env%rad(i)%funccut(nnp_env%n_rad(i))) THEN
581 nnp_env%max_cut = nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
582 END IF
583 ELSE IF (symfnct_type == 3) THEN
584 nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
585 READ (line, *) dummy, ele, symfnct_type, &
586 nnp_env%ang(i)%ele1(nnp_env%n_ang(i)), &
587 nnp_env%ang(i)%ele2(nnp_env%n_ang(i)), &
588 nnp_env%ang(i)%eta(nnp_env%n_ang(i)), &
589 nnp_env%ang(i)%lam(nnp_env%n_ang(i)), &
590 nnp_env%ang(i)%zeta(nnp_env%n_ang(i)), &
591 nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
592 nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)) = &
593 2.0_dp**(1.0_dp - nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
594 IF (nnp_env%max_cut < nnp_env%ang(i)%funccut(nnp_env%n_ang(i))) THEN
595 nnp_env%max_cut = nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
596 END IF
597 ELSE
598 CALL cp_abort(__location__, trim(printtag)// &
599 "| Symmetry function type not supported")
600 END IF
601 END IF
602 END DO
603 first = .false.
604 ELSE
605 IF (first) CALL cp_abort(__location__, trim(printtag)// &
606 "| no symfunction_short specified in NNP_INPUT_FILE")
607 ! no additional symfnct found
608 EXIT
609 END IF
610 END DO
611
612 DO i = 1, nele
613 DO j = 1, nnp_env%n_rad(i)
614 CALL get_ptable_info(nnp_env%rad(i)%ele(j), number=nnp_env%rad(i)%nuc_ele(j))
615 END DO
616 DO j = 1, nnp_env%n_ang(i)
617 CALL get_ptable_info(nnp_env%ang(i)%ele1(j), number=nnp_env%ang(i)%nuc_ele1(j))
618 CALL get_ptable_info(nnp_env%ang(i)%ele2(j), number=nnp_env%ang(i)%nuc_ele2(j))
619 ! sort ele1 and ele2
620 IF (nnp_env%ang(i)%nuc_ele1(j) > nnp_env%ang(i)%nuc_ele2(j)) THEN
621 ele = nnp_env%ang(i)%ele1(j)
622 nnp_env%ang(i)%ele1(j) = nnp_env%ang(i)%ele2(j)
623 nnp_env%ang(i)%ele2(j) = ele
624 nuc_ele = nnp_env%ang(i)%nuc_ele1(j)
625 nnp_env%ang(i)%nuc_ele1(j) = nnp_env%ang(i)%nuc_ele2(j)
626 nnp_env%ang(i)%nuc_ele2(j) = nuc_ele
627 END IF
628 END DO
629 END DO
630 ! Done with input.nn file
631 CALL parser_release(parser)
632
633 ! sort symmetry functions and output information
634 CALL nnp_sort_acsf(nnp_env)
635 CALL nnp_write_acsf(nnp_env, logger%para_env, printtag)
636 CALL nnp_write_arc(nnp_env, logger%para_env, printtag)
637
638 ! read scaling information from file
639 IF (nnp_env%scale_acsf .OR. nnp_env%center_acsf .OR. nnp_env%scale_sigma_acsf) THEN
640 IF (unit_nr > 0) THEN
641 WRITE (unit_nr, *) trim(printtag)//"| Reading scaling information from file: ", trim(file_name)
642 END IF
643 CALL section_vals_val_get(nnp_env%nnp_input, "SCALE_FILE_NAME", &
644 c_val=file_name)
645 CALL parser_create(parser, file_name, para_env=logger%para_env)
646
647 ! Get number of elements in scaling file
648 CALL parser_read_line(parser, 1)
649 k = 0
650 DO WHILE (k < 7)
651 READ (parser%input_line, *, iostat=io) test_array(1:k)
652 IF (io == -1) EXIT
653 k = k + 1
654 END DO
655 k = k - 1
656
657 IF (k == 5 .AND. nnp_env%scale_sigma_acsf) THEN
658 cpabort("Sigma scaling requested, but scaling.data does not contain sigma.")
659 END IF
660
661 CALL parser_reset(parser)
662 DO i = 1, nnp_env%n_ele
663 DO j = 1, nnp_env%n_rad(i)
664 CALL parser_read_line(parser, 1)
665 IF (nnp_env%scale_sigma_acsf) THEN
666 READ (parser%input_line, *) dummy, dummy, &
667 nnp_env%rad(i)%loc_min(j), &
668 nnp_env%rad(i)%loc_max(j), &
669 nnp_env%rad(i)%loc_av(j), &
670 nnp_env%rad(i)%sigma(j)
671 ELSE
672 READ (parser%input_line, *) dummy, dummy, &
673 nnp_env%rad(i)%loc_min(j), &
674 nnp_env%rad(i)%loc_max(j), &
675 nnp_env%rad(i)%loc_av(j)
676 END IF
677 END DO
678 DO j = 1, nnp_env%n_ang(i)
679 CALL parser_read_line(parser, 1)
680 IF (nnp_env%scale_sigma_acsf) THEN
681 READ (parser%input_line, *) dummy, dummy, &
682 nnp_env%ang(i)%loc_min(j), &
683 nnp_env%ang(i)%loc_max(j), &
684 nnp_env%ang(i)%loc_av(j), &
685 nnp_env%ang(i)%sigma(j)
686 ELSE
687 READ (parser%input_line, *) dummy, dummy, &
688 nnp_env%ang(i)%loc_min(j), &
689 nnp_env%ang(i)%loc_max(j), &
690 nnp_env%ang(i)%loc_av(j)
691 END IF
692 END DO
693 END DO
694 CALL parser_release(parser)
695
696 ! Reject degenerate scaling.data at load time. The (loc_max - loc_min)
697 ! and sigma denominators are used unguarded inside nnp_scale_acsf;
698 ! a zero range or zero sigma there produces silent NaN forces.
699 IF (nnp_env%scale_acsf) THEN
700 DO i = 1, nnp_env%n_ele
701 DO j = 1, nnp_env%n_rad(i)
702 IF (nnp_env%rad(i)%loc_max(j) <= nnp_env%rad(i)%loc_min(j)) THEN
703 WRITE (line, '(A,I0,A,I0,A,2(1X,ES13.6))') &
704 "scaling.data: radial sf range non-positive for element ", i, &
705 " sf ", j, " (loc_min, loc_max) =", &
706 nnp_env%rad(i)%loc_min(j), nnp_env%rad(i)%loc_max(j)
707 cpabort(trim(line))
708 END IF
709 END DO
710 DO j = 1, nnp_env%n_ang(i)
711 IF (nnp_env%ang(i)%loc_max(j) <= nnp_env%ang(i)%loc_min(j)) THEN
712 WRITE (line, '(A,I0,A,I0,A,2(1X,ES13.6))') &
713 "scaling.data: angular sf range non-positive for element ", i, &
714 " sf ", j, " (loc_min, loc_max) =", &
715 nnp_env%ang(i)%loc_min(j), nnp_env%ang(i)%loc_max(j)
716 cpabort(trim(line))
717 END IF
718 END DO
719 END DO
720 END IF
721 IF (nnp_env%scale_sigma_acsf) THEN
722 DO i = 1, nnp_env%n_ele
723 DO j = 1, nnp_env%n_rad(i)
724 IF (nnp_env%rad(i)%sigma(j) <= 0.0_dp) THEN
725 WRITE (line, '(A,I0,A,I0,A,1X,ES13.6)') &
726 "scaling.data: radial sf sigma non-positive for element ", i, &
727 " sf ", j, " sigma =", nnp_env%rad(i)%sigma(j)
728 cpabort(trim(line))
729 END IF
730 END DO
731 DO j = 1, nnp_env%n_ang(i)
732 IF (nnp_env%ang(i)%sigma(j) <= 0.0_dp) THEN
733 WRITE (line, '(A,I0,A,I0,A,1X,ES13.6)') &
734 "scaling.data: angular sf sigma non-positive for element ", i, &
735 " sf ", j, " sigma =", nnp_env%ang(i)%sigma(j)
736 cpabort(trim(line))
737 END IF
738 END DO
739 END DO
740 END IF
741 END IF
742
743 CALL nnp_init_acsf_groups(nnp_env)
744
745 ! Allocate per-layer weight tables, then load each committee member
746 ! from disk. Allocation lives at the call site because n_committee /
747 ! n_layer / n_nodes are owned here; the loader only fills them.
748 DO i = 1, nnp_env%n_ele
749 DO j = 2, nnp_env%n_layer
750 ALLOCATE (nnp_env%arc(i)%layer(j)%weights(nnp_env%arc(i)%n_nodes(j - 1), &
751 nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
752 ALLOCATE (nnp_env%arc(i)%layer(j)%bweights(nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
753 END DO
754 END DO
755 CALL nnp_read_committee_weights(nnp_env, model_section, printtag)
756
757 nnp_env%expol = 0
758
759 ! Bias the standard deviation of committee disagreement
760 NULLIFY (bias_section)
761 explicit = .false.
762 !HELIUM NNP does not currently define a bias term
763 bias_section => section_vals_get_subs_vals(nnp_env%nnp_input, "BIAS", can_return_null=.true.)
764 IF (ASSOCIATED(bias_section)) CALL section_vals_get(bias_section, explicit=explicit)
765 nnp_env%bias = .false.
766 IF (explicit) THEN
767 IF (nnp_env%n_committee > 1) THEN
768 IF (unit_nr > 0) THEN
769 WRITE (unit_nr, *) "NNP| Biasing of committee disagreement enabled"
770 END IF
771 nnp_env%bias = .true.
772 ALLOCATE (nnp_env%bias_forces(3, nnp_env%num_atoms))
773 ALLOCATE (nnp_env%bias_e_avrg(nnp_env%n_committee))
774 CALL section_vals_val_get(bias_section, "SIGMA_0", r_val=nnp_env%bias_sigma0)
775 CALL section_vals_val_get(bias_section, "K_B", r_val=nnp_env%bias_kb)
776 nnp_env%bias_e_avrg(:) = 0.0_dp
777 CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", explicit=explicit)
778 nnp_env%bias_align = explicit
779 IF (explicit) THEN
780 NULLIFY (work)
781 CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", r_vals=work)
782 IF (SIZE(work) /= nnp_env%n_committee) THEN
783 cpabort("ALIGN_NNP_ENERGIES size mismatch wrt committee size.")
784 END IF
785 nnp_env%bias_e_avrg(:) = work
786 IF (unit_nr > 0) THEN
787 WRITE (unit_nr, *) trim(printtag)//"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
788 END IF
789 END IF
790 ELSE
791 cpwarn("NNP committee size is 1, BIAS section is ignored.")
792 END IF
793 END IF
794
795 IF (unit_nr > 0) THEN
796 WRITE (unit_nr, *) trim(printtag)//"| NNP force environment initialized"
797 END IF
798
799 CALL timestop(handle)
800
801 END SUBROUTINE nnp_init_model
802
803! **************************************************************************************************
804!> \brief Read committee weights from disk into pre-allocated arc layer tables.
805!>
806!> One file per (committee member, element), named <base>.<nuc>.data. The file
807!> is a flat list of doubles; the order on disk is layer-major over (i->j weights
808!> then biases) so a single counter walks every value into its slot.
809!>
810!> \param nnp_env NNP environment with arc()%layer()%weights/bweights already allocated.
811!> \param model_section &MODEL section block; supplies the per-member WEIGHTS base path.
812!> \param printtag Log-line prefix.
813! **************************************************************************************************
814 SUBROUTINE nnp_read_committee_weights(nnp_env, model_section, printtag)
815 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
816 TYPE(section_vals_type), POINTER :: model_section
817 CHARACTER(LEN=*), INTENT(IN) :: printtag
818
819 CHARACTER(len=default_path_length) :: base_name, file_name
820 INTEGER :: i, i_com, iweight, j, k, l, n_weight, &
821 unit_nr
822 LOGICAL :: at_end
823 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights
824 TYPE(cp_logger_type), POINTER :: logger
825 TYPE(cp_parser_type) :: parser
826
827 NULLIFY (logger)
828 logger => cp_get_default_logger()
829 unit_nr = cp_logger_get_default_unit_nr(logger)
830
831 DO i_com = 1, nnp_env%n_committee
832 CALL section_vals_val_get(model_section, "WEIGHTS", c_val=base_name, i_rep_section=i_com)
833 IF (unit_nr > 0) THEN
834 WRITE (unit_nr, *) trim(printtag)//"| Initializing weights for model: ", i_com
835 END IF
836 DO i = 1, nnp_env%n_ele
837 WRITE (file_name, '(A,I0.3,A)') trim(base_name)//".", nnp_env%nuc_ele(i), ".data"
838 IF (unit_nr > 0) THEN
839 WRITE (unit_nr, *) trim(printtag)//"| Reading weights from file: ", trim(file_name)
840 END IF
841 CALL parser_create(parser, file_name, para_env=logger%para_env)
842 n_weight = 0
843 DO WHILE (.true.)
844 CALL parser_read_line(parser, 1, at_end)
845 IF (at_end) EXIT
846 n_weight = n_weight + 1
847 END DO
848
849 ALLOCATE (weights(n_weight))
850
851 CALL parser_reset(parser)
852 DO j = 1, n_weight
853 CALL parser_read_line(parser, 1)
854 READ (parser%input_line, *) weights(j)
855 END DO
856 CALL parser_release(parser)
857
858 ! sort weights into corresponding arrays
859 iweight = 0
860 DO j = 2, nnp_env%n_layer
861 DO k = 1, nnp_env%arc(i)%n_nodes(j - 1)
862 DO l = 1, nnp_env%arc(i)%n_nodes(j)
863 iweight = iweight + 1
864 nnp_env%arc(i)%layer(j)%weights(k, l, i_com) = weights(iweight)
865 END DO
866 END DO
867
868 DO k = 1, nnp_env%arc(i)%n_nodes(j)
869 iweight = iweight + 1
870 nnp_env%arc(i)%layer(j)%bweights(k, i_com) = weights(iweight)
871 END DO
872 END DO
873
874 DEALLOCATE (weights)
875 END DO
876 END DO
877
878 END SUBROUTINE nnp_read_committee_weights
879
880END MODULE nnp_environment
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public schran2020b
integer, save, public schran2020a
integer, save, public behler2011
integer, save, public behler2007
Handles all functions related to the CELL.
Definition cell_types.F:15
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_read_line(parser, nline, at_end)
Read the next line from a logical unit "unit" (I/O node only). Skip (nline-1) lines and skip also all...
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_reset(parser)
Resets the parser: rewinding the unit and re-initializing all parser structures.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
Initialize a small environment for a particular calculation.
subroutine, public cp_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, use_motion_section, qmmm, qmmm_env, exclusions, elkind)
Creates allocates and fills subsys from given input.
types that represent a subsys, i.e. a part of the system
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public distribution_1d_release(distribution_1d)
releases the given distribution_1d
Distribution methods for atoms, particles, or molecules.
subroutine, public distribute_molecules_1d(atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, force_env_section, prev_molecule_kind_set, prev_local_molecules)
Distribute molecules and particles.
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_path_length
Definition kinds.F:58
Interface to the message passing library MPI.
Define the molecule kind structure types and the corresponding functionality.
subroutine, public write_molecule_kind_set(molecule_kind_set, subsys_section)
Write a moleculeatomic kind set data set to the output unit.
Define the data structure for the molecule information.
Functionality for atom centered symmetry functions for neural network potentials.
Definition nnp_acsf.F:16
subroutine, public nnp_write_acsf(nnp, para_env, printtag)
Print a summary of the active symmetry-function set on the source rank. Emits one line per element li...
Definition nnp_acsf.F:1782
subroutine, public nnp_init_acsf_groups(nnp)
Pack symmetry functions into groups that share input parameters. Builds nnprad(i)symfgrp(:) and nnpan...
Definition nnp_acsf.F:1611
subroutine, public nnp_sort_acsf(nnp)
Sort radial and angular symmetry functions in canonical order. Radial SFs are sorted by eta (ascendin...
Definition nnp_acsf.F:1379
subroutine, public nnp_sort_ele(ele, nuc_ele)
Sort an (ele, nuc_ele) pair of arrays in ascending order of atomic number. Used to canonicalise eleme...
Definition nnp_acsf.F:1340
Data types for neural network potentials.
integer, parameter, public nnp_actfnct_lin
subroutine, public nnp_env_set(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)
Sets various attributes of the nnp environment.
integer, parameter, public nnp_actfnct_cos
integer, parameter, public nnp_actfnct_invsig
integer, parameter, public nnp_actfnct_sig
integer, parameter, public nnp_actfnct_exp
integer, parameter, public nnp_actfnct_softplus
integer, parameter, public nnp_actfnct_quad
integer, parameter, public nnp_actfnct_gaus
integer, parameter, public nnp_actfnct_tanh
Methods dealing with Neural Network potentials.
subroutine, public nnp_init_model(nnp_env, printtag)
Initialize the Neural Network Potential.
subroutine, public nnp_init(nnp_env, root_section, para_env, force_env_section, subsys_section, use_motion_section)
Read and initialize all the information for neural network potentials.
Methods dealing with core routines for artificial neural networks.
Definition nnp_model.F:13
subroutine, public nnp_write_arc(nnp, para_env, printtag)
Write neural network architecture information.
Definition nnp_model.F:47
Define methods related to particle_type.
subroutine, public write_fist_particle_coordinates(particle_set, subsys_section, charges)
Write the atomic coordinates to the output unit.
subroutine, public write_structure_data(particle_set, cell, input_section)
Write structure data requested by a separate structure data input section to the output unit....
subroutine, public write_particle_distances(particle_set, cell, subsys_section)
Write the matrix of the particle distances to the output unit.
Define the data structure for the particle information.
Periodic Table related data definitions.
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.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
structure to store local (to a processor) ordered lists of integers.
stores all the informations relevant to an mpi environment
Main data type collecting all relevant data for neural network potentials.