Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 : !> \date 2020-10-10
12 : ! **************************************************************************************************
13 : MODULE nnp_force
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE cp_log_handling, ONLY: cp_get_default_logger,&
17 : cp_logger_type
18 : USE cp_output_handling, ONLY: cp_p_file,&
19 : cp_print_key_finished_output,&
20 : cp_print_key_should_output,&
21 : cp_print_key_unit_nr
22 : USE cp_subsys_types, ONLY: cp_subsys_get,&
23 : cp_subsys_type
24 : USE cp_units, ONLY: cp_unit_from_cp2k
25 : USE distribution_1d_types, ONLY: distribution_1d_type
26 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: default_path_length,&
30 : default_string_length,&
31 : dp
32 : USE nnp_acsf, ONLY: nnp_calc_acsf
33 : USE nnp_environment_types, ONLY: nnp_env_get,&
34 : nnp_type
35 : USE nnp_model, ONLY: nnp_gradients,&
36 : nnp_predict
37 : USE particle_types, ONLY: particle_type
38 : USE periodic_table, ONLY: get_ptable_info
39 : USE physcon, ONLY: angstrom
40 : USE virial_types, ONLY: virial_type
41 : #include "./base/base_uses.f90"
42 :
43 : IMPLICIT NONE
44 :
45 : PRIVATE
46 :
47 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_force'
49 :
50 : PUBLIC :: nnp_calc_energy_force
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief Calculate the energy and force for a given configuration with the NNP
56 : !> \param nnp ...
57 : !> \param calc_forces ...
58 : !> \date 2020-10-10
59 : !> \author Christoph Schran (christoph.schran@rub.de)
60 : ! **************************************************************************************************
61 308 : SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
62 : TYPE(nnp_type), INTENT(INOUT), POINTER :: nnp
63 : LOGICAL, INTENT(IN) :: calc_forces
64 :
65 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_calc_energy_force'
66 :
67 : INTEGER :: handle, i, i_com, ig, ind, istart, j, k, &
68 : m, mecalc
69 308 : INTEGER, ALLOCATABLE, DIMENSION(:) :: allcalc
70 : LOGICAL :: calc_stress
71 308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: denergydsym
72 308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dsymdxyz, stress
73 308 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
74 : TYPE(cp_logger_type), POINTER :: logger
75 : TYPE(cp_subsys_type), POINTER :: subsys
76 : TYPE(distribution_1d_type), POINTER :: local_particles
77 308 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
78 : TYPE(section_vals_type), POINTER :: print_section
79 : TYPE(virial_type), POINTER :: virial
80 :
81 308 : CALL timeset(routineN, handle)
82 :
83 308 : NULLIFY (particle_set, logger, local_particles, subsys, &
84 308 : atomic_kind_set)
85 308 : logger => cp_get_default_logger()
86 :
87 308 : CPASSERT(ASSOCIATED(nnp))
88 : CALL nnp_env_get(nnp_env=nnp, particle_set=particle_set, &
89 : subsys=subsys, local_particles=local_particles, &
90 308 : atomic_kind_set=atomic_kind_set)
91 :
92 : CALL cp_subsys_get(subsys, &
93 308 : virial=virial)
94 :
95 308 : calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
96 22 : IF (calc_stress .AND. .NOT. calc_forces) CPABORT('Stress cannot be calculated without forces')
97 :
98 : ! Initialize energy and gradient:
99 412874 : nnp%atomic_energy(:, :) = 0.0_dp
100 696234 : IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
101 696234 : IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
102 2596 : IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
103 :
104 : !fill coord array
105 308 : ig = 1
106 924 : DO i = 1, nnp%n_ele
107 110880 : DO j = 1, nnp%num_atoms
108 110572 : IF (nnp%ele(i) == particle_set(j)%atomic_kind%element_symbol) THEN
109 219912 : DO m = 1, 3
110 219912 : nnp%coord(m, ig) = particle_set(j)%r(m)
111 : END DO
112 54978 : nnp%atoms(ig) = nnp%ele(i)
113 54978 : CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
114 54978 : nnp%ele_ind(ig) = i
115 54978 : nnp%sort(ig) = j
116 54978 : nnp%sort_inv(j) = ig
117 54978 : ig = ig + 1
118 : END IF
119 : END DO
120 : END DO
121 :
122 : ! parallization:
123 : mecalc = nnp%num_atoms/logger%para_env%num_pe + &
124 : MIN(MOD(nnp%num_atoms, logger%para_env%num_pe)/ &
125 308 : (logger%para_env%mepos + 1), 1)
126 924 : ALLOCATE (allcalc(logger%para_env%num_pe))
127 924 : allcalc(:) = 0
128 308 : CALL logger%para_env%allgather(mecalc, allcalc)
129 308 : istart = 1
130 462 : DO i = 2, logger%para_env%mepos + 1
131 462 : istart = istart + allcalc(i - 1)
132 : END DO
133 :
134 : ! reset extrapolation status
135 308 : nnp%output_expol = .FALSE.
136 :
137 : ! calc atomic contribution to energy and force
138 27797 : DO i = istart, istart + mecalc - 1
139 :
140 : ! determine index of atom type and offset
141 27489 : ind = nnp%ele_ind(i)
142 :
143 : ! reset input nodes of ele(ind):
144 797181 : nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
145 :
146 : ! compute sym fnct values
147 27489 : IF (calc_forces) THEN
148 : !reset input grads of ele(ind):
149 368445 : nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
150 50820 : ALLOCATE (dsymdxyz(3, nnp%arc(ind)%n_nodes(1), nnp%num_atoms))
151 274955604 : dsymdxyz(:, :, :) = 0.0_dp
152 12705 : IF (calc_stress) THEN
153 6336 : ALLOCATE (stress(3, 3, nnp%arc(ind)%n_nodes(1)))
154 770880 : stress(:, :, :) = 0.0_dp
155 2112 : CALL nnp_calc_acsf(nnp, i, dsymdxyz, stress)
156 : ELSE
157 10593 : CALL nnp_calc_acsf(nnp, i, dsymdxyz)
158 : END IF
159 : ELSE
160 14784 : CALL nnp_calc_acsf(nnp, i)
161 : END IF
162 :
163 232617 : DO i_com = 1, nnp%n_committee
164 : ! predict energy
165 205128 : CALL nnp_predict(nnp%arc(ind), nnp, i_com)
166 : nnp%atomic_energy(i, i_com) = nnp%arc(ind)%layer(nnp%n_layer)%node(1) + &
167 205128 : nnp%atom_energies(ind)
168 :
169 : ! predict forces
170 232617 : IF (calc_forces) THEN
171 260568 : ALLOCATE (denergydsym(nnp%arc(ind)%n_nodes(1)))
172 2518824 : denergydsym(:) = 0.0_dp
173 86856 : CALL nnp_gradients(nnp%arc(ind), nnp, i_com, denergydsym)
174 2518824 : DO j = 1, nnp%arc(ind)%n_nodes(1)
175 467972736 : DO k = 1, nnp%num_atoms
176 1864595040 : DO m = 1, 3
177 1862163072 : nnp%myforce(m, k, i_com) = nnp%myforce(m, k, i_com) - denergydsym(j)*dsymdxyz(m, j, k)
178 : END DO
179 : END DO
180 2518824 : IF (calc_stress) THEN
181 : nnp%committee_stress(:, :, i_com) = nnp%committee_stress(:, :, i_com) - &
182 6150144 : denergydsym(j)*stress(:, :, j)
183 : END IF
184 : END DO
185 86856 : DEALLOCATE (denergydsym)
186 : END IF
187 : END DO
188 :
189 : !deallocate memory
190 27797 : IF (calc_forces) THEN
191 12705 : DEALLOCATE (dsymdxyz)
192 12705 : IF (calc_stress) THEN
193 2112 : DEALLOCATE (stress)
194 : END IF
195 : END IF
196 :
197 : END DO ! loop over num_atoms
198 :
199 : ! calculate energy:
200 308 : CALL logger%para_env%sum(nnp%atomic_energy(:, :))
201 412874 : nnp%committee_energy(:) = SUM(nnp%atomic_energy, 1)
202 2618 : nnp%nnp_potential_energy = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
203 :
204 308 : IF (calc_forces) THEN
205 : ! bring myforce to force array
206 25564 : DO j = 1, nnp%num_atoms
207 101794 : DO k = 1, 3
208 622776 : nnp%committee_forces(k, (nnp%sort(j)), :) = nnp%myforce(k, j, :)
209 : END DO
210 : END DO
211 154 : CALL logger%para_env%sum(nnp%committee_forces)
212 622930 : nnp%nnp_forces(:, :) = SUM(nnp%committee_forces, DIM=3)/REAL(nnp%n_committee, dp)
213 25564 : DO j = 1, nnp%num_atoms
214 178024 : particle_set(j)%f(:) = nnp%nnp_forces(:, j)
215 : END DO
216 : END IF
217 :
218 308 : IF (calc_stress) THEN
219 22 : CALL logger%para_env%sum(nnp%committee_stress)
220 2134 : virial%pv_virial = SUM(nnp%committee_stress, DIM=3)/REAL(nnp%n_committee, dp)
221 : END IF
222 :
223 : ! Bias the standard deviation of committee disagreement
224 308 : IF (nnp%bias) THEN
225 44 : CALL nnp_bias_sigma(nnp, calc_forces)
226 44 : nnp%nnp_potential_energy = nnp%nnp_potential_energy + nnp%bias_energy
227 44 : IF (calc_forces) THEN
228 8492 : DO j = 1, nnp%num_atoms
229 59180 : particle_set(j)%f(:) = particle_set(j)%f(:) + nnp%bias_forces(:, j)
230 : END DO
231 : END IF
232 : ! print properties if requested
233 44 : print_section => section_vals_get_subs_vals(nnp%nnp_input, "BIAS%PRINT")
234 44 : CALL nnp_print_bias(nnp, print_section)
235 : END IF
236 :
237 : ! print properties if requested
238 308 : print_section => section_vals_get_subs_vals(nnp%nnp_input, "PRINT")
239 308 : CALL nnp_print(nnp, print_section)
240 :
241 308 : DEALLOCATE (allcalc)
242 :
243 308 : CALL timestop(handle)
244 :
245 616 : END SUBROUTINE nnp_calc_energy_force
246 :
247 : ! **************************************************************************************************
248 : !> \brief Calculate bias potential and force based on standard deviation of committee disagreement
249 : !> \param nnp ...
250 : !> \param calc_forces ...
251 : !> \date 2020-10-10
252 : !> \author Christoph Schran (christoph.schran@rub.de)
253 : ! **************************************************************************************************
254 44 : SUBROUTINE nnp_bias_sigma(nnp, calc_forces)
255 : TYPE(nnp_type), INTENT(INOUT) :: nnp
256 : LOGICAL, INTENT(IN) :: calc_forces
257 :
258 : CHARACTER(len=*), PARAMETER :: routineN = 'nnp_bias_sigma'
259 :
260 : INTEGER :: handle, i
261 : REAL(KIND=dp) :: avrg, pref, sigma
262 :
263 44 : CALL timeset(routineN, handle)
264 :
265 : ! init
266 44 : sigma = 0.0_dp
267 44 : nnp%bias_energy = 0.0_dp
268 33836 : IF (calc_forces) nnp%bias_forces = 0.0_dp
269 :
270 : ! Subtract reference energy of each committee member, if requested
271 44 : IF (nnp%bias_align) THEN
272 : ! committee energy not used afterward, therefore overwritten
273 396 : nnp%committee_energy(:) = nnp%committee_energy(:) - nnp%bias_e_avrg(:)
274 : END IF
275 :
276 : ! <E> = 1/n sum(E_i)
277 : ! sigma = sqrt(1/n sum((E_i - <E>)**2))
278 : ! = sqrt(1/n sum(dE_i**2))
279 396 : avrg = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
280 396 : DO i = 1, nnp%n_committee
281 396 : sigma = sigma + (nnp%committee_energy(i) - avrg)**2
282 : END DO
283 44 : sigma = SQRT(sigma/REAL(nnp%n_committee, dp))
284 44 : nnp%bias_sigma = sigma
285 :
286 44 : IF (sigma > nnp%bias_sigma0) THEN
287 : ! E_b = 0.5 * kb * (sigma - sigma_0)**2
288 44 : nnp%bias_energy = 0.5_dp*nnp%bias_kb*(sigma - nnp%bias_sigma0)**2
289 :
290 44 : IF (calc_forces) THEN
291 : ! nabla(E_b) = kb*(sigma - sigma_0)*nabla(sigma)
292 : ! nabla(sigma) = 1/sigma * 1/n sum(dE_i* nabla(dE_i))
293 : ! nabla(dE_i) = nabla(E_i) - nabla(<E>)
294 44 : pref = nnp%bias_kb*(1.0_dp - nnp%bias_sigma0/sigma)
295 396 : DO i = 1, nnp%n_committee
296 : nnp%bias_forces(:, :) = nnp%bias_forces(:, :) + &
297 : (nnp%committee_energy(i) - avrg)* &
298 270732 : (nnp%committee_forces(:, :, i) - nnp%nnp_forces(:, :))
299 : END DO
300 44 : pref = pref/REAL(nnp%n_committee, dp)
301 33836 : nnp%bias_forces(:, :) = nnp%bias_forces(:, :)*pref
302 : END IF
303 : END IF
304 :
305 44 : CALL timestop(handle)
306 :
307 44 : END SUBROUTINE nnp_bias_sigma
308 :
309 : ! **************************************************************************************************
310 : !> \brief Print properties according to the requests in input file
311 : !> \param nnp ...
312 : !> \param print_section ...
313 : !> \date 2020-10-10
314 : !> \author Christoph Schran (christoph.schran@rub.de)
315 : ! **************************************************************************************************
316 616 : SUBROUTINE nnp_print(nnp, print_section)
317 : TYPE(nnp_type), INTENT(INOUT) :: nnp
318 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
319 :
320 : INTEGER :: unit_nr
321 : LOGICAL :: explicit, file_is_new
322 : TYPE(cp_logger_type), POINTER :: logger
323 : TYPE(section_vals_type), POINTER :: print_key
324 :
325 308 : NULLIFY (logger, print_key)
326 308 : logger => cp_get_default_logger()
327 :
328 308 : print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
329 308 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
330 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
331 308 : middle_name="nnp-energies", is_new_file=file_is_new)
332 308 : IF (unit_nr > 0) CALL nnp_print_energies(nnp, unit_nr, file_is_new)
333 308 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
334 : END IF
335 :
336 308 : print_key => section_vals_get_subs_vals(print_section, "FORCES")
337 308 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
338 286 : IF (unit_nr > 0) CALL nnp_print_forces(nnp, print_key)
339 : END IF
340 :
341 308 : print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
342 308 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
343 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
344 264 : middle_name="nnp-forces-std")
345 264 : IF (unit_nr > 0) CALL nnp_print_force_sigma(nnp, unit_nr)
346 264 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
347 : END IF
348 :
349 : ! Output structures with extrapolation warning on any processor
350 308 : CALL logger%para_env%sum(nnp%output_expol)
351 308 : IF (nnp%output_expol) THEN
352 22 : print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
353 22 : 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=".xyz", &
355 22 : middle_name="nnp-extrapolation")
356 22 : IF (unit_nr > 0) CALL nnp_print_expol(nnp, unit_nr)
357 22 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
358 : END IF
359 : END IF
360 :
361 308 : print_key => section_vals_get_subs_vals(print_section, "SUM_FORCE")
362 :
363 : CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
364 308 : explicit=explicit)
365 308 : IF (explicit) THEN
366 0 : 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=".dat", &
368 0 : middle_name="nnp-sumforce", is_new_file=file_is_new)
369 0 : IF (unit_nr > 0) CALL nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
370 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
371 : END IF
372 : END IF
373 :
374 308 : END SUBROUTINE nnp_print
375 :
376 : ! **************************************************************************************************
377 : !> \brief Print NNP energies and standard deviation sigma
378 : !> \param nnp ...
379 : !> \param unit_nr ...
380 : !> \param file_is_new ...
381 : !> \date 2020-10-10
382 : !> \author Christoph Schran (christoph.schran@rub.de)
383 : ! **************************************************************************************************
384 154 : SUBROUTINE nnp_print_energies(nnp, unit_nr, file_is_new)
385 : TYPE(nnp_type), INTENT(INOUT) :: nnp
386 : INTEGER, INTENT(IN) :: unit_nr
387 : LOGICAL, INTENT(IN) :: file_is_new
388 :
389 : CHARACTER(len=default_path_length) :: fmt_string
390 : INTEGER :: i
391 : REAL(KIND=dp) :: std
392 :
393 154 : IF (file_is_new) THEN
394 6 : WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
395 6 : WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
396 47 : DO i = 1, nnp%n_committee
397 47 : WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
398 : END DO
399 6 : WRITE (unit_nr, "(A)") ""
400 : END IF
401 :
402 154 : fmt_string = "(2X,2(F20.9))"
403 154 : WRITE (fmt_string, "(A,I3,A)") "(2X", nnp%n_committee + 2, "(F20.9))"
404 206437 : std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
405 154 : std = std/REAL(nnp%n_committee, dp)
406 154 : std = SQRT(std)
407 206437 : WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
408 :
409 154 : END SUBROUTINE nnp_print_energies
410 :
411 : ! **************************************************************************************************
412 : !> \brief Print nnp forces
413 : !> \param nnp ...
414 : !> \param print_key ...
415 : !> \date 2020-10-10
416 : !> \author Christoph Schran (christoph.schran@rub.de)
417 : ! **************************************************************************************************
418 0 : SUBROUTINE nnp_print_forces(nnp, print_key)
419 : TYPE(nnp_type), INTENT(INOUT) :: nnp
420 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_key
421 :
422 : CHARACTER(len=default_path_length) :: fmt_string, middle_name
423 : INTEGER :: i, j, unit_nr
424 : TYPE(cp_logger_type), POINTER :: logger
425 :
426 0 : NULLIFY (logger)
427 0 : logger => cp_get_default_logger()
428 :
429 : ! TODO use write_particle_coordinates from particle_methods.F
430 0 : DO i = 1, nnp%n_committee
431 0 : WRITE (fmt_string, *) i
432 0 : WRITE (middle_name, "(A,A)") "nnp-forces-", ADJUSTL(TRIM(fmt_string))
433 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
434 0 : middle_name=TRIM(middle_name))
435 0 : IF (unit_nr > 0) THEN
436 0 : WRITE (unit_nr, *) nnp%num_atoms
437 0 : WRITE (unit_nr, "(A,1X,A,A,F20.9)") "NNP forces [a.u.] of committee member", &
438 0 : ADJUSTL(TRIM(fmt_string)), "energy [a.u.]=", nnp%committee_energy(i)
439 :
440 0 : fmt_string = "(A4,1X,3F20.10)"
441 0 : DO j = 1, nnp%num_atoms
442 0 : WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(j)), nnp%committee_forces(:, j, i)
443 : END DO
444 : END IF
445 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
446 : END DO
447 :
448 0 : END SUBROUTINE nnp_print_forces
449 :
450 : ! **************************************************************************************************
451 : !> \brief Print standard deviation sigma of NNP forces
452 : !> \param nnp ...
453 : !> \param unit_nr ...
454 : !> \date 2020-10-10
455 : !> \author Christoph Schran (christoph.schran@rub.de)
456 : ! **************************************************************************************************
457 132 : SUBROUTINE nnp_print_force_sigma(nnp, unit_nr)
458 : TYPE(nnp_type), INTENT(INOUT) :: nnp
459 : INTEGER, INTENT(IN) :: unit_nr
460 :
461 : INTEGER :: i, j
462 : REAL(KIND=dp), DIMENSION(3) :: var
463 :
464 132 : IF (unit_nr > 0) THEN
465 132 : WRITE (unit_nr, *) nnp%num_atoms
466 132 : WRITE (unit_nr, "(A,1X,A)") "NNP sigma of forces [a.u.]"
467 :
468 25476 : DO i = 1, nnp%num_atoms
469 25344 : var = 0.0_dp
470 228096 : DO j = 1, nnp%n_committee
471 836352 : var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
472 : END DO
473 101376 : var = var/REAL(nnp%n_committee, dp)
474 101376 : var = SQRT(var)
475 25476 : WRITE (unit_nr, "(A4,1X,3F20.10)") nnp%atoms(nnp%sort_inv(i)), var
476 : END DO
477 : END IF
478 :
479 132 : END SUBROUTINE nnp_print_force_sigma
480 :
481 : ! **************************************************************************************************
482 : !> \brief Print structures with extrapolation warning
483 : !> \param nnp ...
484 : !> \param unit_nr ...
485 : !> \date 2020-10-10
486 : !> \author Christoph Schran (christoph.schran@rub.de)
487 : ! **************************************************************************************************
488 11 : SUBROUTINE nnp_print_expol(nnp, unit_nr)
489 : TYPE(nnp_type), INTENT(INOUT) :: nnp
490 : INTEGER, INTENT(IN) :: unit_nr
491 :
492 : CHARACTER(len=default_path_length) :: fmt_string
493 : INTEGER :: i
494 : REAL(KIND=dp) :: mass, unit_conv
495 : REAL(KIND=dp), DIMENSION(3) :: com
496 :
497 11 : nnp%expol = nnp%expol + 1
498 11 : WRITE (unit_nr, *) nnp%num_atoms
499 11 : WRITE (unit_nr, "(A,1X,I6)") "NNP extrapolation point N =", nnp%expol
500 :
501 : ! move to COM of solute and wrap the box
502 : ! coord not needed afterwards, therefore manipulation ok
503 11 : com = 0.0_dp
504 11 : mass = 0.0_dp
505 44 : DO i = 1, nnp%num_atoms
506 33 : CALL get_ptable_info(nnp%atoms(i), amass=unit_conv)
507 132 : com(:) = com(:) + nnp%coord(:, i)*unit_conv
508 77 : mass = mass + unit_conv
509 : END DO
510 44 : com(:) = com(:)/mass
511 :
512 44 : DO i = 1, nnp%num_atoms
513 143 : nnp%coord(:, i) = nnp%coord(:, i) - com(:)
514 : END DO
515 :
516 : ! write out coordinates
517 11 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
518 11 : fmt_string = "(A4,1X,3F20.10)"
519 44 : DO i = 1, nnp%num_atoms
520 : WRITE (unit_nr, fmt_string) &
521 33 : nnp%atoms(nnp%sort_inv(i)), &
522 33 : nnp%coord(1, nnp%sort_inv(i))*unit_conv, &
523 33 : nnp%coord(2, nnp%sort_inv(i))*unit_conv, &
524 77 : nnp%coord(3, nnp%sort_inv(i))*unit_conv
525 : END DO
526 :
527 11 : END SUBROUTINE nnp_print_expol
528 :
529 : ! **************************************************************************************************
530 : !> \brief Print properties number according the requests in input file
531 : !> \param nnp ...
532 : !> \param print_section ...
533 : !> \date 2020-10-10
534 : !> \author Christoph Schran (christoph.schran@rub.de)
535 : ! **************************************************************************************************
536 44 : SUBROUTINE nnp_print_bias(nnp, print_section)
537 : TYPE(nnp_type), INTENT(INOUT) :: nnp
538 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
539 :
540 : INTEGER :: unit_nr
541 : LOGICAL :: file_is_new
542 : TYPE(cp_logger_type), POINTER :: logger
543 : TYPE(section_vals_type), POINTER :: print_key
544 :
545 44 : NULLIFY (logger, print_key)
546 44 : logger => cp_get_default_logger()
547 :
548 44 : print_key => section_vals_get_subs_vals(print_section, "BIAS_ENERGY")
549 44 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
550 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
551 44 : middle_name="nnp-bias-energy", is_new_file=file_is_new)
552 44 : IF (unit_nr > 0) CALL nnp_print_bias_energy(nnp, unit_nr, file_is_new)
553 44 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
554 : END IF
555 :
556 44 : print_key => section_vals_get_subs_vals(print_section, "BIAS_FORCES")
557 44 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
558 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
559 0 : middle_name="nnp-bias-forces")
560 0 : IF (unit_nr > 0) CALL nnp_print_bias_forces(nnp, unit_nr)
561 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
562 : END IF
563 :
564 44 : END SUBROUTINE nnp_print_bias
565 :
566 : ! **************************************************************************************************
567 : !> \brief Print NNP bias energies
568 : !> \param nnp ...
569 : !> \param unit_nr ...
570 : !> \param file_is_new ...
571 : !> \date 2020-10-10
572 : !> \author Christoph Schran (christoph.schran@rub.de)
573 : ! **************************************************************************************************
574 22 : SUBROUTINE nnp_print_bias_energy(nnp, unit_nr, file_is_new)
575 : TYPE(nnp_type), INTENT(INOUT) :: nnp
576 : INTEGER, INTENT(IN) :: unit_nr
577 : LOGICAL, INTENT(IN) :: file_is_new
578 :
579 : CHARACTER(len=default_path_length) :: fmt_string
580 : INTEGER :: i
581 :
582 22 : IF (file_is_new) THEN
583 1 : WRITE (unit_nr, "(A1)", ADVANCE='no') "#"
584 1 : WRITE (unit_nr, "(2(2X,A19))", ADVANCE='no') "Sigma [a.u.]", "Bias energy [a.u.]"
585 9 : DO i = 1, nnp%n_committee
586 9 : IF (nnp%bias_align) THEN
587 8 : WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "shifted E_NNP", i
588 : ELSE
589 0 : WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "E_NNP", i
590 : END IF
591 : END DO
592 1 : WRITE (unit_nr, "(A)") ""
593 :
594 : END IF
595 :
596 22 : WRITE (fmt_string, "(A,I3,A)") "(2X,", nnp%n_committee + 2, "(F20.9,1X))"
597 22 : WRITE (unit_nr, fmt_string) nnp%bias_sigma, nnp%bias_energy, nnp%committee_energy
598 :
599 22 : END SUBROUTINE nnp_print_bias_energy
600 :
601 : ! **************************************************************************************************
602 : !> \brief Print NNP bias forces
603 : !> \param nnp ...
604 : !> \param unit_nr ...
605 : !> \date 2020-10-10
606 : !> \author Christoph Schran (christoph.schran@rub.de)
607 : ! **************************************************************************************************
608 0 : SUBROUTINE nnp_print_bias_forces(nnp, unit_nr)
609 : TYPE(nnp_type), INTENT(INOUT) :: nnp
610 : INTEGER, INTENT(IN) :: unit_nr
611 :
612 : CHARACTER(len=default_path_length) :: fmt_string
613 : INTEGER :: i
614 :
615 : ! TODO use write_particle_coordinates from particle_methods.F
616 0 : WRITE (unit_nr, *) nnp%num_atoms
617 0 : WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
618 :
619 0 : fmt_string = "(A4,1X,3F20.10)"
620 0 : DO i = 1, nnp%num_atoms
621 0 : WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(i)), nnp%bias_forces(:, i)
622 : END DO
623 :
624 0 : END SUBROUTINE nnp_print_bias_forces
625 :
626 : ! **************************************************************************************************
627 : !> \brief Print NNP summed forces
628 : !> \param nnp ...
629 : !> \param print_section ...
630 : !> \param unit_nr ...
631 : !> \param file_is_new ...
632 : !> \date 2020-10-10
633 : !> \author Christoph Schran (christoph.schran@rub.de)
634 : ! **************************************************************************************************
635 0 : SUBROUTINE nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
636 : TYPE(nnp_type), INTENT(INOUT) :: nnp
637 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
638 : INTEGER, INTENT(IN) :: unit_nr
639 : LOGICAL, INTENT(IN) :: file_is_new
640 :
641 : CHARACTER(len=default_path_length) :: fmt_string
642 : CHARACTER(LEN=default_string_length), &
643 0 : DIMENSION(:), POINTER :: atomlist
644 : INTEGER :: i, ig, j, n
645 : REAL(KIND=dp), DIMENSION(3) :: rvec
646 :
647 0 : NULLIFY (atomlist)
648 0 : IF (file_is_new) THEN
649 0 : WRITE (unit_nr, "(A)") "# Summed forces [a.u.]"
650 : END IF
651 :
652 0 : rvec = 0.0_dp
653 :
654 : ! get atoms to sum over:
655 : CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
656 0 : c_vals=atomlist)
657 0 : IF (ASSOCIATED(atomlist)) THEN
658 0 : n = SIZE(atomlist)
659 0 : DO i = 1, nnp%num_atoms
660 0 : DO j = 1, n
661 0 : ig = nnp%sort_inv(i)
662 0 : IF (TRIM(ADJUSTL(atomlist(j))) == TRIM(ADJUSTL(nnp%atoms(ig)))) THEN
663 0 : rvec(:) = rvec(:) + nnp%nnp_forces(:, i)
664 : END IF
665 : END DO
666 : END DO
667 : END IF
668 :
669 0 : fmt_string = "(3(F20.10,1X))"
670 0 : WRITE (unit_nr, fmt_string) rvec
671 :
672 0 : END SUBROUTINE nnp_print_sumforces
673 :
674 : END MODULE nnp_force
|
