(git:3f0d10b)
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! **************************************************************************************************
15MODULE nnp_force
16
17 USE atomic_kind_types, ONLY: atomic_kind_type
18 USE cp_log_handling, ONLY: cp_get_default_logger,&
19 cp_logger_get_default_io_unit,&
20 cp_logger_type
21 USE cp_output_handling, ONLY: cp_p_file,&
22 cp_print_key_finished_output,&
23 cp_print_key_should_output,&
24 cp_print_key_unit_nr
25 USE cp_subsys_types, ONLY: cp_subsys_get,&
26 cp_subsys_type
27 USE cp_units, ONLY: cp_unit_from_cp2k
28 USE distribution_1d_types, ONLY: distribution_1d_type
29 USE input_constants, ONLY: nnp_gpu_prec_mixed,&
30 nnp_use_gpu_on
31 USE input_section_types, ONLY: section_vals_get_subs_vals,&
32 section_vals_type,&
33 section_vals_val_get
34 USE kinds, ONLY: default_path_length,&
35 default_string_length,&
36 dp
37 USE nnp_acsf, ONLY: nnp_calc_acsf,&
38 nnp_prepare_neighbor_cache
39 USE nnp_environment_types, ONLY: nnp_arc_type,&
40 nnp_env_get,&
41 nnp_type
42 USE nnp_gpu, ONLY: nnp_gpu_energy_force
43 USE nnp_model, ONLY: nnp_gradients,&
44 nnp_predict
45 USE particle_types, ONLY: particle_type
46 USE periodic_table, ONLY: get_ptable_info
47 USE physcon, ONLY: angstrom
48 USE virial_types, ONLY: virial_type
49
50!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
51#include "./base/base_uses.f90"
52
53 IMPLICIT NONE
54
55 PRIVATE
56
57 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
58 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_force'
59
60 PUBLIC :: nnp_calc_energy_force, nnp_scatter_dgdr_to_forces
61
62CONTAINS
63
64! **************************************************************************************************
65!> \brief Calculate the energy and force for a given configuration with the NNP
66!> \param nnp ...
67!> \param calc_forces ...
68!> \date 2020-10-10
69!> \author Christoph Schran (christoph.schran@rub.de)
70!> \author Claudio Malvino (claudiormal@gmail.com)
71! **************************************************************************************************
72 SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
73 TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
74 LOGICAL, INTENT(IN) :: calc_forces
75
76 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_calc_energy_force'
77
78 INTEGER :: handle, handle_loop, i, i_com, ie, ig, &
79 il, ind, istart, j, k, m, &
80 max_input_nodes, mecalc, &
81 n_input_nodes, nthreads, t, tid
82 INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
83 LOGICAL :: calc_stress, run_cpu_path, use_gpu_path, &
84 validate
85 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
86 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: val_energy
87 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: stress, val_force, val_stress
88 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: my_cstress_t, my_force_t
89 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
90 TYPE(cp_logger_type), POINTER :: logger
91 TYPE(cp_subsys_type), POINTER :: subsys
92 TYPE(distribution_1d_type), POINTER :: local_particles
93 TYPE(nnp_arc_type), ALLOCATABLE, DIMENSION(:) :: my_arc
94 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
95 TYPE(section_vals_type), POINTER :: print_section
96 TYPE(virial_type), POINTER :: virial
97
98 CALL timeset(routinen, handle)
99
100 NULLIFY (particle_set, logger, local_particles, subsys, &
101 atomic_kind_set)
102 logger => cp_get_default_logger()
103
104 cpassert(ASSOCIATED(nnp))
105 cpassert(nnp%n_committee >= 1)
106 CALL nnp_env_get(nnp_env=nnp, particle_set=particle_set, &
107 subsys=subsys, local_particles=local_particles, &
108 atomic_kind_set=atomic_kind_set)
109
110 CALL cp_subsys_get(subsys, &
111 virial=virial)
112
113 calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
114 IF (calc_stress .AND. .NOT. calc_forces) THEN
115 cpabort('Stress cannot be calculated without forces')
116 END IF
117
118 nnp%atomic_energy(:, :) = 0.0_dp
119 IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
120 IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
121 IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
122
123 !fill coord array
124 ig = 1
125 DO i = 1, nnp%n_ele
126 DO j = 1, nnp%num_atoms
127 IF (nnp%ele(i) == particle_set(j)%atomic_kind%element_symbol) THEN
128 DO m = 1, 3
129 nnp%coord(m, ig) = particle_set(j)%r(m)
130 END DO
131 nnp%atoms(ig) = nnp%ele(i)
132 CALL get_ptable_info(nnp%atoms(ig) (1:2), number=nnp%nuc_atoms(ig))
133 nnp%ele_ind(ig) = i
134 nnp%sort(ig) = j
135 nnp%sort_inv(j) = ig
136 ig = ig + 1
137 END IF
138 END DO
139 END DO
140
141 ! parallization:
142 mecalc = nnp%num_atoms/logger%para_env%num_pe + &
143 min(mod(nnp%num_atoms, logger%para_env%num_pe)/ &
144 (logger%para_env%mepos + 1), 1)
145 ALLOCATE (allcalc(logger%para_env%num_pe))
146 allcalc(:) = 0
147 CALL logger%para_env%allgather(mecalc, allcalc)
148 istart = 1
149 DO i = 2, logger%para_env%mepos + 1
150 istart = istart + allcalc(i - 1)
151 END DO
152
153 ! reset extrapolation status
154 nnp%output_expol = .false.
155
156 CALL nnp_prepare_neighbor_cache(nnp)
157
158 max_input_nodes = 0
159 DO i = 1, nnp%n_ele
160 max_input_nodes = max(max_input_nodes, nnp%arc(i)%n_nodes(1))
161 END DO
162
163 ! The GPU path runs the whole per-rank step (neighbour walk, descriptors,
164 ! network and force assembly) on the device. It handles energies, forces
165 ! and, when a stress tensor is requested, the analytical virial the
166 ! scatter kernels accumulate alongside the forces. Energy-only evaluations
167 ! always use the CPU loop.
168 use_gpu_path = (nnp%use_gpu == nnp_use_gpu_on) .AND. calc_forces
169 ! VALIDATE runs the CPU reference in shadow mode against the device
170 ! results, so both paths execute on a validating evaluation.
171 validate = use_gpu_path .AND. nnp%validate_gpu
172 run_cpu_path = (.NOT. use_gpu_path) .OR. validate
173
174 ! number of threads that will run the atom loop below. The per-thread
175 ! accumulator slabs belong to the CPU loop: the device kernels scatter
176 ! straight into nnp%myforce and nnp%committee_stress.
177 nthreads = 1
178!$ nthreads = omp_get_max_threads()
179 IF (calc_forces .AND. run_cpu_path) THEN
180 ! dGdr lives in per-element workspace; no global (3, n_sf, num_atoms) slab.
181 ALLOCATE (my_force_t(3, nnp%num_atoms, nnp%n_committee, 0:nthreads - 1))
182 my_force_t(:, :, :, :) = 0.0_dp
183 IF (calc_stress) THEN
184 ALLOCATE (my_cstress_t(3, 3, nnp%n_committee, 0:nthreads - 1))
185 my_cstress_t(:, :, :, :) = 0.0_dp
186 END IF
187 END IF
188
189 ! Per-rank timer around the per-atom work. CP2K aggregates the label as
190 ! avg/max across ranks, so max - avg measures the load imbalance.
191 CALL timeset('nnp_per_rank_atom_loop', handle_loop)
192
193 IF (use_gpu_path) THEN
194 CALL nnp_gpu_energy_force(nnp, istart, mecalc, max_input_nodes, calc_stress)
195 IF (validate) THEN
196 ! Snapshot the device results and clear the accumulators, so the
197 ! CPU reference below rebuilds them from scratch for the comparison.
198 ALLOCATE (val_energy, source=nnp%atomic_energy)
199 ALLOCATE (val_force, source=nnp%myforce)
200 IF (calc_stress) ALLOCATE (val_stress, source=nnp%committee_stress)
201 nnp%atomic_energy(:, :) = 0.0_dp
202 nnp%myforce(:, :, :) = 0.0_dp
203 IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
204 END IF
205 END IF
206 IF (run_cpu_path) THEN
207
208 ! Thread the loop over atomic centres. Each thread works on its own copy
209 ! of the networks (my_arc), its own workspace column (selected by
210 ! omp_get_thread_num inside the ACSF routines) and its own slab of
211 ! my_force_t/my_cstress_t. SCHEDULE(STATIC) hands out contiguous atom
212 ! blocks in thread order, so the ordered sum over threads after the
213 ! region reproduces the serial summation order exactly.
214 !$OMP PARALLEL DEFAULT(NONE) &
215 !$OMP SHARED(calc_forces, calc_stress, istart, max_input_nodes, mecalc, &
216 !$OMP my_cstress_t, my_force_t, nnp) &
217 !$OMP PRIVATE(i, i_com, ind, ie, il, j, n_input_nodes, tid, my_arc, denergydsym, stress)
218
219 tid = 1
220!$ tid = omp_get_thread_num() + 1
221
222 ! Copy the networks for this thread. ALLOCATE with SOURCE= rather than
223 ! plain assignment: allocate-on-assignment to an ALLOCATABLE is what
224 ! -Werror=realloc-lhs rejects in the debug builds. SOURCE= copies the
225 ! shape, the bounds and the values, so the copy is structurally identical
226 ! to the original either way. %layer itself is a POINTER array and is
227 ! allocated and freed explicitly.
228 ALLOCATE (my_arc(nnp%n_ele))
229 DO ie = 1, nnp%n_ele
230 ALLOCATE (my_arc(ie)%n_nodes, source=nnp%arc(ie)%n_nodes)
231 ALLOCATE (my_arc(ie)%layer(nnp%n_layer))
232 DO il = 1, nnp%n_layer
233 IF (ALLOCATED(nnp%arc(ie)%layer(il)%weights)) THEN
234 ALLOCATE (my_arc(ie)%layer(il)%weights, source=nnp%arc(ie)%layer(il)%weights)
235 END IF
236 IF (ALLOCATED(nnp%arc(ie)%layer(il)%bweights)) THEN
237 ALLOCATE (my_arc(ie)%layer(il)%bweights, source=nnp%arc(ie)%layer(il)%bweights)
238 END IF
239 IF (ALLOCATED(nnp%arc(ie)%layer(il)%node)) THEN
240 ALLOCATE (my_arc(ie)%layer(il)%node, source=nnp%arc(ie)%layer(il)%node)
241 END IF
242 IF (ALLOCATED(nnp%arc(ie)%layer(il)%node_grad)) THEN
243 ALLOCATE (my_arc(ie)%layer(il)%node_grad, source=nnp%arc(ie)%layer(il)%node_grad)
244 END IF
245 IF (ALLOCATED(nnp%arc(ie)%layer(il)%tmp_der)) THEN
246 ALLOCATE (my_arc(ie)%layer(il)%tmp_der, source=nnp%arc(ie)%layer(il)%tmp_der)
247 END IF
248 END DO
249 END DO
250
251 IF (calc_forces) THEN
252 ALLOCATE (denergydsym(max_input_nodes))
253 denergydsym(:) = 0.0_dp
254 IF (calc_stress) THEN
255 ALLOCATE (stress(3, 3, max_input_nodes))
256 stress(:, :, :) = 0.0_dp
257 END IF
258 END IF
259
260 ! calc atomic contribution to energy and force
261 !$OMP DO SCHEDULE(STATIC)
262 DO i = istart, istart + mecalc - 1
263
264 ! determine index of atom type and offset
265 ind = nnp%ele_ind(i)
266 n_input_nodes = nnp%arc(ind)%n_nodes(1)
267
268 ! reset input nodes of ele(ind):
269 my_arc(ind)%layer(1)%node(:) = 0.0_dp
270
271 ! compute sym fnct values
272 IF (calc_forces) THEN
273 !reset input grads of ele(ind):
274 my_arc(ind)%layer(1)%node_grad(:) = 0.0_dp
275 IF (calc_stress) THEN
276 stress(:, :, 1:n_input_nodes) = 0.0_dp
277 CALL nnp_calc_acsf(nnp, i, .true., my_arc(ind), stress(:, :, 1:n_input_nodes))
278 ELSE
279 CALL nnp_calc_acsf(nnp, i, .true., my_arc(ind))
280 END IF
281 ELSE
282 CALL nnp_calc_acsf(nnp, i, .false., my_arc(ind))
283 END IF
284
285 DO i_com = 1, nnp%n_committee
286 ! predict energy
287 CALL nnp_predict(my_arc(ind), nnp, i_com)
288 nnp%atomic_energy(i, i_com) = my_arc(ind)%layer(nnp%n_layer)%node(1) + &
289 nnp%atom_energies(ind)
290
291 ! predict forces
292 IF (calc_forces) THEN
293 denergydsym(1:n_input_nodes) = 0.0_dp
294 CALL nnp_gradients(my_arc(ind), nnp, i_com, denergydsym(1:n_input_nodes))
295
296 ! Force scatter into this thread's private accumulator slab; see
297 ! nnp_scatter_dgdr_to_forces below, also reused by the helium-NNP
298 ! coupling in src/motion/helium_interactions.F.
299 CALL nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym(1:n_input_nodes), &
300 my_force_t(:, :, i_com, tid - 1))
301
302 IF (calc_stress) THEN
303 DO j = 1, n_input_nodes
304 my_cstress_t(:, :, i_com, tid - 1) = my_cstress_t(:, :, i_com, tid - 1) - &
305 denergydsym(j)*stress(:, :, j)
306 END DO
307 END IF
308 END IF
309 END DO
310
311 END DO ! loop over num_atoms
312 !$OMP END DO
313
314 IF (calc_forces) THEN
315 DEALLOCATE (denergydsym)
316 IF (calc_stress) DEALLOCATE (stress)
317 END IF
318 ! Free this thread's copy of the networks. %layer is a POINTER array, and
319 ! deallocating a pointer array does NOT release the ALLOCATABLE components
320 ! of its elements, so each one is released explicitly. Without this the
321 ! whole per-thread network copy leaks once per thread per MD step; the
322 ! LeakSanitizer in CP2K's pdbg CI build reports it. (%n_nodes needs no such
323 ! treatment: my_arc is ALLOCATABLE, and deallocating an allocatable does
324 ! cascade to its allocatable components.)
325 DO ie = 1, nnp%n_ele
326 DO il = 1, nnp%n_layer
327 IF (ALLOCATED(my_arc(ie)%layer(il)%weights)) DEALLOCATE (my_arc(ie)%layer(il)%weights)
328 IF (ALLOCATED(my_arc(ie)%layer(il)%bweights)) DEALLOCATE (my_arc(ie)%layer(il)%bweights)
329 IF (ALLOCATED(my_arc(ie)%layer(il)%node)) DEALLOCATE (my_arc(ie)%layer(il)%node)
330 IF (ALLOCATED(my_arc(ie)%layer(il)%node_grad)) DEALLOCATE (my_arc(ie)%layer(il)%node_grad)
331 IF (ALLOCATED(my_arc(ie)%layer(il)%tmp_der)) DEALLOCATE (my_arc(ie)%layer(il)%tmp_der)
332 END DO
333 DEALLOCATE (my_arc(ie)%layer)
334 END DO
335 DEALLOCATE (my_arc)
336 !$OMP END PARALLEL
337
338 ! Sum the per-thread slabs in ascending thread order. Together with the
339 ! static schedule above this keeps the summation order of the serial
340 ! loop, so forces are bit-identical to serial for any thread count.
341 IF (calc_forces) THEN
342 DO t = 0, nthreads - 1
343 nnp%myforce(:, :, :) = nnp%myforce(:, :, :) + my_force_t(:, :, :, t)
344 END DO
345 DEALLOCATE (my_force_t)
346 IF (calc_stress) THEN
347 DO t = 0, nthreads - 1
348 nnp%committee_stress(:, :, :) = nnp%committee_stress(:, :, :) + my_cstress_t(:, :, :, t)
349 END DO
350 DEALLOCATE (my_cstress_t)
351 END IF
352 END IF
353
354 END IF
355
356 CALL timestop(handle_loop)
357
358 ! Compare the CPU reference built above against the device snapshot, then
359 ! continue with the device results.
360 IF (validate) THEN
361 IF (calc_stress) THEN
362 CALL nnp_gpu_validate(nnp, logger, val_energy, val_force, val_stress)
363 ELSE
364 CALL nnp_gpu_validate(nnp, logger, val_energy, val_force)
365 END IF
366 END IF
367
368 ! calculate energy:
369 CALL logger%para_env%sum(nnp%atomic_energy(:, :))
370 nnp%committee_energy(:) = sum(nnp%atomic_energy, 1)
371 nnp%nnp_potential_energy = sum(nnp%committee_energy)/real(nnp%n_committee, dp)
372
373 IF (calc_forces) THEN
374 ! bring myforce to force array
375 DO j = 1, nnp%num_atoms
376 DO k = 1, 3
377 nnp%committee_forces(k, (nnp%sort(j)), :) = nnp%myforce(k, j, :)
378 END DO
379 END DO
380 CALL logger%para_env%sum(nnp%committee_forces)
381 nnp%nnp_forces(:, :) = sum(nnp%committee_forces, dim=3)/real(nnp%n_committee, dp)
382 DO j = 1, nnp%num_atoms
383 particle_set(j)%f(:) = nnp%nnp_forces(:, j)
384 END DO
385 END IF
386
387 IF (calc_stress) THEN
388 CALL logger%para_env%sum(nnp%committee_stress)
389 virial%pv_virial = sum(nnp%committee_stress, dim=3)/real(nnp%n_committee, dp)
390 END IF
391
392 ! Bias the standard deviation of committee disagreement
393 IF (nnp%bias) THEN
394 CALL nnp_bias_sigma(nnp, calc_forces)
395 nnp%nnp_potential_energy = nnp%nnp_potential_energy + nnp%bias_energy
396 IF (calc_forces) THEN
397 DO j = 1, nnp%num_atoms
398 particle_set(j)%f(:) = particle_set(j)%f(:) + nnp%bias_forces(:, j)
399 END DO
400 END IF
401 ! print properties if requested
402 print_section => section_vals_get_subs_vals(nnp%nnp_input, "BIAS%PRINT")
403 CALL nnp_print_bias(nnp, print_section)
404 END IF
405
406 ! print properties if requested
407 print_section => section_vals_get_subs_vals(nnp%nnp_input, "PRINT")
408 CALL nnp_print(nnp, print_section)
409
410 DEALLOCATE (allcalc)
411
412 CALL timestop(handle)
413
414 END SUBROUTINE nnp_calc_energy_force
415
416! **************************************************************************************************
417!> \brief Compare the CPU reference results against the device results, keep
418!> the device results, and abort when the two disagree.
419!>
420!> Runs on the rank-local accumulators, before any parallel reduction: the
421!> per-atom committee energies, the per-rank force contributions, and, when a
422!> stress tensor is requested, the per-rank virial contributions. The metric is
423!> the relative difference floored at one, so small values are compared
424!> absolutely. The atomic additions of the device descriptor accumulation give
425!> its results a scheduling-dependent summation order, which bounds the
426!> agreement of the all-FP64 path; the FP32 angular symmetry functions of the
427!> MIXED path perturb the results much further, so its tolerances are looser.
428!>
429!> \param nnp NNP environment holding the CPU results; overwritten with
430!> the device results on success.
431!> \param logger Default logger, for its parallel environment.
432!> \param gpu_energy Device per-atom committee energies.
433!> \param gpu_force Device per-rank force contributions.
434!> \param gpu_stress Device per-rank virial contributions; present exactly when
435!> the evaluation computed a stress tensor.
436!> \author Dhruv Sharma (ds2173@cam.ac.uk)
437! **************************************************************************************************
438 SUBROUTINE nnp_gpu_validate(nnp, logger, gpu_energy, gpu_force, gpu_stress)
439 TYPE(nnp_type), INTENT(INOUT) :: nnp
440 TYPE(cp_logger_type), POINTER :: logger
441 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: gpu_energy
442 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: gpu_force
443 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
444 OPTIONAL :: gpu_stress
445
446 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_gpu_validate'
447
448 CHARACTER(LEN=2*default_string_length) :: message
449 INTEGER :: handle, unit_nr
450 REAL(kind=dp) :: diff_energy, diff_force, diff_stress, &
451 tol_energy, tol_force
452
453 CALL timeset(routinen, handle)
454
455 IF (nnp%gpu_precision == nnp_gpu_prec_mixed) THEN
456 tol_energy = 1.0e-5_dp
457 tol_force = 1.0e-4_dp
458 ELSE
459 tol_energy = 1.0e-9_dp
460 tol_force = 1.0e-8_dp
461 END IF
462
463 diff_energy = maxval(abs(nnp%atomic_energy - gpu_energy)/ &
464 max(1.0_dp, abs(nnp%atomic_energy)))
465 diff_force = maxval(abs(nnp%myforce - gpu_force)/ &
466 max(1.0_dp, abs(nnp%myforce)))
467 diff_stress = 0.0_dp
468 IF (PRESENT(gpu_stress)) THEN
469 diff_stress = maxval(abs(nnp%committee_stress - gpu_stress)/ &
470 max(1.0_dp, abs(nnp%committee_stress)))
471 END IF
472
473 ! Reduce before deciding, so every rank reaches the same verdict.
474 CALL logger%para_env%max(diff_energy)
475 CALL logger%para_env%max(diff_force)
476 CALL logger%para_env%max(diff_stress)
477
478 IF (diff_energy > tol_energy .OR. diff_force > tol_force .OR. &
479 diff_stress > tol_force) THEN
480 WRITE (message, "(A,3ES11.3,A,2ES10.1)") &
481 "NNP GPU validation failure. Max rel. diff (energies, forces, "// &
482 "stress):", diff_energy, diff_force, diff_stress, &
483 ", tolerances:", tol_energy, tol_force
484 cpabort(trim(message))
485 END IF
486
487 unit_nr = cp_logger_get_default_io_unit(logger)
488 IF (unit_nr > 0) THEN
489 WRITE (unit_nr, "(T2,A,3ES11.3)") &
490 "NNP| Validated GPU force evaluation, max rel. diff: ", &
491 diff_energy, diff_force, diff_stress
492 END IF
493
494 ! Keep the device results as the results of the run, matching an
495 ! evaluation without validation.
496 nnp%atomic_energy(:, :) = gpu_energy(:, :)
497 nnp%myforce(:, :, :) = gpu_force(:, :, :)
498 IF (PRESENT(gpu_stress)) nnp%committee_stress(:, :, :) = gpu_stress(:, :, :)
499
500 CALL timestop(handle)
501
502 END SUBROUTINE nnp_gpu_validate
503
504! **************************************************************************************************
505!> \brief Calculate bias potential and force based on standard deviation of committee disagreement
506!> \param nnp ...
507!> \param calc_forces ...
508!> \date 2020-10-10
509!> \author Christoph Schran (christoph.schran@rub.de)
510! **************************************************************************************************
511 SUBROUTINE nnp_bias_sigma(nnp, calc_forces)
512 TYPE(nnp_type), INTENT(INOUT) :: nnp
513 LOGICAL, INTENT(IN) :: calc_forces
514
515 CHARACTER(len=*), PARAMETER :: routinen = 'nnp_bias_sigma'
516
517 INTEGER :: handle, i
518 REAL(kind=dp) :: avrg, pref, sigma
519
520 CALL timeset(routinen, handle)
521
522 ! init
523 sigma = 0.0_dp
524 nnp%bias_energy = 0.0_dp
525 IF (calc_forces) nnp%bias_forces = 0.0_dp
526
527 ! Subtract reference energy of each committee member, if requested
528 IF (nnp%bias_align) THEN
529 ! committee energy not used afterward, therefore overwritten
530 nnp%committee_energy(:) = nnp%committee_energy(:) - nnp%bias_e_avrg(:)
531 END IF
532
533 ! <E> = 1/n sum(E_i)
534 ! sigma = sqrt(1/n sum((E_i - <E>)**2))
535 ! = sqrt(1/n sum(dE_i**2))
536 avrg = sum(nnp%committee_energy)/real(nnp%n_committee, dp)
537 DO i = 1, nnp%n_committee
538 sigma = sigma + (nnp%committee_energy(i) - avrg)**2
539 END DO
540 sigma = sqrt(sigma/real(nnp%n_committee, dp))
541 nnp%bias_sigma = sigma
542
543 IF (sigma > nnp%bias_sigma0) THEN
544 ! E_b = 0.5 * kb * (sigma - sigma_0)**2
545 nnp%bias_energy = 0.5_dp*nnp%bias_kb*(sigma - nnp%bias_sigma0)**2
546
547 IF (calc_forces) THEN
548 ! nabla(E_b) = kb*(sigma - sigma_0)*nabla(sigma)
549 ! nabla(sigma) = 1/sigma * 1/n sum(dE_i* nabla(dE_i))
550 ! nabla(dE_i) = nabla(E_i) - nabla(<E>)
551 pref = nnp%bias_kb*(1.0_dp - nnp%bias_sigma0/sigma)
552 DO i = 1, nnp%n_committee
553 nnp%bias_forces(:, :) = nnp%bias_forces(:, :) + &
554 (nnp%committee_energy(i) - avrg)* &
555 (nnp%committee_forces(:, :, i) - nnp%nnp_forces(:, :))
556 END DO
557 pref = pref/real(nnp%n_committee, dp)
558 nnp%bias_forces(:, :) = nnp%bias_forces(:, :)*pref
559 END IF
560 END IF
561
562 CALL timestop(handle)
563
564 END SUBROUTINE nnp_bias_sigma
565
566! **************************************************************************************************
567!> \brief Print properties according to the requests in input file
568!> \param nnp ...
569!> \param print_section ...
570!> \date 2020-10-10
571!> \author Christoph Schran (christoph.schran@rub.de)
572! **************************************************************************************************
573 SUBROUTINE nnp_print(nnp, print_section)
574 TYPE(nnp_type), INTENT(INOUT) :: nnp
575 TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
576
577 INTEGER :: unit_nr
578 LOGICAL :: explicit, file_is_new
579 TYPE(cp_logger_type), POINTER :: logger
580 TYPE(section_vals_type), POINTER :: print_key
581
582 NULLIFY (logger, print_key)
583 logger => cp_get_default_logger()
584
585 print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
586 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
587 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
588 middle_name="nnp-energies", is_new_file=file_is_new)
589 IF (unit_nr > 0) CALL nnp_print_energies(nnp, unit_nr, file_is_new)
590 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
591 END IF
592
593 print_key => section_vals_get_subs_vals(print_section, "FORCES")
594 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
595 CALL nnp_print_forces(nnp, print_key)
596 END IF
597
598 print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
599 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
600 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
601 middle_name="nnp-forces-std")
602 IF (unit_nr > 0) CALL nnp_print_force_sigma(nnp, unit_nr)
603 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
604 END IF
605
606 ! Output structures with extrapolation warning on any processor
607 CALL logger%para_env%sum(nnp%output_expol)
608 IF (nnp%output_expol) THEN
609 print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
610 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
611 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
612 middle_name="nnp-extrapolation")
613 IF (unit_nr > 0) CALL nnp_print_expol(nnp, unit_nr)
614 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
615 END IF
616 END IF
617
618 print_key => section_vals_get_subs_vals(print_section, "SUM_FORCE")
619
620 CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
621 explicit=explicit)
622 IF (explicit) THEN
623 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
624 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".dat", &
625 middle_name="nnp-sumforce", is_new_file=file_is_new)
626 IF (unit_nr > 0) CALL nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
627 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
628 END IF
629 END IF
630
631 END SUBROUTINE nnp_print
632
633! **************************************************************************************************
634!> \brief Print NNP energies and standard deviation sigma
635!> \param nnp ...
636!> \param unit_nr ...
637!> \param file_is_new ...
638!> \date 2020-10-10
639!> \author Christoph Schran (christoph.schran@rub.de)
640! **************************************************************************************************
641 SUBROUTINE nnp_print_energies(nnp, unit_nr, file_is_new)
642 TYPE(nnp_type), INTENT(INOUT) :: nnp
643 INTEGER, INTENT(IN) :: unit_nr
644 LOGICAL, INTENT(IN) :: file_is_new
645
646 CHARACTER(LEN=12) :: fmt_string
647 INTEGER :: i
648 REAL(kind=dp) :: std
649
650 IF (file_is_new) THEN
651 WRITE (unit_nr, "(A1,1X,A20)", advance='no') "#", "NNP Average [a.u.],"
652 WRITE (unit_nr, "(A20)", advance='no') "NNP sigma [a.u.]"
653 DO i = 1, nnp%n_committee
654 WRITE (unit_nr, "(A17,I3)", advance='no') "NNP", i
655 END DO
656 WRITE (unit_nr, "(A)") ""
657 END IF
658
659 fmt_string = "(2X, F20.9)"
660 WRITE (unit=fmt_string(5:6), fmt="(I2)") nnp%n_committee + 2
661 std = sum((sum(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
662 std = std/real(nnp%n_committee, dp)
663 std = sqrt(std)
664 WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, sum(nnp%atomic_energy, 1)
665
666 END SUBROUTINE nnp_print_energies
667
668! **************************************************************************************************
669!> \brief Print nnp forces
670!> \param nnp ...
671!> \param print_key ...
672!> \date 2020-10-10
673!> \author Christoph Schran (christoph.schran@rub.de)
674! **************************************************************************************************
675 SUBROUTINE nnp_print_forces(nnp, print_key)
676 TYPE(nnp_type), INTENT(INOUT) :: nnp
677 TYPE(section_vals_type), INTENT(IN), POINTER :: print_key
678
679 CHARACTER(len=default_path_length) :: fmt_string, middle_name
680 INTEGER :: i, j, unit_nr
681 TYPE(cp_logger_type), POINTER :: logger
682
683 NULLIFY (logger)
684 logger => cp_get_default_logger()
685
686 DO i = 1, nnp%n_committee
687 WRITE (fmt_string, *) i
688 WRITE (middle_name, "(A,A)") "nnp-forces-", adjustl(trim(fmt_string))
689 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
690 middle_name=trim(middle_name))
691 IF (unit_nr > 0) THEN
692 WRITE (unit_nr, *) nnp%num_atoms
693 WRITE (unit_nr, "(A,1X,A,A,F20.9)") "NNP forces [a.u.] of committee member", &
694 adjustl(trim(fmt_string)), "energy [a.u.]=", nnp%committee_energy(i)
695
696 fmt_string = "(A4,1X,3F20.10)"
697 DO j = 1, nnp%num_atoms
698 WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(j)), nnp%committee_forces(:, j, i)
699 END DO
700 END IF
701 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
702 END DO
703
704 END SUBROUTINE nnp_print_forces
705
706! **************************************************************************************************
707!> \brief Print standard deviation sigma of NNP forces
708!> \param nnp ...
709!> \param unit_nr ...
710!> \date 2020-10-10
711!> \author Christoph Schran (christoph.schran@rub.de)
712! **************************************************************************************************
713 SUBROUTINE nnp_print_force_sigma(nnp, unit_nr)
714 TYPE(nnp_type), INTENT(INOUT) :: nnp
715 INTEGER, INTENT(IN) :: unit_nr
716
717 INTEGER :: i, j
718 REAL(kind=dp), DIMENSION(3) :: var
719
720 IF (unit_nr > 0) THEN
721 WRITE (unit_nr, *) nnp%num_atoms
722 WRITE (unit_nr, "(A,1X,A)") "NNP sigma of forces [a.u.]"
723
724 DO i = 1, nnp%num_atoms
725 var = 0.0_dp
726 DO j = 1, nnp%n_committee
727 var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
728 END DO
729 var = var/real(nnp%n_committee, dp)
730 var = sqrt(var)
731 WRITE (unit_nr, "(A4,1X,3F20.10)") nnp%atoms(nnp%sort_inv(i)), var
732 END DO
733 END IF
734
735 END SUBROUTINE nnp_print_force_sigma
736
737! **************************************************************************************************
738!> \brief Print structures with extrapolation warning
739!> \param nnp ...
740!> \param unit_nr ...
741!> \date 2020-10-10
742!> \author Christoph Schran (christoph.schran@rub.de)
743! **************************************************************************************************
744 SUBROUTINE nnp_print_expol(nnp, unit_nr)
745 TYPE(nnp_type), INTENT(INOUT) :: nnp
746 INTEGER, INTENT(IN) :: unit_nr
747
748 CHARACTER(len=default_path_length) :: fmt_string
749 INTEGER :: i
750 REAL(kind=dp) :: mass, unit_conv
751 REAL(kind=dp), DIMENSION(3) :: com
752
753 nnp%expol = nnp%expol + 1
754 WRITE (unit_nr, *) nnp%num_atoms
755 WRITE (unit_nr, "(A,1X,I6)") "NNP extrapolation point N =", nnp%expol
756
757 ! move to COM of solute and wrap the box
758 ! coord not needed afterwards, therefore manipulation ok
759 com = 0.0_dp
760 mass = 0.0_dp
761 DO i = 1, nnp%num_atoms
762 CALL get_ptable_info(nnp%atoms(i) (1:2), amass=unit_conv)
763 com(:) = com(:) + nnp%coord(:, i)*unit_conv
764 mass = mass + unit_conv
765 END DO
766 com(:) = com(:)/mass
767
768 DO i = 1, nnp%num_atoms
769 nnp%coord(:, i) = nnp%coord(:, i) - com(:)
770 END DO
771
772 ! write out coordinates
773 unit_conv = cp_unit_from_cp2k(1.0_dp, trim("angstrom"))
774 fmt_string = "(A4,1X,3F20.10)"
775 DO i = 1, nnp%num_atoms
776 WRITE (unit_nr, fmt_string) &
777 nnp%atoms(nnp%sort_inv(i)), &
778 nnp%coord(1, nnp%sort_inv(i))*unit_conv, &
779 nnp%coord(2, nnp%sort_inv(i))*unit_conv, &
780 nnp%coord(3, nnp%sort_inv(i))*unit_conv
781 END DO
782
783 END SUBROUTINE nnp_print_expol
784
785! **************************************************************************************************
786!> \brief Print properties number according the requests in input file
787!> \param nnp ...
788!> \param print_section ...
789!> \date 2020-10-10
790!> \author Christoph Schran (christoph.schran@rub.de)
791! **************************************************************************************************
792 SUBROUTINE nnp_print_bias(nnp, print_section)
793 TYPE(nnp_type), INTENT(INOUT) :: nnp
794 TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
795
796 INTEGER :: unit_nr
797 LOGICAL :: file_is_new
798 TYPE(cp_logger_type), POINTER :: logger
799 TYPE(section_vals_type), POINTER :: print_key
800
801 NULLIFY (logger, print_key)
802 logger => cp_get_default_logger()
803
804 print_key => section_vals_get_subs_vals(print_section, "BIAS_ENERGY")
805 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
806 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
807 middle_name="nnp-bias-energy", is_new_file=file_is_new)
808 IF (unit_nr > 0) CALL nnp_print_bias_energy(nnp, unit_nr, file_is_new)
809 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
810 END IF
811
812 print_key => section_vals_get_subs_vals(print_section, "BIAS_FORCES")
813 IF (btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
814 unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
815 middle_name="nnp-bias-forces")
816 IF (unit_nr > 0) CALL nnp_print_bias_forces(nnp, unit_nr)
817 CALL cp_print_key_finished_output(unit_nr, logger, print_key)
818 END IF
819
820 END SUBROUTINE nnp_print_bias
821
822! **************************************************************************************************
823!> \brief Print NNP bias energies
824!> \param nnp ...
825!> \param unit_nr ...
826!> \param file_is_new ...
827!> \date 2020-10-10
828!> \author Christoph Schran (christoph.schran@rub.de)
829! **************************************************************************************************
830 SUBROUTINE nnp_print_bias_energy(nnp, unit_nr, file_is_new)
831 TYPE(nnp_type), INTENT(INOUT) :: nnp
832 INTEGER, INTENT(IN) :: unit_nr
833 LOGICAL, INTENT(IN) :: file_is_new
834
835 CHARACTER(len=default_path_length) :: fmt_string
836 INTEGER :: i
837
838 IF (file_is_new) THEN
839 WRITE (unit_nr, "(A1)", advance='no') "#"
840 WRITE (unit_nr, "(2(2X,A19))", advance='no') "Sigma [a.u.]", "Bias energy [a.u.]"
841 DO i = 1, nnp%n_committee
842 IF (nnp%bias_align) THEN
843 WRITE (unit_nr, "(2X,A16,I3)", advance='no') "shifted E_NNP", i
844 ELSE
845 WRITE (unit_nr, "(2X,A16,I3)", advance='no') "E_NNP", i
846 END IF
847 END DO
848 WRITE (unit_nr, "(A)") ""
849
850 END IF
851
852 WRITE (fmt_string, "(A,I3,A)") "(2X,", nnp%n_committee + 2, "(F20.9,1X))"
853 WRITE (unit_nr, fmt_string) nnp%bias_sigma, nnp%bias_energy, nnp%committee_energy
854
855 END SUBROUTINE nnp_print_bias_energy
856
857! **************************************************************************************************
858!> \brief Print NNP bias forces
859!> \param nnp ...
860!> \param unit_nr ...
861!> \date 2020-10-10
862!> \author Christoph Schran (christoph.schran@rub.de)
863! **************************************************************************************************
864 SUBROUTINE nnp_print_bias_forces(nnp, unit_nr)
865 TYPE(nnp_type), INTENT(INOUT) :: nnp
866 INTEGER, INTENT(IN) :: unit_nr
867
868 CHARACTER(len=default_path_length) :: fmt_string
869 INTEGER :: i
870
871 WRITE (unit_nr, *) nnp%num_atoms
872 WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
873
874 fmt_string = "(A4,1X,3F20.10)"
875 DO i = 1, nnp%num_atoms
876 WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(i)), nnp%bias_forces(:, i)
877 END DO
878
879 END SUBROUTINE nnp_print_bias_forces
880
881! **************************************************************************************************
882!> \brief Print NNP summed forces
883!> \param nnp ...
884!> \param print_section ...
885!> \param unit_nr ...
886!> \param file_is_new ...
887!> \date 2020-10-10
888!> \author Christoph Schran (christoph.schran@rub.de)
889! **************************************************************************************************
890 SUBROUTINE nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
891 TYPE(nnp_type), INTENT(INOUT) :: nnp
892 TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
893 INTEGER, INTENT(IN) :: unit_nr
894 LOGICAL, INTENT(IN) :: file_is_new
895
896 CHARACTER(len=default_path_length) :: fmt_string
897 CHARACTER(LEN=default_string_length), &
898 DIMENSION(:), POINTER :: atomlist
899 INTEGER :: i, ig, j, n
900 REAL(kind=dp), DIMENSION(3) :: rvec
901
902 NULLIFY (atomlist)
903 IF (file_is_new) THEN
904 WRITE (unit_nr, "(A)") "# Summed forces [a.u.]"
905 END IF
906
907 rvec = 0.0_dp
908
909 ! get atoms to sum over:
910 CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
911 c_vals=atomlist)
912 IF (ASSOCIATED(atomlist)) THEN
913 n = SIZE(atomlist)
914 DO i = 1, nnp%num_atoms
915 DO j = 1, n
916 ig = nnp%sort_inv(i)
917 IF (trim(adjustl(atomlist(j))) == trim(adjustl(nnp%atoms(ig)))) THEN
918 rvec(:) = rvec(:) + nnp%nnp_forces(:, i)
919 END IF
920 END DO
921 END DO
922 END IF
923
924 fmt_string = "(3(F20.10,1X))"
925 WRITE (unit_nr, fmt_string) rvec
926
927 END SUBROUTINE nnp_print_sumforces
928
929! **************************************************************************************************
930!> \brief Scatter the per-neighbour dG/dr arrays held in the nnp_neighbor_workspace
931!> into a per-atom Cartesian force array. Three contributions per central
932!> atom i:
933!> 1) self: self_dGdr(:, m) -> force_xyz(:, i)
934!> 2) radial: walk neighbor%rad(s)%ind(j), read dGdr_rad(s)%data(:, sf, j)
935!> 3) angular: walk neighbor%ang1(s)%ind / ang2(s)%ind, read dGdr_ang_jj/kk
936!> No global (n_sf, num_atoms) slab is touched. Reused by the helium-NNP
937!> coupling in helium_interactions.
938!>
939!> Precondition: nnp_calc_acsf(nnp, i, calc_forces=.TRUE.[, stress]) must
940!> have run for this same atom i immediately before; the per-element
941!> workspace it fills is overwritten on every ACSF call.
942!> \param nnp ...
943!> \param ind central-atom element index
944!> \param i central atom index (absolute)
945!> \param denergydsym dE/dG_k for k = 1..n_input_nodes(ind)
946!> \param force_xyz (3, num_atoms) destination -- accumulated, NOT overwritten
947!> \author Dhruv Sharma (ds2173@cam.ac.uk)
948! **************************************************************************************************
949 SUBROUTINE nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, force_xyz)
950
951 TYPE(nnp_type), INTENT(IN) :: nnp
952 INTEGER, INTENT(IN) :: ind, i
953 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: denergydsym
954 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: force_xyz
955
956 INTEGER :: j, k, k_atom, m, n_ang1_s, n_ang2_s, &
957 n_input_nodes, n_symf_s, off, s, sf, &
958 tid
959 LOGICAL :: homo_grp
960 REAL(kind=dp) :: de
961
962 tid = 1
963!$ tid = omp_get_thread_num() + 1
964 associate(workspace => nnp%neighbor_interface_state%workspace(ind, tid), &
965 neighbor => nnp%neighbor_interface_state%workspace(ind, tid)%neighbor, &
966 self_dgdr => nnp%neighbor_interface_state%workspace(ind, tid)%self_dGdr)
967
968 n_input_nodes = workspace%n_input_nodes
969
970 ! Self contribution
971 DO j = 1, n_input_nodes
972 de = denergydsym(j)
973 force_xyz(1, i) = force_xyz(1, i) - de*self_dgdr(1, j)
974 force_xyz(2, i) = force_xyz(2, i) - de*self_dgdr(2, j)
975 force_xyz(3, i) = force_xyz(3, i) - de*self_dgdr(3, j)
976 END DO
977
978 ! Radial neighbours
979 DO s = 1, nnp%rad(ind)%n_symfgrp
980 n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
981 associate(rad_buf => workspace%dGdr_rad(s)%data)
982 DO j = 1, neighbor%n_rad(s)
983 k_atom = neighbor%rad(s)%ind(j)
984 DO sf = 1, n_symf_s
985 m = nnp%rad(ind)%symfgrp(s)%symf(sf)
986 de = denergydsym(m)
987 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*rad_buf(1, sf, j)
988 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*rad_buf(2, sf, j)
989 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*rad_buf(3, sf, j)
990 END DO
991 END DO
992 END associate
993 END DO
994
995 ! Angular neighbours
996 off = nnp%n_rad(ind)
997 DO s = 1, nnp%ang(ind)%n_symfgrp
998 n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
999 n_ang1_s = neighbor%n_ang1(s)
1000 homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
1001 IF (homo_grp) THEN
1002 n_ang2_s = n_ang1_s ! kk buffer is also indexed in ang1(s)%ind
1003 ELSE
1004 n_ang2_s = neighbor%n_ang2(s)
1005 END IF
1006
1007 associate(jj_buf => workspace%dGdr_ang_jj(s)%data, &
1008 kk_buf => workspace%dGdr_ang_kk(s)%data)
1009 ! jj-side
1010 DO j = 1, n_ang1_s
1011 k_atom = neighbor%ang1(s)%ind(j)
1012 DO sf = 1, n_symf_s
1013 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
1014 de = denergydsym(m)
1015 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*jj_buf(1, sf, j)
1016 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*jj_buf(2, sf, j)
1017 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*jj_buf(3, sf, j)
1018 END DO
1019 END DO
1020
1021 ! kk-side: homo reads ang1(s)%ind, hetero reads ang2(s)%ind
1022 IF (homo_grp) THEN
1023 DO k = 1, n_ang2_s
1024 k_atom = neighbor%ang1(s)%ind(k)
1025 DO sf = 1, n_symf_s
1026 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
1027 de = denergydsym(m)
1028 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
1029 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
1030 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
1031 END DO
1032 END DO
1033 ELSE
1034 DO k = 1, n_ang2_s
1035 k_atom = neighbor%ang2(s)%ind(k)
1036 DO sf = 1, n_symf_s
1037 m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
1038 de = denergydsym(m)
1039 force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
1040 force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
1041 force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
1042 END DO
1043 END DO
1044 END IF
1045 END associate
1046 END DO
1047
1048 END associate
1049
1050 END SUBROUTINE nnp_scatter_dgdr_to_forces
1051
1052END MODULE nnp_force
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_finished_output
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,...
Definition cp_output_handling.F:1156
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:107
cp_output_handling::cp_print_key_should_output
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...
Definition cp_output_handling.F:345
cp_output_handling::cp_print_key_unit_nr
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, file_access)
...
Definition cp_output_handling.F:935
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
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
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
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:1311
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::nnp_gpu_prec_mixed
integer, parameter, public nnp_gpu_prec_mixed
Definition input_constants.F:72
input_constants::nnp_use_gpu_on
integer, parameter, public nnp_use_gpu_on
Definition input_constants.F:65
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
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
Definition input_section_types.F:735
input_section_types::section_vals_val_get
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
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
nnp_acsf
Functionality for atom centered symmetry functions for neural network potentials.
Definition nnp_acsf.F:16
nnp_acsf::nnp_calc_acsf
subroutine, public nnp_calc_acsf(nnp, i, calc_forces, arc, stress)
Calculate atom centered symmetry functions for given atom i.
Definition nnp_acsf.F:89
nnp_acsf::nnp_prepare_neighbor_cache
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:690
nnp_environment_types
Data types for neural network potentials.
Definition nnp_environment_types.F:14
nnp_environment_types::nnp_env_get
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.
Definition nnp_environment_types.F:829
nnp_force
Methods dealing with Neural Network potentials.
Definition nnp_force.F:15
nnp_force::nnp_scatter_dgdr_to_forces
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:950
nnp_force::nnp_calc_energy_force
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:73
nnp_gpu
GPU offload driver for the NNP force evaluation. Hosts the ISO_C_BINDING interfaces to the CUDA launc...
Definition nnp_gpu.F:16
nnp_gpu::nnp_gpu_energy_force
subroutine, public nnp_gpu_energy_force(nnp, istart, mecalc, n_input_max, calc_stress)
Run the whole per-step NNP force evaluation on the device: the cell-list walk, the descriptor build,...
Definition nnp_gpu.F:1318
nnp_model
Methods dealing with core routines for artificial neural networks.
Definition nnp_model.F:13
nnp_model::nnp_predict
subroutine, public nnp_predict(arc, nnp, i_com)
Predict energy by evaluating neural network.
Definition nnp_model.F:82
nnp_model::nnp_gradients
subroutine, public nnp_gradients(arc, nnp, i_com, denergydsym)
Calculate gradients of neural network.
Definition nnp_model.F:169
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
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 periodic_table.F:85
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
nnp_environment_types::nnp_arc_type
Data type for artificial neural networks.
Definition nnp_environment_types.F:468
nnp_environment_types::nnp_type
Main data type collecting all relevant data for neural network potentials.
Definition nnp_environment_types.F:97
particle_types::particle_type
Definition particle_types.F:35
virial_types::virial_type
Definition virial_types.F:26