(git:6d276e9)
Loading...
Searching...
No Matches
neb_io.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief I/O Module for Nudged Elastic Band Calculation
10!> \note
11!> Numerical accuracy for parallel runs:
12!> Each replica starts the SCF run from the one optimized
13!> in a previous run. It may happen then energies and derivatives
14!> of a serial run and a parallel run could be slightly different
15!> 'cause of a different starting density matrix.
16!> Exact results are obtained using:
17!> EXTRAPOLATION USE_GUESS in QS section (Teo 09.2006)
18!> \author Teodoro Laino 10.2006
19! **************************************************************************************************
20MODULE neb_io
21 USE cell_types, ONLY: cell_type
23 USE cp_files, ONLY: close_file,&
39 USE header, ONLY: cp2k_footer
40 USE input_constants, ONLY: band_md_opt,&
41 do_sm,&
43 dump_xmol,&
61 USE kinds, ONLY: default_path_length,&
63 dp
64 USE machine, ONLY: m_flush
66 USE neb_types, ONLY: neb_type,&
71 USE physcon, ONLY: angstrom
73#include "../base/base_uses.f90"
74
75 IMPLICIT NONE
76 PRIVATE
77 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_io'
78
79 PUBLIC :: read_neb_section, &
85
86CONTAINS
87
88! **************************************************************************************************
89!> \brief Read data from the NEB input section
90!> \param neb_env ...
91!> \param neb_section ...
92!> \author Teodoro Laino 09.2006
93! **************************************************************************************************
94 SUBROUTINE read_neb_section(neb_env, neb_section)
95 TYPE(neb_type), POINTER :: neb_env
96 TYPE(section_vals_type), POINTER :: neb_section
97
98 LOGICAL :: explicit
99 TYPE(section_vals_type), POINTER :: wrk_section
100
101 cpassert(ASSOCIATED(neb_env))
102 neb_env%istep = 0
103 CALL section_vals_val_get(neb_section, "BAND_TYPE", i_val=neb_env%id_type)
104 CALL section_vals_val_get(neb_section, "NUMBER_OF_REPLICA", i_val=neb_env%number_of_replica)
105 CALL section_vals_val_get(neb_section, "K_SPRING", r_val=neb_env%K)
106 CALL section_vals_val_get(neb_section, "ROTATE_FRAMES", l_val=neb_env%rotate_frames)
107 CALL section_vals_val_get(neb_section, "ALIGN_FRAMES", l_val=neb_env%align_frames)
108 CALL section_vals_val_get(neb_section, "OPTIMIZE_BAND%OPTIMIZE_END_POINTS", l_val=neb_env%optimize_end_points)
109 ! Climb Image NEB
110 CALL section_vals_val_get(neb_section, "CI_NEB%NSTEPS_IT", i_val=neb_env%nsteps_it)
111 ! Band Optimization Type
112 CALL section_vals_val_get(neb_section, "OPTIMIZE_BAND%OPT_TYPE", i_val=neb_env%opt_type)
113 ! Use colvars
114 CALL section_vals_val_get(neb_section, "USE_COLVARS", l_val=neb_env%use_colvar)
115 CALL section_vals_val_get(neb_section, "POT_TYPE", i_val=neb_env%pot_type)
116 ! Before continuing let's do some consistency check between keywords
117 IF (neb_env%pot_type /= pot_neb_full) THEN
118 ! Requires the use of colvars
119 IF (.NOT. neb_env%use_colvar) THEN
120 CALL cp_abort(__location__, &
121 "A potential energy function based on free energy or minimum energy"// &
122 " was requested without enabling the usage of COLVARS. Both methods"// &
123 " are based on COLVARS definition.")
124 END IF
125 ! Moreover let's check if the proper sections have been defined..
126 SELECT CASE (neb_env%pot_type)
127 CASE (pot_neb_fe)
128 wrk_section => section_vals_get_subs_vals(neb_env%root_section, "MOTION%MD")
129 CALL section_vals_get(wrk_section, explicit=explicit)
130 IF (.NOT. explicit) THEN
131 CALL cp_abort(__location__, &
132 "A free energy BAND (colvars projected) calculation is requested"// &
133 " but NONE MD section was defined in the input.")
134 END IF
135 CASE (pot_neb_me)
136 wrk_section => section_vals_get_subs_vals(neb_env%root_section, "MOTION%GEO_OPT")
137 CALL section_vals_get(wrk_section, explicit=explicit)
138 IF (.NOT. explicit) THEN
139 CALL cp_abort(__location__, &
140 "A minimum energy BAND (colvars projected) calculation is requested"// &
141 " but NONE GEO_OPT section was defined in the input.")
142 END IF
143 END SELECT
144 ELSE
145 IF (neb_env%use_colvar) THEN
146 CALL cp_abort(__location__, &
147 "A band calculation was requested with a full potential energy. USE_COLVAR cannot"// &
148 " be set for this kind of calculation!")
149 END IF
150 END IF
151 ! String Method
152 CALL section_vals_val_get(neb_section, "STRING_METHOD%SMOOTHING", r_val=neb_env%smoothing)
153 CALL section_vals_val_get(neb_section, "STRING_METHOD%SPLINE_ORDER", i_val=neb_env%spline_order)
154 neb_env%reparametrize_frames = .false.
155 IF (neb_env%id_type == do_sm) THEN
156 neb_env%reparametrize_frames = .true.
157 END IF
158 END SUBROUTINE read_neb_section
159
160! **************************************************************************************************
161!> \brief dump final structures after a NEB run
162!> \param neb_env ...
163!> \param energies ...
164!> \param coords ...
165!> \param particle_set ...
166!> \param logger ...
167!> \param output_unit ...
168!> \param converged ...
169!> \par
170!> History
171!> 06.2026 - Created
172!> \author HE Zilong
173!> \version 1.0
174! **************************************************************************************************
175 SUBROUTINE dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged)
176 TYPE(neb_type), POINTER :: neb_env
177 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies
178 TYPE(neb_var_type), POINTER :: coords
179 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
180 TYPE(cp_logger_type), POINTER :: logger
181 INTEGER, INTENT(IN) :: output_unit
182 LOGICAL :: converged
183
184 CHARACTER(len=*), PARAMETER :: routinen = 'dump_neb_final'
185
186 CHARACTER(LEN=1024) :: cell_str, ener_str, lm_str, record, &
187 replica_str, title
188 CHARACTER(LEN=4) :: l_ener
189 CHARACTER(LEN=5) :: pbc_str
190 INTEGER :: irep, iw
191 LOGICAL :: print_kind
192 REAL(kind=dp) :: unit_conv
193 TYPE(cell_type), POINTER :: cell
194 TYPE(section_vals_type), POINTER :: final_band_section
195
196 NULLIFY (final_band_section)
197 final_band_section => section_vals_get_subs_vals(neb_env%neb_section, "FINAL_BAND")
198 CALL force_env_get(neb_env%force_env, cell=cell) ! For now NEB has constant cell
199 pbc_str = "F F F"
200 IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
201 IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
202 IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
203 WRITE (unit=cell_str, fmt="(9(1X,F19.10))") &
204 cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom
205 unit_conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
206
207 ! Print a message to log
208 record = cp_print_key_generate_filename(logger, final_band_section, &
209 extension=".xyz", &
210 my_local=.false.)
211 IF (output_unit > 0) THEN
212 IF (converged) THEN
213 WRITE (unit=output_unit, fmt="(/,T2,A)") &
214 routinen//": Band task converged, writing XYZ trajectory gladly:"
215 ELSE
216 WRITE (unit=output_unit, fmt="(/,T2,A)") &
217 routinen//": Band task not yet converged, writing XYZ trajectory anyway:"
218 END IF
219 WRITE (unit=output_unit, fmt="(T3,A)") trim(record)
220 END IF
221
222 ! Write actual trajectory file
223 iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "FINAL_BAND", &
224 extension=".xyz", file_form="FORMATTED", file_status="REPLACE")
225 CALL section_vals_val_get(neb_env%neb_section, "FINAL_BAND%PRINT_ATOM_KIND", &
226 l_val=print_kind)
227 DO irep = 1, neb_env%number_of_replica
228 l_ener = "(**)"
229 IF (irep > 1) THEN
230 IF (energies(irep) - energies(irep - 1) > 0) THEN
231 l_ener(2:2) = "+"
232 ELSE
233 l_ener(2:2) = "-"
234 END IF
235 END IF
236 IF (irep < neb_env%number_of_replica) THEN
237 IF (energies(irep + 1) - energies(irep) < 0) THEN
238 l_ener(3:3) = "+"
239 ELSE
240 l_ener(3:3) = "-"
241 END IF
242 END IF
243 SELECT CASE (l_ener)
244 CASE ("(++)") ! local maximum
245 WRITE (lm_str, '(A)') "Ener_loc_max=T Ener_loc_min=F"
246 CASE ("(--)") ! local minimum
247 WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=T"
248 CASE DEFAULT
249 WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=F"
250 END SELECT
251 WRITE (unit=replica_str, fmt="(I8)") irep
252 WRITE (unit=ener_str, fmt="(F20.10)") energies(irep)
253 WRITE (unit=title, fmt="(A)") &
254 'Lattice="'//trim(adjustl(cell_str))//'" '// &
255 'Properties=species:S:1:pos:R:3 '// &
256 'pbc="'//pbc_str//'" '// &
257 'Replica='//trim(adjustl(replica_str))//' '// &
258 'Energy='//trim(adjustl(ener_str))//' '// &
259 trim(adjustl(lm_str))
260 IF (iw > 0) THEN
261 ! The iw condition does not hold for certain ranks/processes
262 ! that write to <proj>-BAND<n>.out where n > neb_env%number_of_replica
263 CALL write_particle_coordinates(particle_set, iw, dump_extxyz, "POS", title, &
264 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
265 print_kind=print_kind)
266 CALL m_flush(iw)
267 END IF
268 END DO
269
270 IF (output_unit > 0) THEN
271 WRITE (unit=output_unit, fmt='(/,T2,A)') &
272 routinen//": Done!"
273 END IF
274
275 CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, "FINAL_BAND")
276
277 END SUBROUTINE dump_neb_final
278
279! **************************************************************************************************
280!> \brief dump print info of a NEB run
281!> \param neb_env ...
282!> \param coords ...
283!> \param vels ...
284!> \param forces ...
285!> \param particle_set ...
286!> \param logger ...
287!> \param istep ...
288!> \param energies ...
289!> \param distances ...
290!> \param output_unit ...
291!> \author Teodoro Laino 09.2006
292! **************************************************************************************************
293 SUBROUTINE dump_neb_info(neb_env, coords, vels, forces, particle_set, logger, &
294 istep, energies, distances, output_unit)
295 TYPE(neb_type), POINTER :: neb_env
296 TYPE(neb_var_type), POINTER :: coords
297 TYPE(neb_var_type), OPTIONAL, POINTER :: vels, forces
298 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
299 TYPE(cp_logger_type), POINTER :: logger
300 INTEGER, INTENT(IN) :: istep
301 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies, distances
302 INTEGER, INTENT(IN) :: output_unit
303
304 CHARACTER(len=*), PARAMETER :: routinen = 'dump_neb_info'
305
306 CHARACTER(LEN=20) :: mytype
307 CHARACTER(LEN=4) :: l_ener
308 CHARACTER(LEN=default_string_length) :: line, title, unit_str
309 INTEGER :: crd, ener, frc, handle, i, irep, n_max, &
310 n_min, ndig, ndigl, plt, ttst, vel
311 LOGICAL :: explicit, lval, plot_rel_energy, &
312 print_kind
313 REAL(kind=dp) :: ener_min, ener_range, f_ann, tmp_r1, &
314 unit_conv
315 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ekin, temperatures
316 TYPE(cell_type), POINTER :: cell
317 TYPE(enumeration_type), POINTER :: enum
318 TYPE(keyword_type), POINTER :: keyword
319 TYPE(section_type), POINTER :: section
320 TYPE(section_vals_type), POINTER :: run_info_section, tc_section, vc_section
321
322 CALL timeset(routinen, handle)
323 ndig = ceiling(log10(real(neb_env%number_of_replica + 1, kind=dp)))
324 CALL force_env_get(neb_env%force_env, cell=cell)
325 DO irep = 1, neb_env%number_of_replica
326 ndigl = ceiling(log10(real(irep + 1, kind=dp)))
327 WRITE (line, '(A,'//cp_to_string(ndig)//'("0"),T'//cp_to_string(11 + ndig + 1 - ndigl)//',I0)') "Replica_nr_", irep
328 crd = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "TRAJECTORY", &
329 extension=".xyz", file_form="FORMATTED", middle_name="pos-"//trim(line))
330 IF (PRESENT(vels)) THEN
331 vel = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "VELOCITIES", &
332 extension=".xyz", file_form="FORMATTED", middle_name="vel-"//trim(line))
333 END IF
334 IF (PRESENT(forces)) THEN
335 frc = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "FORCES", &
336 extension=".xyz", file_form="FORMATTED", middle_name="force-"//trim(line))
337 END IF
338 ! Dump Trajectory
339 IF (crd > 0) THEN
340 ! Gather units of measure for output
341 CALL section_vals_val_get(neb_env%motion_print_section, "TRAJECTORY%UNIT", &
342 c_val=unit_str)
343 CALL section_vals_val_get(neb_env%motion_print_section, "TRAJECTORY%PRINT_ATOM_KIND", &
344 l_val=print_kind)
345 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
346 ! This information can be digested by Molden
347 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i =", istep, ", E =", energies(irep)
348 CALL write_particle_coordinates(particle_set, crd, dump_xmol, "POS", title, &
349 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
350 print_kind=print_kind)
351 CALL m_flush(crd)
352 END IF
353 ! Dump Velocities
354 IF (vel > 0 .AND. PRESENT(vels)) THEN
355 ! Gather units of measure for output
356 CALL section_vals_val_get(neb_env%motion_print_section, "VELOCITIES%UNIT", &
357 c_val=unit_str)
358 CALL section_vals_val_get(neb_env%motion_print_section, "VELOCITIES%PRINT_ATOM_KIND", &
359 l_val=print_kind)
360 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
361 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i =", istep, ", E =", energies(irep)
362 CALL write_particle_coordinates(particle_set, vel, dump_xmol, "VEL", title, &
363 cell=cell, array=vels%xyz(:, irep), unit_conv=unit_conv, &
364 print_kind=print_kind)
365 CALL m_flush(vel)
366 END IF
367 ! Dump Forces
368 IF (frc > 0 .AND. PRESENT(forces)) THEN
369 ! Gather units of measure for output
370 CALL section_vals_val_get(neb_env%motion_print_section, "FORCES%UNIT", &
371 c_val=unit_str)
372 CALL section_vals_val_get(neb_env%motion_print_section, "FORCES%PRINT_ATOM_KIND", &
373 l_val=print_kind)
374 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
375 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i =", istep, ", E =", energies(irep)
376 CALL write_particle_coordinates(particle_set, frc, dump_xmol, "FRC", title, &
377 cell=cell, array=forces%xyz(:, irep), unit_conv=unit_conv, &
378 print_kind=print_kind)
379 CALL m_flush(frc)
380 END IF
381 CALL cp_print_key_finished_output(crd, logger, neb_env%motion_print_section, &
382 "TRAJECTORY")
383 IF (PRESENT(vels)) THEN
384 CALL cp_print_key_finished_output(vel, logger, neb_env%motion_print_section, &
385 "VELOCITIES")
386 END IF
387 IF (PRESENT(forces)) THEN
388 CALL cp_print_key_finished_output(frc, logger, neb_env%motion_print_section, &
389 "FORCES")
390 END IF
391 END DO
392 ! NEB summary info on screen
393 IF (output_unit > 0) THEN
394 tc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%TEMP_CONTROL")
395 vc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%VEL_CONTROL")
396 run_info_section => section_vals_get_subs_vals(neb_env%neb_section, "PROGRAM_RUN_INFO")
397 CALL section_vals_val_get(run_info_section, "PLOT_REL_ENERGY", l_val=plot_rel_energy)
398 ALLOCATE (temperatures(neb_env%number_of_replica))
399 ALLOCATE (ekin(neb_env%number_of_replica))
400 CALL get_temperatures(vels, particle_set, temperatures, ekin=ekin)
401 WRITE (output_unit, '(/)', advance="NO")
402 WRITE (output_unit, fmt='(A,A)') ' **************************************', &
403 '*****************************************'
404 NULLIFY (section, keyword, enum)
405 CALL create_band_section(section)
406 keyword => section_get_keyword(section, "BAND_TYPE")
407 CALL keyword_get(keyword, enum=enum)
408 mytype = trim(enum_i2c(enum, neb_env%id_type))
409 WRITE (output_unit, fmt='(A,T61,A)') &
410 ' BAND TYPE =', adjustr(mytype)
411 CALL section_release(section)
412 WRITE (output_unit, fmt='(A,T61,A)') &
413 ' BAND TYPE OPTIMIZATION =', adjustr(neb_env%opt_type_label(1:20))
414 WRITE (output_unit, '( A,T71,I10 )') &
415 ' STEP NUMBER =', istep
416 IF (neb_env%rotate_frames) WRITE (output_unit, '( A,T71,L10 )') &
417 ' RMSD DISTANCE DEFINITION =', neb_env%rotate_frames
418 ! velocity control parameters output
419 CALL section_vals_get(vc_section, explicit=explicit)
420 IF (explicit) THEN
421 CALL section_vals_val_get(vc_section, "PROJ_VELOCITY_VERLET", l_val=lval)
422 IF (lval) WRITE (output_unit, '( A,T71,L10 )') &
423 ' PROJECTED VELOCITY VERLET =', lval
424 CALL section_vals_val_get(vc_section, "SD_LIKE", l_val=lval)
425 IF (lval) WRITE (output_unit, '( A,T71,L10)') &
426 ' STEEPEST DESCENT LIKE =', lval
427 CALL section_vals_val_get(vc_section, "ANNEALING", r_val=f_ann)
428 IF (f_ann /= 1.0_dp) THEN
429 WRITE (output_unit, '( A,T71,F10.5)') &
430 ' ANNEALING FACTOR = ', f_ann
431 END IF
432 END IF
433 ! temperature control parameters output
434 CALL section_vals_get(tc_section, explicit=explicit)
435 IF (explicit) THEN
436 CALL section_vals_val_get(tc_section, "TEMP_TOL_STEPS", i_val=ttst)
437 IF (istep <= ttst) THEN
438 CALL section_vals_val_get(tc_section, "TEMPERATURE", r_val=f_ann)
439 tmp_r1 = cp_unit_from_cp2k(f_ann, "K")
440 WRITE (output_unit, '( A,T71,F10.5)') &
441 ' TEMPERATURE TARGET =', tmp_r1
442 END IF
443 END IF
444 WRITE (output_unit, '( A,T71,I10 )') &
445 ' NUMBER OF NEB REPLICA =', neb_env%number_of_replica
446 ! switch between a longer visual format and a compact data-only print format
447 IF (plot_rel_energy) THEN
448 cpassert(SIZE(distances) == neb_env%number_of_replica - 1)
449 cpassert(SIZE(energies) == neb_env%number_of_replica)
450 cpassert(SIZE(temperatures) == neb_env%number_of_replica)
451 ener_min = minval(energies(:))
452 ener_range = maxval(energies(:)) - ener_min
453 n_max = 0
454 n_min = 0
455 WRITE (output_unit, '(T2,A,T22,A,T35,A,T52,A)') &
456 'REPLICA', 'ENERGY [au]', 'TEMPERATURE [K]', 'o-------------------------> E'
457 DO i = 1, SIZE(distances)
458 plt = floor((energies(i) - ener_min)/ener_range*25)
459 l_ener = "(**)"
460 IF (i > 1) THEN
461 IF (energies(i) - energies(i - 1) > 0) THEN
462 l_ener(2:2) = "+"
463 ELSE
464 l_ener(2:2) = "-"
465 END IF
466 END IF
467 IF (energies(i + 1) - energies(i) < 0) THEN
468 l_ener(3:3) = "+"
469 ELSE
470 l_ener(3:3) = "-"
471 END IF
472 SELECT CASE (l_ener)
473 CASE ("(++)") ! local maximum
474 n_max = n_max + 1
475 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "X"
476 CASE ("(--)") ! local minimum
477 n_min = n_min + 1
478 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "x"
479 CASE DEFAULT
480 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "O"
481 END SELECT
482 WRITE (output_unit, '(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') &
483 i, energies(i), l_ener, temperatures(i), trim(line)
484 WRITE (output_unit, '(T2,A,1X,F16.6,T52,A)') &
485 "DISTANCE = ", distances(i), "|"
486 END DO
487 plt = floor((energies(neb_env%number_of_replica) - ener_min)/ener_range*25)
488 l_ener = "(**)"
489 IF (energies(neb_env%number_of_replica) - energies(neb_env%number_of_replica - 1) > 0) THEN
490 l_ener(2:2) = "+"
491 ELSE
492 l_ener(2:2) = "-"
493 END IF
494 ! The last point would not be local maximum or minimum, as is the first
495 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "O"
496 WRITE (output_unit, '(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') &
497 neb_env%number_of_replica, energies(neb_env%number_of_replica), &
498 l_ener, temperatures(neb_env%number_of_replica), trim(line)
499 WRITE (output_unit, '(T52,A)') "v Nr."
500 WRITE (output_unit, '(T2,A,T44,2(1X,I4))') &
501 "NUMBER OF LOCAL MAXIMA (X) and MINIMA (x):", n_max, n_min
502 ELSE
503 WRITE (output_unit, '( A,T17,4F16.6)') &
504 ' DISTANCES REP =', distances(1:min(4, SIZE(distances)))
505 IF (SIZE(distances) > 4) THEN
506 WRITE (output_unit, '( T17,4F16.6)') distances(5:SIZE(distances))
507 END IF
508 WRITE (output_unit, '( A,T17,4F16.6)') &
509 ' ENERGIES [au] =', energies(1:min(4, SIZE(energies)))
510 IF (SIZE(energies) > 4) THEN
511 WRITE (output_unit, '( T17,4F16.6)') energies(5:SIZE(energies))
512 END IF
513 IF (neb_env%opt_type == band_md_opt) THEN
514 WRITE (output_unit, '( A,T33,4(1X,F11.5))') &
515 ' REPLICA TEMPERATURES (K) =', temperatures(1:min(4, SIZE(temperatures)))
516 DO i = 5, SIZE(temperatures), 4
517 WRITE (output_unit, '( T33,4(1X,F11.5))') &
518 temperatures(i:min(i + 3, SIZE(temperatures)))
519 END DO
520 END IF
521 END IF
522 WRITE (output_unit, '( A,T56,F25.14)') &
523 ' BAND TOTAL ENERGY [au] =', sum(energies(:) + ekin(:)) + &
524 neb_env%spring_energy
525 WRITE (output_unit, fmt='(A,A)') ' **************************************', &
526 '*****************************************'
527 DEALLOCATE (ekin)
528 DEALLOCATE (temperatures)
529 END IF
530 ! Ener file
531 ener = cp_print_key_unit_nr(logger, neb_env%neb_section, "ENERGY", &
532 extension=".ener", file_form="FORMATTED")
533 IF (ener > 0) THEN
534 WRITE (line, '(I0)') 2*neb_env%number_of_replica - 1
535 WRITE (ener, '(I10,'//trim(line)//'(1X,F20.9))') istep, &
536 energies, distances
537 END IF
538 CALL cp_print_key_finished_output(ener, logger, neb_env%neb_section, &
539 "ENERGY")
540
541 ! Dump Restarts
542 CALL cp_add_default_logger(logger)
543 CALL write_restart(force_env=neb_env%force_env, &
544 root_section=neb_env%root_section, &
545 coords=coords, &
546 vels=vels)
548
549 CALL timestop(handle)
550
551 END SUBROUTINE dump_neb_info
552
553! **************************************************************************************************
554!> \brief dump coordinates of a replica NEB
555!> \param particle_set ...
556!> \param coords ...
557!> \param i_rep ...
558!> \param ienum ...
559!> \param iw ...
560!> \param use_colvar ...
561!> \author Teodoro Laino 09.2006
562! **************************************************************************************************
563 SUBROUTINE dump_replica_coordinates(particle_set, coords, i_rep, ienum, iw, use_colvar)
564
565 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
566 TYPE(neb_var_type), POINTER :: coords
567 INTEGER, INTENT(IN) :: i_rep, ienum, iw
568 LOGICAL, INTENT(IN) :: use_colvar
569
570 INTEGER :: iatom, j
571 REAL(kind=dp), DIMENSION(3) :: r
572
573 IF (iw > 0) THEN
574 WRITE (iw, '(/,T2,"NEB|",75("*"))')
575 WRITE (iw, '(T2,"NEB|",1X,A,I0,A)') &
576 "Geometry for Replica Nr. ", ienum, " in Angstrom"
577 DO iatom = 1, SIZE(particle_set)
578 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, coords%xyz(:, i_rep))
579 WRITE (iw, '(T2,"NEB|",1X,A10,5X,3F15.9)') &
580 trim(particle_set(iatom)%atomic_kind%name), r(1:3)*angstrom
581 END DO
582 IF (use_colvar) THEN
583 WRITE (iw, '(/,T2,"NEB|",1X,A10)') "COLLECTIVE VARIABLES:"
584 WRITE (iw, '(T2,"NEB|",16X,3F15.9)') &
585 (coords%int(j, i_rep), j=1, SIZE(coords%int(:, :), 1))
586 END IF
587 WRITE (iw, '(T2,"NEB|",75("*"))')
588 CALL m_flush(iw)
589 END IF
590
591 END SUBROUTINE dump_replica_coordinates
592
593! **************************************************************************************************
594!> \brief Handles the correct file names during a band calculation
595!> \param rep_env ...
596!> \param irep ...
597!> \param n_rep ...
598!> \param istep ...
599!> \author Teodoro Laino 06.2009
600! **************************************************************************************************
601 SUBROUTINE handle_band_file_names(rep_env, irep, n_rep, istep)
602 TYPE(replica_env_type), POINTER :: rep_env
603 INTEGER, INTENT(IN) :: irep, n_rep, istep
604
605 CHARACTER(len=*), PARAMETER :: routinen = 'handle_band_file_names'
606
607 CHARACTER(LEN=default_path_length) :: output_file_path, replica_proj_name
608 INTEGER :: handle, handle2, i, ierr, j, lp, unit_nr
609 TYPE(cp_logger_type), POINTER :: logger, sub_logger
610 TYPE(f_env_type), POINTER :: f_env
611 TYPE(section_vals_type), POINTER :: root_section
612
613 CALL timeset(routinen, handle)
614 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
615 handle=handle2)
616 logger => cp_get_default_logger()
617 CALL force_env_get(f_env%force_env, root_section=root_section)
618 j = irep + (rep_env%local_rep_indices(1) - 1)
619 ! Get replica_project_name
620 replica_proj_name = get_replica_project_name(rep_env, n_rep, j)
621 lp = len_trim(replica_proj_name)
622 CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", &
623 c_val=trim(replica_proj_name))
624 logger%iter_info%project_name = trim(replica_proj_name)
625
626 ! We change the file on which is pointing the global logger and error
627 output_file_path = replica_proj_name(1:lp)//".out"
628 CALL section_vals_val_set(root_section, "GLOBAL%OUTPUT_FILE_NAME", &
629 c_val=trim(output_file_path))
630 IF (logger%default_global_unit_nr > 0) THEN
631 CALL close_file(logger%default_global_unit_nr)
632 CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
633 file_action="WRITE", file_position="APPEND", &
634 unit_number=logger%default_global_unit_nr, &
635 skip_get_unit_number=.true.)
636 WRITE (unit=logger%default_global_unit_nr, fmt="(/,(T2,A79))") &
637 "*******************************************************************************", &
638 "** BAND EVALUATION OF ENERGIES AND FORCES **", &
639 "*******************************************************************************"
640 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,T79,A)") "**", "**"
641 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,T79,A)") "**", "**"
642 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,I5,T41,A,I5,T79,A)") &
643 "** Replica Env Nr. :", rep_env%local_rep_indices(1) - 1, "Replica Band Nr. :", j, "**"
644 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,I5,T79,A)") &
645 "** Band Step Nr. :", istep, "**"
646 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A79)") &
647 "*******************************************************************************"
648 END IF
649
650 ! Handle specific case for mixed_env
651 SELECT CASE (f_env%force_env%in_use)
652 CASE (use_mixed_force)
653 DO i = 1, f_env%force_env%mixed_env%ngroups
654 IF (modulo(i - 1, f_env%force_env%mixed_env%ngroups) == &
655 f_env%force_env%mixed_env%group_distribution(f_env%force_env%mixed_env%para_env%mepos)) THEN
656 sub_logger => f_env%force_env%mixed_env%sub_logger(i)%p
657 sub_logger%iter_info%project_name = replica_proj_name(1:lp)//"-r-"//trim(adjustl(cp_to_string(i)))
658
659 unit_nr = sub_logger%default_global_unit_nr
660 IF (unit_nr > 0) THEN
661 CALL close_file(unit_nr)
662
663 output_file_path = replica_proj_name(1:lp)//"-r-"//trim(adjustl(cp_to_string(i)))//".out"
664 CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
665 file_action="WRITE", file_position="APPEND", &
666 unit_number=unit_nr, skip_get_unit_number=.true.)
667 END IF
668 END IF
669 END DO
670 END SELECT
671
672 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
673 cpassert(ierr == 0)
674 CALL timestop(handle)
675
676 END SUBROUTINE handle_band_file_names
677
678! **************************************************************************************************
679!> \brief Constructs project names for BAND replicas
680!> \param rep_env ...
681!> \param n_rep ...
682!> \param j ...
683!> \return ...
684!> \author Teodoro Laino 06.2009
685! **************************************************************************************************
686 FUNCTION get_replica_project_name(rep_env, n_rep, j) RESULT(replica_proj_name)
687 TYPE(replica_env_type), POINTER :: rep_env
688 INTEGER, INTENT(IN) :: n_rep, j
689 CHARACTER(LEN=default_path_length) :: replica_proj_name
690
691 CHARACTER(LEN=default_string_length) :: padding
692 INTEGER :: i, lp, ndigits
693
694! Setup new replica project name and output file
695
696 replica_proj_name = rep_env%original_project_name
697 ! Find padding
698 ndigits = ceiling(log10(real(n_rep + 1, kind=dp))) - &
699 ceiling(log10(real(j + 1, kind=dp)))
700 padding = ""
701 DO i = 1, ndigits
702 padding(i:i) = "0"
703 END DO
704 lp = len_trim(replica_proj_name)
705 replica_proj_name(lp + 1:len(replica_proj_name)) = "-BAND"// &
706 trim(padding)//adjustl(cp_to_string(j))
707 END FUNCTION get_replica_project_name
708
709! **************************************************************************************************
710!> \brief Print some mapping infos in the replica_env setup output files
711!> i.e. prints in which files one can find information for each band
712!> replica
713!> \param rep_env ...
714!> \param neb_env ...
715!> \author Teodoro Laino 06.2009
716! **************************************************************************************************
717 SUBROUTINE neb_rep_env_map_info(rep_env, neb_env)
718 TYPE(replica_env_type), POINTER :: rep_env
719 TYPE(neb_type), POINTER :: neb_env
720
721 CHARACTER(LEN=default_path_length) :: replica_proj_name
722 INTEGER :: handle2, ierr, irep, n_rep, n_rep_neb, &
723 output_unit
724 TYPE(cp_logger_type), POINTER :: logger
725 TYPE(f_env_type), POINTER :: f_env
726
727 n_rep_neb = neb_env%number_of_replica
728 n_rep = rep_env%nrep
729 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
730 handle=handle2)
731 logger => cp_get_default_logger()
732 output_unit = logger%default_global_unit_nr
733 IF (output_unit > 0) THEN
734 WRITE (unit=output_unit, fmt='(/,(T2,A79))') &
735 "*******************************************************************************", &
736 "** MAPPING OF BAND REPLICA TO REPLICA ENV **", &
737 "*******************************************************************************"
738 WRITE (unit=output_unit, fmt='(T2,A,I6,T32,A,T79,A)') &
739 "** Replica Env Nr.: ", rep_env%local_rep_indices(1) - 1, &
740 "working on the following BAND replicas", "**"
741 WRITE (unit=output_unit, fmt='(T2,A79)') &
742 "** **"
743 END IF
744 DO irep = 1, n_rep_neb, n_rep
745 replica_proj_name = get_replica_project_name(rep_env, n_rep_neb, irep + rep_env%local_rep_indices(1) - 1)
746 IF (output_unit > 0) THEN
747 WRITE (unit=output_unit, fmt='(T2,A,I6,T32,A,T79,A)') &
748 "** Band Replica Nr.: ", irep + rep_env%local_rep_indices(1) - 1, &
749 "Output available on file: "//trim(replica_proj_name)//".out", "**"
750 END IF
751 END DO
752 IF (output_unit > 0) THEN
753 WRITE (unit=output_unit, fmt='(T2,A79)') &
754 "** **", &
755 "*******************************************************************************"
756 WRITE (unit=output_unit, fmt='(/)')
757 END IF
758 ! update runtime info before printing the footer
759 CALL get_runtime_info()
760 ! print footer
761 CALL cp2k_footer(output_unit)
762 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
763 cpassert(ierr == 0)
764 END SUBROUTINE neb_rep_env_map_info
765
766END MODULE neb_io
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Handles all functions related to the CELL.
Definition cell_types.F:15
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:22
subroutine, public get_runtime_info()
...
Definition cp2k_info.F:381
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
various routines to log and control the output. The idea is that decisions about where to log should ...
subroutine, public cp_rm_default_logger()
the cousin of cp_add_default_logger, decrements the stack, so that the default logger is what it has ...
subroutine, public cp_add_default_logger(logger)
adds a default logger. MUST be called before logging occours
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1251
interface to use cp2k as library
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
Interface for the force calculations.
integer, parameter, public use_mixed_force
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
subroutine, public cp2k_footer(iw, wdir)
...
Definition header.F:69
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public pot_neb_me
integer, parameter, public band_md_opt
integer, parameter, public pot_neb_fe
integer, parameter, public dump_xmol
integer, parameter, public dump_extxyz
integer, parameter, public do_sm
integer, parameter, public pot_neb_full
subroutine, public create_band_section(section)
creates the section for a BAND run
Set of routines to dump the restart file of CP2K.
subroutine, public write_restart(md_env, force_env, root_section, coords, vels, pint_env, helium_env)
checks if a restart needs to be written and does so, updating all necessary fields in the input file....
represents an enumeration, i.e. a mapping between integers and strings
character(len=default_string_length) function, public enum_i2c(enum, i)
maps an integer to a string
represents keywords in an input
subroutine, public keyword_get(keyword, names, usage, description, type_of_var, n_var, default_value, lone_keyword_value, repeats, enum, citations)
...
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
recursive type(keyword_type) function, pointer, public section_get_keyword(section, keyword_name)
returns the requested keyword
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
I/O Module for Nudged Elastic Band Calculation.
Definition neb_io.F:20
subroutine, public read_neb_section(neb_env, neb_section)
Read data from the NEB input section.
Definition neb_io.F:95
subroutine, public dump_replica_coordinates(particle_set, coords, i_rep, ienum, iw, use_colvar)
dump coordinates of a replica NEB
Definition neb_io.F:564
subroutine, public neb_rep_env_map_info(rep_env, neb_env)
Print some mapping infos in the replica_env setup output files i.e. prints in which files one can fin...
Definition neb_io.F:718
subroutine, public dump_neb_info(neb_env, coords, vels, forces, particle_set, logger, istep, energies, distances, output_unit)
dump print info of a NEB run
Definition neb_io.F:295
subroutine, public handle_band_file_names(rep_env, irep, n_rep, istep)
Handles the correct file names during a band calculation.
Definition neb_io.F:602
subroutine, public dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged)
dump final structures after a NEB run
Definition neb_io.F:176
Module with utility to perform MD Nudged Elastic Band Calculation.
subroutine, public get_temperatures(vels, particle_set, temperatures, ekin, factor)
Computes temperatures.
Typo for Nudged Elastic Band Calculation.
Definition neb_types.F:20
Define methods related to particle_type.
subroutine, public write_particle_coordinates(particle_set, iunit, output_format, content, title, cell, array, unit_conv, charge_occup, charge_beta, charge_extended, print_kind)
Should be able to write a few formats e.g. xmol, and some binary format (dcd) some format can be used...
Define the data structure for the particle information.
pure real(kind=dp) function, dimension(3), public get_particle_pos_or_vel(iatom, particle_set, vector)
Return the atomic position or velocity of atom iatom in x from a packed vector even if core-shell par...
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
types used to handle many replica of the same system that differ only in atom positions,...
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represent a keyword in the input
represent a section of the input file
keeps replicated information about the replicas