73#include "../base/base_uses.f90"
77 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'neb_io'
101 cpassert(
ASSOCIATED(neb_env))
108 CALL section_vals_val_get(neb_section,
"OPTIMIZE_BAND%OPTIMIZE_END_POINTS", l_val=neb_env%optimize_end_points)
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.")
125 SELECT CASE (neb_env%pot_type)
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.")
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.")
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!")
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.
171 SUBROUTINE dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged)
173 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: energies
177 INTEGER,
INTENT(IN) :: output_unit
180 CHARACTER(len=*),
PARAMETER :: routinen =
'dump_neb_final'
182 CHARACTER(LEN=1024) :: cell_str, ener_str, lm_str, record, &
184 CHARACTER(LEN=4) :: l_ener
185 CHARACTER(LEN=5) :: pbc_str
187 LOGICAL :: print_kind
188 REAL(kind=
dp) :: unit_conv
192 NULLIFY (final_band_section)
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))") &
207 IF (output_unit > 0)
THEN
209 WRITE (unit=output_unit, fmt=
"(/,T2,A)") &
210 routinen//
": Band task converged, writing XYZ trajectory gladly:"
212 WRITE (unit=output_unit, fmt=
"(/,T2,A)") &
213 routinen//
": Band task not yet converged, writing XYZ trajectory anyway:"
215 WRITE (unit=output_unit, fmt=
"(T3,A)") trim(record)
220 extension=
".xyz", file_form=
"FORMATTED", file_status=
"REPLACE")
223 DO irep = 1, neb_env%number_of_replica
226 IF (energies(irep) - energies(irep - 1) > 0)
THEN
232 IF (irep < neb_env%number_of_replica)
THEN
233 IF (energies(irep + 1) - energies(irep) < 0)
THEN
241 WRITE (lm_str,
'(A)')
"Ener_loc_max=T Ener_loc_min=F"
243 WRITE (lm_str,
'(A)')
"Ener_loc_max=F Ener_loc_min=T"
245 WRITE (lm_str,
'(A)')
"Ener_loc_max=F Ener_loc_min=F"
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))
260 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
261 print_kind=print_kind)
266 IF (output_unit > 0) &
267 WRITE (unit=output_unit, fmt=
'(/,T2,A)') &
288 SUBROUTINE dump_neb_info(neb_env, coords, vels, forces, particle_set, logger, &
289 istep, energies, distances, output_unit)
295 INTEGER,
INTENT(IN) :: istep
296 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: energies, distances
297 INTEGER,
INTENT(IN) :: output_unit
299 CHARACTER(len=*),
PARAMETER :: routinen =
'dump_neb_info'
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, &
308 REAL(kind=
dp) :: ener_min, ener_range, f_ann, tmp_r1, &
310 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: ekin, temperatures
317 CALL timeset(routinen, handle)
318 ndig = ceiling(log10(real(neb_env%number_of_replica + 1, kind=
dp)))
320 DO irep = 1, neb_env%number_of_replica
321 ndigl = ceiling(log10(real(irep + 1, kind=
dp)))
324 extension=
".xyz", file_form=
"FORMATTED", middle_name=
"pos-"//trim(line))
325 IF (
PRESENT(vels))
THEN
327 extension=
".xyz", file_form=
"FORMATTED", middle_name=
"vel-"//trim(line))
329 IF (
PRESENT(forces))
THEN
331 extension=
".xyz", file_form=
"FORMATTED", middle_name=
"force-"//trim(line))
342 WRITE (unit=title, fmt=
"(A,I8,A,F20.10)")
" i =", istep,
", E =", energies(irep)
344 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
345 print_kind=print_kind)
349 IF (vel > 0 .AND.
PRESENT(vels))
THEN
356 WRITE (unit=title, fmt=
"(A,I8,A,F20.10)")
" i =", istep,
", E =", energies(irep)
358 cell=cell, array=vels%xyz(:, irep), unit_conv=unit_conv, &
359 print_kind=print_kind)
363 IF (frc > 0 .AND.
PRESENT(forces))
THEN
370 WRITE (unit=title, fmt=
"(A,I8,A,F20.10)")
" i =", istep,
", E =", energies(irep)
372 cell=cell, array=forces%xyz(:, irep), unit_conv=unit_conv, &
373 print_kind=print_kind)
378 IF (
PRESENT(vels))
THEN
382 IF (
PRESENT(forces))
THEN
388 IF (output_unit > 0)
THEN
393 ALLOCATE (temperatures(neb_env%number_of_replica))
394 ALLOCATE (ekin(neb_env%number_of_replica))
396 WRITE (output_unit,
'(/)', advance=
"NO")
397 WRITE (output_unit, fmt=
'(A,A)')
' **************************************', &
398 '*****************************************'
399 NULLIFY (section, keyword, enum)
403 mytype = trim(
enum_i2c(enum, neb_env%id_type))
404 WRITE (output_unit, fmt=
'(A,T61,A)') &
405 ' BAND TYPE =', adjustr(mytype)
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
417 IF (lval)
WRITE (output_unit,
'( A,T71,L10 )') &
418 ' PROJECTED VELOCITY VERLET =', lval
420 IF (lval)
WRITE (output_unit,
'( A,T71,L10)') &
421 ' STEEPEST DESCENT LIKE =', lval
423 IF (f_ann /= 1.0_dp)
THEN
424 WRITE (output_unit,
'( A,T71,F10.5)') &
425 ' ANNEALING FACTOR = ', f_ann
432 IF (istep <= ttst)
THEN
435 WRITE (output_unit,
'( A,T71,F10.5)') &
436 ' TEMPERATURE TARGET =', tmp_r1
439 WRITE (output_unit,
'( A,T71,I10 )') &
440 ' NUMBER OF NEB REPLICA =', neb_env%number_of_replica
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
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)
456 IF (energies(i) - energies(i - 1) > 0)
THEN
462 IF (energies(i + 1) - energies(i) < 0)
THEN
470 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"X"
473 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"x"
475 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"O"
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),
"|"
482 plt = floor((energies(neb_env%number_of_replica) - ener_min)/ener_range*25)
484 IF (energies(neb_env%number_of_replica) - energies(neb_env%number_of_replica - 1) > 0)
THEN
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
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))
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))
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)))
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 '*****************************************'
523 DEALLOCATE (temperatures)
527 extension=
".ener", file_form=
"FORMATTED")
529 WRITE (line,
'(I0)') 2*neb_env%number_of_replica - 1
530 WRITE (ener,
'(I10,'//trim(line)//
'(1X,F20.9))') istep, &
539 root_section=neb_env%root_section, &
544 CALL timestop(handle)
562 INTEGER,
INTENT(IN) :: i_rep, ienum, iw
563 LOGICAL,
INTENT(IN) :: use_colvar
566 REAL(kind=
dp),
DIMENSION(3) :: r
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)
574 WRITE (iw,
'(T2,"NEB|",1X,A10,5X,3F15.9)') &
575 trim(particle_set(iatom)%atomic_kind%name), r(1:3)*
angstrom
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))
582 WRITE (iw,
'(T2,"NEB|",75("*"))')
598 INTEGER,
INTENT(IN) :: irep, n_rep, istep
600 CHARACTER(len=*),
PARAMETER :: routinen =
'handle_band_file_names'
602 CHARACTER(LEN=default_path_length) :: output_file_path, replica_proj_name
603 INTEGER :: handle, handle2, i, ierr, j, lp, unit_nr
608 CALL timeset(routinen, handle)
612 CALL force_env_get(f_env%force_env, root_section=root_section)
613 j = irep + (rep_env%local_rep_indices(1) - 1)
615 replica_proj_name = get_replica_project_name(rep_env, n_rep, j)
616 lp = len_trim(replica_proj_name)
618 c_val=trim(replica_proj_name))
619 logger%iter_info%project_name = trim(replica_proj_name)
622 output_file_path = replica_proj_name(1:lp)//
".out"
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 "*******************************************************************************"
646 SELECT CASE (f_env%force_env%in_use)
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)))
654 unit_nr = sub_logger%default_global_unit_nr
655 IF (unit_nr > 0)
THEN
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.)
669 CALL timestop(handle)
681 FUNCTION get_replica_project_name(rep_env, n_rep, j)
RESULT(replica_proj_name)
683 INTEGER,
INTENT(IN) :: n_rep, j
684 CHARACTER(LEN=default_path_length) :: replica_proj_name
686 CHARACTER(LEN=default_string_length) :: padding
687 INTEGER :: i, lp, ndigits
691 replica_proj_name = rep_env%original_project_name
693 ndigits = ceiling(log10(real(n_rep + 1, kind=
dp))) - &
694 ceiling(log10(real(j + 1, kind=
dp)))
699 lp = len_trim(replica_proj_name)
700 replica_proj_name(lp + 1:len(replica_proj_name)) =
"-BAND"// &
702 END FUNCTION get_replica_project_name
716 CHARACTER(LEN=default_path_length) :: replica_proj_name
717 INTEGER :: handle2, ierr, irep, n_rep, n_rep_neb, &
722 n_rep_neb = neb_env%number_of_replica
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)') &
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",
"**"
747 IF (output_unit > 0)
THEN
748 WRITE (unit=output_unit, fmt=
'(T2,A79)') &
750 "*******************************************************************************"
751 WRITE (unit=output_unit, fmt=
'(/)')
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.
some minimal info about CP2K, including its version and license
subroutine, public get_runtime_info()
...
Utility routines to open and close files. Tracking of preconnections.
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.
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.
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,...
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
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
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
integer, parameter, public default_path_length
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
I/O Module for Nudged Elastic Band Calculation.
subroutine, public read_neb_section(neb_env, neb_section)
Read data from the NEB input section.
subroutine, public dump_replica_coordinates(particle_set, coords, i_rep, ienum, iw, use_colvar)
dump coordinates of a replica NEB
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...
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
subroutine, public handle_band_file_names(rep_env, irep, n_rep, istep)
Handles the correct file names during a band calculation.
subroutine, public dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged)
dump final structures after a NEB run
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.
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:
real(kind=dp), parameter, public angstrom
types used to handle many replica of the same system that differ only in atom positions,...
Type defining parameters related to the simulation cell.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
keeps replicated information about the replicas