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 : !> \author Claudio Malvino (claudiormal@gmail.com)
13 : !> \date 2020-10-10
14 : ! **************************************************************************************************
15 : MODULE 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_type
20 : USE cp_output_handling, ONLY: cp_p_file,&
21 : cp_print_key_finished_output,&
22 : cp_print_key_should_output,&
23 : cp_print_key_unit_nr
24 : USE cp_subsys_types, ONLY: cp_subsys_get,&
25 : cp_subsys_type
26 : USE cp_units, ONLY: cp_unit_from_cp2k
27 : USE distribution_1d_types, ONLY: distribution_1d_type
28 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
29 : section_vals_type,&
30 : section_vals_val_get
31 : USE kinds, ONLY: default_path_length,&
32 : default_string_length,&
33 : dp
34 : USE nnp_acsf, ONLY: nnp_calc_acsf,&
35 : nnp_prepare_neighbor_cache
36 : USE nnp_environment_types, ONLY: nnp_arc_type,&
37 : nnp_env_get,&
38 : nnp_type
39 : USE nnp_model, ONLY: nnp_gradients,&
40 : nnp_predict
41 : USE particle_types, ONLY: particle_type
42 : USE periodic_table, ONLY: get_ptable_info
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 :
56 : PUBLIC :: nnp_calc_energy_force, nnp_scatter_dgdr_to_forces
57 :
58 : CONTAINS
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 320 : 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 320 : INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
79 : LOGICAL :: calc_stress
80 320 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
81 320 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: stress
82 320 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: my_cstress_t, my_force_t
83 320 : 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 320 : TYPE(nnp_arc_type), ALLOCATABLE, DIMENSION(:) :: my_arc
88 320 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
89 : TYPE(section_vals_type), POINTER :: print_section
90 : TYPE(virial_type), POINTER :: virial
91 :
92 320 : CALL timeset(routineN, handle)
93 :
94 320 : NULLIFY (particle_set, logger, local_particles, subsys, &
95 320 : atomic_kind_set)
96 320 : logger => cp_get_default_logger()
97 :
98 320 : CPASSERT(ASSOCIATED(nnp))
99 320 : 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 320 : atomic_kind_set=atomic_kind_set)
103 :
104 : CALL cp_subsys_get(subsys, &
105 320 : virial=virial)
106 :
107 320 : calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
108 22 : IF (calc_stress .AND. .NOT. calc_forces) THEN
109 0 : CPABORT('Stress cannot be calculated without forces')
110 : END IF
111 :
112 486710 : nnp%atomic_energy(:, :) = 0.0_dp
113 991254 : IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
114 991254 : IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
115 2608 : IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
116 :
117 : !fill coord array
118 320 : ig = 1
119 960 : DO i = 1, nnp%n_ele
120 129348 : DO j = 1, nnp%num_atoms
121 129028 : IF (nnp%ele(i) == particle_set(j)%atomic_kind%element_symbol) THEN
122 256776 : DO m = 1, 3
123 256776 : nnp%coord(m, ig) = particle_set(j)%r(m)
124 : END DO
125 64194 : nnp%atoms(ig) = nnp%ele(i)
126 64194 : CALL get_ptable_info(nnp%atoms(ig) (1:2), number=nnp%nuc_atoms(ig))
127 64194 : nnp%ele_ind(ig) = i
128 64194 : nnp%sort(ig) = j
129 64194 : nnp%sort_inv(j) = ig
130 64194 : 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 320 : (logger%para_env%mepos + 1), 1)
139 960 : ALLOCATE (allcalc(logger%para_env%num_pe))
140 320 : allcalc(:) = 0
141 320 : CALL logger%para_env%allgather(mecalc, allcalc)
142 320 : istart = 1
143 480 : DO i = 2, logger%para_env%mepos + 1
144 480 : istart = istart + allcalc(i - 1)
145 : END DO
146 :
147 : ! reset extrapolation status
148 320 : nnp%output_expol = .FALSE.
149 :
150 320 : CALL nnp_prepare_neighbor_cache(nnp)
151 :
152 320 : max_input_nodes = 0
153 960 : DO i = 1, nnp%n_ele
154 960 : 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 320 : nthreads = 1
159 320 : !$ nthreads = omp_get_max_threads()
160 320 : IF (calc_forces) THEN
161 : ! dGdr lives in per-element workspace; no global (3, n_sf, num_atoms) slab.
162 830 : ALLOCATE (my_force_t(3, nnp%num_atoms, nnp%n_committee, 0:nthreads - 1))
163 166 : my_force_t(:, :, :, :) = 0.0_dp
164 166 : IF (calc_stress) THEN
165 88 : ALLOCATE (my_cstress_t(3, 3, nnp%n_committee, 0:nthreads - 1))
166 22 : 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 320 : 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 320 : !$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 320 : IF (calc_forces) THEN
308 332 : DO t = 0, nthreads - 1
309 991266 : nnp%myforce(:, :, :) = nnp%myforce(:, :, :) + my_force_t(:, :, :, t)
310 : END DO
311 166 : DEALLOCATE (my_force_t)
312 166 : IF (calc_stress) THEN
313 44 : DO t = 0, nthreads - 1
314 2332 : nnp%committee_stress(:, :, :) = nnp%committee_stress(:, :, :) + my_cstress_t(:, :, :, t)
315 : END DO
316 22 : DEALLOCATE (my_cstress_t)
317 : END IF
318 : END IF
319 :
320 320 : CALL timestop(handle_loop)
321 :
322 : ! calculate energy:
323 320 : CALL logger%para_env%sum(nnp%atomic_energy(:, :))
324 486710 : nnp%committee_energy(:) = SUM(nnp%atomic_energy, 1)
325 2726 : nnp%nnp_potential_energy = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
326 :
327 320 : IF (calc_forces) THEN
328 : ! bring myforce to force array
329 34792 : DO j = 1, nnp%num_atoms
330 138670 : DO k = 1, 3
331 880824 : nnp%committee_forces(k, (nnp%sort(j)), :) = nnp%myforce(k, j, :)
332 : END DO
333 : END DO
334 166 : CALL logger%para_env%sum(nnp%committee_forces)
335 880990 : nnp%nnp_forces(:, :) = SUM(nnp%committee_forces, DIM=3)/REAL(nnp%n_committee, dp)
336 34792 : DO j = 1, nnp%num_atoms
337 242548 : particle_set(j)%f(:) = nnp%nnp_forces(:, j)
338 : END DO
339 : END IF
340 :
341 320 : IF (calc_stress) THEN
342 22 : CALL logger%para_env%sum(nnp%committee_stress)
343 2134 : 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 320 : IF (nnp%bias) THEN
348 44 : CALL nnp_bias_sigma(nnp, calc_forces)
349 44 : nnp%nnp_potential_energy = nnp%nnp_potential_energy + nnp%bias_energy
350 44 : IF (calc_forces) THEN
351 8492 : DO j = 1, nnp%num_atoms
352 59180 : particle_set(j)%f(:) = particle_set(j)%f(:) + nnp%bias_forces(:, j)
353 : END DO
354 : END IF
355 : ! print properties if requested
356 44 : print_section => section_vals_get_subs_vals(nnp%nnp_input, "BIAS%PRINT")
357 44 : CALL nnp_print_bias(nnp, print_section)
358 : END IF
359 :
360 : ! print properties if requested
361 320 : print_section => section_vals_get_subs_vals(nnp%nnp_input, "PRINT")
362 320 : CALL nnp_print(nnp, print_section)
363 :
364 320 : DEALLOCATE (allcalc)
365 :
366 320 : CALL timestop(handle)
367 :
368 960 : 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 44 : 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 44 : CALL timeset(routineN, handle)
387 :
388 : ! init
389 44 : sigma = 0.0_dp
390 44 : nnp%bias_energy = 0.0_dp
391 33836 : IF (calc_forces) nnp%bias_forces = 0.0_dp
392 :
393 : ! Subtract reference energy of each committee member, if requested
394 44 : IF (nnp%bias_align) THEN
395 : ! committee energy not used afterward, therefore overwritten
396 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 396 : avrg = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
403 396 : DO i = 1, nnp%n_committee
404 396 : sigma = sigma + (nnp%committee_energy(i) - avrg)**2
405 : END DO
406 44 : sigma = SQRT(sigma/REAL(nnp%n_committee, dp))
407 44 : nnp%bias_sigma = sigma
408 :
409 44 : IF (sigma > nnp%bias_sigma0) THEN
410 : ! E_b = 0.5 * kb * (sigma - sigma_0)**2
411 44 : nnp%bias_energy = 0.5_dp*nnp%bias_kb*(sigma - nnp%bias_sigma0)**2
412 :
413 44 : 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 44 : pref = nnp%bias_kb*(1.0_dp - nnp%bias_sigma0/sigma)
418 396 : DO i = 1, nnp%n_committee
419 : nnp%bias_forces(:, :) = nnp%bias_forces(:, :) + &
420 : (nnp%committee_energy(i) - avrg)* &
421 270732 : (nnp%committee_forces(:, :, i) - nnp%nnp_forces(:, :))
422 : END DO
423 44 : pref = pref/REAL(nnp%n_committee, dp)
424 33836 : nnp%bias_forces(:, :) = nnp%bias_forces(:, :)*pref
425 : END IF
426 : END IF
427 :
428 44 : CALL timestop(handle)
429 :
430 44 : 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 640 : 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 320 : NULLIFY (logger, print_key)
449 320 : logger => cp_get_default_logger()
450 :
451 320 : print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
452 320 : 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 310 : middle_name="nnp-energies", is_new_file=file_is_new)
455 310 : IF (unit_nr > 0) CALL nnp_print_energies(nnp, unit_nr, file_is_new)
456 310 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
457 : END IF
458 :
459 320 : print_key => section_vals_get_subs_vals(print_section, "FORCES")
460 320 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
461 286 : CALL nnp_print_forces(nnp, print_key)
462 : END IF
463 :
464 320 : print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
465 320 : 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 264 : middle_name="nnp-forces-std")
468 264 : IF (unit_nr > 0) CALL nnp_print_force_sigma(nnp, unit_nr)
469 264 : 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 320 : CALL logger%para_env%sum(nnp%output_expol)
474 320 : IF (nnp%output_expol) THEN
475 22 : print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
476 22 : 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 22 : middle_name="nnp-extrapolation")
479 22 : IF (unit_nr > 0) CALL nnp_print_expol(nnp, unit_nr)
480 22 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
481 : END IF
482 : END IF
483 :
484 320 : 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 320 : explicit=explicit)
488 320 : IF (explicit) THEN
489 0 : 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 0 : middle_name="nnp-sumforce", is_new_file=file_is_new)
492 0 : IF (unit_nr > 0) CALL nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
493 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
494 : END IF
495 : END IF
496 :
497 320 : 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 155 : 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 155 : IF (file_is_new) THEN
517 7 : WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
518 7 : WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
519 56 : DO i = 1, nnp%n_committee
520 56 : WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
521 : END DO
522 7 : WRITE (unit_nr, "(A)") ""
523 : END IF
524 :
525 155 : fmt_string = "(2X, F20.9)"
526 155 : WRITE (UNIT=fmt_string(5:6), FMT="(I2)") nnp%n_committee + 2
527 212590 : std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
528 155 : std = std/REAL(nnp%n_committee, dp)
529 155 : std = SQRT(std)
530 212590 : WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
531 :
532 155 : 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 286 : 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 286 : NULLIFY (logger)
550 286 : logger => cp_get_default_logger()
551 :
552 2420 : DO i = 1, nnp%n_committee
553 2134 : WRITE (fmt_string, *) i
554 2134 : 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 2134 : middle_name=TRIM(middle_name))
557 2134 : IF (unit_nr > 0) THEN
558 1067 : WRITE (unit_nr, *) nnp%num_atoms
559 1067 : WRITE (unit_nr, "(A,1X,A,A,F20.9)") "NNP forces [a.u.] of committee member", &
560 2134 : ADJUSTL(TRIM(fmt_string)), "energy [a.u.]=", nnp%committee_energy(i)
561 :
562 1067 : fmt_string = "(A4,1X,3F20.10)"
563 205931 : DO j = 1, nnp%num_atoms
564 205931 : WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(j)), nnp%committee_forces(:, j, i)
565 : END DO
566 : END IF
567 2420 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
568 : END DO
569 :
570 286 : 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 132 : 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 132 : IF (unit_nr > 0) THEN
587 132 : WRITE (unit_nr, *) nnp%num_atoms
588 132 : WRITE (unit_nr, "(A,1X,A)") "NNP sigma of forces [a.u.]"
589 :
590 25476 : DO i = 1, nnp%num_atoms
591 25344 : var = 0.0_dp
592 228096 : DO j = 1, nnp%n_committee
593 836352 : var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
594 : END DO
595 101376 : var = var/REAL(nnp%n_committee, dp)
596 101376 : var = SQRT(var)
597 25476 : WRITE (unit_nr, "(A4,1X,3F20.10)") nnp%atoms(nnp%sort_inv(i)), var
598 : END DO
599 : END IF
600 :
601 132 : 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 11 : 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 11 : nnp%expol = nnp%expol + 1
620 11 : WRITE (unit_nr, *) nnp%num_atoms
621 11 : 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 11 : com = 0.0_dp
626 11 : mass = 0.0_dp
627 44 : DO i = 1, nnp%num_atoms
628 33 : CALL get_ptable_info(nnp%atoms(i) (1:2), amass=unit_conv)
629 132 : com(:) = com(:) + nnp%coord(:, i)*unit_conv
630 77 : mass = mass + unit_conv
631 : END DO
632 44 : com(:) = com(:)/mass
633 :
634 44 : DO i = 1, nnp%num_atoms
635 143 : nnp%coord(:, i) = nnp%coord(:, i) - com(:)
636 : END DO
637 :
638 : ! write out coordinates
639 11 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
640 11 : fmt_string = "(A4,1X,3F20.10)"
641 44 : DO i = 1, nnp%num_atoms
642 : WRITE (unit_nr, fmt_string) &
643 33 : nnp%atoms(nnp%sort_inv(i)), &
644 33 : nnp%coord(1, nnp%sort_inv(i))*unit_conv, &
645 33 : nnp%coord(2, nnp%sort_inv(i))*unit_conv, &
646 77 : nnp%coord(3, nnp%sort_inv(i))*unit_conv
647 : END DO
648 :
649 11 : 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 44 : 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 44 : NULLIFY (logger, print_key)
668 44 : logger => cp_get_default_logger()
669 :
670 44 : print_key => section_vals_get_subs_vals(print_section, "BIAS_ENERGY")
671 44 : 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 44 : middle_name="nnp-bias-energy", is_new_file=file_is_new)
674 44 : IF (unit_nr > 0) CALL nnp_print_bias_energy(nnp, unit_nr, file_is_new)
675 44 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
676 : END IF
677 :
678 44 : print_key => section_vals_get_subs_vals(print_section, "BIAS_FORCES")
679 44 : 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 0 : middle_name="nnp-bias-forces")
682 0 : IF (unit_nr > 0) CALL nnp_print_bias_forces(nnp, unit_nr)
683 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
684 : END IF
685 :
686 44 : 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 22 : 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 22 : IF (file_is_new) THEN
705 1 : WRITE (unit_nr, "(A1)", ADVANCE='no') "#"
706 1 : WRITE (unit_nr, "(2(2X,A19))", ADVANCE='no') "Sigma [a.u.]", "Bias energy [a.u.]"
707 9 : DO i = 1, nnp%n_committee
708 9 : IF (nnp%bias_align) THEN
709 8 : WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "shifted E_NNP", i
710 : ELSE
711 0 : WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "E_NNP", i
712 : END IF
713 : END DO
714 1 : WRITE (unit_nr, "(A)") ""
715 :
716 : END IF
717 :
718 22 : WRITE (fmt_string, "(A,I3,A)") "(2X,", nnp%n_committee + 2, "(F20.9,1X))"
719 22 : WRITE (unit_nr, fmt_string) nnp%bias_sigma, nnp%bias_energy, nnp%committee_energy
720 :
721 22 : 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 0 : 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 0 : WRITE (unit_nr, *) nnp%num_atoms
738 0 : WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
739 :
740 0 : fmt_string = "(A4,1X,3F20.10)"
741 0 : DO i = 1, nnp%num_atoms
742 0 : WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(i)), nnp%bias_forces(:, i)
743 : END DO
744 :
745 0 : 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 0 : 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 0 : DIMENSION(:), POINTER :: atomlist
765 : INTEGER :: i, ig, j, n
766 : REAL(KIND=dp), DIMENSION(3) :: rvec
767 :
768 0 : NULLIFY (atomlist)
769 0 : IF (file_is_new) THEN
770 0 : WRITE (unit_nr, "(A)") "# Summed forces [a.u.]"
771 : END IF
772 :
773 0 : rvec = 0.0_dp
774 :
775 : ! get atoms to sum over:
776 : CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
777 0 : c_vals=atomlist)
778 0 : IF (ASSOCIATED(atomlist)) THEN
779 0 : n = SIZE(atomlist)
780 0 : DO i = 1, nnp%num_atoms
781 0 : DO j = 1, n
782 0 : ig = nnp%sort_inv(i)
783 0 : IF (TRIM(ADJUSTL(atomlist(j))) == TRIM(ADJUSTL(nnp%atoms(ig)))) THEN
784 0 : rvec(:) = rvec(:) + nnp%nnp_forces(:, i)
785 : END IF
786 : END DO
787 : END DO
788 : END IF
789 :
790 0 : fmt_string = "(3(F20.10,1X))"
791 0 : WRITE (unit_nr, fmt_string) rvec
792 :
793 0 : 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 130120 : 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 130120 : tid = 1
829 130120 : !$ 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 130120 : n_input_nodes = workspace%n_input_nodes
835 :
836 : ! Self contribution
837 3725480 : DO j = 1, n_input_nodes
838 3595360 : de = denergydsym(j)
839 3595360 : force_xyz(1, i) = force_xyz(1, i) - de*self_dGdr(1, j)
840 3595360 : force_xyz(2, i) = force_xyz(2, i) - de*self_dGdr(2, j)
841 3725480 : force_xyz(3, i) = force_xyz(3, i) - de*self_dGdr(3, j)
842 : END DO
843 :
844 : ! Radial neighbours
845 393560 : DO s = 1, nnp%rad(ind)%n_symfgrp
846 263440 : n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
847 130120 : ASSOCIATE (rad_buf => workspace%dGdr_rad(s)%data)
848 13413642 : DO j = 1, neighbor%n_rad(s)
849 13150202 : k_atom = neighbor%rad(s)%ind(j)
850 118591990 : DO sf = 1, n_symf_s
851 105178348 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
852 105178348 : de = denergydsym(m)
853 105178348 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*rad_buf(1, sf, j)
854 105178348 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*rad_buf(2, sf, j)
855 118328550 : 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 130120 : off = nnp%n_rad(ind)
863 561720 : DO s = 1, nnp%ang(ind)%n_symfgrp
864 301480 : n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
865 301480 : n_ang1_s = neighbor%n_ang1(s)
866 301480 : homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
867 301480 : IF (homo_grp) THEN
868 : n_ang2_s = n_ang1_s ! kk buffer is also indexed in ang1(s)%ind
869 : ELSE
870 133320 : n_ang2_s = neighbor%n_ang2(s)
871 : END IF
872 :
873 : ASSOCIATE (jj_buf => workspace%dGdr_ang_jj(s)%data, &
874 130120 : kk_buf => workspace%dGdr_ang_kk(s)%data)
875 : ! jj-side
876 16399965 : DO j = 1, n_ang1_s
877 16098485 : k_atom = neighbor%ang1(s)%ind(j)
878 104160489 : DO sf = 1, n_symf_s
879 87760524 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
880 87760524 : de = denergydsym(m)
881 87760524 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*jj_buf(1, sf, j)
882 87760524 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*jj_buf(2, sf, j)
883 103859009 : 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 602960 : IF (homo_grp) THEN
889 7483470 : DO k = 1, n_ang2_s
890 7315310 : k_atom = neighbor%ang1(s)%ind(k)
891 42632218 : DO sf = 1, n_symf_s
892 35148748 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
893 35148748 : de = denergydsym(m)
894 35148748 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
895 35148748 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
896 42464058 : 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 4503547 : DO k = 1, n_ang2_s
901 4370227 : k_atom = neighbor%ang2(s)%ind(k)
902 30817724 : DO sf = 1, n_symf_s
903 26314177 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
904 26314177 : de = denergydsym(m)
905 26314177 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
906 26314177 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
907 30684404 : 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 130120 : END SUBROUTINE nnp_scatter_dgdr_to_forces
917 :
918 : END MODULE nnp_force
|
