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