Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Methods dealing with Neural Network potentials
10 : !> \author Christoph Schran (christoph.schran@rub.de)
11 : !> \author Dhruv Sharma (ds2173@cam.ac.uk)
12 : !> \date 2020-10-10
13 : ! **************************************************************************************************
14 : MODULE nnp_force
15 :
16 : USE atomic_kind_types, ONLY: atomic_kind_type
17 : USE cp_log_handling, ONLY: cp_get_default_logger,&
18 : cp_logger_type
19 : USE cp_output_handling, ONLY: cp_p_file,&
20 : cp_print_key_finished_output,&
21 : cp_print_key_should_output,&
22 : cp_print_key_unit_nr
23 : USE cp_subsys_types, ONLY: cp_subsys_get,&
24 : cp_subsys_type
25 : USE cp_units, ONLY: cp_unit_from_cp2k
26 : USE distribution_1d_types, ONLY: distribution_1d_type
27 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
28 : section_vals_type,&
29 : section_vals_val_get
30 : USE kinds, ONLY: default_path_length,&
31 : default_string_length,&
32 : dp
33 : USE nnp_acsf, ONLY: nnp_calc_acsf,&
34 : nnp_prepare_neighbor_cache
35 : USE nnp_environment_types, ONLY: nnp_env_get,&
36 : nnp_type
37 : USE nnp_model, ONLY: nnp_gradients,&
38 : nnp_predict
39 : USE particle_types, ONLY: particle_type
40 : USE periodic_table, ONLY: get_ptable_info
41 : USE physcon, ONLY: angstrom
42 : USE virial_types, ONLY: virial_type
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : PRIVATE
48 :
49 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
50 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_force'
51 :
52 : PUBLIC :: nnp_calc_energy_force, nnp_scatter_dgdr_to_forces
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Calculate the energy and force for a given configuration with the NNP
58 : !> \param nnp ...
59 : !> \param calc_forces ...
60 : !> \date 2020-10-10
61 : !> \author Christoph Schran (christoph.schran@rub.de)
62 : ! **************************************************************************************************
63 320 : SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
64 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
65 : LOGICAL, INTENT(IN) :: calc_forces
66 :
67 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_calc_energy_force'
68 :
69 : INTEGER :: handle, handle_loop, i, i_com, ig, ind, &
70 : istart, j, k, m, max_input_nodes, &
71 : mecalc, n_input_nodes
72 320 : INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
73 : LOGICAL :: calc_stress
74 320 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
75 320 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: stress
76 320 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
77 : TYPE(cp_logger_type), POINTER :: logger
78 : TYPE(cp_subsys_type), POINTER :: subsys
79 : TYPE(distribution_1d_type), POINTER :: local_particles
80 320 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
81 : TYPE(section_vals_type), POINTER :: print_section
82 : TYPE(virial_type), POINTER :: virial
83 :
84 320 : CALL timeset(routineN, handle)
85 :
86 320 : NULLIFY (particle_set, logger, local_particles, subsys, &
87 320 : atomic_kind_set)
88 320 : logger => cp_get_default_logger()
89 :
90 320 : CPASSERT(ASSOCIATED(nnp))
91 320 : CPASSERT(nnp%n_committee >= 1)
92 : CALL nnp_env_get(nnp_env=nnp, particle_set=particle_set, &
93 : subsys=subsys, local_particles=local_particles, &
94 320 : atomic_kind_set=atomic_kind_set)
95 :
96 : CALL cp_subsys_get(subsys, &
97 320 : virial=virial)
98 :
99 320 : calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
100 22 : IF (calc_stress .AND. .NOT. calc_forces) THEN
101 0 : CPABORT('Stress cannot be calculated without forces')
102 : END IF
103 :
104 486710 : nnp%atomic_energy(:, :) = 0.0_dp
105 991254 : IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
106 991254 : IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
107 2608 : IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
108 :
109 : !fill coord array
110 320 : ig = 1
111 960 : DO i = 1, nnp%n_ele
112 129348 : DO j = 1, nnp%num_atoms
113 129028 : IF (nnp%ele(i) == particle_set(j)%atomic_kind%element_symbol) THEN
114 256776 : DO m = 1, 3
115 256776 : nnp%coord(m, ig) = particle_set(j)%r(m)
116 : END DO
117 64194 : nnp%atoms(ig) = nnp%ele(i)
118 64194 : CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
119 64194 : nnp%ele_ind(ig) = i
120 64194 : nnp%sort(ig) = j
121 64194 : nnp%sort_inv(j) = ig
122 64194 : ig = ig + 1
123 : END IF
124 : END DO
125 : END DO
126 :
127 : ! parallization:
128 : mecalc = nnp%num_atoms/logger%para_env%num_pe + &
129 : MIN(MOD(nnp%num_atoms, logger%para_env%num_pe)/ &
130 320 : (logger%para_env%mepos + 1), 1)
131 960 : ALLOCATE (allcalc(logger%para_env%num_pe))
132 320 : allcalc(:) = 0
133 320 : CALL logger%para_env%allgather(mecalc, allcalc)
134 320 : istart = 1
135 480 : DO i = 2, logger%para_env%mepos + 1
136 480 : istart = istart + allcalc(i - 1)
137 : END DO
138 :
139 : ! reset extrapolation status
140 320 : nnp%output_expol = .FALSE.
141 :
142 320 : CALL nnp_prepare_neighbor_cache(nnp)
143 :
144 320 : max_input_nodes = 0
145 320 : IF (calc_forces) THEN
146 498 : DO i = 1, nnp%n_ele
147 498 : max_input_nodes = MAX(max_input_nodes, nnp%arc(i)%n_nodes(1))
148 : END DO
149 : ! dGdr lives in per-element workspace; no global (3, n_sf, num_atoms) slab.
150 166 : IF (calc_stress) THEN
151 66 : ALLOCATE (stress(3, 3, max_input_nodes))
152 22 : stress(:, :, :) = 0.0_dp
153 : END IF
154 498 : ALLOCATE (denergydsym(max_input_nodes))
155 166 : denergydsym(:) = 0.0_dp
156 : END IF
157 :
158 : ! Per-rank timer around the per-atom work. CP2K aggregates the label as
159 : ! avg/max across ranks, so max - avg measures the load imbalance.
160 320 : CALL timeset('nnp_per_rank_atom_loop', handle_loop)
161 :
162 : ! calc atomic contribution to energy and force
163 32417 : DO i = istart, istart + mecalc - 1
164 :
165 : ! determine index of atom type and offset
166 32097 : ind = nnp%ele_ind(i)
167 32097 : n_input_nodes = nnp%arc(ind)%n_nodes(1)
168 :
169 : ! reset input nodes of ele(ind):
170 930813 : nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
171 :
172 : ! compute sym fnct values
173 32097 : IF (calc_forces) THEN
174 : !reset input grads of ele(ind):
175 502077 : nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
176 17313 : IF (calc_stress) THEN
177 770880 : stress(:, :, 1:n_input_nodes) = 0.0_dp
178 2112 : CALL nnp_calc_acsf(nnp, i, .TRUE., stress(:, :, 1:n_input_nodes))
179 : ELSE
180 15201 : CALL nnp_calc_acsf(nnp, i, .TRUE.)
181 : END IF
182 : ELSE
183 14784 : CALL nnp_calc_acsf(nnp, i, .FALSE.)
184 : END IF
185 :
186 274409 : DO i_com = 1, nnp%n_committee
187 : ! predict energy
188 241992 : CALL nnp_predict(nnp%arc(ind), nnp, i_com)
189 : nnp%atomic_energy(i, i_com) = nnp%arc(ind)%layer(nnp%n_layer)%node(1) + &
190 241992 : nnp%atom_energies(ind)
191 :
192 : ! predict forces
193 274089 : IF (calc_forces) THEN
194 3587880 : denergydsym(1:n_input_nodes) = 0.0_dp
195 123720 : CALL nnp_gradients(nnp%arc(ind), nnp, i_com, denergydsym(1:n_input_nodes))
196 :
197 : ! Force scatter over the per-element workspace dGdr arrays; see
198 : ! nnp_scatter_dgdr_to_forces below, also reused by the helium-NNP
199 : ! coupling in src/motion/helium_interactions.F.
200 : CALL nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym(1:n_input_nodes), &
201 123720 : nnp%myforce(:, :, i_com))
202 :
203 123720 : IF (calc_stress) THEN
204 489984 : DO j = 1, n_input_nodes
205 : nnp%committee_stress(:, :, i_com) = nnp%committee_stress(:, :, i_com) - &
206 6167040 : denergydsym(j)*stress(:, :, j)
207 : END DO
208 : END IF
209 : END IF
210 : END DO
211 :
212 : END DO ! loop over num_atoms
213 :
214 320 : CALL timestop(handle_loop)
215 :
216 320 : IF (calc_forces) THEN
217 166 : DEALLOCATE (denergydsym)
218 166 : IF (calc_stress) THEN
219 22 : DEALLOCATE (stress)
220 : END IF
221 : END IF
222 :
223 : ! calculate energy:
224 320 : CALL logger%para_env%sum(nnp%atomic_energy(:, :))
225 486710 : nnp%committee_energy(:) = SUM(nnp%atomic_energy, 1)
226 2726 : nnp%nnp_potential_energy = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
227 :
228 320 : IF (calc_forces) THEN
229 : ! bring myforce to force array
230 34792 : DO j = 1, nnp%num_atoms
231 138670 : DO k = 1, 3
232 880824 : nnp%committee_forces(k, (nnp%sort(j)), :) = nnp%myforce(k, j, :)
233 : END DO
234 : END DO
235 166 : CALL logger%para_env%sum(nnp%committee_forces)
236 880990 : nnp%nnp_forces(:, :) = SUM(nnp%committee_forces, DIM=3)/REAL(nnp%n_committee, dp)
237 34792 : DO j = 1, nnp%num_atoms
238 242548 : particle_set(j)%f(:) = nnp%nnp_forces(:, j)
239 : END DO
240 : END IF
241 :
242 320 : IF (calc_stress) THEN
243 22 : CALL logger%para_env%sum(nnp%committee_stress)
244 2134 : virial%pv_virial = SUM(nnp%committee_stress, DIM=3)/REAL(nnp%n_committee, dp)
245 : END IF
246 :
247 : ! Bias the standard deviation of committee disagreement
248 320 : IF (nnp%bias) THEN
249 44 : CALL nnp_bias_sigma(nnp, calc_forces)
250 44 : nnp%nnp_potential_energy = nnp%nnp_potential_energy + nnp%bias_energy
251 44 : IF (calc_forces) THEN
252 8492 : DO j = 1, nnp%num_atoms
253 59180 : particle_set(j)%f(:) = particle_set(j)%f(:) + nnp%bias_forces(:, j)
254 : END DO
255 : END IF
256 : ! print properties if requested
257 44 : print_section => section_vals_get_subs_vals(nnp%nnp_input, "BIAS%PRINT")
258 44 : CALL nnp_print_bias(nnp, print_section)
259 : END IF
260 :
261 : ! print properties if requested
262 320 : print_section => section_vals_get_subs_vals(nnp%nnp_input, "PRINT")
263 320 : CALL nnp_print(nnp, print_section)
264 :
265 320 : DEALLOCATE (allcalc)
266 :
267 320 : CALL timestop(handle)
268 :
269 960 : END SUBROUTINE nnp_calc_energy_force
270 :
271 : ! **************************************************************************************************
272 : !> \brief Calculate bias potential and force based on standard deviation of committee disagreement
273 : !> \param nnp ...
274 : !> \param calc_forces ...
275 : !> \date 2020-10-10
276 : !> \author Christoph Schran (christoph.schran@rub.de)
277 : ! **************************************************************************************************
278 44 : SUBROUTINE nnp_bias_sigma(nnp, calc_forces)
279 : TYPE(nnp_type), INTENT(INOUT) :: nnp
280 : LOGICAL, INTENT(IN) :: calc_forces
281 :
282 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_bias_sigma'
283 :
284 : INTEGER :: handle, i
285 : REAL(KIND=dp) :: avrg, pref, sigma
286 :
287 44 : CALL timeset(routineN, handle)
288 :
289 : ! init
290 44 : sigma = 0.0_dp
291 44 : nnp%bias_energy = 0.0_dp
292 33836 : IF (calc_forces) nnp%bias_forces = 0.0_dp
293 :
294 : ! Subtract reference energy of each committee member, if requested
295 44 : IF (nnp%bias_align) THEN
296 : ! committee energy not used afterward, therefore overwritten
297 396 : nnp%committee_energy(:) = nnp%committee_energy(:) - nnp%bias_e_avrg(:)
298 : END IF
299 :
300 : ! <E> = 1/n sum(E_i)
301 : ! sigma = sqrt(1/n sum((E_i - <E>)**2))
302 : ! = sqrt(1/n sum(dE_i**2))
303 396 : avrg = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
304 396 : DO i = 1, nnp%n_committee
305 396 : sigma = sigma + (nnp%committee_energy(i) - avrg)**2
306 : END DO
307 44 : sigma = SQRT(sigma/REAL(nnp%n_committee, dp))
308 44 : nnp%bias_sigma = sigma
309 :
310 44 : IF (sigma > nnp%bias_sigma0) THEN
311 : ! E_b = 0.5 * kb * (sigma - sigma_0)**2
312 44 : nnp%bias_energy = 0.5_dp*nnp%bias_kb*(sigma - nnp%bias_sigma0)**2
313 :
314 44 : IF (calc_forces) THEN
315 : ! nabla(E_b) = kb*(sigma - sigma_0)*nabla(sigma)
316 : ! nabla(sigma) = 1/sigma * 1/n sum(dE_i* nabla(dE_i))
317 : ! nabla(dE_i) = nabla(E_i) - nabla(<E>)
318 44 : pref = nnp%bias_kb*(1.0_dp - nnp%bias_sigma0/sigma)
319 396 : DO i = 1, nnp%n_committee
320 : nnp%bias_forces(:, :) = nnp%bias_forces(:, :) + &
321 : (nnp%committee_energy(i) - avrg)* &
322 270732 : (nnp%committee_forces(:, :, i) - nnp%nnp_forces(:, :))
323 : END DO
324 44 : pref = pref/REAL(nnp%n_committee, dp)
325 33836 : nnp%bias_forces(:, :) = nnp%bias_forces(:, :)*pref
326 : END IF
327 : END IF
328 :
329 44 : CALL timestop(handle)
330 :
331 44 : END SUBROUTINE nnp_bias_sigma
332 :
333 : ! **************************************************************************************************
334 : !> \brief Print properties according to the requests in input file
335 : !> \param nnp ...
336 : !> \param print_section ...
337 : !> \date 2020-10-10
338 : !> \author Christoph Schran (christoph.schran@rub.de)
339 : ! **************************************************************************************************
340 640 : SUBROUTINE nnp_print(nnp, print_section)
341 : TYPE(nnp_type), INTENT(INOUT) :: nnp
342 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
343 :
344 : INTEGER :: unit_nr
345 : LOGICAL :: explicit, file_is_new
346 : TYPE(cp_logger_type), POINTER :: logger
347 : TYPE(section_vals_type), POINTER :: print_key
348 :
349 320 : NULLIFY (logger, print_key)
350 320 : logger => cp_get_default_logger()
351 :
352 320 : print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
353 320 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
354 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
355 310 : middle_name="nnp-energies", is_new_file=file_is_new)
356 310 : IF (unit_nr > 0) CALL nnp_print_energies(nnp, unit_nr, file_is_new)
357 310 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
358 : END IF
359 :
360 320 : print_key => section_vals_get_subs_vals(print_section, "FORCES")
361 320 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
362 286 : CALL nnp_print_forces(nnp, print_key)
363 : END IF
364 :
365 320 : print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
366 320 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
367 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
368 264 : middle_name="nnp-forces-std")
369 264 : IF (unit_nr > 0) CALL nnp_print_force_sigma(nnp, unit_nr)
370 264 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
371 : END IF
372 :
373 : ! Output structures with extrapolation warning on any processor
374 320 : CALL logger%para_env%sum(nnp%output_expol)
375 320 : IF (nnp%output_expol) THEN
376 22 : print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
377 22 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
378 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
379 22 : middle_name="nnp-extrapolation")
380 22 : IF (unit_nr > 0) CALL nnp_print_expol(nnp, unit_nr)
381 22 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
382 : END IF
383 : END IF
384 :
385 320 : print_key => section_vals_get_subs_vals(print_section, "SUM_FORCE")
386 :
387 : CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
388 320 : explicit=explicit)
389 320 : IF (explicit) THEN
390 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
391 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".dat", &
392 0 : middle_name="nnp-sumforce", is_new_file=file_is_new)
393 0 : IF (unit_nr > 0) CALL nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
394 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
395 : END IF
396 : END IF
397 :
398 320 : END SUBROUTINE nnp_print
399 :
400 : ! **************************************************************************************************
401 : !> \brief Print NNP energies and standard deviation sigma
402 : !> \param nnp ...
403 : !> \param unit_nr ...
404 : !> \param file_is_new ...
405 : !> \date 2020-10-10
406 : !> \author Christoph Schran (christoph.schran@rub.de)
407 : ! **************************************************************************************************
408 155 : SUBROUTINE nnp_print_energies(nnp, unit_nr, file_is_new)
409 : TYPE(nnp_type), INTENT(INOUT) :: nnp
410 : INTEGER, INTENT(IN) :: unit_nr
411 : LOGICAL, INTENT(IN) :: file_is_new
412 :
413 : CHARACTER(LEN=12) :: fmt_string
414 : INTEGER :: i
415 : REAL(KIND=dp) :: std
416 :
417 155 : IF (file_is_new) THEN
418 7 : WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
419 7 : WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
420 56 : DO i = 1, nnp%n_committee
421 56 : WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
422 : END DO
423 7 : WRITE (unit_nr, "(A)") ""
424 : END IF
425 :
426 155 : fmt_string = "(2X, F20.9)"
427 155 : WRITE (UNIT=fmt_string(5:6), FMT="(I2)") nnp%n_committee + 2
428 212590 : std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
429 155 : std = std/REAL(nnp%n_committee, dp)
430 155 : std = SQRT(std)
431 212590 : WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
432 :
433 155 : END SUBROUTINE nnp_print_energies
434 :
435 : ! **************************************************************************************************
436 : !> \brief Print nnp forces
437 : !> \param nnp ...
438 : !> \param print_key ...
439 : !> \date 2020-10-10
440 : !> \author Christoph Schran (christoph.schran@rub.de)
441 : ! **************************************************************************************************
442 286 : SUBROUTINE nnp_print_forces(nnp, print_key)
443 : TYPE(nnp_type), INTENT(INOUT) :: nnp
444 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_key
445 :
446 : CHARACTER(len=default_path_length) :: fmt_string, middle_name
447 : INTEGER :: i, j, unit_nr
448 : TYPE(cp_logger_type), POINTER :: logger
449 :
450 286 : NULLIFY (logger)
451 286 : logger => cp_get_default_logger()
452 :
453 2420 : DO i = 1, nnp%n_committee
454 2134 : WRITE (fmt_string, *) i
455 2134 : WRITE (middle_name, "(A,A)") "nnp-forces-", ADJUSTL(TRIM(fmt_string))
456 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
457 2134 : middle_name=TRIM(middle_name))
458 2134 : IF (unit_nr > 0) THEN
459 1067 : WRITE (unit_nr, *) nnp%num_atoms
460 1067 : WRITE (unit_nr, "(A,1X,A,A,F20.9)") "NNP forces [a.u.] of committee member", &
461 2134 : ADJUSTL(TRIM(fmt_string)), "energy [a.u.]=", nnp%committee_energy(i)
462 :
463 1067 : fmt_string = "(A4,1X,3F20.10)"
464 205931 : DO j = 1, nnp%num_atoms
465 205931 : WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(j)), nnp%committee_forces(:, j, i)
466 : END DO
467 : END IF
468 2420 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
469 : END DO
470 :
471 286 : END SUBROUTINE nnp_print_forces
472 :
473 : ! **************************************************************************************************
474 : !> \brief Print standard deviation sigma of NNP forces
475 : !> \param nnp ...
476 : !> \param unit_nr ...
477 : !> \date 2020-10-10
478 : !> \author Christoph Schran (christoph.schran@rub.de)
479 : ! **************************************************************************************************
480 132 : SUBROUTINE nnp_print_force_sigma(nnp, unit_nr)
481 : TYPE(nnp_type), INTENT(INOUT) :: nnp
482 : INTEGER, INTENT(IN) :: unit_nr
483 :
484 : INTEGER :: i, j
485 : REAL(KIND=dp), DIMENSION(3) :: var
486 :
487 132 : IF (unit_nr > 0) THEN
488 132 : WRITE (unit_nr, *) nnp%num_atoms
489 132 : WRITE (unit_nr, "(A,1X,A)") "NNP sigma of forces [a.u.]"
490 :
491 25476 : DO i = 1, nnp%num_atoms
492 25344 : var = 0.0_dp
493 228096 : DO j = 1, nnp%n_committee
494 836352 : var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
495 : END DO
496 101376 : var = var/REAL(nnp%n_committee, dp)
497 101376 : var = SQRT(var)
498 25476 : WRITE (unit_nr, "(A4,1X,3F20.10)") nnp%atoms(nnp%sort_inv(i)), var
499 : END DO
500 : END IF
501 :
502 132 : END SUBROUTINE nnp_print_force_sigma
503 :
504 : ! **************************************************************************************************
505 : !> \brief Print structures with extrapolation warning
506 : !> \param nnp ...
507 : !> \param unit_nr ...
508 : !> \date 2020-10-10
509 : !> \author Christoph Schran (christoph.schran@rub.de)
510 : ! **************************************************************************************************
511 11 : SUBROUTINE nnp_print_expol(nnp, unit_nr)
512 : TYPE(nnp_type), INTENT(INOUT) :: nnp
513 : INTEGER, INTENT(IN) :: unit_nr
514 :
515 : CHARACTER(len=default_path_length) :: fmt_string
516 : INTEGER :: i
517 : REAL(KIND=dp) :: mass, unit_conv
518 : REAL(KIND=dp), DIMENSION(3) :: com
519 :
520 11 : nnp%expol = nnp%expol + 1
521 11 : WRITE (unit_nr, *) nnp%num_atoms
522 11 : WRITE (unit_nr, "(A,1X,I6)") "NNP extrapolation point N =", nnp%expol
523 :
524 : ! move to COM of solute and wrap the box
525 : ! coord not needed afterwards, therefore manipulation ok
526 11 : com = 0.0_dp
527 11 : mass = 0.0_dp
528 44 : DO i = 1, nnp%num_atoms
529 33 : CALL get_ptable_info(nnp%atoms(i), amass=unit_conv)
530 132 : com(:) = com(:) + nnp%coord(:, i)*unit_conv
531 77 : mass = mass + unit_conv
532 : END DO
533 44 : com(:) = com(:)/mass
534 :
535 44 : DO i = 1, nnp%num_atoms
536 143 : nnp%coord(:, i) = nnp%coord(:, i) - com(:)
537 : END DO
538 :
539 : ! write out coordinates
540 11 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
541 11 : fmt_string = "(A4,1X,3F20.10)"
542 44 : DO i = 1, nnp%num_atoms
543 : WRITE (unit_nr, fmt_string) &
544 33 : nnp%atoms(nnp%sort_inv(i)), &
545 33 : nnp%coord(1, nnp%sort_inv(i))*unit_conv, &
546 33 : nnp%coord(2, nnp%sort_inv(i))*unit_conv, &
547 77 : nnp%coord(3, nnp%sort_inv(i))*unit_conv
548 : END DO
549 :
550 11 : END SUBROUTINE nnp_print_expol
551 :
552 : ! **************************************************************************************************
553 : !> \brief Print properties number according the requests in input file
554 : !> \param nnp ...
555 : !> \param print_section ...
556 : !> \date 2020-10-10
557 : !> \author Christoph Schran (christoph.schran@rub.de)
558 : ! **************************************************************************************************
559 44 : SUBROUTINE nnp_print_bias(nnp, print_section)
560 : TYPE(nnp_type), INTENT(INOUT) :: nnp
561 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
562 :
563 : INTEGER :: unit_nr
564 : LOGICAL :: file_is_new
565 : TYPE(cp_logger_type), POINTER :: logger
566 : TYPE(section_vals_type), POINTER :: print_key
567 :
568 44 : NULLIFY (logger, print_key)
569 44 : logger => cp_get_default_logger()
570 :
571 44 : print_key => section_vals_get_subs_vals(print_section, "BIAS_ENERGY")
572 44 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
573 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
574 44 : middle_name="nnp-bias-energy", is_new_file=file_is_new)
575 44 : IF (unit_nr > 0) CALL nnp_print_bias_energy(nnp, unit_nr, file_is_new)
576 44 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
577 : END IF
578 :
579 44 : print_key => section_vals_get_subs_vals(print_section, "BIAS_FORCES")
580 44 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
581 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
582 0 : middle_name="nnp-bias-forces")
583 0 : IF (unit_nr > 0) CALL nnp_print_bias_forces(nnp, unit_nr)
584 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
585 : END IF
586 :
587 44 : END SUBROUTINE nnp_print_bias
588 :
589 : ! **************************************************************************************************
590 : !> \brief Print NNP bias energies
591 : !> \param nnp ...
592 : !> \param unit_nr ...
593 : !> \param file_is_new ...
594 : !> \date 2020-10-10
595 : !> \author Christoph Schran (christoph.schran@rub.de)
596 : ! **************************************************************************************************
597 22 : SUBROUTINE nnp_print_bias_energy(nnp, unit_nr, file_is_new)
598 : TYPE(nnp_type), INTENT(INOUT) :: nnp
599 : INTEGER, INTENT(IN) :: unit_nr
600 : LOGICAL, INTENT(IN) :: file_is_new
601 :
602 : CHARACTER(len=default_path_length) :: fmt_string
603 : INTEGER :: i
604 :
605 22 : IF (file_is_new) THEN
606 1 : WRITE (unit_nr, "(A1)", ADVANCE='no') "#"
607 1 : WRITE (unit_nr, "(2(2X,A19))", ADVANCE='no') "Sigma [a.u.]", "Bias energy [a.u.]"
608 9 : DO i = 1, nnp%n_committee
609 9 : IF (nnp%bias_align) THEN
610 8 : WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "shifted E_NNP", i
611 : ELSE
612 0 : WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "E_NNP", i
613 : END IF
614 : END DO
615 1 : WRITE (unit_nr, "(A)") ""
616 :
617 : END IF
618 :
619 22 : WRITE (fmt_string, "(A,I3,A)") "(2X,", nnp%n_committee + 2, "(F20.9,1X))"
620 22 : WRITE (unit_nr, fmt_string) nnp%bias_sigma, nnp%bias_energy, nnp%committee_energy
621 :
622 22 : END SUBROUTINE nnp_print_bias_energy
623 :
624 : ! **************************************************************************************************
625 : !> \brief Print NNP bias forces
626 : !> \param nnp ...
627 : !> \param unit_nr ...
628 : !> \date 2020-10-10
629 : !> \author Christoph Schran (christoph.schran@rub.de)
630 : ! **************************************************************************************************
631 0 : SUBROUTINE nnp_print_bias_forces(nnp, unit_nr)
632 : TYPE(nnp_type), INTENT(INOUT) :: nnp
633 : INTEGER, INTENT(IN) :: unit_nr
634 :
635 : CHARACTER(len=default_path_length) :: fmt_string
636 : INTEGER :: i
637 :
638 0 : WRITE (unit_nr, *) nnp%num_atoms
639 0 : WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
640 :
641 0 : fmt_string = "(A4,1X,3F20.10)"
642 0 : DO i = 1, nnp%num_atoms
643 0 : WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(i)), nnp%bias_forces(:, i)
644 : END DO
645 :
646 0 : END SUBROUTINE nnp_print_bias_forces
647 :
648 : ! **************************************************************************************************
649 : !> \brief Print NNP summed forces
650 : !> \param nnp ...
651 : !> \param print_section ...
652 : !> \param unit_nr ...
653 : !> \param file_is_new ...
654 : !> \date 2020-10-10
655 : !> \author Christoph Schran (christoph.schran@rub.de)
656 : ! **************************************************************************************************
657 0 : SUBROUTINE nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
658 : TYPE(nnp_type), INTENT(INOUT) :: nnp
659 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
660 : INTEGER, INTENT(IN) :: unit_nr
661 : LOGICAL, INTENT(IN) :: file_is_new
662 :
663 : CHARACTER(len=default_path_length) :: fmt_string
664 : CHARACTER(LEN=default_string_length), &
665 0 : DIMENSION(:), POINTER :: atomlist
666 : INTEGER :: i, ig, j, n
667 : REAL(KIND=dp), DIMENSION(3) :: rvec
668 :
669 0 : NULLIFY (atomlist)
670 0 : IF (file_is_new) THEN
671 0 : WRITE (unit_nr, "(A)") "# Summed forces [a.u.]"
672 : END IF
673 :
674 0 : rvec = 0.0_dp
675 :
676 : ! get atoms to sum over:
677 : CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
678 0 : c_vals=atomlist)
679 0 : IF (ASSOCIATED(atomlist)) THEN
680 0 : n = SIZE(atomlist)
681 0 : DO i = 1, nnp%num_atoms
682 0 : DO j = 1, n
683 0 : ig = nnp%sort_inv(i)
684 0 : IF (TRIM(ADJUSTL(atomlist(j))) == TRIM(ADJUSTL(nnp%atoms(ig)))) THEN
685 0 : rvec(:) = rvec(:) + nnp%nnp_forces(:, i)
686 : END IF
687 : END DO
688 : END DO
689 : END IF
690 :
691 0 : fmt_string = "(3(F20.10,1X))"
692 0 : WRITE (unit_nr, fmt_string) rvec
693 :
694 0 : END SUBROUTINE nnp_print_sumforces
695 :
696 : ! **************************************************************************************************
697 : !> \brief Scatter the per-neighbour dG/dr arrays held in the nnp_neighbor_workspace
698 : !> into a per-atom Cartesian force array. Three contributions per central
699 : !> atom i:
700 : !> 1) self: self_dGdr(:, m) -> force_xyz(:, i)
701 : !> 2) radial: walk neighbor%rad(s)%ind(j), read dGdr_rad(s)%data(:, sf, j)
702 : !> 3) angular: walk neighbor%ang1(s)%ind / ang2(s)%ind, read dGdr_ang_jj/kk
703 : !> No global (n_sf, num_atoms) slab is touched. Reused by the helium-NNP
704 : !> coupling in helium_interactions.
705 : !>
706 : !> Precondition: nnp_calc_acsf(nnp, i, calc_forces=.TRUE.[, stress]) must
707 : !> have run for this same atom i immediately before; the per-element
708 : !> workspace it fills is overwritten on every ACSF call.
709 : !> \param nnp ...
710 : !> \param ind central-atom element index
711 : !> \param i central atom index (absolute)
712 : !> \param denergydsym dE/dG_k for k = 1..n_input_nodes(ind)
713 : !> \param force_xyz (3, num_atoms) destination -- accumulated, NOT overwritten
714 : !> \author Dhruv Sharma (ds2173@cam.ac.uk)
715 : ! **************************************************************************************************
716 130120 : SUBROUTINE nnp_scatter_dgdr_to_forces(nnp, ind, i, denergydsym, force_xyz)
717 :
718 : TYPE(nnp_type), INTENT(IN) :: nnp
719 : INTEGER, INTENT(IN) :: ind, i
720 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: denergydsym
721 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: force_xyz
722 :
723 : INTEGER :: j, k, k_atom, m, n_ang1_s, n_ang2_s, &
724 : n_input_nodes, n_symf_s, off, s, sf
725 : LOGICAL :: homo_grp
726 : REAL(KIND=dp) :: de
727 :
728 : ASSOCIATE (workspace => nnp%neighbor_interface_state%workspace(ind), &
729 : neighbor => nnp%neighbor_interface_state%workspace(ind)%neighbor, &
730 : self_dGdr => nnp%neighbor_interface_state%workspace(ind)%self_dGdr)
731 :
732 130120 : n_input_nodes = workspace%n_input_nodes
733 :
734 : ! Self contribution
735 3725480 : DO j = 1, n_input_nodes
736 3595360 : de = denergydsym(j)
737 3595360 : force_xyz(1, i) = force_xyz(1, i) - de*self_dGdr(1, j)
738 3595360 : force_xyz(2, i) = force_xyz(2, i) - de*self_dGdr(2, j)
739 3725480 : force_xyz(3, i) = force_xyz(3, i) - de*self_dGdr(3, j)
740 : END DO
741 :
742 : ! Radial neighbours
743 393560 : DO s = 1, nnp%rad(ind)%n_symfgrp
744 263440 : n_symf_s = nnp%rad(ind)%symfgrp(s)%n_symf
745 130120 : ASSOCIATE (rad_buf => workspace%dGdr_rad(s)%data)
746 13413642 : DO j = 1, neighbor%n_rad(s)
747 13150202 : k_atom = neighbor%rad(s)%ind(j)
748 118591990 : DO sf = 1, n_symf_s
749 105178348 : m = nnp%rad(ind)%symfgrp(s)%symf(sf)
750 105178348 : de = denergydsym(m)
751 105178348 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*rad_buf(1, sf, j)
752 105178348 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*rad_buf(2, sf, j)
753 118328550 : force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*rad_buf(3, sf, j)
754 : END DO
755 : END DO
756 : END ASSOCIATE
757 : END DO
758 :
759 : ! Angular neighbours
760 130120 : off = nnp%n_rad(ind)
761 561720 : DO s = 1, nnp%ang(ind)%n_symfgrp
762 301480 : n_symf_s = nnp%ang(ind)%symfgrp(s)%n_symf
763 301480 : n_ang1_s = neighbor%n_ang1(s)
764 301480 : homo_grp = (nnp%ang(ind)%symfgrp(s)%ele(1) == nnp%ang(ind)%symfgrp(s)%ele(2))
765 301480 : IF (homo_grp) THEN
766 : n_ang2_s = n_ang1_s ! kk buffer is also indexed in ang1(s)%ind
767 : ELSE
768 133320 : n_ang2_s = neighbor%n_ang2(s)
769 : END IF
770 :
771 : ASSOCIATE (jj_buf => workspace%dGdr_ang_jj(s)%data, &
772 130120 : kk_buf => workspace%dGdr_ang_kk(s)%data)
773 : ! jj-side
774 16399965 : DO j = 1, n_ang1_s
775 16098485 : k_atom = neighbor%ang1(s)%ind(j)
776 104160489 : DO sf = 1, n_symf_s
777 87760524 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
778 87760524 : de = denergydsym(m)
779 87760524 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*jj_buf(1, sf, j)
780 87760524 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*jj_buf(2, sf, j)
781 103859009 : force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*jj_buf(3, sf, j)
782 : END DO
783 : END DO
784 :
785 : ! kk-side: homo reads ang1(s)%ind, hetero reads ang2(s)%ind
786 602960 : IF (homo_grp) THEN
787 7483470 : DO k = 1, n_ang2_s
788 7315310 : k_atom = neighbor%ang1(s)%ind(k)
789 42632218 : DO sf = 1, n_symf_s
790 35148748 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
791 35148748 : de = denergydsym(m)
792 35148748 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
793 35148748 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
794 42464058 : force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
795 : END DO
796 : END DO
797 : ELSE
798 4503547 : DO k = 1, n_ang2_s
799 4370227 : k_atom = neighbor%ang2(s)%ind(k)
800 30817724 : DO sf = 1, n_symf_s
801 26314177 : m = off + nnp%ang(ind)%symfgrp(s)%symf(sf)
802 26314177 : de = denergydsym(m)
803 26314177 : force_xyz(1, k_atom) = force_xyz(1, k_atom) - de*kk_buf(1, sf, k)
804 26314177 : force_xyz(2, k_atom) = force_xyz(2, k_atom) - de*kk_buf(2, sf, k)
805 30684404 : force_xyz(3, k_atom) = force_xyz(3, k_atom) - de*kk_buf(3, sf, k)
806 : END DO
807 : END DO
808 : END IF
809 : END ASSOCIATE
810 : END DO
811 :
812 : END ASSOCIATE
813 :
814 130120 : END SUBROUTINE nnp_scatter_dgdr_to_forces
815 :
816 : END MODULE nnp_force
|
