(git:4dc5213)
Loading...
Searching...
No Matches
nnp_force.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Methods dealing with Neural Network potentials
10!> \author Christoph Schran (christoph.schran@rub.de)
11!> \author Dhruv Sharma (ds2173@cam.ac.uk)
12!> \author Claudio Malvino (claudiormal@gmail.com)
13!> \date 2020-10-10
14! **************************************************************************************************
16
20 USE cp_output_handling, ONLY: cp_p_file,&
31 USE kinds, ONLY: default_path_length,&
33 dp
34 USE nnp_acsf, ONLY: nnp_calc_acsf,&
39 USE nnp_model, ONLY: nnp_gradients,&
43 USE physcon, ONLY: angstrom
44 USE virial_types, ONLY: virial_type
45
46!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
47#include "./base/base_uses.f90"
48
49 IMPLICIT NONE
50
51 PRIVATE
52
53 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_force'
55
57
58CONTAINS
59
60! **************************************************************************************************
61!> \brief Calculate the energy and force for a given configuration with the NNP
62!> \param nnp ...
63!> \param calc_forces ...
64!> \date 2020-10-10
65!> \author Christoph Schran (christoph.schran@rub.de)
66!> \author Claudio Malvino (claudiormal@gmail.com)
67! **************************************************************************************************
68 SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
69 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
70 LOGICAL, INTENT(IN) :: calc_forces
71
72 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_calc_energy_force'
73
74 INTEGER :: handle, handle_loop, i, i_com, ie, ig, &
75 il, ind, istart, j, k, m, &
76 max_input_nodes, mecalc, &
77 n_input_nodes, nthreads, t, tid
78 INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
79 LOGICAL :: calc_stress
80 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
81 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: stress
82 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: my_cstress_t, my_force_t
83 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
84 TYPE(cp_logger_type), POINTER :: logger
85 TYPE(cp_subsys_type), POINTER :: subsys
86 TYPE(distribution_1d_type), POINTER :: local_particles
87 TYPE(nnp_arc_type), ALLOCATABLE, DIMENSION(:) :: my_arc
88 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
89 TYPE(section_vals_type), POINTER :: print_section
90 TYPE(virial_type), POINTER :: virial
91
92 CALL timeset(routinen, handle)
93
94 NULLIFY (particle_set, logger, local_particles, subsys, &
95 atomic_kind_set)
96 logger => cp_get_default_logger()
97
98 cpassert(ASSOCIATED(nnp))
99 cpassert(nnp%n_committee >= 1)
100 CALL nnp_env_get(nnp_env=nnp, particle_set=particle_set, &
101 subsys=subsys, local_particles=local_particles, &
102 atomic_kind_set=atomic_kind_set)
103
104 CALL cp_subsys_get(subsys, &
105 virial=virial)
106
107 calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
108 IF (calc_stress .AND. .NOT. calc_forces) THEN
109 cpabort('Stress cannot be calculated without forces')
110 END IF
111
112 nnp%atomic_energy(:, :) = 0.0_dp
113 IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
114 IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
115 IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
116
117 !fill coord array
118 ig = 1
119 DO i = 1, nnp%n_ele
120 DO j = 1, nnp%num_atoms
121 IF (nnp%ele(i) == particle_set(j)%atomic_kind%element_symbol) THEN
122 DO m = 1, 3
123 nnp%coord(m, ig) = particle_set(j)%r(m)
124 END DO
125 nnp%atoms(ig) = nnp%ele(i)
126 CALL get_ptable_info(nnp%atoms(ig) (1:2), number=nnp%nuc_atoms(ig))
127 nnp%ele_ind(ig) = i
128 nnp%sort(ig) = j
129 nnp%sort_inv(j) = ig
130 ig = ig + 1
131 END IF
132 END DO
133 END DO
134
135 ! parallization:
136 mecalc = nnp%num_atoms/logger%para_env%num_pe + &
137 min(mod(nnp%num_atoms, logger%para_env%num_pe)/ &
138 (logger%para_env%mepos + 1), 1)
139 ALLOCATE (allcalc(logger%para_env%num_pe))
140 allcalc(:) = 0
141 CALL logger%para_env%allgather(mecalc, allcalc)
142 istart = 1
143 DO i = 2, logger%para_env%mepos + 1
144 istart = istart + allcalc(i - 1)
145 END DO
146
147 ! reset extrapolation status
148 nnp%output_expol = .false.
149
151
152 max_input_nodes = 0
153 DO i = 1, nnp%n_ele
154 max_input_nodes = max(max_input_nodes, nnp%arc(i)%n_nodes(1))
155 END DO
156
157 ! number of threads that will run the atom loop below
158 nthreads = 1
159!$ nthreads = omp_get_max_threads()
160 IF (calc_forces) THEN
161 ! dGdr lives in per-element workspace; no global (3, n_sf, num_atoms) slab.
162 ALLOCATE (my_force_t(3, nnp%num_atoms, nnp%n_committee, 0:nthreads - 1))
163 my_force_t(:, :, :, :) = 0.0_dp
164 IF (calc_stress) THEN
165 ALLOCATE (my_cstress_t(3, 3, nnp%n_committee, 0:nthreads - 1))
166 my_cstress_t(:, :, :, :) = 0.0_dp
167 END IF
168 END IF
169
170 ! Per-rank timer around the per-atom work. CP2K aggregates the label as
171 ! avg/max across ranks, so max - avg measures the load imbalance.
172 CALL timeset('nnp_per_rank_atom_loop', handle_loop)
173
174 ! Thread the loop over atomic centres. Each thread works on its own copy
175 ! of the networks (my_arc), its own workspace column (selected by
176 ! omp_get_thread_num inside the ACSF routines) and its own slab of
177 ! my_force_t/my_cstress_t. SCHEDULE(STATIC) hands out contiguous atom
178 ! blocks in thread order, so the ordered sum over threads after the
179 ! region reproduces the serial summation order exactly.
180 !$OMP PARALLEL DEFAULT(NONE) &
181 !$OMP SHARED(calc_forces, calc_stress, istart, max_input_nodes, mecalc, &
182 !$OMP my_cstress_t, my_force_t, nnp) &
183 !$OMP PRIVATE(i, i_com, ind, ie, il, j, n_input_nodes, tid, my_arc, denergydsym, stress)
184
185 tid = 1
186!$ tid = omp_get_thread_num() + 1
187
188 ! Copy the networks for this thread. ALLOCATE with SOURCE= rather than
189 ! plain assignment: allocate-on-assignment to an ALLOCATABLE is what
190 ! -Werror=realloc-lhs rejects in the debug builds. SOURCE= copies the
191 ! shape, the bounds and the values, so the copy is structurally identical
192 ! to the original either way. %layer itself is a POINTER array and is
193 ! allocated and freed explicitly.
194 ALLOCATE (my_arc(nnp%n_ele))
195 DO ie = 1, nnp%n_ele
196 ALLOCATE (my_arc(ie)%n_nodes, source=nnp%arc(ie)%n_nodes)
197 ALLOCATE (my_arc(ie)%layer(nnp%n_layer))
198 DO il = 1, nnp%n_layer
199 IF (ALLOCATED(nnp%arc(ie)%layer(il)%weights)) THEN
200 ALLOCATE (my_arc(ie)%layer(il)%weights, source=nnp%arc(ie)%layer(il)%weights)
201 END IF
202 IF (ALLOCATED(nnp%arc(ie)%layer(il)%bweights)) THEN
203 ALLOCATE (my_arc(ie)%layer(il)%bweights, source=nnp%arc(ie)%layer(il)%bweights)
204 END IF
205 IF (ALLOCATED(nnp%arc(ie)%layer(il)%node)) THEN
206 ALLOCATE (my_arc(ie)%layer(il)%node, source=nnp%arc(ie)%layer(il)%node)
207 END IF
208 IF (ALLOCATED(nnp%arc(ie)%layer(il)%node_grad)) THEN
209 ALLOCATE (my_arc(ie)%layer(il)%node_grad, source=nnp%arc(ie)%layer(il)%node_grad)
210 END IF
211 IF (ALLOCATED(nnp%arc(ie)%layer(il)%tmp_der)) THEN
212 ALLOCATE (my_arc(ie)%layer(il)%tmp_der, source=nnp%arc(ie)%layer(il)%tmp_der)
213 END IF
214 END DO
215 END DO
216
217 IF (calc_forces) THEN
218 ALLOCATE (denergydsym(max_input_nodes))
219 denergydsym(:) = 0.0_dp
220 IF (calc_stress) THEN
221 ALLOCATE (stress(3, 3, max_input_nodes))
222 stress(:, :, :) = 0.0_dp
223 END IF
224 END IF
225
226 ! calc atomic contribution to energy and force
227 !$OMP DO SCHEDULE(STATIC)
228 DO i = istart, istart + mecalc - 1
229
230 ! determine index of atom type and offset
231 ind = nnp%ele_ind(i)
232 n_input_nodes = nnp%arc(ind)%n_nodes(1)
233
234 ! reset input nodes of ele(ind):
235 my_arc(ind)%layer(1)%node(:) = 0.0_dp
236
237 ! compute sym fnct values
238 IF (calc_forces) THEN
239 !reset input grads of ele(ind):
240 my_arc(ind)%layer(1)%node_grad(:) = 0.0_dp
241 IF (calc_stress) THEN
242 stress(:, :, 1:n_input_nodes) = 0.0_dp
243 CALL nnp_calc_acsf(nnp, i, .true., my_arc(ind), stress(:, :, 1:n_input_nodes))
244 ELSE
245 CALL nnp_calc_acsf(nnp, i, .true., my_arc(ind))
246 END IF
247 ELSE
248 CALL nnp_calc_acsf(nnp, i, .false., my_arc(ind))
249 END IF
250
251 DO i_com = 1, nnp%n_committee
252 ! predict energy
253 CALL nnp_predict(my_arc(ind), nnp, i_com)
254 nnp%atomic_energy(i, i_com) = my_arc(ind)%layer(nnp%n_layer)%node(1) + &
255 nnp%atom_energies(ind)
256
257 ! predict forces
258 IF (calc_forces) THEN
259 denergydsym(1:n_input_nodes) = 0.0_dp
260 CALL nnp_gradients(my_arc(ind), nnp, i_com, denergydsym(1:n_input_nodes))
261
262 ! Force scatter into this thread's private accumulator slab; see
263 ! nnp_scatter_dgdr_to_forces below, also reused by the helium-NNP
264 ! coupling in src/motion/helium_interactions.F.
265 CALL nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym(1:n_input_nodes), &
266 my_force_t(:, :, i_com, tid - 1))
267
268 IF (calc_stress) THEN
269 DO j = 1, n_input_nodes
270 my_cstress_t(:, :, i_com, tid - 1) = my_cstress_t(:, :, i_com, tid - 1) - &
271 denergydsym(j)*stress(:, :, j)
272 END DO
273 END IF
274 END IF
275 END DO
276
277 END DO ! loop over num_atoms
278 !$OMP END DO
279
280 IF (calc_forces) THEN
281 DEALLOCATE (denergydsym)
282 IF (calc_stress) DEALLOCATE (stress)
283 END IF
284 ! Free this thread's copy of the networks. %layer is a POINTER array, and
285 ! deallocating a pointer array does NOT release the ALLOCATABLE components
286 ! of its elements, so each one is released explicitly. Without this the
287 ! whole per-thread network copy leaks once per thread per MD step; the
288 ! LeakSanitizer in CP2K's pdbg CI build reports it. (%n_nodes needs no such
289 ! treatment: my_arc is ALLOCATABLE, and deallocating an allocatable does
290 ! cascade to its allocatable components.)
291 DO ie = 1, nnp%n_ele
292 DO il = 1, nnp%n_layer
293 IF (ALLOCATED(my_arc(ie)%layer(il)%weights)) DEALLOCATE (my_arc(ie)%layer(il)%weights)
294 IF (ALLOCATED(my_arc(ie)%layer(il)%bweights)) DEALLOCATE (my_arc(ie)%layer(il)%bweights)
295 IF (ALLOCATED(my_arc(ie)%layer(il)%node)) DEALLOCATE (my_arc(ie)%layer(il)%node)
296 IF (ALLOCATED(my_arc(ie)%layer(il)%node_grad)) DEALLOCATE (my_arc(ie)%layer(il)%node_grad)
297 IF (ALLOCATED(my_arc(ie)%layer(il)%tmp_der)) DEALLOCATE (my_arc(ie)%layer(il)%tmp_der)
298 END DO
299 DEALLOCATE (my_arc(ie)%layer)
300 END DO
301 DEALLOCATE (my_arc)
302 !$OMP END PARALLEL
303
304 ! Sum the per-thread slabs in ascending thread order. Together with the
305 ! static schedule above this keeps the summation order of the serial
306 ! loop, so forces are bit-identical to serial for any thread count.
307 IF (calc_forces) THEN
308 DO t = 0, nthreads - 1
309 nnp%myforce(:, :, :) = nnp%myforce(:, :, :) + my_force_t(:, :, :, t)
310 END DO
311 DEALLOCATE (my_force_t)
312 IF (calc_stress) THEN
313 DO t = 0, nthreads - 1
314 nnp%committee_stress(:, :, :) = nnp%committee_stress(:, :, :) + my_cstress_t(:, :, :, t)
315 END DO
316 DEALLOCATE (my_cstress_t)
317 END IF
318 END IF
319
320 CALL timestop(handle_loop)
321
322 ! calculate energy:
323 CALL logger%para_env%sum(nnp%atomic_energy(:, :))
324 nnp%committee_energy(:) = sum(nnp%atomic_energy, 1)
325 nnp%nnp_potential_energy = sum(nnp%committee_energy)/real(nnp%n_committee, dp)
326
327 IF (calc_forces) THEN
328 ! bring myforce to force array
329 DO j = 1, nnp%num_atoms
330 DO k = 1, 3
331 nnp%committee_forces(k, (nnp%sort(j)), :) = nnp%myforce(k, j, :)
332 END DO
333 END DO
334 CALL logger%para_env%sum(nnp%committee_forces)
335 nnp%nnp_forces(:, :) = sum(nnp%committee_forces, dim=3)/real(nnp%n_committee, dp)
336 DO j = 1, nnp%num_atoms
337 particle_set(j)%f(:) = nnp%nnp_forces(:, j)
338 END DO
339 END IF
340
341 IF (calc_stress) THEN
342 CALL logger%para_env%sum(nnp%committee_stress)
343 virial%pv_virial = sum(nnp%committee_stress, dim=3)/real(nnp%n_committee, dp)
344 END IF
345
346 ! Bias the standard deviation of committee disagreement
347 IF (nnp%bias) THEN
348 CALL nnp_bias_sigma(nnp, calc_forces)
349 nnp%nnp_potential_energy = nnp%nnp_potential_energy + nnp%bias_energy
350 IF (calc_forces) THEN
351 DO j = 1, nnp%num_atoms
352 particle_set(j)%f(:) = particle_set(j)%f(:) + nnp%bias_forces(:, j)
353 END DO
354 END IF
355 ! print properties if requested
356 print_section => section_vals_get_subs_vals(nnp%nnp_input, "BIAS%PRINT")
357 CALL nnp_print_bias(nnp, print_section)
358 END IF
359
360 ! print properties if requested
361 print_section => section_vals_get_subs_vals(nnp%nnp_input, "PRINT")
362 CALL nnp_print(nnp, print_section)
363
364 DEALLOCATE (allcalc)
365
366 CALL timestop(handle)
367
368 END SUBROUTINE nnp_calc_energy_force
369
370! **************************************************************************************************
371!> \brief Calculate bias potential and force based on standard deviation of committee disagreement
372!> \param nnp ...
373!> \param calc_forces ...
374!> \date 2020-10-10
375!> \author Christoph Schran (christoph.schran@rub.de)
376! **************************************************************************************************
377 SUBROUTINE nnp_bias_sigma(nnp, calc_forces)
378 TYPE(nnp_type), INTENT(INOUT) :: nnp
379 LOGICAL, INTENT(IN) :: calc_forces
380
381 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_bias_sigma'
382
383 INTEGER :: handle, i
384 REAL(kind=dp) :: avrg, pref, sigma
385
386 CALL timeset(routinen, handle)
387
388 ! init
389 sigma = 0.0_dp
390 nnp%bias_energy = 0.0_dp
391 IF (calc_forces) nnp%bias_forces = 0.0_dp
392
393 ! Subtract reference energy of each committee member, if requested
394 IF (nnp%bias_align) THEN
395 ! committee energy not used afterward, therefore overwritten
396 nnp%committee_energy(:) = nnp%committee_energy(:) - nnp%bias_e_avrg(:)
397 END IF
398
399 ! <E> = 1/n sum(E_i)
400 ! sigma = sqrt(1/n sum((E_i - <E>)**2))
401 ! = sqrt(1/n sum(dE_i**2))
402 avrg = sum(nnp%committee_energy)/real(nnp%n_committee, dp)
403 DO i = 1, nnp%n_committee
404 sigma = sigma + (nnp%committee_energy(i) - avrg)**2
405 END DO
406 sigma = sqrt(sigma/real(nnp%n_committee, dp))
407 nnp%bias_sigma = sigma
408
409 IF (sigma > nnp%bias_sigma0) THEN
410 ! E_b = 0.5 * kb * (sigma - sigma_0)**2
411 nnp%bias_energy = 0.5_dp*nnp%bias_kb*(sigma - nnp%bias_sigma0)**2
412
413 IF (calc_forces) THEN
414 ! nabla(E_b) = kb*(sigma - sigma_0)*nabla(sigma)
415 ! nabla(sigma) = 1/sigma * 1/n sum(dE_i* nabla(dE_i))
416 ! nabla(dE_i) = nabla(E_i) - nabla(<E>)
417 pref = nnp%bias_kb*(1.0_dp - nnp%bias_sigma0/sigma)
418 DO i = 1, nnp%n_committee
419 nnp%bias_forces(:, :) = nnp%bias_forces(:, :) + &
420 (nnp%committee_energy(i) - avrg)* &
421 (nnp%committee_forces(:, :, i) - nnp%nnp_forces(:, :))
422 END DO
423 pref = pref/real(nnp%n_committee, dp)
424 nnp%bias_forces(:, :) = nnp%bias_forces(:, :)*pref
425 END IF
426 END IF
427
428 CALL timestop(handle)
429
430 END SUBROUTINE nnp_bias_sigma
431
432! **************************************************************************************************
433!> \brief Print properties according to the requests in input file
434!> \param nnp ...
435!> \param print_section ...
436!> \date 2020-10-10
437!> \author Christoph Schran (christoph.schran@rub.de)
438! **************************************************************************************************
439 SUBROUTINE nnp_print(nnp, print_section)
440 TYPE(nnp_type), INTENT(INOUT) :: nnp
441 TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
442
443 INTEGER :: unit_nr
444 LOGICAL :: explicit, file_is_new
445 TYPE(cp_logger_type), POINTER :: logger
446 TYPE(section_vals_type), POINTER :: print_key
447
448 NULLIFY (logger, print_key)
449 logger => cp_get_default_logger()
450
451 print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
452 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
453 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
454 middle_name="nnp-energies", is_new_file=file_is_new)
455 IF (unit_nr > 0) CALL nnp_print_energies(nnp, unit_nr, file_is_new)
456 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
457 END IF
458
459 print_key => section_vals_get_subs_vals(print_section, "FORCES")
460 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
461 CALL nnp_print_forces(nnp, print_key)
462 END IF
463
464 print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
465 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
466 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
467 middle_name="nnp-forces-std")
468 IF (unit_nr > 0) CALL nnp_print_force_sigma(nnp, unit_nr)
469 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
470 END IF
471
472 ! Output structures with extrapolation warning on any processor
473 CALL logger%para_env%sum(nnp%output_expol)
474 IF (nnp%output_expol) THEN
475 print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
476 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
477 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
478 middle_name="nnp-extrapolation")
479 IF (unit_nr > 0) CALL nnp_print_expol(nnp, unit_nr)
480 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
481 END IF
482 END IF
483
484 print_key => section_vals_get_subs_vals(print_section, "SUM_FORCE")
485
486 CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
487 explicit=explicit)
488 IF (explicit) THEN
489 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
490 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".dat", &
491 middle_name="nnp-sumforce", is_new_file=file_is_new)
492 IF (unit_nr > 0) CALL nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
493 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
494 END IF
495 END IF
496
497 END SUBROUTINE nnp_print
498
499! **************************************************************************************************
500!> \brief Print NNP energies and standard deviation sigma
501!> \param nnp ...
502!> \param unit_nr ...
503!> \param file_is_new ...
504!> \date 2020-10-10
505!> \author Christoph Schran (christoph.schran@rub.de)
506! **************************************************************************************************
507 SUBROUTINE nnp_print_energies(nnp, unit_nr, file_is_new)
508 TYPE(nnp_type), INTENT(INOUT) :: nnp
509 INTEGER, INTENT(IN) :: unit_nr
510 LOGICAL, INTENT(IN) :: file_is_new
511
512 CHARACTER(LEN=12) :: fmt_string
513 INTEGER :: i
514 REAL(kind=dp) :: std
515
516 IF (file_is_new) THEN
517 WRITE (unit_nr, "(A1,1X,A20)", advance='no') "#", "NNP Average [a.u.],"
518 WRITE (unit_nr, "(A20)", advance='no') "NNP sigma [a.u.]"
519 DO i = 1, nnp%n_committee
520 WRITE (unit_nr, "(A17,I3)", advance='no') "NNP", i
521 END DO
522 WRITE (unit_nr, "(A)") ""
523 END IF
524
525 fmt_string = "(2X, F20.9)"
526 WRITE (unit=fmt_string(5:6), fmt="(I2)") nnp%n_committee + 2
527 std = sum((sum(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
528 std = std/real(nnp%n_committee, dp)
529 std = sqrt(std)
530 WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, sum(nnp%atomic_energy, 1)
531
532 END SUBROUTINE nnp_print_energies
533
534! **************************************************************************************************
535!> \brief Print nnp forces
536!> \param nnp ...
537!> \param print_key ...
538!> \date 2020-10-10
539!> \author Christoph Schran (christoph.schran@rub.de)
540! **************************************************************************************************
541 SUBROUTINE nnp_print_forces(nnp, print_key)
542 TYPE(nnp_type), INTENT(INOUT) :: nnp
543 TYPE(section_vals_type), INTENT(IN), POINTER :: print_key
544
545 CHARACTER(len=default_path_length) :: fmt_string, middle_name
546 INTEGER :: i, j, unit_nr
547 TYPE(cp_logger_type), POINTER :: logger
548
549 NULLIFY (logger)
550 logger => cp_get_default_logger()
551
552 DO i = 1, nnp%n_committee
553 WRITE (fmt_string, *) i
554 WRITE (middle_name, "(A,A)") "nnp-forces-", adjustl(trim(fmt_string))
555 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
556 middle_name=trim(middle_name))
557 IF (unit_nr > 0) THEN
558 WRITE (unit_nr, *) nnp%num_atoms
559 WRITE (unit_nr, "(A,1X,A,A,F20.9)") "NNP forces [a.u.] of committee member", &
560 adjustl(trim(fmt_string)), "energy [a.u.]=", nnp%committee_energy(i)
561
562 fmt_string = "(A4,1X,3F20.10)"
563 DO j = 1, nnp%num_atoms
564 WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(j)), nnp%committee_forces(:, j, i)
565 END DO
566 END IF
567 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
568 END DO
569
570 END SUBROUTINE nnp_print_forces
571
572! **************************************************************************************************
573!> \brief Print standard deviation sigma of NNP forces
574!> \param nnp ...
575!> \param unit_nr ...
576!> \date 2020-10-10
577!> \author Christoph Schran (christoph.schran@rub.de)
578! **************************************************************************************************
579 SUBROUTINE nnp_print_force_sigma(nnp, unit_nr)
580 TYPE(nnp_type), INTENT(INOUT) :: nnp
581 INTEGER, INTENT(IN) :: unit_nr
582
583 INTEGER :: i, j
584 REAL(kind=dp), DIMENSION(3) :: var
585
586 IF (unit_nr > 0) THEN
587 WRITE (unit_nr, *) nnp%num_atoms
588 WRITE (unit_nr, "(A,1X,A)") "NNP sigma of forces [a.u.]"
589
590 DO i = 1, nnp%num_atoms
591 var = 0.0_dp
592 DO j = 1, nnp%n_committee
593 var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
594 END DO
595 var = var/real(nnp%n_committee, dp)
596 var = sqrt(var)
597 WRITE (unit_nr, "(A4,1X,3F20.10)") nnp%atoms(nnp%sort_inv(i)), var
598 END DO
599 END IF
600
601 END SUBROUTINE nnp_print_force_sigma
602
603! **************************************************************************************************
604!> \brief Print structures with extrapolation warning
605!> \param nnp ...
606!> \param unit_nr ...
607!> \date 2020-10-10
608!> \author Christoph Schran (christoph.schran@rub.de)
609! **************************************************************************************************
610 SUBROUTINE nnp_print_expol(nnp, unit_nr)
611 TYPE(nnp_type), INTENT(INOUT) :: nnp
612 INTEGER, INTENT(IN) :: unit_nr
613
614 CHARACTER(len=default_path_length) :: fmt_string
615 INTEGER :: i
616 REAL(kind=dp) :: mass, unit_conv
617 REAL(kind=dp), DIMENSION(3) :: com
618
619 nnp%expol = nnp%expol + 1
620 WRITE (unit_nr, *) nnp%num_atoms
621 WRITE (unit_nr, "(A,1X,I6)") "NNP extrapolation point N =", nnp%expol
622
623 ! move to COM of solute and wrap the box
624 ! coord not needed afterwards, therefore manipulation ok
625 com = 0.0_dp
626 mass = 0.0_dp
627 DO i = 1, nnp%num_atoms
628 CALL get_ptable_info(nnp%atoms(i) (1:2), amass=unit_conv)
629 com(:) = com(:) + nnp%coord(:, i)*unit_conv
630 mass = mass + unit_conv
631 END DO
632 com(:) = com(:)/mass
633
634 DO i = 1, nnp%num_atoms
635 nnp%coord(:, i) = nnp%coord(:, i) - com(:)
636 END DO
637
638 ! write out coordinates
639 unit_conv = cp_unit_from_cp2k(1.0_dp, trim("angstrom"))
640 fmt_string = "(A4,1X,3F20.10)"
641 DO i = 1, nnp%num_atoms
642 WRITE (unit_nr, fmt_string) &
643 nnp%atoms(nnp%sort_inv(i)), &
644 nnp%coord(1, nnp%sort_inv(i))*unit_conv, &
645 nnp%coord(2, nnp%sort_inv(i))*unit_conv, &
646 nnp%coord(3, nnp%sort_inv(i))*unit_conv
647 END DO
648
649 END SUBROUTINE nnp_print_expol
650
651! **************************************************************************************************
652!> \brief Print properties number according the requests in input file
653!> \param nnp ...
654!> \param print_section ...
655!> \date 2020-10-10
656!> \author Christoph Schran (christoph.schran@rub.de)
657! **************************************************************************************************
658 SUBROUTINE nnp_print_bias(nnp, print_section)
659 TYPE(nnp_type), INTENT(INOUT) :: nnp
660 TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
661
662 INTEGER :: unit_nr
663 LOGICAL :: file_is_new
664 TYPE(cp_logger_type), POINTER :: logger
665 TYPE(section_vals_type), POINTER :: print_key
666
667 NULLIFY (logger, print_key)
668 logger => cp_get_default_logger()
669
670 print_key => section_vals_get_subs_vals(print_section, "BIAS_ENERGY")
671 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
672 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
673 middle_name="nnp-bias-energy", is_new_file=file_is_new)
674 IF (unit_nr > 0) CALL nnp_print_bias_energy(nnp, unit_nr, file_is_new)
675 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
676 END IF
677
678 print_key => section_vals_get_subs_vals(print_section, "BIAS_FORCES")
679 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
680 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
681 middle_name="nnp-bias-forces")
682 IF (unit_nr > 0) CALL nnp_print_bias_forces(nnp, unit_nr)
683 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
684 END IF
685
686 END SUBROUTINE nnp_print_bias
687
688! **************************************************************************************************
689!> \brief Print NNP bias energies
690!> \param nnp ...
691!> \param unit_nr ...
692!> \param file_is_new ...
693!> \date 2020-10-10
694!> \author Christoph Schran (christoph.schran@rub.de)
695! **************************************************************************************************
696 SUBROUTINE nnp_print_bias_energy(nnp, unit_nr, file_is_new)
697 TYPE(nnp_type), INTENT(INOUT) :: nnp
698 INTEGER, INTENT(IN) :: unit_nr
699 LOGICAL, INTENT(IN) :: file_is_new
700
701 CHARACTER(len=default_path_length) :: fmt_string
702 INTEGER :: i
703
704 IF (file_is_new) THEN
705 WRITE (unit_nr, "(A1)", advance='no') "#"
706 WRITE (unit_nr, "(2(2X,A19))", advance='no') "Sigma [a.u.]", "Bias energy [a.u.]"
707 DO i = 1, nnp%n_committee
708 IF (nnp%bias_align) THEN
709 WRITE (unit_nr, "(2X,A16,I3)", advance='no') "shifted E_NNP", i
710 ELSE
711 WRITE (unit_nr, "(2X,A16,I3)", advance='no') "E_NNP", i
712 END IF
713 END DO
714 WRITE (unit_nr, "(A)") ""
715
716 END IF
717
718 WRITE (fmt_string, "(A,I3,A)") "(2X,", nnp%n_committee + 2, "(F20.9,1X))"
719 WRITE (unit_nr, fmt_string) nnp%bias_sigma, nnp%bias_energy, nnp%committee_energy
720
721 END SUBROUTINE nnp_print_bias_energy
722
723! **************************************************************************************************
724!> \brief Print NNP bias forces
725!> \param nnp ...
726!> \param unit_nr ...
727!> \date 2020-10-10
728!> \author Christoph Schran (christoph.schran@rub.de)
729! **************************************************************************************************
730 SUBROUTINE nnp_print_bias_forces(nnp, unit_nr)
731 TYPE(nnp_type), INTENT(INOUT) :: nnp
732 INTEGER, INTENT(IN) :: unit_nr
733
734 CHARACTER(len=default_path_length) :: fmt_string
735 INTEGER :: i
736
737 WRITE (unit_nr, *) nnp%num_atoms
738 WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
739
740 fmt_string = "(A4,1X,3F20.10)"
741 DO i = 1, nnp%num_atoms
742 WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(i)), nnp%bias_forces(:, i)
743 END DO
744
745 END SUBROUTINE nnp_print_bias_forces
746
747! **************************************************************************************************
748!> \brief Print NNP summed forces
749!> \param nnp ...
750!> \param print_section ...
751!> \param unit_nr ...
752!> \param file_is_new ...
753!> \date 2020-10-10
754!> \author Christoph Schran (christoph.schran@rub.de)
755! **************************************************************************************************
756 SUBROUTINE nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
757 TYPE(nnp_type), INTENT(INOUT) :: nnp
758 TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
759 INTEGER, INTENT(IN) :: unit_nr
760 LOGICAL, INTENT(IN) :: file_is_new
761
762 CHARACTER(len=default_path_length) :: fmt_string
763 CHARACTER(LEN=default_string_length), &
764 DIMENSION(:), POINTER :: atomlist
765 INTEGER :: i, ig, j, n
766 REAL(kind=dp), DIMENSION(3) :: rvec
767
768 NULLIFY (atomlist)
769 IF (file_is_new) THEN
770 WRITE (unit_nr, "(A)") "# Summed forces [a.u.]"
771 END IF
772
773 rvec = 0.0_dp
774
775 ! get atoms to sum over:
776 CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
777 c_vals=atomlist)
778 IF (ASSOCIATED(atomlist)) THEN
779 n = SIZE(atomlist)
780 DO i = 1, nnp%num_atoms
781 DO j = 1, n
782 ig = nnp%sort_inv(i)
783 IF (trim(adjustl(atomlist(j))) == trim(adjustl(nnp%atoms(ig)))) THEN
784 rvec(:) = rvec(:) + nnp%nnp_forces(:, i)
785 END IF
786 END DO
787 END DO
788 END IF
789
790 fmt_string = "(3(F20.10,1X))"
791 WRITE (unit_nr, fmt_string) rvec
792
793 END SUBROUTINE nnp_print_sumforces
794
795! **************************************************************************************************
796!> \brief Scatter the per-neighbour dG/dr arrays held in the nnp_neighbor_workspace
797!> into a per-atom Cartesian force array. Three contributions per central
798!> atom i:
799!> 1) self: self_dGdr(:, m) -> force_xyz(:, i)
800!> 2) radial: walk neighbor%rad(s)%ind(j), read dGdr_rad(s)%data(:, sf, j)
801!> 3) angular: walk neighbor%ang1(s)%ind / ang2(s)%ind, read dGdr_ang_jj/kk
802!> No global (n_sf, num_atoms) slab is touched. Reused by the helium-NNP
803!> coupling in helium_interactions.
804!>
805!> Precondition: nnp_calc_acsf(nnp, i, calc_forces=.TRUE.[, stress]) must
806!> have run for this same atom i immediately before; the per-element
807!> workspace it fills is overwritten on every ACSF call.
808!> \param nnp ...
809!> \param ind central-atom element index
810!> \param i central atom index (absolute)
811!> \param denergydsym dE/dG_k for k = 1..n_input_nodes(ind)
812!> \param force_xyz (3, num_atoms) destination -- accumulated, NOT overwritten
813!> \author Dhruv Sharma (ds2173@cam.ac.uk)
814! **************************************************************************************************
815 SUBROUTINE nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, force_xyz)
816
817 TYPE(nnp_type), INTENT(IN) :: nnp
818 INTEGER, INTENT(IN) :: ind, i
819 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: denergydsym
820 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force_xyz
821
822 INTEGER :: j, k, k_atom, m, n_ang1_s, n_ang2_s, &
823 n_input_nodes, n_symf_s, off, s, sf, &
824 tid
825 LOGICAL :: homo_grp
826 REAL(kind=dp) :: de
827
828 tid = 1
829!$ tid = omp_get_thread_num() + 1
830 associate(workspace => nnp%neighbor_interface_state%workspace(ind, tid), &
831 neighbor => nnp%neighbor_interface_state%workspace(ind, tid)%neighbor, &
832 self_dgdr => nnp%neighbor_interface_state%workspace(ind, tid)%self_dGdr)
833
834 n_input_nodes = workspace%n_input_nodes
835
836 ! Self contribution
837 DO j = 1, n_input_nodes
838 de = denergydsym(j)
839 force_xyz(1, i) = force_xyz(1, i) - de*self_dgdr(1, j)
840 force_xyz(2, i) = force_xyz(2, i) - de*self_dgdr(2, j)
841 force_xyz(3, i) = force_xyz(3, i) - de*self_dgdr(3, j)
842 END DO
843
844 ! Radial neighbours
845 DO s = 1, nnp%rad(ind)%n_symfgrp
846 n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
847 associate(rad_buf => workspace%dGdr_rad(s)%data)
848 DO j = 1, neighbor%n_rad(s)
849 k_atom = neighbor%rad(s)%ind(j)
850 DO sf = 1, n_symf_s
851 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
852 de = denergydsym(m)
853 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*rad_buf(1, sf, j)
854 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*rad_buf(2, sf, j)
855 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*rad_buf(3, sf, j)
856 END DO
857 END DO
858 END associate
859 END DO
860
861 ! Angular neighbours
862 off = nnp%n_rad(ind)
863 DO s = 1, nnp%ang(ind)%n_symfgrp
864 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
865 n_ang1_s = neighbor%n_ang1(s)
866 homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
867 IF (homo_grp) THEN
868 n_ang2_s = n_ang1_s ! kk buffer is also indexed in ang1(s)%ind
869 ELSE
870 n_ang2_s = neighbor%n_ang2(s)
871 END IF
872
873 associate(jj_buf => workspace%dGdr_ang_jj(s)%data, &
874 kk_buf => workspace%dGdr_ang_kk(s)%data)
875 ! jj-side
876 DO j = 1, n_ang1_s
877 k_atom = neighbor%ang1(s)%ind(j)
878 DO sf = 1, n_symf_s
879 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
880 de = denergydsym(m)
881 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*jj_buf(1, sf, j)
882 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*jj_buf(2, sf, j)
883 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*jj_buf(3, sf, j)
884 END DO
885 END DO
886
887 ! kk-side: homo reads ang1(s)%ind, hetero reads ang2(s)%ind
888 IF (homo_grp) THEN
889 DO k = 1, n_ang2_s
890 k_atom = neighbor%ang1(s)%ind(k)
891 DO sf = 1, n_symf_s
892 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
893 de = denergydsym(m)
894 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
895 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
896 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
897 END DO
898 END DO
899 ELSE
900 DO k = 1, n_ang2_s
901 k_atom = neighbor%ang2(s)%ind(k)
902 DO sf = 1, n_symf_s
903 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
904 de = denergydsym(m)
905 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
906 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
907 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
908 END DO
909 END DO
910 END IF
911 END associate
912 END DO
913
914 END associate
915
916 END SUBROUTINE nnp_scatter_dgdr_to_forces
917
918END MODULE nnp_force
Define the atomic kind types and their sub types.
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1251
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
Functionality for atom centered symmetry functions for neural network potentials.
Definition nnp_acsf.F:16
subroutine, public nnp_calc_acsf(nnp, i, calc_forces, arc, stress)
Calculate atom centered symmetry functions for given atom i.
Definition nnp_acsf.F:88
subroutine, public nnp_prepare_neighbor_cache(nnp)
Prepare or update the linked-cell / Verlet cache for the current geometry. Lazily allocates nnpcell_l...
Definition nnp_acsf.F:689
Data types for neural network potentials.
subroutine, public nnp_env_get(nnp_env, nnp_forces, subsys, atomic_kind_set, particle_set, local_particles, molecule_kind_set, molecule_set, local_molecules, nnp_input, force_env_input, cell, cell_ref, use_ref_cell, nnp_potential_energy, virial)
Returns various attributes of the nnp environment.
Methods dealing with Neural Network potentials.
Definition nnp_force.F:15
subroutine, public nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, force_xyz)
Scatter the per-neighbour dG/dr arrays held in the nnp_neighbor_workspace into a per-atom Cartesian f...
Definition nnp_force.F:816
subroutine, public nnp_calc_energy_force(nnp, calc_forces)
Calculate the energy and force for a given configuration with the NNP.
Definition nnp_force.F:69
Methods dealing with core routines for artificial neural networks.
Definition nnp_model.F:13
subroutine, public nnp_predict(arc, nnp, i_com)
Predict energy by evaluating neural network.
Definition nnp_model.F:82
subroutine, public nnp_gradients(arc, nnp, i_com, denergydsym)
Calculate gradients of neural network.
Definition nnp_model.F:169
Define the data structure for the particle information.
Periodic Table related data definitions.
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
Provides all information about an atomic kind.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represents a system: atoms, molecules, their pos,vel,...
structure to store local (to a processor) ordered lists of integers.
Data type for artificial neural networks.
Main data type collecting all relevant data for neural network potentials.