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 Methods dealing with Neural Network potentials
10 : !> \author Christoph Schran (christoph.schran@rub.de)
11 : !> \date 2020-10-10
12 : ! **************************************************************************************************
13 : MODULE nnp_environment
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE bibliography, ONLY: Behler2007,&
17 : Behler2011,&
18 : Schran2020a,&
19 : Schran2020b,&
20 : cite_reference
21 : USE cell_methods, ONLY: read_cell,&
22 : write_cell
23 : USE cell_types, ONLY: cell_release,&
24 : cell_type,&
25 : get_cell
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_get_default_unit_nr,&
28 : cp_logger_type
29 : USE cp_parser_methods, ONLY: parser_read_line,&
30 : parser_search_string
31 : USE cp_parser_types, ONLY: cp_parser_type,&
32 : parser_create,&
33 : parser_release,&
34 : parser_reset
35 : USE cp_subsys_methods, ONLY: cp_subsys_create
36 : USE cp_subsys_types, ONLY: cp_subsys_set,&
37 : cp_subsys_type
38 : USE distribution_1d_types, ONLY: distribution_1d_release,&
39 : distribution_1d_type
40 : USE distribution_methods, ONLY: distribute_molecules_1d
41 : USE input_section_types, ONLY: section_vals_get,&
42 : section_vals_get_subs_vals,&
43 : section_vals_type,&
44 : section_vals_val_get
45 : USE kinds, ONLY: default_path_length,&
46 : dp
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE molecule_kind_types, ONLY: molecule_kind_type,&
49 : write_molecule_kind_set
50 : USE molecule_types, ONLY: molecule_type
51 : USE nnp_acsf, ONLY: nnp_init_acsf_groups,&
52 : nnp_sort_acsf,&
53 : nnp_sort_ele,&
54 : nnp_write_acsf
55 : USE nnp_environment_types, ONLY: &
56 : nnp_actfnct_cos, nnp_actfnct_exp, nnp_actfnct_gaus, nnp_actfnct_invsig, nnp_actfnct_lin, &
57 : nnp_actfnct_quad, nnp_actfnct_sig, nnp_actfnct_softplus, nnp_actfnct_tanh, nnp_env_set, &
58 : nnp_type
59 : USE nnp_model, ONLY: nnp_write_arc
60 : USE particle_methods, ONLY: write_fist_particle_coordinates,&
61 : write_particle_distances,&
62 : write_structure_data
63 : USE particle_types, ONLY: particle_type
64 : USE periodic_table, ONLY: get_ptable_info
65 : #include "./base/base_uses.f90"
66 :
67 : IMPLICIT NONE
68 :
69 : PRIVATE
70 :
71 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
72 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_environment'
73 :
74 : PUBLIC :: nnp_init
75 : PUBLIC :: nnp_init_model
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief Read and initialize all the information for neural network potentials
81 : !> \param nnp_env ...
82 : !> \param root_section ...
83 : !> \param para_env ...
84 : !> \param force_env_section ...
85 : !> \param subsys_section ...
86 : !> \param use_motion_section ...
87 : !> \date 2020-10-10
88 : !> \author Christoph Schran (christoph.schran@rub.de)
89 : ! **************************************************************************************************
90 28 : SUBROUTINE nnp_init(nnp_env, root_section, para_env, force_env_section, subsys_section, &
91 : use_motion_section)
92 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
93 : TYPE(section_vals_type), INTENT(IN), POINTER :: root_section
94 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
95 : TYPE(section_vals_type), INTENT(INOUT), POINTER :: force_env_section, subsys_section
96 : LOGICAL, INTENT(IN) :: use_motion_section
97 :
98 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_init'
99 :
100 : INTEGER :: handle
101 : LOGICAL :: explicit, use_ref_cell
102 : REAL(KIND=dp), DIMENSION(3) :: abc
103 : TYPE(cell_type), POINTER :: cell, cell_ref
104 : TYPE(cp_subsys_type), POINTER :: subsys
105 : TYPE(section_vals_type), POINTER :: cell_section, nnp_section
106 :
107 14 : CALL timeset(routineN, handle)
108 14 : CALL cite_reference(Behler2007)
109 14 : CALL cite_reference(Behler2011)
110 14 : CALL cite_reference(Schran2020a)
111 14 : CALL cite_reference(Schran2020b)
112 :
113 14 : CPASSERT(ASSOCIATED(nnp_env))
114 :
115 14 : NULLIFY (cell_section, nnp_section, cell, cell_ref, subsys)
116 :
117 14 : IF (.NOT. ASSOCIATED(subsys_section)) THEN
118 0 : subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
119 : END IF
120 14 : cell_section => section_vals_get_subs_vals(subsys_section, "CELL")
121 14 : nnp_section => section_vals_get_subs_vals(force_env_section, "NNP")
122 14 : CALL section_vals_get(nnp_section, explicit=explicit)
123 14 : IF (.NOT. explicit) THEN
124 0 : CPWARN("NNP section not explicitly stated. Using default file names.")
125 : END IF
126 :
127 : CALL nnp_env_set(nnp_env=nnp_env, nnp_input=nnp_section, &
128 14 : force_env_input=force_env_section)
129 :
130 : CALL read_cell(cell=cell, cell_ref=cell_ref, use_ref_cell=use_ref_cell, cell_section=cell_section, &
131 14 : para_env=para_env)
132 14 : CALL get_cell(cell=cell, abc=abc)
133 14 : CALL write_cell(cell=cell, subsys_section=subsys_section)
134 :
135 : CALL cp_subsys_create(subsys, para_env, root_section, &
136 : force_env_section=force_env_section, subsys_section=subsys_section, &
137 14 : use_motion_section=use_motion_section)
138 :
139 : CALL nnp_init_subsys(nnp_env=nnp_env, subsys=subsys, cell=cell, &
140 : cell_ref=cell_ref, use_ref_cell=use_ref_cell, &
141 14 : subsys_section=subsys_section)
142 :
143 14 : CALL cell_release(cell)
144 14 : CALL cell_release(cell_ref)
145 :
146 14 : CALL timestop(handle)
147 :
148 14 : END SUBROUTINE nnp_init
149 :
150 : ! **************************************************************************************************
151 : !> \brief Read and initialize all the information for neural network potentials
152 : !> \param nnp_env ...
153 : !> \param subsys ...
154 : !> \param cell ...
155 : !> \param cell_ref ...
156 : !> \param use_ref_cell ...
157 : !> \param subsys_section ...
158 : !> \date 2020-10-10
159 : !> \author Christoph Schran (christoph.schran@rub.de)
160 : ! **************************************************************************************************
161 14 : SUBROUTINE nnp_init_subsys(nnp_env, subsys, cell, cell_ref, use_ref_cell, subsys_section)
162 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
163 : TYPE(cp_subsys_type), INTENT(IN), POINTER :: subsys
164 : TYPE(cell_type), INTENT(INOUT), POINTER :: cell, cell_ref
165 : LOGICAL, INTENT(IN) :: use_ref_cell
166 : TYPE(section_vals_type), INTENT(IN), POINTER :: subsys_section
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_init_subsys'
169 :
170 : INTEGER :: handle, natom
171 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
172 : TYPE(distribution_1d_type), POINTER :: local_molecules, local_particles
173 14 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
174 14 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
175 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
176 :
177 14 : CALL timeset(routineN, handle)
178 :
179 : NULLIFY (atomic_kind_set, molecule_kind_set, particle_set, molecule_set, &
180 14 : local_molecules, local_particles)
181 :
182 14 : particle_set => subsys%particles%els
183 14 : atomic_kind_set => subsys%atomic_kinds%els
184 14 : molecule_kind_set => subsys%molecule_kinds%els
185 14 : molecule_set => subsys%molecules%els
186 :
187 : !Print the molecule kind set
188 14 : CALL write_molecule_kind_set(molecule_kind_set, subsys_section)
189 :
190 : !Print the atomic coordinates
191 14 : CALL write_fist_particle_coordinates(particle_set, subsys_section)
192 : CALL write_particle_distances(particle_set, cell=cell, &
193 14 : subsys_section=subsys_section)
194 : CALL write_structure_data(particle_set, cell=cell, &
195 14 : input_section=subsys_section)
196 :
197 : !Distribute molecules and atoms using the new data structures
198 : CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
199 : particle_set=particle_set, &
200 : local_particles=local_particles, &
201 : molecule_kind_set=molecule_kind_set, &
202 : molecule_set=molecule_set, &
203 : local_molecules=local_molecules, &
204 14 : force_env_section=nnp_env%force_env_input)
205 :
206 14 : natom = SIZE(particle_set)
207 :
208 42 : ALLOCATE (nnp_env%nnp_forces(3, natom))
209 :
210 9254 : nnp_env%nnp_forces(:, :) = 0.0_dp
211 :
212 14 : nnp_env%nnp_potential_energy = 0.0_dp
213 :
214 : ! Set up arrays for calculation:
215 14 : nnp_env%num_atoms = natom
216 42 : ALLOCATE (nnp_env%ele_ind(natom))
217 42 : ALLOCATE (nnp_env%nuc_atoms(natom))
218 42 : ALLOCATE (nnp_env%coord(3, natom))
219 42 : ALLOCATE (nnp_env%atoms(natom))
220 42 : ALLOCATE (nnp_env%sort(natom))
221 42 : ALLOCATE (nnp_env%sort_inv(natom))
222 :
223 14 : CALL cp_subsys_set(subsys, cell=cell)
224 :
225 : CALL nnp_env_set(nnp_env=nnp_env, subsys=subsys, &
226 : cell_ref=cell_ref, use_ref_cell=use_ref_cell, &
227 : local_molecules=local_molecules, &
228 14 : local_particles=local_particles)
229 :
230 14 : CALL distribution_1d_release(local_particles)
231 14 : CALL distribution_1d_release(local_molecules)
232 :
233 14 : CALL nnp_init_model(nnp_env=nnp_env, printtag="NNP")
234 :
235 14 : CALL timestop(handle)
236 :
237 14 : END SUBROUTINE nnp_init_subsys
238 :
239 : ! **************************************************************************************************
240 : !> \brief Initialize the Neural Network Potential
241 : !> \param nnp_env ...
242 : !> \param printtag ...
243 : !> \date 2020-10-10
244 : !> \author Christoph Schran (christoph.schran@rub.de)
245 : ! **************************************************************************************************
246 15 : SUBROUTINE nnp_init_model(nnp_env, printtag)
247 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp_env
248 : CHARACTER(LEN=*), INTENT(IN) :: printtag
249 :
250 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_init_model'
251 : INTEGER, PARAMETER :: def_str_len = 256, &
252 : default_path_length = 256
253 :
254 15 : CHARACTER(len=1), ALLOCATABLE, DIMENSION(:) :: cactfnct
255 : CHARACTER(len=2) :: ele
256 : CHARACTER(len=def_str_len) :: dummy, line
257 : CHARACTER(len=default_path_length) :: base_name, file_name
258 : INTEGER :: handle, i, i_com, io, iweight, j, k, l, &
259 : n_weight, nele, nuc_ele, symfnct_type, &
260 : unit_nr
261 : LOGICAL :: at_end, atom_e_found, explicit, first, &
262 : found
263 : REAL(KIND=dp) :: energy
264 15 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weights
265 : REAL(KIND=dp), DIMENSION(7) :: test_array
266 15 : REAL(KIND=dp), DIMENSION(:), POINTER :: work
267 : TYPE(cp_logger_type), POINTER :: logger
268 : TYPE(cp_parser_type) :: parser
269 : TYPE(section_vals_type), POINTER :: bias_section, model_section
270 :
271 15 : CALL timeset(routineN, handle)
272 :
273 15 : NULLIFY (logger)
274 :
275 15 : logger => cp_get_default_logger()
276 :
277 15 : IF (logger%para_env%is_source()) THEN
278 8 : unit_nr = cp_logger_get_default_unit_nr(logger)
279 8 : WRITE (unit_nr, *) ""
280 8 : WRITE (unit_nr, *) TRIM(printtag)//"| Neural Network Potential Force Environment"
281 : END IF
282 :
283 15 : model_section => section_vals_get_subs_vals(nnp_env%nnp_input, "MODEL")
284 15 : CALL section_vals_get(model_section, n_repetition=nnp_env%n_committee)
285 60 : ALLOCATE (nnp_env%atomic_energy(nnp_env%num_atoms, nnp_env%n_committee))
286 45 : ALLOCATE (nnp_env%committee_energy(nnp_env%n_committee))
287 60 : ALLOCATE (nnp_env%myforce(3, nnp_env%num_atoms, nnp_env%n_committee))
288 60 : ALLOCATE (nnp_env%committee_forces(3, nnp_env%num_atoms, nnp_env%n_committee))
289 45 : ALLOCATE (nnp_env%committee_stress(3, 3, nnp_env%n_committee))
290 :
291 15 : CALL section_vals_val_get(nnp_env%nnp_input, "NNP_INPUT_FILE_NAME", c_val=file_name)
292 15 : CALL parser_create(parser, file_name, para_env=logger%para_env)
293 :
294 : ! read number of elements and cut_type and check for scale and center
295 15 : nnp_env%scale_acsf = .FALSE.
296 15 : nnp_env%scale_sigma_acsf = .FALSE.
297 : ! Defaults for scale min and max:
298 15 : nnp_env%scmin = 0.0_dp
299 15 : nnp_env%scmax = 1.0_dp
300 15 : nnp_env%center_acsf = .FALSE.
301 15 : nnp_env%normnodes = .FALSE.
302 15 : nnp_env%n_hlayer = 0
303 :
304 15 : IF (logger%para_env%is_source()) THEN
305 8 : unit_nr = cp_logger_get_default_unit_nr(logger)
306 8 : WRITE (unit_nr, *) TRIM(printtag)//"| Reading NNP input from file: ", TRIM(file_name)
307 : END IF
308 :
309 : CALL parser_search_string(parser, "number_of_elements", .TRUE., found, line, &
310 15 : search_from_begin_of_file=.TRUE.)
311 15 : IF (found) THEN
312 15 : READ (line, *) dummy, nnp_env%n_ele
313 : ELSE
314 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
315 0 : "| number of elements missing in NNP_INPUT_FILE")
316 : END IF
317 :
318 : CALL parser_search_string(parser, "scale_symmetry_functions_sigma", .TRUE., found, &
319 15 : search_from_begin_of_file=.TRUE.)
320 15 : nnp_env%scale_sigma_acsf = found
321 :
322 : CALL parser_search_string(parser, "scale_symmetry_functions", .TRUE., found, &
323 15 : search_from_begin_of_file=.TRUE.)
324 15 : nnp_env%scale_acsf = found
325 :
326 : ! Test if there are two keywords of this:
327 15 : CALL parser_search_string(parser, "scale_symmetry_functions", .TRUE., found)
328 15 : IF (found .AND. nnp_env%scale_sigma_acsf) THEN
329 0 : CPWARN('Two scaling keywords in the input, we will ignore sigma scaling in this case')
330 0 : nnp_env%scale_sigma_acsf = .FALSE.
331 15 : ELSE IF (.NOT. found .AND. nnp_env%scale_sigma_acsf) THEN
332 0 : nnp_env%scale_acsf = .FALSE.
333 : END IF
334 :
335 : CALL parser_search_string(parser, "scale_min_short_atomic", .TRUE., found, line, &
336 15 : search_from_begin_of_file=.TRUE.)
337 15 : IF (found) READ (line, *) dummy, nnp_env%scmin
338 :
339 : CALL parser_search_string(parser, "scale_max_short_atomic", .TRUE., found, line, &
340 15 : search_from_begin_of_file=.TRUE.)
341 15 : IF (found) READ (line, *) dummy, nnp_env%scmax
342 :
343 : CALL parser_search_string(parser, "center_symmetry_functions", .TRUE., found, &
344 15 : search_from_begin_of_file=.TRUE.)
345 15 : nnp_env%center_acsf = found
346 : ! n2p2 overwrites sigma scaling, if centering is requested:
347 15 : IF (nnp_env%scale_sigma_acsf .AND. nnp_env%center_acsf) THEN
348 0 : nnp_env%scale_sigma_acsf = .FALSE.
349 : END IF
350 : ! Print warning if centering and scaling is requested:
351 15 : IF (nnp_env%center_acsf .AND. nnp_env%scale_acsf) THEN
352 15 : 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
353 : CALL cp_warn(__LOCATION__, &
354 : "Centering and scaling of symmetry functions requested while scale_min_short_atomic != 0 and/or "// &
355 : "scale_max_short_atomic != 1. Make sure that scaling and centering of symmetry functions in CP2K "// &
356 : "is consistent with your training code. "// &
357 0 : "In CP2K: G* = (G - ave(G)) / (max(G) - min(G)) * (Smax - Smin) + Smin")
358 : END IF
359 : END IF
360 :
361 : CALL parser_search_string(parser, "normalize_nodes", .TRUE., found, &
362 15 : search_from_begin_of_file=.TRUE.)
363 15 : nnp_env%normnodes = found
364 :
365 : CALL parser_search_string(parser, "cutoff_type", .TRUE., found, line, &
366 15 : search_from_begin_of_file=.TRUE.)
367 15 : IF (found) THEN
368 15 : READ (line, *) dummy, nnp_env%cut_type
369 : ELSE
370 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
371 0 : "| no cutoff type specified in NNP_INPUT_FILE")
372 : END IF
373 :
374 : CALL parser_search_string(parser, "global_hidden_layers_short", .TRUE., found, line, &
375 15 : search_from_begin_of_file=.TRUE.)
376 15 : IF (found) THEN
377 15 : READ (line, *) dummy, nnp_env%n_hlayer
378 : ELSE
379 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
380 0 : "| number of hidden layers missing in NNP_INPUT_FILE")
381 : END IF
382 15 : nnp_env%n_layer = nnp_env%n_hlayer + 2
383 :
384 15 : nele = nnp_env%n_ele
385 76 : ALLOCATE (nnp_env%rad(nele))
386 76 : ALLOCATE (nnp_env%ang(nele))
387 45 : ALLOCATE (nnp_env%n_rad(nele))
388 30 : ALLOCATE (nnp_env%n_ang(nele))
389 45 : ALLOCATE (nnp_env%actfnct(nnp_env%n_hlayer + 1))
390 45 : ALLOCATE (cactfnct(nnp_env%n_hlayer + 1))
391 30 : ALLOCATE (nnp_env%ele(nele))
392 30 : ALLOCATE (nnp_env%nuc_ele(nele))
393 76 : ALLOCATE (nnp_env%arc(nele))
394 46 : DO i = 1, nele
395 217 : ALLOCATE (nnp_env%arc(i)%layer(nnp_env%n_layer))
396 108 : ALLOCATE (nnp_env%arc(i)%n_nodes(nnp_env%n_layer))
397 : END DO
398 45 : ALLOCATE (nnp_env%n_hnodes(nnp_env%n_hlayer))
399 45 : ALLOCATE (nnp_env%atom_energies(nele))
400 46 : nnp_env%atom_energies = 0.0_dp
401 :
402 : ! read elements, broadcast and sort
403 15 : CALL parser_reset(parser)
404 : DO
405 30 : CALL parser_search_string(parser, "elements", .TRUE., found, line)
406 30 : IF (found) THEN
407 30 : READ (line, *) dummy
408 30 : IF (TRIM(ADJUSTL(dummy)) == "elements") THEN
409 15 : READ (line, *) dummy, nnp_env%ele(:)
410 15 : CALL nnp_sort_ele(nnp_env%ele, nnp_env%nuc_ele)
411 15 : EXIT
412 : END IF
413 : ELSE
414 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
415 0 : "| elements not specified in NNP_INPUT_FILE")
416 : END IF
417 : END DO
418 :
419 : CALL parser_search_string(parser, "remove_atom_energies", .TRUE., atom_e_found, &
420 15 : search_from_begin_of_file=.TRUE.)
421 :
422 15 : IF (atom_e_found) THEN
423 0 : CALL parser_reset(parser)
424 0 : i = 0
425 : DO
426 0 : CALL parser_search_string(parser, "atom_energy", .TRUE., found, line)
427 0 : IF (found) THEN
428 0 : READ (line, *) dummy, ele, energy
429 0 : DO j = 1, nele
430 0 : IF (nnp_env%ele(j) == TRIM(ele)) THEN
431 0 : i = i + 1
432 0 : nnp_env%atom_energies(j) = energy
433 : END IF
434 : END DO
435 0 : IF (i == nele) EXIT
436 : ELSE
437 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
438 0 : "| atom energies are not specified")
439 : END IF
440 : END DO
441 : END IF
442 :
443 : CALL parser_search_string(parser, "global_nodes_short", .TRUE., found, line, &
444 15 : search_from_begin_of_file=.TRUE.)
445 15 : IF (found) THEN
446 15 : READ (line, *) dummy, nnp_env%n_hnodes(:)
447 : ELSE
448 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
449 0 : "NNP| global_nodes_short not specified in NNP_INPUT_FILE")
450 : END IF
451 :
452 : CALL parser_search_string(parser, "global_activation_short", .TRUE., found, line, &
453 15 : search_from_begin_of_file=.TRUE.)
454 15 : IF (found) THEN
455 15 : READ (line, *) dummy, cactfnct(:)
456 : ELSE
457 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
458 0 : "| global_activation_short not specified in NNP_INPUT_FILE")
459 : END IF
460 :
461 60 : DO i = 1, nnp_env%n_hlayer + 1
462 15 : SELECT CASE (cactfnct(i))
463 : CASE ("t")
464 30 : nnp_env%actfnct(i) = nnp_actfnct_tanh
465 : CASE ("g")
466 0 : nnp_env%actfnct(i) = nnp_actfnct_gaus
467 : CASE ("l")
468 15 : nnp_env%actfnct(i) = nnp_actfnct_lin
469 : CASE ("c")
470 0 : nnp_env%actfnct(i) = nnp_actfnct_cos
471 : CASE ("s")
472 0 : nnp_env%actfnct(i) = nnp_actfnct_sig
473 : CASE ("S")
474 0 : nnp_env%actfnct(i) = nnp_actfnct_invsig
475 : CASE ("e")
476 0 : nnp_env%actfnct(i) = nnp_actfnct_exp
477 : CASE ("p")
478 0 : nnp_env%actfnct(i) = nnp_actfnct_softplus
479 : CASE ("h")
480 0 : nnp_env%actfnct(i) = nnp_actfnct_quad
481 : CASE DEFAULT
482 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
483 45 : "| Activation function unkown")
484 : END SELECT
485 : END DO
486 :
487 : ! determine n_rad and n_ang
488 46 : DO i = 1, nele
489 31 : nnp_env%n_rad(i) = 0
490 46 : nnp_env%n_ang(i) = 0
491 : END DO
492 :
493 : ! count symfunctions
494 15 : CALL parser_reset(parser)
495 15 : first = .TRUE.
496 : DO
497 871 : CALL parser_search_string(parser, "symfunction_short", .TRUE., found, line)
498 871 : IF (found) THEN
499 856 : READ (line, *) dummy, ele, symfnct_type
500 2626 : DO i = 1, nele
501 2626 : IF (TRIM(ele) .EQ. nnp_env%ele(i)) THEN
502 856 : IF (symfnct_type .EQ. 2) THEN
503 486 : nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
504 370 : ELSE IF (symfnct_type .EQ. 3) THEN
505 370 : nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
506 : ELSE
507 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
508 0 : "| Symmetry function type not supported")
509 : END IF
510 : END IF
511 : END DO
512 : first = .FALSE.
513 : ELSE
514 15 : IF (first) CALL cp_abort(__LOCATION__, TRIM(printtag)// &
515 0 : "| no symfunction_short specified in NNP_INPUT_FILE")
516 : ! no additional symfnct found
517 : EXIT
518 : END IF
519 : END DO
520 :
521 46 : DO i = 1, nele
522 93 : ALLOCATE (nnp_env%rad(i)%y(nnp_env%n_rad(i)))
523 93 : ALLOCATE (nnp_env%rad(i)%funccut(nnp_env%n_rad(i)))
524 93 : ALLOCATE (nnp_env%rad(i)%eta(nnp_env%n_rad(i)))
525 93 : ALLOCATE (nnp_env%rad(i)%rs(nnp_env%n_rad(i)))
526 93 : ALLOCATE (nnp_env%rad(i)%loc_min(nnp_env%n_rad(i)))
527 93 : ALLOCATE (nnp_env%rad(i)%loc_max(nnp_env%n_rad(i)))
528 93 : ALLOCATE (nnp_env%rad(i)%loc_av(nnp_env%n_rad(i)))
529 93 : ALLOCATE (nnp_env%rad(i)%sigma(nnp_env%n_rad(i)))
530 62 : ALLOCATE (nnp_env%rad(i)%ele(nnp_env%n_rad(i)))
531 93 : ALLOCATE (nnp_env%rad(i)%nuc_ele(nnp_env%n_rad(i)))
532 517 : nnp_env%rad(i)%funccut = 0.0_dp
533 517 : nnp_env%rad(i)%eta = 0.0_dp
534 517 : nnp_env%rad(i)%rs = 0.0_dp
535 517 : nnp_env%rad(i)%ele = 'X'
536 517 : nnp_env%rad(i)%nuc_ele = 0
537 :
538 93 : ALLOCATE (nnp_env%ang(i)%y(nnp_env%n_ang(i)))
539 93 : ALLOCATE (nnp_env%ang(i)%funccut(nnp_env%n_ang(i)))
540 93 : ALLOCATE (nnp_env%ang(i)%eta(nnp_env%n_ang(i)))
541 93 : ALLOCATE (nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
542 93 : ALLOCATE (nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)))
543 93 : ALLOCATE (nnp_env%ang(i)%lam(nnp_env%n_ang(i)))
544 93 : ALLOCATE (nnp_env%ang(i)%loc_min(nnp_env%n_ang(i)))
545 93 : ALLOCATE (nnp_env%ang(i)%loc_max(nnp_env%n_ang(i)))
546 93 : ALLOCATE (nnp_env%ang(i)%loc_av(nnp_env%n_ang(i)))
547 93 : ALLOCATE (nnp_env%ang(i)%sigma(nnp_env%n_ang(i)))
548 62 : ALLOCATE (nnp_env%ang(i)%ele1(nnp_env%n_ang(i)))
549 62 : ALLOCATE (nnp_env%ang(i)%ele2(nnp_env%n_ang(i)))
550 93 : ALLOCATE (nnp_env%ang(i)%nuc_ele1(nnp_env%n_ang(i)))
551 93 : ALLOCATE (nnp_env%ang(i)%nuc_ele2(nnp_env%n_ang(i)))
552 401 : nnp_env%ang(i)%funccut = 0.0_dp
553 401 : nnp_env%ang(i)%eta = 0.0_dp
554 401 : nnp_env%ang(i)%zeta = 0.0_dp
555 401 : nnp_env%ang(i)%prefzeta = 1.0_dp
556 401 : nnp_env%ang(i)%lam = 0.0_dp
557 401 : nnp_env%ang(i)%ele1 = 'X'
558 401 : nnp_env%ang(i)%ele2 = 'X'
559 401 : nnp_env%ang(i)%nuc_ele1 = 0
560 401 : nnp_env%ang(i)%nuc_ele2 = 0
561 :
562 : ! set number of nodes
563 31 : nnp_env%arc(i)%n_nodes(1) = nnp_env%n_rad(i) + nnp_env%n_ang(i)
564 93 : nnp_env%arc(i)%n_nodes(2:nnp_env%n_layer - 1) = nnp_env%n_hnodes
565 31 : nnp_env%arc(i)%n_nodes(nnp_env%n_layer) = 1
566 170 : DO j = 1, nnp_env%n_layer
567 372 : ALLOCATE (nnp_env%arc(i)%layer(j)%node(nnp_env%arc(i)%n_nodes(j)))
568 372 : ALLOCATE (nnp_env%arc(i)%layer(j)%node_grad(nnp_env%arc(i)%n_nodes(j)))
569 527 : ALLOCATE (nnp_env%arc(i)%layer(j)%tmp_der(nnp_env%arc(i)%n_nodes(1), nnp_env%arc(i)%n_nodes(j)))
570 : END DO
571 : END DO
572 :
573 : ! read, bcast and sort symfnct parameters
574 46 : DO i = 1, nele
575 31 : nnp_env%n_rad(i) = 0
576 46 : nnp_env%n_ang(i) = 0
577 : END DO
578 15 : CALL parser_reset(parser)
579 15 : first = .TRUE.
580 15 : nnp_env%max_cut = 0.0_dp
581 : DO
582 871 : CALL parser_search_string(parser, "symfunction_short", .TRUE., found, line)
583 871 : IF (found) THEN
584 856 : READ (line, *) dummy, ele, symfnct_type
585 2626 : DO i = 1, nele
586 2626 : IF (TRIM(ele) .EQ. nnp_env%ele(i)) THEN
587 856 : IF (symfnct_type .EQ. 2) THEN
588 486 : nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
589 486 : READ (line, *) dummy, ele, symfnct_type, &
590 486 : nnp_env%rad(i)%ele(nnp_env%n_rad(i)), &
591 486 : nnp_env%rad(i)%eta(nnp_env%n_rad(i)), &
592 486 : nnp_env%rad(i)%rs(nnp_env%n_rad(i)), &
593 972 : nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
594 486 : IF (nnp_env%max_cut < nnp_env%rad(i)%funccut(nnp_env%n_rad(i))) THEN
595 16 : nnp_env%max_cut = nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
596 : END IF
597 370 : ELSE IF (symfnct_type .EQ. 3) THEN
598 370 : nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
599 370 : READ (line, *) dummy, ele, symfnct_type, &
600 370 : nnp_env%ang(i)%ele1(nnp_env%n_ang(i)), &
601 370 : nnp_env%ang(i)%ele2(nnp_env%n_ang(i)), &
602 370 : nnp_env%ang(i)%eta(nnp_env%n_ang(i)), &
603 370 : nnp_env%ang(i)%lam(nnp_env%n_ang(i)), &
604 370 : nnp_env%ang(i)%zeta(nnp_env%n_ang(i)), &
605 740 : nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
606 : nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)) = &
607 370 : 2.0_dp**(1.0_dp - nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
608 370 : IF (nnp_env%max_cut < nnp_env%ang(i)%funccut(nnp_env%n_ang(i))) THEN
609 0 : nnp_env%max_cut = nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
610 : END IF
611 : ELSE
612 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
613 0 : "| Symmetry function type not supported")
614 : END IF
615 : END IF
616 : END DO
617 : first = .FALSE.
618 : ELSE
619 15 : IF (first) CALL cp_abort(__LOCATION__, TRIM(printtag)// &
620 0 : "| no symfunction_short specified in NNP_INPUT_FILE")
621 : ! no additional symfnct found
622 : EXIT
623 : END IF
624 : END DO
625 :
626 46 : DO i = 1, nele
627 517 : DO j = 1, nnp_env%n_rad(i)
628 517 : CALL get_ptable_info(nnp_env%rad(i)%ele(j), number=nnp_env%rad(i)%nuc_ele(j))
629 : END DO
630 416 : DO j = 1, nnp_env%n_ang(i)
631 370 : CALL get_ptable_info(nnp_env%ang(i)%ele1(j), number=nnp_env%ang(i)%nuc_ele1(j))
632 370 : CALL get_ptable_info(nnp_env%ang(i)%ele2(j), number=nnp_env%ang(i)%nuc_ele2(j))
633 : ! sort ele1 and ele2
634 401 : IF (nnp_env%ang(i)%nuc_ele1(j) .GT. nnp_env%ang(i)%nuc_ele2(j)) THEN
635 164 : ele = nnp_env%ang(i)%ele1(j)
636 164 : nnp_env%ang(i)%ele1(j) = nnp_env%ang(i)%ele2(j)
637 164 : nnp_env%ang(i)%ele2(j) = ele
638 164 : nuc_ele = nnp_env%ang(i)%nuc_ele1(j)
639 164 : nnp_env%ang(i)%nuc_ele1(j) = nnp_env%ang(i)%nuc_ele2(j)
640 164 : nnp_env%ang(i)%nuc_ele2(j) = nuc_ele
641 : END IF
642 : END DO
643 : END DO
644 : ! Done with input.nn file
645 15 : CALL parser_release(parser)
646 :
647 : ! sort symmetry functions and output information
648 15 : CALL nnp_sort_acsf(nnp_env)
649 15 : CALL nnp_write_acsf(nnp_env, logger%para_env, printtag)
650 15 : CALL nnp_write_arc(nnp_env, logger%para_env, printtag)
651 :
652 : ! read scaling information from file
653 15 : IF (nnp_env%scale_acsf .OR. nnp_env%center_acsf .OR. nnp_env%scale_sigma_acsf) THEN
654 15 : IF (logger%para_env%is_source()) THEN
655 8 : WRITE (unit_nr, *) TRIM(printtag)//"| Reading scaling information from file: ", TRIM(file_name)
656 : END IF
657 : CALL section_vals_val_get(nnp_env%nnp_input, "SCALE_FILE_NAME", &
658 15 : c_val=file_name)
659 15 : CALL parser_create(parser, file_name, para_env=logger%para_env)
660 :
661 : ! Get number of elements in scaling file
662 15 : CALL parser_read_line(parser, 1)
663 15 : k = 0
664 119 : DO WHILE (k < 7)
665 105 : READ (parser%input_line, *, IOSTAT=io) test_array(1:k)
666 105 : IF (io == -1) EXIT
667 105 : k = k + 1
668 : END DO
669 15 : k = k - 1
670 :
671 15 : IF (k == 5 .AND. nnp_env%scale_sigma_acsf) THEN
672 0 : CPABORT("Sigma scaling requested, but scaling.data does not contain sigma.")
673 : END IF
674 :
675 15 : CALL parser_reset(parser)
676 46 : DO i = 1, nnp_env%n_ele
677 517 : DO j = 1, nnp_env%n_rad(i)
678 486 : CALL parser_read_line(parser, 1)
679 517 : IF (nnp_env%scale_sigma_acsf) THEN
680 0 : READ (parser%input_line, *) dummy, dummy, &
681 0 : nnp_env%rad(i)%loc_min(j), &
682 0 : nnp_env%rad(i)%loc_max(j), &
683 0 : nnp_env%rad(i)%loc_av(j), &
684 0 : nnp_env%rad(i)%sigma(j)
685 : ELSE
686 486 : READ (parser%input_line, *) dummy, dummy, &
687 486 : nnp_env%rad(i)%loc_min(j), &
688 486 : nnp_env%rad(i)%loc_max(j), &
689 972 : nnp_env%rad(i)%loc_av(j)
690 : END IF
691 : END DO
692 416 : DO j = 1, nnp_env%n_ang(i)
693 370 : CALL parser_read_line(parser, 1)
694 401 : IF (nnp_env%scale_sigma_acsf) THEN
695 0 : READ (parser%input_line, *) dummy, dummy, &
696 0 : nnp_env%ang(i)%loc_min(j), &
697 0 : nnp_env%ang(i)%loc_max(j), &
698 0 : nnp_env%ang(i)%loc_av(j), &
699 0 : nnp_env%ang(i)%sigma(j)
700 : ELSE
701 370 : READ (parser%input_line, *) dummy, dummy, &
702 370 : nnp_env%ang(i)%loc_min(j), &
703 370 : nnp_env%ang(i)%loc_max(j), &
704 740 : nnp_env%ang(i)%loc_av(j)
705 : END IF
706 : END DO
707 : END DO
708 15 : CALL parser_release(parser)
709 : END IF
710 :
711 15 : CALL nnp_init_acsf_groups(nnp_env)
712 :
713 : ! read weights from file
714 46 : DO i = 1, nnp_env%n_ele
715 139 : DO j = 2, nnp_env%n_layer
716 0 : ALLOCATE (nnp_env%arc(i)%layer(j)%weights(nnp_env%arc(i)%n_nodes(j - 1), &
717 465 : nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
718 403 : ALLOCATE (nnp_env%arc(i)%layer(j)%bweights(nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
719 : END DO
720 : END DO
721 114 : DO i_com = 1, nnp_env%n_committee
722 99 : CALL section_vals_val_get(model_section, "WEIGHTS", c_val=base_name, i_rep_section=i_com)
723 99 : IF (logger%para_env%is_source()) THEN
724 50 : WRITE (unit_nr, *) TRIM(printtag)//"| Initializing weights for model: ", i_com
725 : END IF
726 313 : DO i = 1, nnp_env%n_ele
727 199 : WRITE (file_name, '(A,I0.3,A)') TRIM(base_name)//".", nnp_env%nuc_ele(i), ".data"
728 199 : IF (logger%para_env%is_source()) THEN
729 101 : WRITE (unit_nr, *) TRIM(printtag)//"| Reading weights from file: ", TRIM(file_name)
730 : END IF
731 199 : CALL parser_create(parser, file_name, para_env=logger%para_env)
732 199 : n_weight = 0
733 205629 : DO WHILE (.TRUE.)
734 205828 : CALL parser_read_line(parser, 1, at_end)
735 205828 : IF (at_end) EXIT
736 205629 : n_weight = n_weight + 1
737 : END DO
738 :
739 597 : ALLOCATE (weights(n_weight))
740 :
741 199 : CALL parser_reset(parser)
742 205828 : DO j = 1, n_weight
743 205629 : CALL parser_read_line(parser, 1)
744 205828 : READ (parser%input_line, *) weights(j)
745 : END DO
746 199 : CALL parser_release(parser)
747 :
748 : ! sort weights into corresponding arrays
749 199 : iweight = 0
750 796 : DO j = 2, nnp_env%n_layer
751 14231 : DO k = 1, nnp_env%arc(i)%n_nodes(j - 1)
752 211671 : DO l = 1, nnp_env%arc(i)%n_nodes(j)
753 197440 : iweight = iweight + 1
754 211074 : nnp_env%arc(i)%layer(j)%weights(k, l, i_com) = weights(iweight)
755 : END DO
756 : END DO
757 :
758 8985 : DO k = 1, nnp_env%arc(i)%n_nodes(j)
759 8189 : iweight = iweight + 1
760 8786 : nnp_env%arc(i)%layer(j)%bweights(k, i_com) = weights(iweight)
761 : END DO
762 : END DO
763 :
764 497 : DEALLOCATE (weights)
765 : END DO
766 : END DO
767 :
768 : !Initialize extrapolation counter
769 15 : nnp_env%expol = 0
770 :
771 : ! Bias the standard deviation of committee disagreement
772 15 : NULLIFY (bias_section)
773 15 : explicit = .FALSE.
774 : !HELIUM NNP does atm not allow for bias (not even defined)
775 15 : bias_section => section_vals_get_subs_vals(nnp_env%nnp_input, "BIAS", can_return_null=.TRUE.)
776 15 : IF (ASSOCIATED(bias_section)) CALL section_vals_get(bias_section, explicit=explicit)
777 15 : nnp_env%bias = .FALSE.
778 15 : IF (explicit) THEN
779 4 : IF (nnp_env%n_committee > 1) THEN
780 4 : IF (logger%para_env%is_source()) THEN
781 2 : WRITE (unit_nr, *) "NNP| Biasing of committee disagreement enabled"
782 : END IF
783 4 : nnp_env%bias = .TRUE.
784 12 : ALLOCATE (nnp_env%bias_forces(3, nnp_env%num_atoms))
785 12 : ALLOCATE (nnp_env%bias_e_avrg(nnp_env%n_committee))
786 4 : CALL section_vals_val_get(bias_section, "SIGMA_0", r_val=nnp_env%bias_sigma0)
787 4 : CALL section_vals_val_get(bias_section, "K_B", r_val=nnp_env%bias_kb)
788 36 : nnp_env%bias_e_avrg(:) = 0.0_dp
789 4 : CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", explicit=explicit)
790 4 : nnp_env%bias_align = explicit
791 4 : IF (explicit) THEN
792 4 : NULLIFY (work)
793 4 : CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", r_vals=work)
794 4 : IF (SIZE(work) .NE. nnp_env%n_committee) THEN
795 0 : CPABORT("ALIGN_NNP_ENERGIES size mismatch wrt committee size.")
796 : END IF
797 68 : nnp_env%bias_e_avrg(:) = work
798 4 : IF (logger%para_env%is_source()) THEN
799 2 : WRITE (unit_nr, *) TRIM(printtag)//"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
800 : END IF
801 : END IF
802 : ELSE
803 0 : CPWARN("NNP committee size is 1, BIAS section is ignored.")
804 : END IF
805 : END IF
806 :
807 15 : IF (logger%para_env%is_source()) THEN
808 8 : WRITE (unit_nr, *) TRIM(printtag)//"| NNP force environment initialized"
809 : END IF
810 :
811 15 : CALL timestop(handle)
812 :
813 75 : END SUBROUTINE nnp_init_model
814 :
815 : END MODULE nnp_environment
|
