(git:871dbd5)
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) &
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 ! Moreover let's check if the proper sections have been defined..
125 SELECT CASE (neb_env%pot_type)
126 CASE (pot_neb_fe)
127 wrk_section => section_vals_get_subs_vals(neb_env%root_section, "MOTION%MD")
128 CALL section_vals_get(wrk_section, explicit=explicit)
129 IF (.NOT. explicit) &
130 CALL cp_abort(__location__, &
131 "A free energy BAND (colvars projected) calculation is requested"// &
132 " but NONE MD section was defined in the input.")
133 CASE (pot_neb_me)
134 wrk_section => section_vals_get_subs_vals(neb_env%root_section, "MOTION%GEO_OPT")
135 CALL section_vals_get(wrk_section, explicit=explicit)
136 IF (.NOT. explicit) &
137 CALL cp_abort(__location__, &
138 "A minimum energy BAND (colvars projected) calculation is requested"// &
139 " but NONE GEO_OPT section was defined in the input.")
140 END SELECT
141 ELSE
142 IF (neb_env%use_colvar) &
143 CALL cp_abort(__location__, &
144 "A band calculation was requested with a full potential energy. USE_COLVAR cannot"// &
145 " be set for this kind of calculation!")
146 END IF
147 ! String Method
148 CALL section_vals_val_get(neb_section, "STRING_METHOD%SMOOTHING", r_val=neb_env%smoothing)
149 CALL section_vals_val_get(neb_section, "STRING_METHOD%SPLINE_ORDER", i_val=neb_env%spline_order)
150 neb_env%reparametrize_frames = .false.
151 IF (neb_env%id_type == do_sm) THEN
152 neb_env%reparametrize_frames = .true.
153 END IF
154 END SUBROUTINE read_neb_section
155
156! **************************************************************************************************
157!> \brief dump final structures after a NEB run
158!> \param neb_env ...
159!> \param energies ...
160!> \param coords ...
161!> \param particle_set ...
162!> \param logger ...
163!> \param output_unit ...
164!> \param converged ...
165!> \par
166!> History
167!> 06.2026 - Created
168!> \author HE Zilong
169!> \version 1.0
170! **************************************************************************************************
171 SUBROUTINE dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged)
172 TYPE(neb_type), POINTER :: neb_env
173 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies
174 TYPE(neb_var_type), POINTER :: coords
175 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
176 TYPE(cp_logger_type), POINTER :: logger
177 INTEGER, INTENT(IN) :: output_unit
178 LOGICAL :: converged
179
180 CHARACTER(len=*), PARAMETER :: routinen = 'dump_neb_final'
181
182 CHARACTER(LEN=1024) :: cell_str, ener_str, lm_str, record, &
183 replica_str, title
184 CHARACTER(LEN=4) :: l_ener
185 CHARACTER(LEN=5) :: pbc_str
186 INTEGER :: irep, iw
187 LOGICAL :: print_kind
188 REAL(kind=dp) :: unit_conv
189 TYPE(cell_type), POINTER :: cell
190 TYPE(section_vals_type), POINTER :: final_band_section
191
192 NULLIFY (final_band_section)
193 final_band_section => section_vals_get_subs_vals(neb_env%neb_section, "FINAL_BAND")
194 CALL force_env_get(neb_env%force_env, cell=cell) ! For now NEB has constant cell
195 pbc_str = "F F F"
196 IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
197 IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
198 IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
199 WRITE (unit=cell_str, fmt="(9(1X,F19.10))") &
200 cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom
201 unit_conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
202
203 ! Print a message to log
204 record = cp_print_key_generate_filename(logger, final_band_section, &
205 extension=".xyz", &
206 my_local=.false.)
207 IF (output_unit > 0) THEN
208 IF (converged) THEN
209 WRITE (unit=output_unit, fmt="(/,T2,A)") &
210 routinen//": Band task converged, writing XYZ trajectory gladly:"
211 ELSE
212 WRITE (unit=output_unit, fmt="(/,T2,A)") &
213 routinen//": Band task not yet converged, writing XYZ trajectory anyway:"
214 END IF
215 WRITE (unit=output_unit, fmt="(T3,A)") trim(record)
216 END IF
217
218 ! Write actual trajectory file
219 iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "FINAL_BAND", &
220 extension=".xyz", file_form="FORMATTED", file_status="REPLACE")
221 CALL section_vals_val_get(neb_env%neb_section, "FINAL_BAND%PRINT_ATOM_KIND", &
222 l_val=print_kind)
223 DO irep = 1, neb_env%number_of_replica
224 l_ener = "(**)"
225 IF (irep > 1) THEN
226 IF (energies(irep) - energies(irep - 1) > 0) THEN
227 l_ener(2:2) = "+"
228 ELSE
229 l_ener(2:2) = "-"
230 END IF
231 END IF
232 IF (irep < neb_env%number_of_replica) THEN
233 IF (energies(irep + 1) - energies(irep) < 0) THEN
234 l_ener(3:3) = "+"
235 ELSE
236 l_ener(3:3) = "-"
237 END IF
238 END IF
239 SELECT CASE (l_ener)
240 CASE ("(++)") ! local maximum
241 WRITE (lm_str, '(A)') "Ener_loc_max=T Ener_loc_min=F"
242 CASE ("(--)") ! local minimum
243 WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=T"
244 CASE DEFAULT
245 WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=F"
246 END SELECT
247 WRITE (unit=replica_str, fmt="(I8)") irep
248 WRITE (unit=ener_str, fmt="(F20.10)") energies(irep)
249 WRITE (unit=title, fmt="(A)") &
250 'Lattice="'//trim(adjustl(cell_str))//'" '// &
251 'Properties=species:S:1:pos:R:3 '// &
252 'pbc="'//pbc_str//'" '// &
253 'Replica='//trim(adjustl(replica_str))//' '// &
254 'Energy='//trim(adjustl(ener_str))//' '// &
255 trim(adjustl(lm_str))
256 IF (iw > 0) THEN
257 ! The iw condition does not hold for certain ranks/processes
258 ! that write to <proj>-BAND<n>.out where n > neb_env%number_of_replica
259 CALL write_particle_coordinates(particle_set, iw, dump_extxyz, "POS", title, &
260 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
261 print_kind=print_kind)
262 CALL m_flush(iw)
263 END IF
264 END DO
265
266 IF (output_unit > 0) &
267 WRITE (unit=output_unit, fmt='(/,T2,A)') &
268 routinen//": Done!"
269
270 CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, "FINAL_BAND")
271
272 END SUBROUTINE dump_neb_final
273
274! **************************************************************************************************
275!> \brief dump print info of a NEB run
276!> \param neb_env ...
277!> \param coords ...
278!> \param vels ...
279!> \param forces ...
280!> \param particle_set ...
281!> \param logger ...
282!> \param istep ...
283!> \param energies ...
284!> \param distances ...
285!> \param output_unit ...
286!> \author Teodoro Laino 09.2006
287! **************************************************************************************************
288 SUBROUTINE dump_neb_info(neb_env, coords, vels, forces, particle_set, logger, &
289 istep, energies, distances, output_unit)
290 TYPE(neb_type), POINTER :: neb_env
291 TYPE(neb_var_type), POINTER :: coords
292 TYPE(neb_var_type), OPTIONAL, POINTER :: vels, forces
293 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
294 TYPE(cp_logger_type), POINTER :: logger
295 INTEGER, INTENT(IN) :: istep
296 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: energies, distances
297 INTEGER, INTENT(IN) :: output_unit
298
299 CHARACTER(len=*), PARAMETER :: routinen = 'dump_neb_info'
300
301 CHARACTER(LEN=20) :: mytype
302 CHARACTER(LEN=4) :: l_ener
303 CHARACTER(LEN=default_string_length) :: line, title, unit_str
304 INTEGER :: crd, ener, frc, handle, i, irep, n_max, &
305 n_min, ndig, ndigl, plt, ttst, vel
306 LOGICAL :: explicit, lval, plot_rel_energy, &
307 print_kind
308 REAL(kind=dp) :: ener_min, ener_range, f_ann, tmp_r1, &
309 unit_conv
310 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ekin, temperatures
311 TYPE(cell_type), POINTER :: cell
312 TYPE(enumeration_type), POINTER :: enum
313 TYPE(keyword_type), POINTER :: keyword
314 TYPE(section_type), POINTER :: section
315 TYPE(section_vals_type), POINTER :: run_info_section, tc_section, vc_section
316
317 CALL timeset(routinen, handle)
318 ndig = ceiling(log10(real(neb_env%number_of_replica + 1, kind=dp)))
319 CALL force_env_get(neb_env%force_env, cell=cell)
320 DO irep = 1, neb_env%number_of_replica
321 ndigl = ceiling(log10(real(irep + 1, kind=dp)))
322 WRITE (line, '(A,'//cp_to_string(ndig)//'("0"),T'//cp_to_string(11 + ndig + 1 - ndigl)//',I0)') "Replica_nr_", irep
323 crd = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "TRAJECTORY", &
324 extension=".xyz", file_form="FORMATTED", middle_name="pos-"//trim(line))
325 IF (PRESENT(vels)) THEN
326 vel = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "VELOCITIES", &
327 extension=".xyz", file_form="FORMATTED", middle_name="vel-"//trim(line))
328 END IF
329 IF (PRESENT(forces)) THEN
330 frc = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "FORCES", &
331 extension=".xyz", file_form="FORMATTED", middle_name="force-"//trim(line))
332 END IF
333 ! Dump Trajectory
334 IF (crd > 0) THEN
335 ! Gather units of measure for output
336 CALL section_vals_val_get(neb_env%motion_print_section, "TRAJECTORY%UNIT", &
337 c_val=unit_str)
338 CALL section_vals_val_get(neb_env%motion_print_section, "TRAJECTORY%PRINT_ATOM_KIND", &
339 l_val=print_kind)
340 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
341 ! This information can be digested by Molden
342 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i =", istep, ", E =", energies(irep)
343 CALL write_particle_coordinates(particle_set, crd, dump_xmol, "POS", title, &
344 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
345 print_kind=print_kind)
346 CALL m_flush(crd)
347 END IF
348 ! Dump Velocities
349 IF (vel > 0 .AND. PRESENT(vels)) THEN
350 ! Gather units of measure for output
351 CALL section_vals_val_get(neb_env%motion_print_section, "VELOCITIES%UNIT", &
352 c_val=unit_str)
353 CALL section_vals_val_get(neb_env%motion_print_section, "VELOCITIES%PRINT_ATOM_KIND", &
354 l_val=print_kind)
355 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
356 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i =", istep, ", E =", energies(irep)
357 CALL write_particle_coordinates(particle_set, vel, dump_xmol, "VEL", title, &
358 cell=cell, array=vels%xyz(:, irep), unit_conv=unit_conv, &
359 print_kind=print_kind)
360 CALL m_flush(vel)
361 END IF
362 ! Dump Forces
363 IF (frc > 0 .AND. PRESENT(forces)) THEN
364 ! Gather units of measure for output
365 CALL section_vals_val_get(neb_env%motion_print_section, "FORCES%UNIT", &
366 c_val=unit_str)
367 CALL section_vals_val_get(neb_env%motion_print_section, "FORCES%PRINT_ATOM_KIND", &
368 l_val=print_kind)
369 unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
370 WRITE (unit=title, fmt="(A,I8,A,F20.10)") " i =", istep, ", E =", energies(irep)
371 CALL write_particle_coordinates(particle_set, frc, dump_xmol, "FRC", title, &
372 cell=cell, array=forces%xyz(:, irep), unit_conv=unit_conv, &
373 print_kind=print_kind)
374 CALL m_flush(frc)
375 END IF
376 CALL cp_print_key_finished_output(crd, logger, neb_env%motion_print_section, &
377 "TRAJECTORY")
378 IF (PRESENT(vels)) THEN
379 CALL cp_print_key_finished_output(vel, logger, neb_env%motion_print_section, &
380 "VELOCITIES")
381 END IF
382 IF (PRESENT(forces)) THEN
383 CALL cp_print_key_finished_output(frc, logger, neb_env%motion_print_section, &
384 "FORCES")
385 END IF
386 END DO
387 ! NEB summary info on screen
388 IF (output_unit > 0) THEN
389 tc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%TEMP_CONTROL")
390 vc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%VEL_CONTROL")
391 run_info_section => section_vals_get_subs_vals(neb_env%neb_section, "PROGRAM_RUN_INFO")
392 CALL section_vals_val_get(run_info_section, "PLOT_REL_ENERGY", l_val=plot_rel_energy)
393 ALLOCATE (temperatures(neb_env%number_of_replica))
394 ALLOCATE (ekin(neb_env%number_of_replica))
395 CALL get_temperatures(vels, particle_set, temperatures, ekin=ekin)
396 WRITE (output_unit, '(/)', advance="NO")
397 WRITE (output_unit, fmt='(A,A)') ' **************************************', &
398 '*****************************************'
399 NULLIFY (section, keyword, enum)
400 CALL create_band_section(section)
401 keyword => section_get_keyword(section, "BAND_TYPE")
402 CALL keyword_get(keyword, enum=enum)
403 mytype = trim(enum_i2c(enum, neb_env%id_type))
404 WRITE (output_unit, fmt='(A,T61,A)') &
405 ' BAND TYPE =', adjustr(mytype)
406 CALL section_release(section)
407 WRITE (output_unit, fmt='(A,T61,A)') &
408 ' BAND TYPE OPTIMIZATION =', adjustr(neb_env%opt_type_label(1:20))
409 WRITE (output_unit, '( A,T71,I10 )') &
410 ' STEP NUMBER =', istep
411 IF (neb_env%rotate_frames) WRITE (output_unit, '( A,T71,L10 )') &
412 ' RMSD DISTANCE DEFINITION =', neb_env%rotate_frames
413 ! velocity control parameters output
414 CALL section_vals_get(vc_section, explicit=explicit)
415 IF (explicit) THEN
416 CALL section_vals_val_get(vc_section, "PROJ_VELOCITY_VERLET", l_val=lval)
417 IF (lval) WRITE (output_unit, '( A,T71,L10 )') &
418 ' PROJECTED VELOCITY VERLET =', lval
419 CALL section_vals_val_get(vc_section, "SD_LIKE", l_val=lval)
420 IF (lval) WRITE (output_unit, '( A,T71,L10)') &
421 ' STEEPEST DESCENT LIKE =', lval
422 CALL section_vals_val_get(vc_section, "ANNEALING", r_val=f_ann)
423 IF (f_ann /= 1.0_dp) THEN
424 WRITE (output_unit, '( A,T71,F10.5)') &
425 ' ANNEALING FACTOR = ', f_ann
426 END IF
427 END IF
428 ! temperature control parameters output
429 CALL section_vals_get(tc_section, explicit=explicit)
430 IF (explicit) THEN
431 CALL section_vals_val_get(tc_section, "TEMP_TOL_STEPS", i_val=ttst)
432 IF (istep <= ttst) THEN
433 CALL section_vals_val_get(tc_section, "TEMPERATURE", r_val=f_ann)
434 tmp_r1 = cp_unit_from_cp2k(f_ann, "K")
435 WRITE (output_unit, '( A,T71,F10.5)') &
436 ' TEMPERATURE TARGET =', tmp_r1
437 END IF
438 END IF
439 WRITE (output_unit, '( A,T71,I10 )') &
440 ' NUMBER OF NEB REPLICA =', neb_env%number_of_replica
441 ! switch between a longer visual format and a compact data-only print format
442 IF (plot_rel_energy) THEN
443 cpassert(SIZE(distances) == neb_env%number_of_replica - 1)
444 cpassert(SIZE(energies) == neb_env%number_of_replica)
445 cpassert(SIZE(temperatures) == neb_env%number_of_replica)
446 ener_min = minval(energies(:))
447 ener_range = maxval(energies(:)) - ener_min
448 n_max = 0
449 n_min = 0
450 WRITE (output_unit, '(T2,A,T22,A,T35,A,T52,A)') &
451 'REPLICA', 'ENERGY [au]', 'TEMPERATURE [K]', 'o-------------------------> E'
452 DO i = 1, SIZE(distances)
453 plt = floor((energies(i) - ener_min)/ener_range*25)
454 l_ener = "(**)"
455 IF (i > 1) THEN
456 IF (energies(i) - energies(i - 1) > 0) THEN
457 l_ener(2:2) = "+"
458 ELSE
459 l_ener(2:2) = "-"
460 END IF
461 END IF
462 IF (energies(i + 1) - energies(i) < 0) THEN
463 l_ener(3:3) = "+"
464 ELSE
465 l_ener(3:3) = "-"
466 END IF
467 SELECT CASE (l_ener)
468 CASE ("(++)") ! local maximum
469 n_max = n_max + 1
470 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "X"
471 CASE ("(--)") ! local minimum
472 n_min = n_min + 1
473 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "x"
474 CASE DEFAULT
475 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "O"
476 END SELECT
477 WRITE (output_unit, '(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') &
478 i, energies(i), l_ener, temperatures(i), trim(line)
479 WRITE (output_unit, '(T2,A,1X,F16.6,T52,A)') &
480 "DISTANCE = ", distances(i), "|"
481 END DO
482 plt = floor((energies(neb_env%number_of_replica) - ener_min)/ener_range*25)
483 l_ener = "(**)"
484 IF (energies(neb_env%number_of_replica) - energies(neb_env%number_of_replica - 1) > 0) THEN
485 l_ener(2:2) = "+"
486 ELSE
487 l_ener(2:2) = "-"
488 END IF
489 ! The last point would not be local maximum or minimum, as is the first
490 WRITE (line, '(A,A,A)') "|", repeat(" ", plt), "O"
491 WRITE (output_unit, '(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') &
492 neb_env%number_of_replica, energies(neb_env%number_of_replica), &
493 l_ener, temperatures(neb_env%number_of_replica), trim(line)
494 WRITE (output_unit, '(T52,A)') "v Nr."
495 WRITE (output_unit, '(T2,A,T44,2(1X,I4))') &
496 "NUMBER OF LOCAL MAXIMA (X) and MINIMA (x):", n_max, n_min
497 ELSE
498 WRITE (output_unit, '( A,T17,4F16.6)') &
499 ' DISTANCES REP =', distances(1:min(4, SIZE(distances)))
500 IF (SIZE(distances) > 4) THEN
501 WRITE (output_unit, '( T17,4F16.6)') distances(5:SIZE(distances))
502 END IF
503 WRITE (output_unit, '( A,T17,4F16.6)') &
504 ' ENERGIES [au] =', energies(1:min(4, SIZE(energies)))
505 IF (SIZE(energies) > 4) THEN
506 WRITE (output_unit, '( T17,4F16.6)') energies(5:SIZE(energies))
507 END IF
508 IF (neb_env%opt_type == band_md_opt) THEN
509 WRITE (output_unit, '( A,T33,4(1X,F11.5))') &
510 ' REPLICA TEMPERATURES (K) =', temperatures(1:min(4, SIZE(temperatures)))
511 DO i = 5, SIZE(temperatures), 4
512 WRITE (output_unit, '( T33,4(1X,F11.5))') &
513 temperatures(i:min(i + 3, SIZE(temperatures)))
514 END DO
515 END IF
516 END IF
517 WRITE (output_unit, '( A,T56,F25.14)') &
518 ' BAND TOTAL ENERGY [au] =', sum(energies(:) + ekin(:)) + &
519 neb_env%spring_energy
520 WRITE (output_unit, fmt='(A,A)') ' **************************************', &
521 '*****************************************'
522 DEALLOCATE (ekin)
523 DEALLOCATE (temperatures)
524 END IF
525 ! Ener file
526 ener = cp_print_key_unit_nr(logger, neb_env%neb_section, "ENERGY", &
527 extension=".ener", file_form="FORMATTED")
528 IF (ener > 0) THEN
529 WRITE (line, '(I0)') 2*neb_env%number_of_replica - 1
530 WRITE (ener, '(I10,'//trim(line)//'(1X,F20.9))') istep, &
531 energies, distances
532 END IF
533 CALL cp_print_key_finished_output(ener, logger, neb_env%neb_section, &
534 "ENERGY")
535
536 ! Dump Restarts
537 CALL cp_add_default_logger(logger)
538 CALL write_restart(force_env=neb_env%force_env, &
539 root_section=neb_env%root_section, &
540 coords=coords, &
541 vels=vels)
543
544 CALL timestop(handle)
545
546 END SUBROUTINE dump_neb_info
547
548! **************************************************************************************************
549!> \brief dump coordinates of a replica NEB
550!> \param particle_set ...
551!> \param coords ...
552!> \param i_rep ...
553!> \param ienum ...
554!> \param iw ...
555!> \param use_colvar ...
556!> \author Teodoro Laino 09.2006
557! **************************************************************************************************
558 SUBROUTINE dump_replica_coordinates(particle_set, coords, i_rep, ienum, iw, use_colvar)
559
560 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
561 TYPE(neb_var_type), POINTER :: coords
562 INTEGER, INTENT(IN) :: i_rep, ienum, iw
563 LOGICAL, INTENT(IN) :: use_colvar
564
565 INTEGER :: iatom, j
566 REAL(kind=dp), DIMENSION(3) :: r
567
568 IF (iw > 0) THEN
569 WRITE (iw, '(/,T2,"NEB|",75("*"))')
570 WRITE (iw, '(T2,"NEB|",1X,A,I0,A)') &
571 "Geometry for Replica Nr. ", ienum, " in Angstrom"
572 DO iatom = 1, SIZE(particle_set)
573 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, coords%xyz(:, i_rep))
574 WRITE (iw, '(T2,"NEB|",1X,A10,5X,3F15.9)') &
575 trim(particle_set(iatom)%atomic_kind%name), r(1:3)*angstrom
576 END DO
577 IF (use_colvar) THEN
578 WRITE (iw, '(/,T2,"NEB|",1X,A10)') "COLLECTIVE VARIABLES:"
579 WRITE (iw, '(T2,"NEB|",16X,3F15.9)') &
580 (coords%int(j, i_rep), j=1, SIZE(coords%int(:, :), 1))
581 END IF
582 WRITE (iw, '(T2,"NEB|",75("*"))')
583 CALL m_flush(iw)
584 END IF
585
586 END SUBROUTINE dump_replica_coordinates
587
588! **************************************************************************************************
589!> \brief Handles the correct file names during a band calculation
590!> \param rep_env ...
591!> \param irep ...
592!> \param n_rep ...
593!> \param istep ...
594!> \author Teodoro Laino 06.2009
595! **************************************************************************************************
596 SUBROUTINE handle_band_file_names(rep_env, irep, n_rep, istep)
597 TYPE(replica_env_type), POINTER :: rep_env
598 INTEGER, INTENT(IN) :: irep, n_rep, istep
599
600 CHARACTER(len=*), PARAMETER :: routinen = 'handle_band_file_names'
601
602 CHARACTER(LEN=default_path_length) :: output_file_path, replica_proj_name
603 INTEGER :: handle, handle2, i, ierr, j, lp, unit_nr
604 TYPE(cp_logger_type), POINTER :: logger, sub_logger
605 TYPE(f_env_type), POINTER :: f_env
606 TYPE(section_vals_type), POINTER :: root_section
607
608 CALL timeset(routinen, handle)
609 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
610 handle=handle2)
611 logger => cp_get_default_logger()
612 CALL force_env_get(f_env%force_env, root_section=root_section)
613 j = irep + (rep_env%local_rep_indices(1) - 1)
614 ! Get replica_project_name
615 replica_proj_name = get_replica_project_name(rep_env, n_rep, j)
616 lp = len_trim(replica_proj_name)
617 CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", &
618 c_val=trim(replica_proj_name))
619 logger%iter_info%project_name = trim(replica_proj_name)
620
621 ! We change the file on which is pointing the global logger and error
622 output_file_path = replica_proj_name(1:lp)//".out"
623 CALL section_vals_val_set(root_section, "GLOBAL%OUTPUT_FILE_NAME", &
624 c_val=trim(output_file_path))
625 IF (logger%default_global_unit_nr > 0) THEN
626 CALL close_file(logger%default_global_unit_nr)
627 CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
628 file_action="WRITE", file_position="APPEND", &
629 unit_number=logger%default_global_unit_nr, &
630 skip_get_unit_number=.true.)
631 WRITE (unit=logger%default_global_unit_nr, fmt="(/,(T2,A79))") &
632 "*******************************************************************************", &
633 "** BAND EVALUATION OF ENERGIES AND FORCES **", &
634 "*******************************************************************************"
635 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,T79,A)") "**", "**"
636 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,T79,A)") "**", "**"
637 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,I5,T41,A,I5,T79,A)") &
638 "** Replica Env Nr. :", rep_env%local_rep_indices(1) - 1, "Replica Band Nr. :", j, "**"
639 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A,I5,T79,A)") &
640 "** Band Step Nr. :", istep, "**"
641 WRITE (unit=logger%default_global_unit_nr, fmt="(T2,A79)") &
642 "*******************************************************************************"
643 END IF
644
645 ! Handle specific case for mixed_env
646 SELECT CASE (f_env%force_env%in_use)
647 CASE (use_mixed_force)
648 DO i = 1, f_env%force_env%mixed_env%ngroups
649 IF (modulo(i - 1, f_env%force_env%mixed_env%ngroups) == &
650 f_env%force_env%mixed_env%group_distribution(f_env%force_env%mixed_env%para_env%mepos)) THEN
651 sub_logger => f_env%force_env%mixed_env%sub_logger(i)%p
652 sub_logger%iter_info%project_name = replica_proj_name(1:lp)//"-r-"//trim(adjustl(cp_to_string(i)))
653
654 unit_nr = sub_logger%default_global_unit_nr
655 IF (unit_nr > 0) THEN
656 CALL close_file(unit_nr)
657
658 output_file_path = replica_proj_name(1:lp)//"-r-"//trim(adjustl(cp_to_string(i)))//".out"
659 CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
660 file_action="WRITE", file_position="APPEND", &
661 unit_number=unit_nr, skip_get_unit_number=.true.)
662 END IF
663 END IF
664 END DO
665 END SELECT
666
667 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
668 cpassert(ierr == 0)
669 CALL timestop(handle)
670
671 END SUBROUTINE handle_band_file_names
672
673! **************************************************************************************************
674!> \brief Constructs project names for BAND replicas
675!> \param rep_env ...
676!> \param n_rep ...
677!> \param j ...
678!> \return ...
679!> \author Teodoro Laino 06.2009
680! **************************************************************************************************
681 FUNCTION get_replica_project_name(rep_env, n_rep, j) RESULT(replica_proj_name)
682 TYPE(replica_env_type), POINTER :: rep_env
683 INTEGER, INTENT(IN) :: n_rep, j
684 CHARACTER(LEN=default_path_length) :: replica_proj_name
685
686 CHARACTER(LEN=default_string_length) :: padding
687 INTEGER :: i, lp, ndigits
688
689! Setup new replica project name and output file
690
691 replica_proj_name = rep_env%original_project_name
692 ! Find padding
693 ndigits = ceiling(log10(real(n_rep + 1, kind=dp))) - &
694 ceiling(log10(real(j + 1, kind=dp)))
695 padding = ""
696 DO i = 1, ndigits
697 padding(i:i) = "0"
698 END DO
699 lp = len_trim(replica_proj_name)
700 replica_proj_name(lp + 1:len(replica_proj_name)) = "-BAND"// &
701 trim(padding)//adjustl(cp_to_string(j))
702 END FUNCTION get_replica_project_name
703
704! **************************************************************************************************
705!> \brief Print some mapping infos in the replica_env setup output files
706!> i.e. prints in which files one can find information for each band
707!> replica
708!> \param rep_env ...
709!> \param neb_env ...
710!> \author Teodoro Laino 06.2009
711! **************************************************************************************************
712 SUBROUTINE neb_rep_env_map_info(rep_env, neb_env)
713 TYPE(replica_env_type), POINTER :: rep_env
714 TYPE(neb_type), POINTER :: neb_env
715
716 CHARACTER(LEN=default_path_length) :: replica_proj_name
717 INTEGER :: handle2, ierr, irep, n_rep, n_rep_neb, &
718 output_unit
719 TYPE(cp_logger_type), POINTER :: logger
720 TYPE(f_env_type), POINTER :: f_env
721
722 n_rep_neb = neb_env%number_of_replica
723 n_rep = rep_env%nrep
724 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
725 handle=handle2)
726 logger => cp_get_default_logger()
727 output_unit = logger%default_global_unit_nr
728 IF (output_unit > 0) THEN
729 WRITE (unit=output_unit, fmt='(/,(T2,A79))') &
730 "*******************************************************************************", &
731 "** MAPPING OF BAND REPLICA TO REPLICA ENV **", &
732 "*******************************************************************************"
733 WRITE (unit=output_unit, fmt='(T2,A,I6,T32,A,T79,A)') &
734 "** Replica Env Nr.: ", rep_env%local_rep_indices(1) - 1, &
735 "working on the following BAND replicas", "**"
736 WRITE (unit=output_unit, fmt='(T2,A79)') &
737 "** **"
738 END IF
739 DO irep = 1, n_rep_neb, n_rep
740 replica_proj_name = get_replica_project_name(rep_env, n_rep_neb, irep + rep_env%local_rep_indices(1) - 1)
741 IF (output_unit > 0) THEN
742 WRITE (unit=output_unit, fmt='(T2,A,I6,T32,A,T79,A)') &
743 "** Band Replica Nr.: ", irep + rep_env%local_rep_indices(1) - 1, &
744 "Output available on file: "//trim(replica_proj_name)//".out", "**"
745 END IF
746 END DO
747 IF (output_unit > 0) THEN
748 WRITE (unit=output_unit, fmt='(T2,A79)') &
749 "** **", &
750 "*******************************************************************************"
751 WRITE (unit=output_unit, fmt='(/)')
752 END IF
753 ! update runtime info before printing the footer
754 CALL get_runtime_info()
755 ! print footer
756 CALL cp2k_footer(output_unit)
757 CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
758 cpassert(ierr == 0)
759 END SUBROUTINE neb_rep_env_map_info
760
761END 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:375
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:1239
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:559
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:713
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:290
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:597
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:172
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