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)
THEN
120 CALL cp_abort(__location__, &
121 "A potential energy function based on free energy or minimum energy"// &
122 " was requested without enabling the usage of COLVARS. Both methods"// &
123 " are based on COLVARS definition.")
126 SELECT CASE (neb_env%pot_type)
130 IF (.NOT. explicit)
THEN
131 CALL cp_abort(__location__, &
132 "A free energy BAND (colvars projected) calculation is requested"// &
133 " but NONE MD section was defined in the input.")
138 IF (.NOT. explicit)
THEN
139 CALL cp_abort(__location__, &
140 "A minimum energy BAND (colvars projected) calculation is requested"// &
141 " but NONE GEO_OPT section was defined in the input.")
145 IF (neb_env%use_colvar)
THEN
146 CALL cp_abort(__location__, &
147 "A band calculation was requested with a full potential energy. USE_COLVAR cannot"// &
148 " be set for this kind of calculation!")
153 CALL section_vals_val_get(neb_section,
"STRING_METHOD%SPLINE_ORDER", i_val=neb_env%spline_order)
154 neb_env%reparametrize_frames = .false.
155 IF (neb_env%id_type ==
do_sm)
THEN
156 neb_env%reparametrize_frames = .true.
175 SUBROUTINE dump_neb_final(neb_env, energies, coords, particle_set, logger, output_unit, converged)
177 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: energies
181 INTEGER,
INTENT(IN) :: output_unit
184 CHARACTER(len=*),
PARAMETER :: routinen =
'dump_neb_final'
186 CHARACTER(LEN=1024) :: cell_str, ener_str, lm_str, record, &
188 CHARACTER(LEN=4) :: l_ener
189 CHARACTER(LEN=5) :: pbc_str
191 LOGICAL :: print_kind
192 REAL(kind=
dp) :: unit_conv
196 NULLIFY (final_band_section)
200 IF (cell%perd(1) == 1) pbc_str(1:1) =
"T"
201 IF (cell%perd(2) == 1) pbc_str(3:3) =
"T"
202 IF (cell%perd(3) == 1) pbc_str(5:5) =
"T"
203 WRITE (unit=cell_str, fmt=
"(9(1X,F19.10))") &
211 IF (output_unit > 0)
THEN
213 WRITE (unit=output_unit, fmt=
"(/,T2,A)") &
214 routinen//
": Band task converged, writing XYZ trajectory gladly:"
216 WRITE (unit=output_unit, fmt=
"(/,T2,A)") &
217 routinen//
": Band task not yet converged, writing XYZ trajectory anyway:"
219 WRITE (unit=output_unit, fmt=
"(T3,A)") trim(record)
224 extension=
".xyz", file_form=
"FORMATTED", file_status=
"REPLACE")
227 DO irep = 1, neb_env%number_of_replica
230 IF (energies(irep) - energies(irep - 1) > 0)
THEN
236 IF (irep < neb_env%number_of_replica)
THEN
237 IF (energies(irep + 1) - energies(irep) < 0)
THEN
245 WRITE (lm_str,
'(A)')
"Ener_loc_max=T Ener_loc_min=F"
247 WRITE (lm_str,
'(A)')
"Ener_loc_max=F Ener_loc_min=T"
249 WRITE (lm_str,
'(A)')
"Ener_loc_max=F Ener_loc_min=F"
251 WRITE (unit=replica_str, fmt=
"(I8)") irep
252 WRITE (unit=ener_str, fmt=
"(F20.10)") energies(irep)
253 WRITE (unit=title, fmt=
"(A)") &
254 'Lattice="'//trim(adjustl(cell_str))//
'" '// &
255 'Properties=species:S:1:pos:R:3 '// &
256 'pbc="'//pbc_str//
'" '// &
257 'Replica='//trim(adjustl(replica_str))//
' '// &
258 'Energy='//trim(adjustl(ener_str))//
' '// &
259 trim(adjustl(lm_str))
264 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
265 print_kind=print_kind)
270 IF (output_unit > 0)
THEN
271 WRITE (unit=output_unit, fmt=
'(/,T2,A)') &
293 SUBROUTINE dump_neb_info(neb_env, coords, vels, forces, particle_set, logger, &
294 istep, energies, distances, output_unit)
300 INTEGER,
INTENT(IN) :: istep
301 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: energies, distances
302 INTEGER,
INTENT(IN) :: output_unit
304 CHARACTER(len=*),
PARAMETER :: routinen =
'dump_neb_info'
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, &
313 REAL(kind=
dp) :: ener_min, ener_range, f_ann, tmp_r1, &
315 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: ekin, temperatures
322 CALL timeset(routinen, handle)
323 ndig = ceiling(log10(real(neb_env%number_of_replica + 1, kind=
dp)))
325 DO irep = 1, neb_env%number_of_replica
326 ndigl = ceiling(log10(real(irep + 1, kind=
dp)))
329 extension=
".xyz", file_form=
"FORMATTED", middle_name=
"pos-"//trim(line))
330 IF (
PRESENT(vels))
THEN
332 extension=
".xyz", file_form=
"FORMATTED", middle_name=
"vel-"//trim(line))
334 IF (
PRESENT(forces))
THEN
336 extension=
".xyz", file_form=
"FORMATTED", middle_name=
"force-"//trim(line))
347 WRITE (unit=title, fmt=
"(A,I8,A,F20.10)")
" i =", istep,
", E =", energies(irep)
349 cell=cell, array=coords%xyz(:, irep), unit_conv=unit_conv, &
350 print_kind=print_kind)
354 IF (vel > 0 .AND.
PRESENT(vels))
THEN
361 WRITE (unit=title, fmt=
"(A,I8,A,F20.10)")
" i =", istep,
", E =", energies(irep)
363 cell=cell, array=vels%xyz(:, irep), unit_conv=unit_conv, &
364 print_kind=print_kind)
368 IF (frc > 0 .AND.
PRESENT(forces))
THEN
375 WRITE (unit=title, fmt=
"(A,I8,A,F20.10)")
" i =", istep,
", E =", energies(irep)
377 cell=cell, array=forces%xyz(:, irep), unit_conv=unit_conv, &
378 print_kind=print_kind)
383 IF (
PRESENT(vels))
THEN
387 IF (
PRESENT(forces))
THEN
393 IF (output_unit > 0)
THEN
398 ALLOCATE (temperatures(neb_env%number_of_replica))
399 ALLOCATE (ekin(neb_env%number_of_replica))
401 WRITE (output_unit,
'(/)', advance=
"NO")
402 WRITE (output_unit, fmt=
'(A,A)')
' **************************************', &
403 '*****************************************'
404 NULLIFY (section, keyword, enum)
408 mytype = trim(
enum_i2c(enum, neb_env%id_type))
409 WRITE (output_unit, fmt=
'(A,T61,A)') &
410 ' BAND TYPE =', adjustr(mytype)
412 WRITE (output_unit, fmt=
'(A,T61,A)') &
413 ' BAND TYPE OPTIMIZATION =', adjustr(neb_env%opt_type_label(1:20))
414 WRITE (output_unit,
'( A,T71,I10 )') &
415 ' STEP NUMBER =', istep
416 IF (neb_env%rotate_frames)
WRITE (output_unit,
'( A,T71,L10 )') &
417 ' RMSD DISTANCE DEFINITION =', neb_env%rotate_frames
422 IF (lval)
WRITE (output_unit,
'( A,T71,L10 )') &
423 ' PROJECTED VELOCITY VERLET =', lval
425 IF (lval)
WRITE (output_unit,
'( A,T71,L10)') &
426 ' STEEPEST DESCENT LIKE =', lval
428 IF (f_ann /= 1.0_dp)
THEN
429 WRITE (output_unit,
'( A,T71,F10.5)') &
430 ' ANNEALING FACTOR = ', f_ann
437 IF (istep <= ttst)
THEN
440 WRITE (output_unit,
'( A,T71,F10.5)') &
441 ' TEMPERATURE TARGET =', tmp_r1
444 WRITE (output_unit,
'( A,T71,I10 )') &
445 ' NUMBER OF NEB REPLICA =', neb_env%number_of_replica
447 IF (plot_rel_energy)
THEN
448 cpassert(
SIZE(distances) == neb_env%number_of_replica - 1)
449 cpassert(
SIZE(energies) == neb_env%number_of_replica)
450 cpassert(
SIZE(temperatures) == neb_env%number_of_replica)
451 ener_min = minval(energies(:))
452 ener_range = maxval(energies(:)) - ener_min
455 WRITE (output_unit,
'(T2,A,T22,A,T35,A,T52,A)') &
456 'REPLICA',
'ENERGY [au]',
'TEMPERATURE [K]',
'o-------------------------> E'
457 DO i = 1,
SIZE(distances)
458 plt = floor((energies(i) - ener_min)/ener_range*25)
461 IF (energies(i) - energies(i - 1) > 0)
THEN
467 IF (energies(i + 1) - energies(i) < 0)
THEN
475 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"X"
478 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"x"
480 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"O"
482 WRITE (output_unit,
'(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') &
483 i, energies(i), l_ener, temperatures(i), trim(line)
484 WRITE (output_unit,
'(T2,A,1X,F16.6,T52,A)') &
485 "DISTANCE = ", distances(i),
"|"
487 plt = floor((energies(neb_env%number_of_replica) - ener_min)/ener_range*25)
489 IF (energies(neb_env%number_of_replica) - energies(neb_env%number_of_replica - 1) > 0)
THEN
495 WRITE (line,
'(A,A,A)')
"|", repeat(
" ", plt),
"O"
496 WRITE (output_unit,
'(T2,I7,T10,F18.8,1X,A,T34,F16.6,T52,A)') &
497 neb_env%number_of_replica, energies(neb_env%number_of_replica), &
498 l_ener, temperatures(neb_env%number_of_replica), trim(line)
499 WRITE (output_unit,
'(T52,A)')
"v Nr."
500 WRITE (output_unit,
'(T2,A,T44,2(1X,I4))') &
501 "NUMBER OF LOCAL MAXIMA (X) and MINIMA (x):", n_max, n_min
503 WRITE (output_unit,
'( A,T17,4F16.6)') &
504 ' DISTANCES REP =', distances(1:min(4,
SIZE(distances)))
505 IF (
SIZE(distances) > 4)
THEN
506 WRITE (output_unit,
'( T17,4F16.6)') distances(5:
SIZE(distances))
508 WRITE (output_unit,
'( A,T17,4F16.6)') &
509 ' ENERGIES [au] =', energies(1:min(4,
SIZE(energies)))
510 IF (
SIZE(energies) > 4)
THEN
511 WRITE (output_unit,
'( T17,4F16.6)') energies(5:
SIZE(energies))
514 WRITE (output_unit,
'( A,T33,4(1X,F11.5))') &
515 ' REPLICA TEMPERATURES (K) =', temperatures(1:min(4,
SIZE(temperatures)))
516 DO i = 5,
SIZE(temperatures), 4
517 WRITE (output_unit,
'( T33,4(1X,F11.5))') &
518 temperatures(i:min(i + 3,
SIZE(temperatures)))
522 WRITE (output_unit,
'( A,T56,F25.14)') &
523 ' BAND TOTAL ENERGY [au] =', sum(energies(:) + ekin(:)) + &
524 neb_env%spring_energy
525 WRITE (output_unit, fmt=
'(A,A)')
' **************************************', &
526 '*****************************************'
528 DEALLOCATE (temperatures)
532 extension=
".ener", file_form=
"FORMATTED")
534 WRITE (line,
'(I0)') 2*neb_env%number_of_replica - 1
535 WRITE (ener,
'(I10,'//trim(line)//
'(1X,F20.9))') istep, &
544 root_section=neb_env%root_section, &
549 IF (output_unit > 0)
CALL m_flush(output_unit)
550 CALL timestop(handle)
568 INTEGER,
INTENT(IN) :: i_rep, ienum, iw
569 LOGICAL,
INTENT(IN) :: use_colvar
572 REAL(kind=
dp),
DIMENSION(3) :: r
575 WRITE (iw,
'(/,T2,"NEB|",75("*"))')
576 WRITE (iw,
'(T2,"NEB|",1X,A,I0,A)') &
577 "Geometry for Replica Nr. ", ienum,
" in Angstrom"
578 DO iatom = 1,
SIZE(particle_set)
580 WRITE (iw,
'(T2,"NEB|",1X,A10,5X,3F15.9)') &
581 trim(particle_set(iatom)%atomic_kind%name), r(1:3)*
angstrom
584 WRITE (iw,
'(/,T2,"NEB|",1X,A10)')
"COLLECTIVE VARIABLES:"
585 WRITE (iw,
'(T2,"NEB|",16X,3F15.9)') &
586 (coords%int(j, i_rep), j=1,
SIZE(coords%int(:, :), 1))
588 WRITE (iw,
'(T2,"NEB|",75("*"))')
604 INTEGER,
INTENT(IN) :: irep, n_rep, istep
606 CHARACTER(len=*),
PARAMETER :: routinen =
'handle_band_file_names'
608 CHARACTER(LEN=default_path_length) :: output_file_path, replica_proj_name
609 INTEGER :: handle, handle2, i, ierr, j, lp, unit_nr
614 CALL timeset(routinen, handle)
618 CALL force_env_get(f_env%force_env, root_section=root_section)
619 j = irep + (rep_env%local_rep_indices(1) - 1)
621 replica_proj_name = get_replica_project_name(rep_env, n_rep, j)
622 lp = len_trim(replica_proj_name)
624 c_val=trim(replica_proj_name))
625 logger%iter_info%project_name = trim(replica_proj_name)
628 output_file_path = replica_proj_name(1:lp)//
".out"
630 c_val=trim(output_file_path))
631 IF (logger%default_global_unit_nr > 0)
THEN
632 CALL close_file(unit_number=logger%default_global_unit_nr)
635 CALL open_file(file_name=output_file_path, &
636 file_status=
"UNKNOWN", &
637 file_action=
"WRITE", &
638 file_position=
"APPEND", &
639 unit_number=logger%default_global_unit_nr)
640 WRITE (unit=logger%default_global_unit_nr, fmt=
"(/,(T2,A79))") &
641 "*******************************************************************************", &
642 "** BAND EVALUATION OF ENERGIES AND FORCES **", &
643 "*******************************************************************************"
644 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A,T79,A)")
"**",
"**"
645 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A,T79,A)")
"**",
"**"
646 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A,I5,T41,A,I5,T79,A)") &
647 "** Replica Env Nr. :", rep_env%local_rep_indices(1) - 1,
"Replica Band Nr. :", j,
"**"
648 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A,I5,T79,A)") &
649 "** Band Step Nr. :", istep,
"**"
650 WRITE (unit=logger%default_global_unit_nr, fmt=
"(T2,A79)") &
651 "*******************************************************************************"
655 SELECT CASE (f_env%force_env%in_use)
657 DO i = 1, f_env%force_env%mixed_env%ngroups
658 IF (
modulo(i - 1, f_env%force_env%mixed_env%ngroups) == &
659 f_env%force_env%mixed_env%group_distribution(f_env%force_env%mixed_env%para_env%mepos))
THEN
660 sub_logger => f_env%force_env%mixed_env%sub_logger(i)%p
661 sub_logger%iter_info%project_name = replica_proj_name(1:lp)//
"-r-"//trim(adjustl(
cp_to_string(i)))
663 unit_nr = sub_logger%default_global_unit_nr
664 IF (unit_nr > 0)
THEN
666 output_file_path = replica_proj_name(1:lp)//
"-r-"//trim(adjustl(
cp_to_string(i)))//
".out"
667 CALL open_file(file_name=output_file_path, &
668 file_status=
"UNKNOWN", &
669 file_action=
"WRITE", &
670 file_position=
"APPEND", &
672 sub_logger%default_global_unit_nr = unit_nr
680 CALL timestop(handle)
692 FUNCTION get_replica_project_name(rep_env, n_rep, j)
RESULT(replica_proj_name)
694 INTEGER,
INTENT(IN) :: n_rep, j
695 CHARACTER(LEN=default_path_length) :: replica_proj_name
697 CHARACTER(LEN=default_string_length) :: padding
698 INTEGER :: i, lp, ndigits
702 replica_proj_name = rep_env%original_project_name
704 ndigits = ceiling(log10(real(n_rep + 1, kind=
dp))) - &
705 ceiling(log10(real(j + 1, kind=
dp)))
710 lp = len_trim(replica_proj_name)
711 replica_proj_name(lp + 1:len(replica_proj_name)) =
"-BAND"// &
713 END FUNCTION get_replica_project_name
727 CHARACTER(LEN=default_path_length) :: replica_proj_name
728 INTEGER :: handle2, ierr, irep, n_rep, n_rep_neb, &
733 n_rep_neb = neb_env%number_of_replica
738 output_unit = logger%default_global_unit_nr
739 IF (output_unit > 0)
THEN
740 WRITE (unit=output_unit, fmt=
'(/,(T2,A79))') &
741 "*******************************************************************************", &
742 "** MAPPING OF BAND REPLICA TO REPLICA ENV **", &
743 "*******************************************************************************"
744 WRITE (unit=output_unit, fmt=
'(T2,A,I6,T32,A,T79,A)') &
745 "** Replica Env Nr.: ", rep_env%local_rep_indices(1) - 1, &
746 "working on the following BAND replicas",
"**"
747 WRITE (unit=output_unit, fmt=
'(T2,A79)') &
750 DO irep = 1, n_rep_neb, n_rep
751 replica_proj_name = get_replica_project_name(rep_env, n_rep_neb, irep + rep_env%local_rep_indices(1) - 1)
752 IF (output_unit > 0)
THEN
753 WRITE (unit=output_unit, fmt=
'(T2,A,I6,T32,A,T79,A)') &
754 "** Band Replica Nr.: ", irep + rep_env%local_rep_indices(1) - 1, &
755 "Output available on file: "//trim(replica_proj_name)//
".out",
"**"
758 IF (output_unit > 0)
THEN
759 WRITE (unit=output_unit, fmt=
'(T2,A79)') &
761 "*******************************************************************************"
762 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.
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
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