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 :
351 : CALL parser_search_string(parser, "normalize_nodes", .TRUE., found, &
352 15 : search_from_begin_of_file=.TRUE.)
353 15 : nnp_env%normnodes = found
354 :
355 : CALL parser_search_string(parser, "cutoff_type", .TRUE., found, line, &
356 15 : search_from_begin_of_file=.TRUE.)
357 15 : IF (found) THEN
358 15 : READ (line, *) dummy, nnp_env%cut_type
359 : ELSE
360 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
361 0 : "| no cutoff type specified in NNP_INPUT_FILE")
362 : END IF
363 :
364 : CALL parser_search_string(parser, "global_hidden_layers_short", .TRUE., found, line, &
365 15 : search_from_begin_of_file=.TRUE.)
366 15 : IF (found) THEN
367 15 : READ (line, *) dummy, nnp_env%n_hlayer
368 : ELSE
369 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
370 0 : "| number of hidden layers missing in NNP_INPUT_FILE")
371 : END IF
372 15 : nnp_env%n_layer = nnp_env%n_hlayer + 2
373 :
374 15 : nele = nnp_env%n_ele
375 76 : ALLOCATE (nnp_env%rad(nele))
376 76 : ALLOCATE (nnp_env%ang(nele))
377 45 : ALLOCATE (nnp_env%n_rad(nele))
378 45 : ALLOCATE (nnp_env%n_ang(nele))
379 45 : ALLOCATE (nnp_env%actfnct(nnp_env%n_hlayer + 1))
380 45 : ALLOCATE (cactfnct(nnp_env%n_hlayer + 1))
381 30 : ALLOCATE (nnp_env%ele(nele))
382 45 : ALLOCATE (nnp_env%nuc_ele(nele))
383 76 : ALLOCATE (nnp_env%arc(nele))
384 46 : DO i = 1, nele
385 217 : ALLOCATE (nnp_env%arc(i)%layer(nnp_env%n_layer))
386 108 : ALLOCATE (nnp_env%arc(i)%n_nodes(nnp_env%n_layer))
387 : END DO
388 45 : ALLOCATE (nnp_env%n_hnodes(nnp_env%n_hlayer))
389 45 : ALLOCATE (nnp_env%atom_energies(nele))
390 46 : nnp_env%atom_energies = 0.0_dp
391 :
392 : ! read elements, broadcast and sort
393 15 : CALL parser_reset(parser)
394 : DO
395 30 : CALL parser_search_string(parser, "elements", .TRUE., found, line)
396 30 : IF (found) THEN
397 30 : READ (line, *) dummy
398 30 : IF (TRIM(ADJUSTL(dummy)) == "elements") THEN
399 15 : READ (line, *) dummy, nnp_env%ele(:)
400 15 : CALL nnp_sort_ele(nnp_env%ele, nnp_env%nuc_ele)
401 15 : EXIT
402 : END IF
403 : ELSE
404 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
405 0 : "| elements not specified in NNP_INPUT_FILE")
406 : END IF
407 : END DO
408 :
409 : CALL parser_search_string(parser, "remove_atom_energies", .TRUE., atom_e_found, &
410 15 : search_from_begin_of_file=.TRUE.)
411 :
412 15 : IF (atom_e_found) THEN
413 0 : CALL parser_reset(parser)
414 0 : i = 0
415 : DO
416 0 : CALL parser_search_string(parser, "atom_energy", .TRUE., found, line)
417 0 : IF (found) THEN
418 0 : READ (line, *) dummy, ele, energy
419 0 : DO j = 1, nele
420 0 : IF (nnp_env%ele(j) == TRIM(ele)) THEN
421 0 : i = i + 1
422 0 : nnp_env%atom_energies(j) = energy
423 : END IF
424 : END DO
425 0 : IF (i == nele) EXIT
426 : ELSE
427 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
428 0 : "| atom energies are not specified")
429 : END IF
430 : END DO
431 : END IF
432 :
433 : CALL parser_search_string(parser, "global_nodes_short", .TRUE., found, line, &
434 15 : search_from_begin_of_file=.TRUE.)
435 15 : IF (found) THEN
436 15 : READ (line, *) dummy, nnp_env%n_hnodes(:)
437 : ELSE
438 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
439 0 : "NNP| global_nodes_short not specified in NNP_INPUT_FILE")
440 : END IF
441 :
442 : CALL parser_search_string(parser, "global_activation_short", .TRUE., found, line, &
443 15 : search_from_begin_of_file=.TRUE.)
444 15 : IF (found) THEN
445 15 : READ (line, *) dummy, cactfnct(:)
446 : ELSE
447 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
448 0 : "| global_activation_short not specified in NNP_INPUT_FILE")
449 : END IF
450 :
451 60 : DO i = 1, nnp_env%n_hlayer + 1
452 15 : SELECT CASE (cactfnct(i))
453 : CASE ("t")
454 30 : nnp_env%actfnct(i) = nnp_actfnct_tanh
455 : CASE ("g")
456 0 : nnp_env%actfnct(i) = nnp_actfnct_gaus
457 : CASE ("l")
458 15 : nnp_env%actfnct(i) = nnp_actfnct_lin
459 : CASE ("c")
460 0 : nnp_env%actfnct(i) = nnp_actfnct_cos
461 : CASE ("s")
462 0 : nnp_env%actfnct(i) = nnp_actfnct_sig
463 : CASE ("S")
464 0 : nnp_env%actfnct(i) = nnp_actfnct_invsig
465 : CASE ("e")
466 0 : nnp_env%actfnct(i) = nnp_actfnct_exp
467 : CASE ("p")
468 0 : nnp_env%actfnct(i) = nnp_actfnct_softplus
469 : CASE ("h")
470 0 : nnp_env%actfnct(i) = nnp_actfnct_quad
471 : CASE DEFAULT
472 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
473 45 : "| Activation function unkown")
474 : END SELECT
475 : END DO
476 :
477 : ! determine n_rad and n_ang
478 46 : DO i = 1, nele
479 31 : nnp_env%n_rad(i) = 0
480 46 : nnp_env%n_ang(i) = 0
481 : END DO
482 :
483 : ! count symfunctions
484 15 : CALL parser_reset(parser)
485 15 : first = .TRUE.
486 : DO
487 871 : CALL parser_search_string(parser, "symfunction_short", .TRUE., found, line)
488 871 : IF (found) THEN
489 856 : READ (line, *) dummy, ele, symfnct_type
490 2626 : DO i = 1, nele
491 2626 : IF (TRIM(ele) .EQ. nnp_env%ele(i)) THEN
492 856 : IF (symfnct_type .EQ. 2) THEN
493 486 : nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
494 370 : ELSE IF (symfnct_type .EQ. 3) THEN
495 370 : nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
496 : ELSE
497 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
498 0 : "| Symmetry function type not supported")
499 : END IF
500 : END IF
501 : END DO
502 : first = .FALSE.
503 : ELSE
504 15 : IF (first) CALL cp_abort(__LOCATION__, TRIM(printtag)// &
505 0 : "| no symfunction_short specified in NNP_INPUT_FILE")
506 : ! no additional symfnct found
507 : EXIT
508 : END IF
509 : END DO
510 :
511 46 : DO i = 1, nele
512 93 : ALLOCATE (nnp_env%rad(i)%y(nnp_env%n_rad(i)))
513 93 : ALLOCATE (nnp_env%rad(i)%funccut(nnp_env%n_rad(i)))
514 93 : ALLOCATE (nnp_env%rad(i)%eta(nnp_env%n_rad(i)))
515 93 : ALLOCATE (nnp_env%rad(i)%rs(nnp_env%n_rad(i)))
516 93 : ALLOCATE (nnp_env%rad(i)%loc_min(nnp_env%n_rad(i)))
517 93 : ALLOCATE (nnp_env%rad(i)%loc_max(nnp_env%n_rad(i)))
518 93 : ALLOCATE (nnp_env%rad(i)%loc_av(nnp_env%n_rad(i)))
519 93 : ALLOCATE (nnp_env%rad(i)%sigma(nnp_env%n_rad(i)))
520 62 : ALLOCATE (nnp_env%rad(i)%ele(nnp_env%n_rad(i)))
521 93 : ALLOCATE (nnp_env%rad(i)%nuc_ele(nnp_env%n_rad(i)))
522 517 : nnp_env%rad(i)%funccut = 0.0_dp
523 517 : nnp_env%rad(i)%eta = 0.0_dp
524 517 : nnp_env%rad(i)%rs = 0.0_dp
525 517 : nnp_env%rad(i)%ele = 'X'
526 517 : nnp_env%rad(i)%nuc_ele = 0
527 :
528 93 : ALLOCATE (nnp_env%ang(i)%y(nnp_env%n_ang(i)))
529 93 : ALLOCATE (nnp_env%ang(i)%funccut(nnp_env%n_ang(i)))
530 93 : ALLOCATE (nnp_env%ang(i)%eta(nnp_env%n_ang(i)))
531 93 : ALLOCATE (nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
532 93 : ALLOCATE (nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)))
533 93 : ALLOCATE (nnp_env%ang(i)%lam(nnp_env%n_ang(i)))
534 93 : ALLOCATE (nnp_env%ang(i)%loc_min(nnp_env%n_ang(i)))
535 93 : ALLOCATE (nnp_env%ang(i)%loc_max(nnp_env%n_ang(i)))
536 93 : ALLOCATE (nnp_env%ang(i)%loc_av(nnp_env%n_ang(i)))
537 93 : ALLOCATE (nnp_env%ang(i)%sigma(nnp_env%n_ang(i)))
538 62 : ALLOCATE (nnp_env%ang(i)%ele1(nnp_env%n_ang(i)))
539 62 : ALLOCATE (nnp_env%ang(i)%ele2(nnp_env%n_ang(i)))
540 93 : ALLOCATE (nnp_env%ang(i)%nuc_ele1(nnp_env%n_ang(i)))
541 93 : ALLOCATE (nnp_env%ang(i)%nuc_ele2(nnp_env%n_ang(i)))
542 401 : nnp_env%ang(i)%funccut = 0.0_dp
543 401 : nnp_env%ang(i)%eta = 0.0_dp
544 401 : nnp_env%ang(i)%zeta = 0.0_dp
545 401 : nnp_env%ang(i)%prefzeta = 1.0_dp
546 401 : nnp_env%ang(i)%lam = 0.0_dp
547 401 : nnp_env%ang(i)%ele1 = 'X'
548 401 : nnp_env%ang(i)%ele2 = 'X'
549 401 : nnp_env%ang(i)%nuc_ele1 = 0
550 401 : nnp_env%ang(i)%nuc_ele2 = 0
551 :
552 : ! set number of nodes
553 31 : nnp_env%arc(i)%n_nodes(1) = nnp_env%n_rad(i) + nnp_env%n_ang(i)
554 93 : nnp_env%arc(i)%n_nodes(2:nnp_env%n_layer - 1) = nnp_env%n_hnodes
555 31 : nnp_env%arc(i)%n_nodes(nnp_env%n_layer) = 1
556 170 : DO j = 1, nnp_env%n_layer
557 372 : ALLOCATE (nnp_env%arc(i)%layer(j)%node(nnp_env%arc(i)%n_nodes(j)))
558 372 : ALLOCATE (nnp_env%arc(i)%layer(j)%node_grad(nnp_env%arc(i)%n_nodes(j)))
559 527 : ALLOCATE (nnp_env%arc(i)%layer(j)%tmp_der(nnp_env%arc(i)%n_nodes(1), nnp_env%arc(i)%n_nodes(j)))
560 : END DO
561 : END DO
562 :
563 : ! read, bcast and sort symfnct parameters
564 46 : DO i = 1, nele
565 31 : nnp_env%n_rad(i) = 0
566 46 : nnp_env%n_ang(i) = 0
567 : END DO
568 15 : CALL parser_reset(parser)
569 15 : first = .TRUE.
570 15 : nnp_env%max_cut = 0.0_dp
571 : DO
572 871 : CALL parser_search_string(parser, "symfunction_short", .TRUE., found, line)
573 871 : IF (found) THEN
574 856 : READ (line, *) dummy, ele, symfnct_type
575 2626 : DO i = 1, nele
576 2626 : IF (TRIM(ele) .EQ. nnp_env%ele(i)) THEN
577 856 : IF (symfnct_type .EQ. 2) THEN
578 486 : nnp_env%n_rad(i) = nnp_env%n_rad(i) + 1
579 486 : READ (line, *) dummy, ele, symfnct_type, &
580 486 : nnp_env%rad(i)%ele(nnp_env%n_rad(i)), &
581 486 : nnp_env%rad(i)%eta(nnp_env%n_rad(i)), &
582 486 : nnp_env%rad(i)%rs(nnp_env%n_rad(i)), &
583 972 : nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
584 486 : IF (nnp_env%max_cut < nnp_env%rad(i)%funccut(nnp_env%n_rad(i))) THEN
585 16 : nnp_env%max_cut = nnp_env%rad(i)%funccut(nnp_env%n_rad(i))
586 : END IF
587 370 : ELSE IF (symfnct_type .EQ. 3) THEN
588 370 : nnp_env%n_ang(i) = nnp_env%n_ang(i) + 1
589 370 : READ (line, *) dummy, ele, symfnct_type, &
590 370 : nnp_env%ang(i)%ele1(nnp_env%n_ang(i)), &
591 370 : nnp_env%ang(i)%ele2(nnp_env%n_ang(i)), &
592 370 : nnp_env%ang(i)%eta(nnp_env%n_ang(i)), &
593 370 : nnp_env%ang(i)%lam(nnp_env%n_ang(i)), &
594 370 : nnp_env%ang(i)%zeta(nnp_env%n_ang(i)), &
595 740 : nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
596 : nnp_env%ang(i)%prefzeta(nnp_env%n_ang(i)) = &
597 370 : 2.0_dp**(1.0_dp - nnp_env%ang(i)%zeta(nnp_env%n_ang(i)))
598 370 : IF (nnp_env%max_cut < nnp_env%ang(i)%funccut(nnp_env%n_ang(i))) THEN
599 0 : nnp_env%max_cut = nnp_env%ang(i)%funccut(nnp_env%n_ang(i))
600 : END IF
601 : ELSE
602 : CALL cp_abort(__LOCATION__, TRIM(printtag)// &
603 0 : "| Symmetry function type not supported")
604 : END IF
605 : END IF
606 : END DO
607 : first = .FALSE.
608 : ELSE
609 15 : IF (first) CALL cp_abort(__LOCATION__, TRIM(printtag)// &
610 0 : "| no symfunction_short specified in NNP_INPUT_FILE")
611 : ! no additional symfnct found
612 : EXIT
613 : END IF
614 : END DO
615 :
616 46 : DO i = 1, nele
617 517 : DO j = 1, nnp_env%n_rad(i)
618 517 : CALL get_ptable_info(nnp_env%rad(i)%ele(j), number=nnp_env%rad(i)%nuc_ele(j))
619 : END DO
620 416 : DO j = 1, nnp_env%n_ang(i)
621 370 : CALL get_ptable_info(nnp_env%ang(i)%ele1(j), number=nnp_env%ang(i)%nuc_ele1(j))
622 370 : CALL get_ptable_info(nnp_env%ang(i)%ele2(j), number=nnp_env%ang(i)%nuc_ele2(j))
623 : ! sort ele1 and ele2
624 401 : IF (nnp_env%ang(i)%nuc_ele1(j) .GT. nnp_env%ang(i)%nuc_ele2(j)) THEN
625 164 : ele = nnp_env%ang(i)%ele1(j)
626 164 : nnp_env%ang(i)%ele1(j) = nnp_env%ang(i)%ele2(j)
627 164 : nnp_env%ang(i)%ele2(j) = ele
628 164 : nuc_ele = nnp_env%ang(i)%nuc_ele1(j)
629 164 : nnp_env%ang(i)%nuc_ele1(j) = nnp_env%ang(i)%nuc_ele2(j)
630 164 : nnp_env%ang(i)%nuc_ele2(j) = nuc_ele
631 : END IF
632 : END DO
633 : END DO
634 : ! Done with input.nn file
635 15 : CALL parser_release(parser)
636 :
637 : ! sort symmetry functions and output information
638 15 : CALL nnp_sort_acsf(nnp_env)
639 15 : CALL nnp_write_acsf(nnp_env, logger%para_env, printtag)
640 15 : CALL nnp_write_arc(nnp_env, logger%para_env, printtag)
641 :
642 : ! read scaling information from file
643 15 : IF (nnp_env%scale_acsf .OR. nnp_env%center_acsf .OR. nnp_env%scale_sigma_acsf) THEN
644 15 : IF (logger%para_env%is_source()) THEN
645 8 : WRITE (unit_nr, *) TRIM(printtag)//"| Reading scaling information from file: ", TRIM(file_name)
646 : END IF
647 : CALL section_vals_val_get(nnp_env%nnp_input, "SCALE_FILE_NAME", &
648 15 : c_val=file_name)
649 15 : CALL parser_create(parser, file_name, para_env=logger%para_env)
650 :
651 : ! Get number of elements in scaling file
652 15 : CALL parser_read_line(parser, 1)
653 15 : k = 0
654 119 : DO WHILE (k < 7)
655 105 : READ (parser%input_line, *, IOSTAT=io) test_array(1:k)
656 105 : IF (io == -1) EXIT
657 105 : k = k + 1
658 : END DO
659 15 : k = k - 1
660 :
661 15 : IF (k == 5 .AND. nnp_env%scale_sigma_acsf) THEN
662 0 : CPABORT("Sigma scaling requested, but scaling.data does not contain sigma.")
663 : END IF
664 :
665 15 : CALL parser_reset(parser)
666 46 : DO i = 1, nnp_env%n_ele
667 517 : DO j = 1, nnp_env%n_rad(i)
668 486 : CALL parser_read_line(parser, 1)
669 517 : IF (nnp_env%scale_sigma_acsf) THEN
670 0 : READ (parser%input_line, *) dummy, dummy, &
671 0 : nnp_env%rad(i)%loc_min(j), &
672 0 : nnp_env%rad(i)%loc_max(j), &
673 0 : nnp_env%rad(i)%loc_av(j), &
674 0 : nnp_env%rad(i)%sigma(j)
675 : ELSE
676 486 : READ (parser%input_line, *) dummy, dummy, &
677 486 : nnp_env%rad(i)%loc_min(j), &
678 486 : nnp_env%rad(i)%loc_max(j), &
679 972 : nnp_env%rad(i)%loc_av(j)
680 : END IF
681 : END DO
682 416 : DO j = 1, nnp_env%n_ang(i)
683 370 : CALL parser_read_line(parser, 1)
684 401 : IF (nnp_env%scale_sigma_acsf) THEN
685 0 : READ (parser%input_line, *) dummy, dummy, &
686 0 : nnp_env%ang(i)%loc_min(j), &
687 0 : nnp_env%ang(i)%loc_max(j), &
688 0 : nnp_env%ang(i)%loc_av(j), &
689 0 : nnp_env%ang(i)%sigma(j)
690 : ELSE
691 370 : READ (parser%input_line, *) dummy, dummy, &
692 370 : nnp_env%ang(i)%loc_min(j), &
693 370 : nnp_env%ang(i)%loc_max(j), &
694 740 : nnp_env%ang(i)%loc_av(j)
695 : END IF
696 : END DO
697 : END DO
698 15 : CALL parser_release(parser)
699 : END IF
700 :
701 15 : CALL nnp_init_acsf_groups(nnp_env)
702 :
703 : ! read weights from file
704 46 : DO i = 1, nnp_env%n_ele
705 139 : DO j = 2, nnp_env%n_layer
706 0 : ALLOCATE (nnp_env%arc(i)%layer(j)%weights(nnp_env%arc(i)%n_nodes(j - 1), &
707 465 : nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
708 403 : ALLOCATE (nnp_env%arc(i)%layer(j)%bweights(nnp_env%arc(i)%n_nodes(j), nnp_env%n_committee))
709 : END DO
710 : END DO
711 114 : DO i_com = 1, nnp_env%n_committee
712 99 : CALL section_vals_val_get(model_section, "WEIGHTS", c_val=base_name, i_rep_section=i_com)
713 99 : IF (logger%para_env%is_source()) THEN
714 50 : WRITE (unit_nr, *) TRIM(printtag)//"| Initializing weights for model: ", i_com
715 : END IF
716 313 : DO i = 1, nnp_env%n_ele
717 199 : WRITE (file_name, '(A,I0.3,A)') TRIM(base_name)//".", nnp_env%nuc_ele(i), ".data"
718 199 : IF (logger%para_env%is_source()) THEN
719 101 : WRITE (unit_nr, *) TRIM(printtag)//"| Reading weights from file: ", TRIM(file_name)
720 : END IF
721 199 : CALL parser_create(parser, file_name, para_env=logger%para_env)
722 199 : n_weight = 0
723 205629 : DO WHILE (.TRUE.)
724 205828 : CALL parser_read_line(parser, 1, at_end)
725 205828 : IF (at_end) EXIT
726 205629 : n_weight = n_weight + 1
727 : END DO
728 :
729 597 : ALLOCATE (weights(n_weight))
730 :
731 199 : CALL parser_reset(parser)
732 205828 : DO j = 1, n_weight
733 205629 : CALL parser_read_line(parser, 1)
734 205828 : READ (parser%input_line, *) weights(j)
735 : END DO
736 199 : CALL parser_release(parser)
737 :
738 : ! sort weights into corresponding arrays
739 199 : iweight = 0
740 796 : DO j = 2, nnp_env%n_layer
741 14231 : DO k = 1, nnp_env%arc(i)%n_nodes(j - 1)
742 211671 : DO l = 1, nnp_env%arc(i)%n_nodes(j)
743 197440 : iweight = iweight + 1
744 211074 : nnp_env%arc(i)%layer(j)%weights(k, l, i_com) = weights(iweight)
745 : END DO
746 : END DO
747 :
748 8985 : DO k = 1, nnp_env%arc(i)%n_nodes(j)
749 8189 : iweight = iweight + 1
750 8786 : nnp_env%arc(i)%layer(j)%bweights(k, i_com) = weights(iweight)
751 : END DO
752 : END DO
753 :
754 298 : DEALLOCATE (weights)
755 : END DO
756 : END DO
757 :
758 : !Initialize extrapolation counter
759 15 : nnp_env%expol = 0
760 :
761 : ! Bias the standard deviation of committee disagreement
762 15 : NULLIFY (bias_section)
763 15 : explicit = .FALSE.
764 : !HELIUM NNP does atm not allow for bias (not even defined)
765 15 : bias_section => section_vals_get_subs_vals(nnp_env%nnp_input, "BIAS", can_return_null=.TRUE.)
766 15 : IF (ASSOCIATED(bias_section)) CALL section_vals_get(bias_section, explicit=explicit)
767 15 : nnp_env%bias = .FALSE.
768 15 : IF (explicit) THEN
769 4 : IF (nnp_env%n_committee > 1) THEN
770 4 : IF (logger%para_env%is_source()) THEN
771 2 : WRITE (unit_nr, *) "NNP| Biasing of committee disagreement enabled"
772 : END IF
773 4 : nnp_env%bias = .TRUE.
774 12 : ALLOCATE (nnp_env%bias_forces(3, nnp_env%num_atoms))
775 12 : ALLOCATE (nnp_env%bias_e_avrg(nnp_env%n_committee))
776 4 : CALL section_vals_val_get(bias_section, "SIGMA_0", r_val=nnp_env%bias_sigma0)
777 4 : CALL section_vals_val_get(bias_section, "K_B", r_val=nnp_env%bias_kb)
778 36 : nnp_env%bias_e_avrg(:) = 0.0_dp
779 4 : CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", explicit=explicit)
780 4 : nnp_env%bias_align = explicit
781 4 : IF (explicit) THEN
782 4 : NULLIFY (work)
783 4 : CALL section_vals_val_get(bias_section, "ALIGN_NNP_ENERGIES", r_vals=work)
784 4 : IF (SIZE(work) .NE. nnp_env%n_committee) THEN
785 0 : CPABORT("ALIGN_NNP_ENERGIES size mismatch wrt committee size.")
786 : END IF
787 68 : nnp_env%bias_e_avrg(:) = work
788 4 : IF (logger%para_env%is_source()) THEN
789 2 : WRITE (unit_nr, *) TRIM(printtag)//"| Biasing is aligned by shifting the energy prediction of the C-NNP members"
790 : END IF
791 : END IF
792 : ELSE
793 0 : CPWARN("NNP committee size is 1, BIAS section is ignored.")
794 : END IF
795 : END IF
796 :
797 15 : IF (logger%para_env%is_source()) THEN
798 8 : WRITE (unit_nr, *) TRIM(printtag)//"| NNP force environment initialized"
799 : END IF
800 :
801 15 : CALL timestop(handle)
802 :
803 60 : END SUBROUTINE nnp_init_model
804 :
805 : END MODULE nnp_environment
|
