From 3d922ff5b668364be98cf18198f79ecffa29ea85 Mon Sep 17 00:00:00 2001 From: HE Zilong <159878975+he-zilong@users.noreply.github.com> Date: Fri, 12 Jun 2026 17:41:51 +0800 Subject: [PATCH] NEB: revamp output (#5382) --- src/input_cp2k_neb.F | 46 ++++++-- src/motion/neb_io.F | 226 +++++++++++++++++++++++++++++++++++---- src/motion/neb_methods.F | 32 ++++-- src/motion/neb_utils.F | 61 ++++++++--- 4 files changed, 319 insertions(+), 46 deletions(-) diff --git a/src/input_cp2k_neb.F b/src/input_cp2k_neb.F index c626dc541b..a84792b1bc 100644 --- a/src/input_cp2k_neb.F +++ b/src/input_cp2k_neb.F @@ -109,7 +109,13 @@ CONTAINS CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="NUMBER_OF_REPLICA", & - description="Specify the number of Replica to use in the BAND", & + description="Specify the number of Replica to use in the BAND. This may "// & + "be equal to or larger than the number of defined &REPLICA sections. If "// & + "larger, the rest of missing replica will automatically be interpolated "// & + "in an iterative bisection procedure: on each step, the largest distance "// & + "between adjacent replica is found and a new replica is inserted there by "// & + "taking the average of adjacent replica; this is repeated until getting "// & + "requested number of replica.", & default_i_val=10) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -263,9 +269,9 @@ CONTAINS CALL section_release(subsection) ! Print key section - CALL cp_print_key_section_create(print_key, __LOCATION__, "program_run_info", & + CALL cp_print_key_section_create(print_key, __LOCATION__, "PROGRAM_RUN_INFO", & description="Controls the printing basic info about the BAND run", & - print_level=medium_print_level, add_last=add_last_numeric, filename="__STD_OUT__") + print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__") CALL keyword_create(keyword, __LOCATION__, name="INITIAL_CONFIGURATION_INFO", & description="Print information for the setup of the initial configuration.", & @@ -274,16 +280,26 @@ CONTAINS CALL section_add_keyword(print_key, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="PLOT_REL_ENERGY", & + description="If a simple plot of relative energy of each replica is shown "// & + "alongside data vertically. This makes the output format less compact than "// & + "default, but creates a visual aid for ease of monitoring with some special "// & + "marks for the local maxima (X) or minima (x) in the plot.", & + usage="PLOT_REL_ENERGY ", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + CALL section_add_subsection(section, print_key) CALL section_release(print_key) - CALL cp_print_key_section_create(print_key, __LOCATION__, "convergence_info", & + CALL cp_print_key_section_create(print_key, __LOCATION__, "CONVERGENCE_INFO", & description="Controls the printing of the convergence criteria during a BAND run", & - print_level=medium_print_level, add_last=add_last_numeric, filename="__STD_OUT__") + print_level=low_print_level, add_last=add_last_numeric, filename="__STD_OUT__") CALL section_add_subsection(section, print_key) CALL section_release(print_key) - CALL cp_print_key_section_create(print_key, __LOCATION__, "replica_info", & + CALL cp_print_key_section_create(print_key, __LOCATION__, "REPLICA_INFO", & description="Controls the printing of each replica info during a BAND run", & print_level=medium_print_level, add_last=add_last_numeric, filename="__STD_OUT__") CALL section_add_subsection(section, print_key) @@ -296,6 +312,24 @@ CONTAINS CALL section_add_subsection(section, print_key) CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, "FINAL_BAND", & + description="Controls the printing of the final structures as an "// & + "XYZ trajectory file after a BAND run, regardless of convergence "// & + "status. The unit is angstrom for coordinates and cell vectors.", & + print_level=low_print_level, common_iter_levels=1, & + filename="FINAL") + + CALL keyword_create(keyword, __LOCATION__, name="PRINT_ATOM_KIND", & + description="Write the atom kind given in the subsys section instead "// & + "of the element symbol in the XYZ trajectory file.", & + usage="PRINT_ATOM_KIND {LOGICAL}", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(print_key, keyword) + CALL keyword_release(keyword) + + CALL section_add_subsection(section, print_key) + CALL section_release(print_key) + CALL cp_print_key_section_create(print_key, __LOCATION__, "BANNER", & description="Controls the printing of the BAND banner", & print_level=low_print_level, common_iter_levels=1, & diff --git a/src/motion/neb_io.F b/src/motion/neb_io.F index 1483233a16..0e23a7bbb3 100644 --- a/src/motion/neb_io.F +++ b/src/motion/neb_io.F @@ -28,6 +28,7 @@ MODULE neb_io cp_rm_default_logger,& cp_to_string USE cp_output_handling, ONLY: cp_print_key_finished_output,& + cp_print_key_generate_filename,& cp_print_key_unit_nr USE cp_units, ONLY: cp_unit_from_cp2k USE f77_interface, ONLY: f_env_add_defaults,& @@ -38,6 +39,7 @@ MODULE neb_io USE header, ONLY: cp2k_footer USE input_constants, ONLY: band_md_opt,& do_sm,& + dump_extxyz,& dump_xmol,& pot_neb_fe,& pot_neb_full,& @@ -75,6 +77,7 @@ MODULE neb_io CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'neb_io' PUBLIC :: read_neb_section, & + dump_neb_final, & dump_neb_info, & dump_replica_coordinates, & handle_band_file_names, & @@ -150,6 +153,124 @@ CONTAINS END IF END SUBROUTINE read_neb_section +! ************************************************************************************************** +!> \brief dump final structures after a NEB run +!> \param neb_env ... +!> \param energies ... +!> \param coords ... +!> \param particle_set ... +!> \param logger ... +!> \param output_unit ... +!> \param converged ... +!> \par +!> History +!> 06.2026 - Created +!> \author HE Zilong +!> \version 1.0 +! ************************************************************************************************** + SUBROUTINE dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged) + TYPE(neb_type), POINTER :: neb_env + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: energies + TYPE(neb_var_type), POINTER :: coords + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(cp_logger_type), POINTER :: logger + INTEGER, INTENT(IN) :: output_unit + LOGICAL :: converged + + CHARACTER(len=*), PARAMETER :: routineN = 'dump_neb_final' + + CHARACTER(LEN=1024) :: cell_str, ener_str, lm_str, record, & + replica_str, title + CHARACTER(LEN=4) :: l_ener + CHARACTER(LEN=5) :: pbc_str + INTEGER :: irep, iw + LOGICAL :: print_kind + REAL(KIND=dp) :: unit_conv + TYPE(cell_type), POINTER :: cell + TYPE(section_vals_type), POINTER :: final_band_section + + NULLIFY (final_band_section) + final_band_section => section_vals_get_subs_vals(neb_env%neb_section, "FINAL_BAND") + CALL force_env_get(neb_env%force_env, cell=cell) ! For now NEB has constant cell + pbc_str = "F F F" + IF (cell%perd(1) == 1) pbc_str(1:1) = "T" + IF (cell%perd(2) == 1) pbc_str(3:3) = "T" + IF (cell%perd(3) == 1) pbc_str(5:5) = "T" + WRITE (UNIT=cell_str, FMT="(9(1X,F19.10))") & + cell%hmat(:, 1)*angstrom, cell%hmat(:, 2)*angstrom, cell%hmat(:, 3)*angstrom + unit_conv = cp_unit_from_cp2k(1.0_dp, "angstrom") + + ! Print a message to log + record = cp_print_key_generate_filename(logger, final_band_section, & + extension=".xyz", & + my_local=.FALSE.) + IF (output_unit > 0) THEN + IF (converged) THEN + WRITE (UNIT=output_unit, FMT="(/,T2,A)") & + routineN//": Band task converged, writing XYZ trajectory gladly:" + ELSE + WRITE (UNIT=output_unit, FMT="(/,T2,A)") & + routineN//": Band task not yet converged, writing XYZ trajectory anyway:" + END IF + WRITE (UNIT=output_unit, FMT="(T3,A)") TRIM(record) + END IF + + ! Write actual trajectory file + iw = cp_print_key_unit_nr(logger, neb_env%neb_section, "FINAL_BAND", & + extension=".xyz", file_form="FORMATTED", file_status="REPLACE") + CALL section_vals_val_get(neb_env%neb_section, "FINAL_BAND%PRINT_ATOM_KIND", & + l_val=print_kind) + DO irep = 1, neb_env%number_of_replica + l_ener = "(**)" + IF (irep > 1) THEN + IF (energies(irep) - energies(irep - 1) > 0) THEN + l_ener(2:2) = "+" + ELSE + l_ener(2:2) = "-" + END IF + END IF + IF (irep < neb_env%number_of_replica) THEN + IF (energies(irep + 1) - energies(irep) < 0) THEN + l_ener(3:3) = "+" + ELSE + l_ener(3:3) = "-" + END IF + END IF + SELECT CASE (l_ener) + CASE ("(++)") ! local maximum + WRITE (lm_str, '(A)') "Ener_loc_max=T Ener_loc_min=F" + CASE ("(--)") ! local minimum + WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=T" + CASE DEFAULT + WRITE (lm_str, '(A)') "Ener_loc_max=F Ener_loc_min=F" + END SELECT + WRITE (UNIT=replica_str, FMT="(I8)") irep + WRITE (UNIT=ener_str, FMT="(F20.10)") energies(irep) + WRITE (UNIT=title, FMT="(A)") & + 'Lattice="'//TRIM(ADJUSTL(cell_str))//'" '// & + 'Properties=species:S:1:pos:R:3 '// & + 'pbc="'//pbc_str//'" '// & + 'Replica='//TRIM(ADJUSTL(replica_str))//' '// & + 'Energy='//TRIM(ADJUSTL(ener_str))//' '// & + TRIM(ADJUSTL(lm_str)) + IF (iw > 0) THEN + ! The iw condition does not hold for certain ranks/processes + ! that write to -BAND.out where n > neb_env%number_of_replica + CALL write_particle_coordinates(particle_set, iw, dump_extxyz, "POS", title, & + cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, & + print_kind=print_kind) + CALL m_flush(iw) + END IF + END DO + + IF (output_unit > 0) & + WRITE (UNIT=output_unit, FMT='(/,T2,A)') & + routineN//": Done!" + + CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, "FINAL_BAND") + + END SUBROUTINE dump_neb_final + ! ************************************************************************************************** !> \brief dump print info of a NEB run !> \param neb_env ... @@ -178,17 +299,20 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'dump_neb_info' CHARACTER(LEN=20) :: mytype + CHARACTER(LEN=4) :: l_ener CHARACTER(LEN=default_string_length) :: line, title, unit_str - INTEGER :: crd, ener, frc, handle, i, irep, ndig, & - ndigl, ttst, vel - LOGICAL :: explicit, lval, print_kind - REAL(KIND=dp) :: f_ann, tmp_r1, unit_conv + INTEGER :: crd, ener, frc, handle, i, irep, n_max, & + n_min, ndig, ndigl, plt, ttst, vel + LOGICAL :: explicit, lval, plot_rel_energy, & + print_kind + REAL(KIND=dp) :: ener_min, ener_range, f_ann, tmp_r1, & + unit_conv REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ekin, temperatures TYPE(cell_type), POINTER :: cell TYPE(enumeration_type), POINTER :: enum TYPE(keyword_type), POINTER :: keyword TYPE(section_type), POINTER :: section - TYPE(section_vals_type), POINTER :: tc_section, vc_section + TYPE(section_vals_type), POINTER :: run_info_section, tc_section, vc_section CALL timeset(routineN, handle) ndig = CEILING(LOG10(REAL(neb_env%number_of_replica + 1, KIND=dp))) @@ -264,6 +388,8 @@ CONTAINS IF (output_unit > 0) THEN tc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%TEMP_CONTROL") vc_section => section_vals_get_subs_vals(neb_env%neb_section, "OPTIMIZE_BAND%MD%VEL_CONTROL") + run_info_section => section_vals_get_subs_vals(neb_env%neb_section, "PROGRAM_RUN_INFO") + CALL section_vals_val_get(run_info_section, "PLOT_REL_ENERGY", l_val=plot_rel_energy) ALLOCATE (temperatures(neb_env%number_of_replica)) ALLOCATE (ekin(neb_env%number_of_replica)) CALL get_temperatures(vels, particle_set, temperatures, ekin=ekin) @@ -312,23 +438,81 @@ CONTAINS END IF WRITE (output_unit, '( A,T71,I10 )') & ' NUMBER OF NEB REPLICA =', neb_env%number_of_replica - WRITE (output_unit, '( A,T17,4F16.6)') & - ' DISTANCES REP =', distances(1:MIN(4, SIZE(distances))) - IF (SIZE(distances) > 4) THEN - WRITE (output_unit, '( T17,4F16.6)') distances(5:SIZE(distances)) - END IF - WRITE (output_unit, '( A,T17,4F16.6)') & - ' ENERGIES [au] =', energies(1:MIN(4, SIZE(energies))) - IF (SIZE(energies) > 4) THEN - WRITE (output_unit, '( T17,4F16.6)') energies(5:SIZE(energies)) - END IF - IF (neb_env%opt_type == band_md_opt) THEN - WRITE (output_unit, '( A,T33,4(1X,F11.5))') & - ' REPLICA TEMPERATURES (K) =', temperatures(1:MIN(4, SIZE(temperatures))) - DO i = 5, SIZE(temperatures), 4 - WRITE (output_unit, '( T33,4(1X,F11.5))') & - temperatures(i:MIN(i + 3, SIZE(temperatures))) + ! switch between a longer visual format and a compact data-only print format + IF (plot_rel_energy) THEN + CPASSERT(SIZE(distances) == neb_env%number_of_replica - 1) + CPASSERT(SIZE(energies) == neb_env%number_of_replica) + CPASSERT(SIZE(temperatures) == neb_env%number_of_replica) + ener_min = MINVAL(energies(:)) + ener_range = MAXVAL(energies(:)) - ener_min + n_max = 0 + n_min = 0 + WRITE (output_unit, '(T2,A,T22,A,T35,A,T52,A)') & + 'REPLICA', 'ENERGY [au]', 'TEMPERATURE [K]', 'o-------------------------> E' + DO i = 1, SIZE(distances) + plt = FLOOR((energies(i) - ener_min)/ener_range*25) + l_ener = "(**)" + IF (i > 1) THEN + IF (energies(i) - energies(i - 1) > 0) THEN + l_ener(2:2) = "+" + ELSE + l_ener(2:2) = "-" + END IF + END IF + IF (energies(i + 1) - energies(i) < 0) THEN + l_ener(3:3) = "+" + ELSE + l_ener(3:3) = "-" + END IF + SELECT CASE (l_ener) + CASE ("(++)") ! local maximum + n_max = n_max + 1 + WRITE (line, '(A,A,A)') "|", REPEAT(" ", plt), "X" + CASE ("(--)") ! local minimum + n_min = n_min + 1 + WRITE (line, '(A,A,A)') "|", REPEAT(" ", plt), "x" + CASE DEFAULT + WRITE (line, '(A,A,A)') "|", REPEAT(" ", plt), "O" + END SELECT + WRITE (output_unit, '(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') & + i, energies(i), l_ener, temperatures(i), TRIM(line) + WRITE (output_unit, '(T2,A,1X,F16.6,T52,A)') & + "DISTANCE = ", distances(i), "|" END DO + plt = FLOOR((energies(neb_env%number_of_replica) - ener_min)/ener_range*25) + l_ener = "(**)" + IF (energies(neb_env%number_of_replica) - energies(neb_env%number_of_replica - 1) > 0) THEN + l_ener(2:2) = "+" + ELSE + l_ener(2:2) = "-" + END IF + ! The last point would not be local maximum or minimum, as is the first + WRITE (line, '(A,A,A)') "|", REPEAT(" ", plt), "O" + WRITE (output_unit, '(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') & + neb_env%number_of_replica, energies(neb_env%number_of_replica), & + l_ener, temperatures(neb_env%number_of_replica), TRIM(line) + WRITE (output_unit, '(T52,A)') "v Nr." + WRITE (output_unit, '(T2,A,T44,2(1X,I4))') & + "NUMBER OF LOCAL MAXIMA (X) and MINIMA (x):", n_max, n_min + ELSE + WRITE (output_unit, '( A,T17,4F16.6)') & + ' DISTANCES REP =', distances(1:MIN(4, SIZE(distances))) + IF (SIZE(distances) > 4) THEN + WRITE (output_unit, '( T17,4F16.6)') distances(5:SIZE(distances)) + END IF + WRITE (output_unit, '( A,T17,4F16.6)') & + ' ENERGIES [au] =', energies(1:MIN(4, SIZE(energies))) + IF (SIZE(energies) > 4) THEN + WRITE (output_unit, '( T17,4F16.6)') energies(5:SIZE(energies)) + END IF + IF (neb_env%opt_type == band_md_opt) THEN + WRITE (output_unit, '( A,T33,4(1X,F11.5))') & + ' REPLICA TEMPERATURES (K) =', temperatures(1:MIN(4, SIZE(temperatures))) + DO i = 5, SIZE(temperatures), 4 + WRITE (output_unit, '( T33,4(1X,F11.5))') & + temperatures(i:MIN(i + 3, SIZE(temperatures))) + END DO + END IF END IF WRITE (output_unit, '( A,T56,F25.14)') & ' BAND TOTAL ENERGY [au] =', SUM(energies(:) + ekin(:)) + & diff --git a/src/motion/neb_methods.F b/src/motion/neb_methods.F index 9631fbe689..9af2b50a23 100644 --- a/src/motion/neb_methods.F +++ b/src/motion/neb_methods.F @@ -48,7 +48,8 @@ MODULE neb_methods section_vals_val_get USE kinds, ONLY: dp USE message_passing, ONLY: mp_para_env_type - USE neb_io, ONLY: dump_neb_info,& + USE neb_io, ONLY: dump_neb_final,& + dump_neb_info,& neb_rep_env_map_info,& read_neb_section USE neb_md_utils, ONLY: control_vels_a,& @@ -334,10 +335,18 @@ CONTAINS energies=energies, & distances=distances, & output_unit=output_unit) - converged = check_convergence(neb_env, Dcoords, forces) - IF (converged) EXIT md_opt_loop + converged = check_convergence(neb_env, Dcoords, forces, logger) + IF (converged) THEN + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79) + WRITE (UNIT=output_unit, FMT="(T2,A,T32,A,T78,A)") & + "***", "BAND TASK COMPLETED", "***" + WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79) + END IF + EXIT md_opt_loop + END IF END DO md_opt_loop - + CALL dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged) DEALLOCATE (mass) DEALLOCATE (energies) DEALLOCATE (distances) @@ -494,15 +503,26 @@ CONTAINS vels=vels, & output_unit=output_unit) - converged = check_convergence(neb_env, sline, forces) - IF (converged) EXIT + converged = check_convergence(neb_env, sline, forces, logger) + IF (converged) THEN + IF (output_unit > 0) THEN + WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79) + WRITE (UNIT=output_unit, FMT="(T2,A,T32,A,T78,A)") & + "***", "BAND TASK COMPLETED", "***" + WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79) + END IF + EXIT + END IF END DO + CALL dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged) DEALLOCATE (energies) DEALLOCATE (distances) DEALLOCATE (err) DEALLOCATE (crr) DEALLOCATE (set_err) CALL neb_var_release(sline) + CALL cp_print_key_finished_output(iw, logger, neb_env%neb_section, & + "REPLICA_INFO") CALL timestop(handle) END SUBROUTINE neb_diis diff --git a/src/motion/neb_utils.F b/src/motion/neb_utils.F index 15b0659016..8a0c83ea53 100644 --- a/src/motion/neb_utils.F +++ b/src/motion/neb_utils.F @@ -148,8 +148,9 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'build_replica_coords' CHARACTER(LEN=default_path_length) :: filename - INTEGER :: handle, i_rep, iatom, ic, input_nr_replica, is, ivar, j, jtarg, n_rep, natom, & + INTEGER :: handle, i_rep, iatom, ic, input_nr_replica, is, ivar, j, jtarg, k, n_rep, natom, & neb_nr_replica, nr_replica_to_interpolate, nval, nvar, shell_index + INTEGER, ALLOCATABLE, DIMENSION(:) :: rep_map LOGICAL :: check, explicit, skip_vel_section REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: distance REAL(KIND=dp), DIMENSION(3) :: r @@ -181,10 +182,14 @@ CONTAINS END DO ! Setup cartesian coordinates and COLVAR (if requested) coords%xyz(:, :) = 0.0_dp + ! Mapping between input replica and actual replica + ALLOCATE (rep_map(neb_nr_replica)) + rep_map(:) = 0 DO i_rep = 1, input_nr_replica coord_section => section_vals_get_subs_vals(replica_section, "COORD", & i_rep_section=i_rep) CALL section_vals_get(coord_section, explicit=explicit) + rep_map(i_rep) = i_rep ! Cartesian Coordinates IF (explicit) THEN CALL section_vals_val_get(coord_section, "_DEFAULT_KEYWORD_", & @@ -293,12 +298,23 @@ CONTAINS END IF END DO ! i_rep ALLOCATE (distance(neb_nr_replica - 1)) + IF (iw > 0) THEN + WRITE (iw, '(T2,A)') 'NEB| Mapping between input and requested replica so far' + WRITE (iw, '(T2,A)') 'NEB| 1, 2, ... = input replica in sections order,' + WRITE (iw, '(T2,A)') 'NEB| -1, -2, ... = added replica in insertion order (if any),' + WRITE (iw, '(T2,A)') 'NEB| 0 = yet to be completed replica (if any)' + DO j = 1, neb_nr_replica, 8 + WRITE (iw, '(T2,A,T9,8(1X,I8))') 'NEB|', rep_map(j:MIN(j + 7, neb_nr_replica)) + END DO + END IF IF (input_nr_replica < neb_nr_replica) THEN ! Interpolate missing replicas nr_replica_to_interpolate = neb_nr_replica - input_nr_replica + k = 0 distance = 0.0_dp IF (iw > 0) THEN - WRITE (iw, '(T2,A,I0,A)') 'NEB| Interpolating ', nr_replica_to_interpolate, ' missing Replica.' + WRITE (iw, '(T2,A,I0,A)') 'NEB| Interpolating ', nr_replica_to_interpolate, ' missing replica '// & + 'by stepwise bisection of distance.' END IF DO WHILE (nr_replica_to_interpolate > 0) ! Compute distances between known images to find the interval @@ -309,12 +325,15 @@ CONTAINS END DO jtarg = MAXLOC(distance(1:input_nr_replica), 1) IF (iw > 0) THEN - WRITE (iw, '(T2,3(A,I0),A)') 'NEB| Interpolating Nr. ', & - nr_replica_to_interpolate, ' missing Replica between Replica Nr. ', & + WRITE (iw, '(/,T2,3(A,I0),A)') 'NEB| Interpolating Nr. ', & + nr_replica_to_interpolate, ' missing Replica; next between Replica Nr. ', & jtarg, ' and ', jtarg + 1, '.' END IF input_nr_replica = input_nr_replica + 1 nr_replica_to_interpolate = nr_replica_to_interpolate - 1 + k = k + 1 + rep_map(jtarg + 2:input_nr_replica) = rep_map(jtarg + 1:input_nr_replica - 1) + IF (jtarg + 1 <= neb_nr_replica .AND. jtarg + 1 >= 1) rep_map(jtarg + 1) = -k ! Interpolation is a simple bisection in XYZ coords%xyz(:, jtarg + 2:input_nr_replica) = coords%xyz(:, jtarg + 1:input_nr_replica - 1) coords%xyz(:, jtarg + 1) = (coords%xyz(:, jtarg) + coords%xyz(:, jtarg + 2))/2.0_dp @@ -332,6 +351,12 @@ CONTAINS input_nr_replica, iw, neb_env%use_colvar) CALL neb_initialize_velocity(vels%wrk, neb_section, particle_set, & jtarg + 1, iw, globenv, neb_env) + IF (iw > 0) THEN + WRITE (iw, '(T2,A)') 'NEB| Mapping between input and requested replica so far' + DO j = 1, neb_nr_replica, 8 + WRITE (iw, '(T2,A,T9,8(1X,I8))') 'NEB|', rep_map(j:MIN(j + 7, neb_nr_replica)) + END DO + END IF END DO END IF vels%wrk(:, 1) = 0.0_dp @@ -344,7 +369,17 @@ CONTAINS CALL neb_replica_distance(particle_set, coords, j, j + 1, distance(j), iw, & rotate=neb_env%align_frames) END DO + ! If there are still large distances, hint that more replica may be requested + IF (MAXVAL(distance)/MINVAL(distance) > 1.5_dp) THEN + IF (iw > 0) THEN + WRITE (iw, '(/,T2,A)') 'NEB| After interpolating replica, the maximum distance between adjacent replica' + WRITE (iw, '(T2,A)') 'NEB| is still larger than 1.5 times the minimum distance. In order to make the' + WRITE (iw, '(T2,A)') 'NEB| distribution of replica more uniform, consider providing more structures' + WRITE (iw, '(T2,A)') 'NEB| with &BAND/&REPLICA sections or increasing &BAND/NUMBER_OF_REPLICA value.' + END IF + END IF DEALLOCATE (distance) + DEALLOCATE (rep_map) CALL timestop(handle) @@ -987,12 +1022,14 @@ CONTAINS !> \param neb_env ... !> \param Dcoords ... !> \param forces ... +!> \param logger ... !> \return ... !> \author Teodoro Laino 10.2006 ! ************************************************************************************************** - FUNCTION check_convergence(neb_env, Dcoords, forces) RESULT(converged) + FUNCTION check_convergence(neb_env, Dcoords, forces, logger) RESULT(converged) TYPE(neb_type), POINTER :: neb_env TYPE(neb_var_type), POINTER :: Dcoords, forces + TYPE(cp_logger_type), POINTER :: logger LOGICAL :: converged CHARACTER(LEN=3), DIMENSION(4) :: labels @@ -1000,11 +1037,9 @@ CONTAINS REAL(KIND=dp) :: max_dr, max_force, my_max_dr, & my_max_force, my_rms_dr, my_rms_force, & rms_dr, rms_force - TYPE(cp_logger_type), POINTER :: logger TYPE(section_vals_type), POINTER :: cc_section - NULLIFY (logger, cc_section) - logger => cp_get_default_logger() + NULLIFY (cc_section) cc_section => section_vals_get_subs_vals(neb_env%neb_section, "CONVERGENCE_CONTROL") CALL section_vals_val_get(cc_section, "MAX_DR", r_val=max_dr) CALL section_vals_val_get(cc_section, "MAX_FORCE", r_val=max_force) @@ -1028,11 +1063,11 @@ CONTAINS ! Print convergence info WRITE (iw, FMT='(A,A)') ' **************************************', & '*****************************************' - WRITE (iw, FMT='(1X,A,2X,F8.5,5X,"[",F8.5,"]",1X,T76,"(",A,")")') & - 'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), & - 'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), & - 'RMS FORCE =', my_rms_force, rms_force, labels(4), & - 'MAX FORCE =', my_max_force, max_force, labels(2) + WRITE (iw, FMT='(1X,A,2X,F16.10,5X,"[",F16.10,1X,"]",T76,"(",A,")")') & + 'RMS DISPLACEMENT =', my_rms_dr, rms_dr, labels(3), & + 'MAX DISPLACEMENT =', my_max_dr, max_dr, labels(1), & + 'RMS FORCE =', my_rms_force, rms_force, labels(4), & + 'MAX FORCE =', my_max_force, max_force, labels(2) WRITE (iw, FMT='(A,A)') ' **************************************', & '*****************************************' END IF