Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
20 : MODULE neb_io
21 : USE cell_types, ONLY: cell_type
22 : USE cp2k_info, ONLY: get_runtime_info
23 : USE cp_files, ONLY: close_file,&
24 : open_file
25 : USE cp_log_handling, ONLY: cp_add_default_logger,&
26 : cp_get_default_logger,&
27 : cp_logger_type,&
28 : cp_rm_default_logger,&
29 : cp_to_string
30 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
31 : cp_print_key_generate_filename,&
32 : cp_print_key_unit_nr
33 : USE cp_units, ONLY: cp_unit_from_cp2k
34 : USE f77_interface, ONLY: f_env_add_defaults,&
35 : f_env_rm_defaults,&
36 : f_env_type
37 : USE force_env_types, ONLY: force_env_get,&
38 : use_mixed_force
39 : USE header, ONLY: cp2k_footer
40 : USE input_constants, ONLY: band_md_opt,&
41 : do_sm,&
42 : dump_extxyz,&
43 : dump_xmol,&
44 : pot_neb_fe,&
45 : pot_neb_full,&
46 : pot_neb_me
47 : USE input_cp2k_neb, ONLY: create_band_section
48 : USE input_cp2k_restarts, ONLY: write_restart
49 : USE input_enumeration_types, ONLY: enum_i2c,&
50 : enumeration_type
51 : USE input_keyword_types, ONLY: keyword_get,&
52 : keyword_type
53 : USE input_section_types, ONLY: section_get_keyword,&
54 : section_release,&
55 : section_type,&
56 : section_vals_get,&
57 : section_vals_get_subs_vals,&
58 : section_vals_type,&
59 : section_vals_val_get,&
60 : section_vals_val_set
61 : USE kinds, ONLY: default_path_length,&
62 : default_string_length,&
63 : dp
64 : USE machine, ONLY: m_flush
65 : USE neb_md_utils, ONLY: get_temperatures
66 : USE neb_types, ONLY: neb_type,&
67 : neb_var_type
68 : USE particle_methods, ONLY: write_particle_coordinates
69 : USE particle_types, ONLY: get_particle_pos_or_vel,&
70 : particle_type
71 : USE physcon, ONLY: angstrom
72 : USE replica_types, ONLY: replica_env_type
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, &
80 : dump_neb_final, &
81 : dump_neb_info, &
82 : dump_replica_coordinates, &
83 : handle_band_file_names, &
84 : neb_rep_env_map_info
85 :
86 : CONTAINS
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 34 : 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 34 : CPASSERT(ASSOCIATED(neb_env))
102 34 : neb_env%istep = 0
103 34 : CALL section_vals_val_get(neb_section, "BAND_TYPE", i_val=neb_env%id_type)
104 34 : CALL section_vals_val_get(neb_section, "NUMBER_OF_REPLICA", i_val=neb_env%number_of_replica)
105 34 : CALL section_vals_val_get(neb_section, "K_SPRING", r_val=neb_env%K)
106 34 : CALL section_vals_val_get(neb_section, "ROTATE_FRAMES", l_val=neb_env%rotate_frames)
107 34 : CALL section_vals_val_get(neb_section, "ALIGN_FRAMES", l_val=neb_env%align_frames)
108 34 : CALL section_vals_val_get(neb_section, "OPTIMIZE_BAND%OPTIMIZE_END_POINTS", l_val=neb_env%optimize_end_points)
109 : ! Climb Image NEB
110 34 : CALL section_vals_val_get(neb_section, "CI_NEB%NSTEPS_IT", i_val=neb_env%nsteps_it)
111 : ! Band Optimization Type
112 34 : CALL section_vals_val_get(neb_section, "OPTIMIZE_BAND%OPT_TYPE", i_val=neb_env%opt_type)
113 : ! Use colvars
114 34 : CALL section_vals_val_get(neb_section, "USE_COLVARS", l_val=neb_env%use_colvar)
115 34 : 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 34 : IF (neb_env%pot_type /= pot_neb_full) THEN
118 : ! Requires the use of colvars
119 4 : 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 0 : " are based on COLVARS definition.")
124 : END IF
125 : ! Moreover let's check if the proper sections have been defined..
126 4 : SELECT CASE (neb_env%pot_type)
127 : CASE (pot_neb_fe)
128 0 : wrk_section => section_vals_get_subs_vals(neb_env%root_section, "MOTION%MD")
129 0 : CALL section_vals_get(wrk_section, explicit=explicit)
130 0 : IF (.NOT. explicit) THEN
131 : CALL cp_abort(__LOCATION__, &
132 : "A free energy BAND (colvars projected) calculation is requested"// &
133 0 : " but NONE MD section was defined in the input.")
134 : END IF
135 : CASE (pot_neb_me)
136 4 : wrk_section => section_vals_get_subs_vals(neb_env%root_section, "MOTION%GEO_OPT")
137 4 : CALL section_vals_get(wrk_section, explicit=explicit)
138 8 : IF (.NOT. explicit) THEN
139 : CALL cp_abort(__LOCATION__, &
140 : "A minimum energy BAND (colvars projected) calculation is requested"// &
141 0 : " but NONE GEO_OPT section was defined in the input.")
142 : END IF
143 : END SELECT
144 : ELSE
145 30 : 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 0 : " be set for this kind of calculation!")
149 : END IF
150 : END IF
151 : ! String Method
152 34 : CALL section_vals_val_get(neb_section, "STRING_METHOD%SMOOTHING", r_val=neb_env%smoothing)
153 34 : CALL section_vals_val_get(neb_section, "STRING_METHOD%SPLINE_ORDER", i_val=neb_env%spline_order)
154 34 : neb_env%reparametrize_frames = .FALSE.
155 34 : IF (neb_env%id_type == do_sm) THEN
156 2 : neb_env%reparametrize_frames = .TRUE.
157 : END IF
158 34 : 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 34 : 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 34 : NULLIFY (final_band_section)
197 34 : final_band_section => section_vals_get_subs_vals(neb_env%neb_section, "FINAL_BAND")
198 34 : CALL force_env_get(neb_env%force_env, cell=cell) ! For now NEB has constant cell
199 34 : pbc_str = "F F F"
200 34 : IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
201 34 : IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
202 34 : IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
203 : WRITE (UNIT=cell_str, FMT="(9(1X,F19.10))") &
204 340 : cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom
205 34 : 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 34 : my_local=.FALSE.)
211 34 : IF (output_unit > 0) THEN
212 17 : IF (converged) THEN
213 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
214 2 : routineN//": Band task converged, writing XYZ trajectory gladly:"
215 : ELSE
216 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
217 15 : routineN//": Band task not yet converged, writing XYZ trajectory anyway:"
218 : END IF
219 17 : 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 34 : extension=".xyz", file_form="FORMATTED", file_status="REPLACE")
225 : CALL section_vals_val_get(neb_env%neb_section, "FINAL_BAND%PRINT_ATOM_KIND", &
226 34 : l_val=print_kind)
227 246 : DO irep = 1, neb_env%number_of_replica
228 212 : l_ener = "(**)"
229 212 : IF (irep > 1) THEN
230 178 : IF (energies(irep) - energies(irep - 1) > 0) THEN
231 106 : l_ener(2:2) = "+"
232 : ELSE
233 72 : l_ener(2:2) = "-"
234 : END IF
235 : END IF
236 212 : IF (irep < neb_env%number_of_replica) THEN
237 178 : IF (energies(irep + 1) - energies(irep) < 0) THEN
238 72 : l_ener(3:3) = "+"
239 : ELSE
240 106 : l_ener(3:3) = "-"
241 : END IF
242 : END IF
243 46 : SELECT CASE (l_ener)
244 : CASE ("(++)") ! local maximum
245 46 : WRITE (lm_str, '(A)') "Ener_loc_max=T Ener_loc_min=F"
246 : CASE ("(--)") ! local minimum
247 34 : WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=T"
248 : CASE DEFAULT
249 212 : WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=F"
250 : END SELECT
251 212 : WRITE (UNIT=replica_str, FMT="(I8)") irep
252 212 : 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 212 : TRIM(ADJUSTL(lm_str))
260 246 : 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 106 : print_kind=print_kind)
266 106 : CALL m_flush(iw)
267 : END IF
268 : END DO
269 :
270 34 : IF (output_unit > 0) THEN
271 : WRITE (UNIT=output_unit, FMT='(/,T2,A)') &
272 17 : routineN//": Done!"
273 : END IF
274 :
275 34 : CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, "FINAL_BAND")
276 :
277 34 : 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 578 : SUBROUTINE dump_neb_info(neb_env, coords, vels, forces, particle_set, logger, &
294 578 : 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 578 : 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 578 : CALL timeset(routineN, handle)
323 578 : ndig = CEILING(LOG10(REAL(neb_env%number_of_replica + 1, KIND=dp)))
324 578 : CALL force_env_get(neb_env%force_env, cell=cell)
325 4152 : DO irep = 1, neb_env%number_of_replica
326 3574 : ndigl = CEILING(LOG10(REAL(irep + 1, KIND=dp)))
327 3574 : 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 3574 : extension=".xyz", file_form="FORMATTED", middle_name="pos-"//TRIM(line))
330 3574 : IF (PRESENT(vels)) THEN
331 : vel = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "VELOCITIES", &
332 3574 : extension=".xyz", file_form="FORMATTED", middle_name="vel-"//TRIM(line))
333 : END IF
334 3574 : IF (PRESENT(forces)) THEN
335 : frc = cp_print_key_unit_nr(logger, neb_env%motion_print_section, "FORCES", &
336 3574 : extension=".xyz", file_form="FORMATTED", middle_name="force-"//TRIM(line))
337 : END IF
338 : ! Dump Trajectory
339 3574 : 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 1565 : c_val=unit_str)
343 : CALL section_vals_val_get(neb_env%motion_print_section, "TRAJECTORY%PRINT_ATOM_KIND", &
344 1565 : l_val=print_kind)
345 1565 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
346 : ! This information can be digested by Molden
347 1565 : 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 1565 : print_kind=print_kind)
351 1565 : CALL m_flush(crd)
352 : END IF
353 : ! Dump Velocities
354 3574 : 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 0 : c_val=unit_str)
358 : CALL section_vals_val_get(neb_env%motion_print_section, "VELOCITIES%PRINT_ATOM_KIND", &
359 0 : l_val=print_kind)
360 0 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
361 0 : 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 0 : print_kind=print_kind)
365 0 : CALL m_flush(vel)
366 : END IF
367 : ! Dump Forces
368 3574 : 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 0 : c_val=unit_str)
372 : CALL section_vals_val_get(neb_env%motion_print_section, "FORCES%PRINT_ATOM_KIND", &
373 0 : l_val=print_kind)
374 0 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
375 0 : 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 0 : print_kind=print_kind)
379 0 : CALL m_flush(frc)
380 : END IF
381 : CALL cp_print_key_finished_output(crd, logger, neb_env%motion_print_section, &
382 3574 : "TRAJECTORY")
383 3574 : IF (PRESENT(vels)) THEN
384 : CALL cp_print_key_finished_output(vel, logger, neb_env%motion_print_section, &
385 3574 : "VELOCITIES")
386 : END IF
387 4152 : IF (PRESENT(forces)) THEN
388 : CALL cp_print_key_finished_output(frc, logger, neb_env%motion_print_section, &
389 3574 : "FORCES")
390 : END IF
391 : END DO
392 : ! NEB summary info on screen
393 578 : IF (output_unit > 0) THEN
394 289 : tc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%TEMP_CONTROL")
395 289 : vc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%VEL_CONTROL")
396 289 : run_info_section => section_vals_get_subs_vals(neb_env%neb_section, "PROGRAM_RUN_INFO")
397 289 : CALL section_vals_val_get(run_info_section, "PLOT_REL_ENERGY", l_val=plot_rel_energy)
398 867 : ALLOCATE (temperatures(neb_env%number_of_replica))
399 578 : ALLOCATE (ekin(neb_env%number_of_replica))
400 289 : CALL get_temperatures(vels, particle_set, temperatures, ekin=ekin)
401 289 : WRITE (output_unit, '(/)', ADVANCE="NO")
402 289 : WRITE (output_unit, FMT='(A,A)') ' **************************************', &
403 578 : '*****************************************'
404 289 : NULLIFY (section, keyword, enum)
405 289 : CALL create_band_section(section)
406 289 : keyword => section_get_keyword(section, "BAND_TYPE")
407 289 : CALL keyword_get(keyword, enum=enum)
408 289 : mytype = TRIM(enum_i2c(enum, neb_env%id_type))
409 : WRITE (output_unit, FMT='(A,T61,A)') &
410 289 : ' BAND TYPE =', ADJUSTR(mytype)
411 289 : CALL section_release(section)
412 : WRITE (output_unit, FMT='(A,T61,A)') &
413 289 : ' BAND TYPE OPTIMIZATION =', ADJUSTR(neb_env%opt_type_label(1:20))
414 : WRITE (output_unit, '( A,T71,I10 )') &
415 289 : ' STEP NUMBER =', istep
416 289 : IF (neb_env%rotate_frames) WRITE (output_unit, '( A,T71,L10 )') &
417 80 : ' RMSD DISTANCE DEFINITION =', neb_env%rotate_frames
418 : ! velocity control parameters output
419 289 : CALL section_vals_get(vc_section, explicit=explicit)
420 289 : IF (explicit) THEN
421 88 : CALL section_vals_val_get(vc_section, "PROJ_VELOCITY_VERLET", l_val=lval)
422 88 : IF (lval) WRITE (output_unit, '( A,T71,L10 )') &
423 77 : ' PROJECTED VELOCITY VERLET =', lval
424 88 : CALL section_vals_val_get(vc_section, "SD_LIKE", l_val=lval)
425 88 : IF (lval) WRITE (output_unit, '( A,T71,L10)') &
426 0 : ' STEEPEST DESCENT LIKE =', lval
427 88 : CALL section_vals_val_get(vc_section, "ANNEALING", r_val=f_ann)
428 88 : IF (f_ann /= 1.0_dp) THEN
429 : WRITE (output_unit, '( A,T71,F10.5)') &
430 88 : ' ANNEALING FACTOR = ', f_ann
431 : END IF
432 : END IF
433 : ! temperature control parameters output
434 289 : CALL section_vals_get(tc_section, explicit=explicit)
435 289 : IF (explicit) THEN
436 32 : CALL section_vals_val_get(tc_section, "TEMP_TOL_STEPS", i_val=ttst)
437 32 : IF (istep <= ttst) THEN
438 22 : CALL section_vals_val_get(tc_section, "TEMPERATURE", r_val=f_ann)
439 22 : tmp_r1 = cp_unit_from_cp2k(f_ann, "K")
440 : WRITE (output_unit, '( A,T71,F10.5)') &
441 22 : ' TEMPERATURE TARGET =', tmp_r1
442 : END IF
443 : END IF
444 : WRITE (output_unit, '( A,T71,I10 )') &
445 289 : ' NUMBER OF NEB REPLICA =', neb_env%number_of_replica
446 : ! switch between a longer visual format and a compact data-only print format
447 289 : IF (plot_rel_energy) THEN
448 0 : CPASSERT(SIZE(distances) == neb_env%number_of_replica - 1)
449 0 : CPASSERT(SIZE(energies) == neb_env%number_of_replica)
450 0 : CPASSERT(SIZE(temperatures) == neb_env%number_of_replica)
451 0 : ener_min = MINVAL(energies(:))
452 0 : ener_range = MAXVAL(energies(:)) - ener_min
453 0 : n_max = 0
454 0 : n_min = 0
455 : WRITE (output_unit, '(T2,A,T22,A,T35,A,T52,A)') &
456 0 : 'REPLICA', 'ENERGY [au]', 'TEMPERATURE [K]', 'o-------------------------> E'
457 0 : DO i = 1, SIZE(distances)
458 0 : plt = FLOOR((energies(i) - ener_min)/ener_range*25)
459 0 : l_ener = "(**)"
460 0 : IF (i > 1) THEN
461 0 : IF (energies(i) - energies(i - 1) > 0) THEN
462 0 : l_ener(2:2) = "+"
463 : ELSE
464 0 : l_ener(2:2) = "-"
465 : END IF
466 : END IF
467 0 : IF (energies(i + 1) - energies(i) < 0) THEN
468 0 : l_ener(3:3) = "+"
469 : ELSE
470 0 : l_ener(3:3) = "-"
471 : END IF
472 0 : SELECT CASE (l_ener)
473 : CASE ("(++)") ! local maximum
474 0 : n_max = n_max + 1
475 0 : WRITE (line, '(A,A,A)') "|", REPEAT(" ", plt), "X"
476 : CASE ("(--)") ! local minimum
477 0 : n_min = n_min + 1
478 0 : WRITE (line, '(A,A,A)') "|", REPEAT(" ", plt), "x"
479 : CASE DEFAULT
480 0 : 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 0 : i, energies(i), l_ener, temperatures(i), TRIM(line)
484 : WRITE (output_unit, '(T2,A,1X,F16.6,T52,A)') &
485 0 : "DISTANCE = ", distances(i), "|"
486 : END DO
487 0 : plt = FLOOR((energies(neb_env%number_of_replica) - ener_min)/ener_range*25)
488 0 : l_ener = "(**)"
489 0 : IF (energies(neb_env%number_of_replica) - energies(neb_env%number_of_replica - 1) > 0) THEN
490 0 : l_ener(2:2) = "+"
491 : ELSE
492 0 : l_ener(2:2) = "-"
493 : END IF
494 : ! The last point would not be local maximum or minimum, as is the first
495 0 : 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 0 : neb_env%number_of_replica, energies(neb_env%number_of_replica), &
498 0 : l_ener, temperatures(neb_env%number_of_replica), TRIM(line)
499 0 : WRITE (output_unit, '(T52,A)') "v Nr."
500 : WRITE (output_unit, '(T2,A,T44,2(1X,I4))') &
501 0 : "NUMBER OF LOCAL MAXIMA (X) and MINIMA (x):", n_max, n_min
502 : ELSE
503 : WRITE (output_unit, '( A,T17,4F16.6)') &
504 289 : ' DISTANCES REP =', distances(1:MIN(4, SIZE(distances)))
505 289 : IF (SIZE(distances) > 4) THEN
506 74 : WRITE (output_unit, '( T17,4F16.6)') distances(5:SIZE(distances))
507 : END IF
508 : WRITE (output_unit, '( A,T17,4F16.6)') &
509 289 : ' ENERGIES [au] =', energies(1:MIN(4, SIZE(energies)))
510 289 : IF (SIZE(energies) > 4) THEN
511 198 : WRITE (output_unit, '( T17,4F16.6)') energies(5:SIZE(energies))
512 : END IF
513 289 : IF (neb_env%opt_type == band_md_opt) THEN
514 : WRITE (output_unit, '( A,T33,4(1X,F11.5))') &
515 88 : ' REPLICA TEMPERATURES (K) =', temperatures(1:MIN(4, SIZE(temperatures)))
516 187 : DO i = 5, SIZE(temperatures), 4
517 : WRITE (output_unit, '( T33,4(1X,F11.5))') &
518 187 : 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 289 : ' BAND TOTAL ENERGY [au] =', SUM(energies(:) + ekin(:)) + &
524 2365 : neb_env%spring_energy
525 289 : WRITE (output_unit, FMT='(A,A)') ' **************************************', &
526 578 : '*****************************************'
527 289 : DEALLOCATE (ekin)
528 1156 : DEALLOCATE (temperatures)
529 : END IF
530 : ! Ener file
531 : ener = cp_print_key_unit_nr(logger, neb_env%neb_section, "ENERGY", &
532 578 : extension=".ener", file_form="FORMATTED")
533 578 : IF (ener > 0) THEN
534 289 : WRITE (line, '(I0)') 2*neb_env%number_of_replica - 1
535 289 : WRITE (ener, '(I10,'//TRIM(line)//'(1X,F20.9))') istep, &
536 578 : energies, distances
537 : END IF
538 : CALL cp_print_key_finished_output(ener, logger, neb_env%neb_section, &
539 578 : "ENERGY")
540 :
541 : ! Dump Restarts
542 578 : 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 578 : vels=vels)
547 578 : CALL cp_rm_default_logger()
548 :
549 578 : CALL timestop(handle)
550 :
551 578 : 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 212 : 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 212 : IF (iw > 0) THEN
574 18 : WRITE (iw, '(/,T2,"NEB|",75("*"))')
575 : WRITE (iw, '(T2,"NEB|",1X,A,I0,A)') &
576 18 : "Geometry for Replica Nr. ", ienum, " in Angstrom"
577 948 : DO iatom = 1, SIZE(particle_set)
578 930 : 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 4668 : TRIM(particle_set(iatom)%atomic_kind%name), r(1:3)*angstrom
581 : END DO
582 18 : IF (use_colvar) THEN
583 10 : WRITE (iw, '(/,T2,"NEB|",1X,A10)') "COLLECTIVE VARIABLES:"
584 : WRITE (iw, '(T2,"NEB|",16X,3F15.9)') &
585 20 : (coords%int(j, i_rep), j=1, SIZE(coords%int(:, :), 1))
586 : END IF
587 18 : WRITE (iw, '(T2,"NEB|",75("*"))')
588 18 : CALL m_flush(iw)
589 : END IF
590 :
591 212 : 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 8376 : 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 2792 : CALL timeset(routineN, handle)
614 : CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
615 2792 : handle=handle2)
616 2792 : logger => cp_get_default_logger()
617 2792 : CALL force_env_get(f_env%force_env, root_section=root_section)
618 2792 : j = irep + (rep_env%local_rep_indices(1) - 1)
619 : ! Get replica_project_name
620 2792 : replica_proj_name = get_replica_project_name(rep_env, n_rep, j)
621 2792 : lp = LEN_TRIM(replica_proj_name)
622 : CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", &
623 2792 : c_val=TRIM(replica_proj_name))
624 2792 : 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 2792 : output_file_path = replica_proj_name(1:lp)//".out"
628 : CALL section_vals_val_set(root_section, "GLOBAL%OUTPUT_FILE_NAME", &
629 2792 : c_val=TRIM(output_file_path))
630 2792 : IF (logger%default_global_unit_nr > 0) THEN
631 2777 : 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 2777 : skip_get_unit_number=.TRUE.)
636 : WRITE (UNIT=logger%default_global_unit_nr, FMT="(/,(T2,A79))") &
637 2777 : "*******************************************************************************", &
638 2777 : "** BAND EVALUATION OF ENERGIES AND FORCES **", &
639 5554 : "*******************************************************************************"
640 2777 : WRITE (UNIT=logger%default_global_unit_nr, FMT="(T2,A,T79,A)") "**", "**"
641 2777 : 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 2777 : "** 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 2777 : "** Band Step Nr. :", istep, "**"
646 : WRITE (UNIT=logger%default_global_unit_nr, FMT="(T2,A79)") &
647 2777 : "*******************************************************************************"
648 : END IF
649 :
650 : ! Handle specific case for mixed_env
651 2822 : SELECT CASE (f_env%force_env%in_use)
652 : CASE (use_mixed_force)
653 2852 : DO i = 1, f_env%force_env%mixed_env%ngroups
654 60 : IF (MODULO(i - 1, f_env%force_env%mixed_env%ngroups) == &
655 30 : f_env%force_env%mixed_env%group_distribution(f_env%force_env%mixed_env%para_env%mepos)) THEN
656 30 : sub_logger => f_env%force_env%mixed_env%sub_logger(i)%p
657 30 : sub_logger%iter_info%project_name = replica_proj_name(1:lp)//"-r-"//TRIM(ADJUSTL(cp_to_string(i)))
658 :
659 30 : unit_nr = sub_logger%default_global_unit_nr
660 30 : IF (unit_nr > 0) THEN
661 30 : CALL close_file(unit_nr)
662 :
663 30 : 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 30 : unit_number=unit_nr, skip_get_unit_number=.TRUE.)
667 : END IF
668 : END IF
669 : END DO
670 : END SELECT
671 :
672 2792 : CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
673 2792 : CPASSERT(ierr == 0)
674 2792 : CALL timestop(handle)
675 :
676 2792 : 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 2916 : 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 2916 : replica_proj_name = rep_env%original_project_name
697 : ! Find padding
698 : ndigits = CEILING(LOG10(REAL(n_rep + 1, KIND=dp))) - &
699 2916 : CEILING(LOG10(REAL(j + 1, KIND=dp)))
700 2916 : padding = ""
701 3618 : DO i = 1, ndigits
702 3618 : padding(i:i) = "0"
703 : END DO
704 2916 : lp = LEN_TRIM(replica_proj_name)
705 : replica_proj_name(lp + 1:LEN(replica_proj_name)) = "-BAND"// &
706 2916 : TRIM(padding)//ADJUSTL(cp_to_string(j))
707 2916 : 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 68 : 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 34 : n_rep_neb = neb_env%number_of_replica
728 34 : n_rep = rep_env%nrep
729 : CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env, &
730 34 : handle=handle2)
731 34 : logger => cp_get_default_logger()
732 34 : output_unit = logger%default_global_unit_nr
733 34 : IF (output_unit > 0) THEN
734 : WRITE (UNIT=output_unit, FMT='(/,(T2,A79))') &
735 33 : "*******************************************************************************", &
736 33 : "** MAPPING OF BAND REPLICA TO REPLICA ENV **", &
737 66 : "*******************************************************************************"
738 : WRITE (UNIT=output_unit, FMT='(T2,A,I6,T32,A,T79,A)') &
739 33 : "** Replica Env Nr.: ", rep_env%local_rep_indices(1) - 1, &
740 66 : "working on the following BAND replicas", "**"
741 : WRITE (UNIT=output_unit, FMT='(T2,A79)') &
742 33 : "** **"
743 : END IF
744 158 : DO irep = 1, n_rep_neb, n_rep
745 124 : replica_proj_name = get_replica_project_name(rep_env, n_rep_neb, irep + rep_env%local_rep_indices(1) - 1)
746 158 : IF (output_unit > 0) THEN
747 : WRITE (UNIT=output_unit, FMT='(T2,A,I6,T32,A,T79,A)') &
748 119 : "** Band Replica Nr.: ", irep + rep_env%local_rep_indices(1) - 1, &
749 238 : "Output available on file: "//TRIM(replica_proj_name)//".out", "**"
750 : END IF
751 : END DO
752 34 : IF (output_unit > 0) THEN
753 : WRITE (UNIT=output_unit, FMT='(T2,A79)') &
754 33 : "** **", &
755 66 : "*******************************************************************************"
756 33 : WRITE (UNIT=output_unit, FMT='(/)')
757 : END IF
758 : ! update runtime info before printing the footer
759 34 : CALL get_runtime_info()
760 : ! print footer
761 34 : CALL cp2k_footer(output_unit)
762 34 : CALL f_env_rm_defaults(f_env=f_env, ierr=ierr, handle=handle2)
763 34 : CPASSERT(ierr == 0)
764 34 : END SUBROUTINE neb_rep_env_map_info
765 :
766 : END MODULE neb_io
|
