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